Nonlinear Hall effect in topological Dirac semimetals in parallel magnetic field

TL;DR

This paper theoretically investigates the nonlinear Hall effect in 2D topological Dirac semimetals under in-plane magnetic fields, deriving quantum kinetic equations and proposing experimental tests in various materials.

Contribution

The authors develop a quantum kinetic framework to analyze second-harmonic responses and nonlinear Hall effects in topological Dirac semimetals with in-plane magnetic fields.

Findings

Derived the quantum kinetic equation for the Wigner distribution function.

Analyzed Berry curvature dipole and field-induced current contributions.

Proposed experimental measurements of Hall resistivity dependence on magnetic field.

Abstract

We compute the second-harmonic response of two-dimensional topological Dirac semimetals subjected to an external in-plane magnetic field. The quantum kinetic equation for the Wigner distribution function is derived and then solved to evaluate the second-order electric-field contributions to the current density. Both the Berry curvature dipole and the field-induced terms in the current are analyzed across a broad range of model parameters. We propose that our theory can be tested experimentally by measuring the dependence of the anomalous Hall resistivity on the in-plane magnetic field in the surface states of the topological insulator SnTe, in WTe$_2$ and WSe$_2$ monolayers, as well as in the Kondo lattice material Ce$_3$Bi$_4$Pd$_3$ at very low temperatures.

Tables

Table 1. Table 1: The transformation properties of the Pauli matrices parametrizing the Hamiltonian and the momentum counted relative to ± 𝑲 \pm\bm{K} under the m x m_{x} mirror operation and the time-reversal operation, 𝒯 \mathcal{T} . As both operations exchange the valley if the original Pauli matrices parametrize the H 𝑲 H_{\bm{K}} the transformed ones parametrize H − 𝑲 H_{-\bm{K}} . Thanks to our choice of the basis, Eq. ( 102 ) the Pauli matrices transform trivially in the above sense.

	${\sigma }_x$	${\sigma }_y$	${\sigma }_z$	$k_x$	$k_y$	$B_x$	$B_y$
$m_x$	${\sigma }_x$	${\sigma }_y$	${\sigma }_z$	$-k_x$	$k_y$	$B_x$	$-B_y$
$𝒯$	${\sigma }_x$	$-{\sigma }_y$	${\sigma }_z$	$-k_x$	$-k_y$	$-B_x$	$-B_y$

Equations

$$\begin{array}{rl}{\hat{H}}_{\pm K}(\mathbf{k})& =\pm \alpha k_x{\hat{\sigma }}_0\pm v{k}_x{\hat{\sigma }}_x+v{k}_y{\hat{\sigma }}_y+m{\hat{\sigma }}_z.\end{array}$$

$${\epsilon }_{\mathbf{k}s}^{(\eta )}=s\alpha k_x+(-1{)}^{\eta }\sqrt{m^2+v^2{k}^2}.$$

$$\begin{array}{rl}|u_{s1}(\mathbf{k})⟩& =\frac{1}{\sqrt{2}}\left(\begin{array}{rl}& -s\sqrt{1-\frac{m}{b_{\mathbf{k}}}}{e}^{-is{\phi }_{\mathbf{k}}}\\ & \sqrt{1+\frac{m}{b_{\mathbf{k}}}}\end{array}\right),\\ |u_{s2}(\mathbf{k})⟩& =\frac{1}{\sqrt{2}}\left(\begin{array}{rl}& s\sqrt{1+\frac{m}{b_{\mathbf{k}}}}{e}^{-is{\phi }_{\mathbf{k}}}\\ & \sqrt{1-\frac{m}{b_{\mathbf{k}}}}\end{array}\right).\end{array}$$

$${\hat{v}}_x^{(s)}=\left[\begin{array}{cc}s\alpha -\frac{v^2{k}_x}{b_{\mathbf{k}}}& \frac{v}{k}\left(\frac{m{k}_x}{b_{\mathbf{k}}}-is{k}_y\right)\\ \frac{v}{k}\left(\frac{m{k}_x}{b_{\mathbf{k}}}+is{k}_y\right)& s\alpha +\frac{v^2{k}_x}{b_{\mathbf{k}}}\end{array}\right].$$

$${\hat{v}}_y^{(s)}=\left[\begin{array}{cc}-\frac{v^2{k}_y}{b_{\mathbf{k}}}& \frac{v}{k}\left(\frac{m{k}_y}{b_{\mathbf{k}}}+is{k}_x\right)\\ \frac{v}{k}\left(\frac{m{k}_y}{b_{\mathbf{k}}}-is{k}_x\right)& \frac{v^2{k}_y}{b_{\mathbf{k}}}\end{array}\right].$$

$${\Omega }_z^{(\eta )}(\mathbf{k})=\left(\frac{sm{v}^2}{2b_{\mathbf{k}}^3}\right)(-1{)}^{1+\eta }.$$

$$w_{\alpha \beta }(\vec{k},\vec{r};ϵ,t)$$

$$\times ⟨{\psi }_{\beta }^{†}\left(\vec{r}+\frac{\vec{\rho }}{2},t+\frac{\eta }{2}\right){\psi }_{\alpha }\left(\vec{r}-\frac{\vec{\rho }}{2},t-\frac{\eta }{2}\right)⟩.$$

$${\hat{w}}_{\mathbf{k}ϵ}^{(0)}=f(ϵ)\left(\begin{array}{cc}\delta \left(ϵ-{ϵ}_{\mathbf{k}s}^{(1)}\right)& 0\\ 0& \delta \left(ϵ-{ϵ}_{\mathbf{k}s}^{(2)}\right)\end{array}\right),$$

$$\displaystyle\left(\frac{\partial}{\partial t}+\frac{1}{\tau}\right)\hat{w}_{\mathbf{k}\epsilon}(t)+i\left\{{\mathbf{b}}_{\mathbf{k}}\cdot{\mbox{$\sigma$}},\hat{w}_{\mathbf{k}\epsilon}(t)\right\}_{-}$$

$$\displaystyle=-\alpha\tilde{\nabla}_{x}\hat{w}_{\mathbf{k}\epsilon}(t)-\frac{v}{2}\left\{\hat{\mbox{$\sigma$}},\tilde{\mbox{$\nabla$}}\hat{w}_{\mathbf{k}\epsilon}(t)\right\}_{+}.$$

$$\tilde{\mbox{$\nabla$}}\hat{w}_{\mathbf{k}\epsilon}(t)=\frac{e{\mathbf{E}}(t)}{\omega}\left(\hat{w}_{\mathbf{k}\epsilon+\frac{\omega}{2}}(t)-\hat{w}_{\mathbf{k}\epsilon-\frac{\omega}{2}}(t)\right).$$

$$\mathbf{E}(t)=\vec{\mathcal{E}}{e}^{-i\omega t}+c.c.,$$

$${\left[{\hat{w}}_{\mathbf{k}ϵ}^{(1)}\right]}_{aa}$$

$${\left[{\hat{w}}_{\mathbf{k}ϵ}^{(1)}\right]}_{a\overline{a}}$$

$$\mathbf{j}(t)=-e\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}\underset{-\infty }{\overset{\infty }{\int }}\text{Tr}[\hat{\mathbf{v}}{\hat{w}}_{\mathbf{k}ϵ}(t)]dϵ,$$

$${\mathcal{D}}_i\equiv \sum _{s=\pm }\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}{\Omega }_z(\mathbf{k})\frac{\partial f({ϵ}_{\mathbf{k}s})}{\partial k_i},$$

$$j_x^{(2\omega )}={\chi }_{xxy}(\omega ){\mathcal{E}}_x{\mathcal{E}}_y.$$

$$[{\hat{v}}_y{]}_{21}[{\hat{v}}_x{]}_{12}=-\frac{v^2}{b_{\mathbf{k}}}\left(ism+\frac{v^2{k}_x{k}_y}{b_{\mathbf{k}}}\right),$$

$$[{\hat{v}}_y{]}_{12}[{\hat{v}}_x{]}_{21}=\frac{v^2}{b_{\mathbf{k}}}\left(ism-\frac{v^2{k}_x{k}_y}{b_{\mathbf{k}}}\right).$$

$${\chi }_{xxy}$$

$$\times \underset{-\infty }{\overset{\infty }{\int }}\left(\frac{v{k}_x}{b_{\mathbf{k}}}\right)\mathcal{W}(\mathbf{k}ϵ;\omega )dϵ,$$

$$\begin{array}{rl}\mathcal{W}(\mathbf{k}ϵ;\omega )& =\sum _{a=1}^2(-1{)}^a{\left[{\hat{w}}_{\mathbf{k}ϵ+\omega }^{(0)}+{\hat{w}}_{\mathbf{k}ϵ-\omega }^{(0)}-2{\hat{w}}_{\mathbf{k}ϵ}^{(0)}\right]}_{aa}\end{array}.$$

$${\chi }_{xxy}$$

$$\times \left[f({ϵ}_{\mathbf{k}s}+\omega )+f({ϵ}_{\mathbf{k}s}-\omega )-2f({ϵ}_{\mathbf{k}s})\right].$$

$${\chi }_{xxy}\approx \sum _{s=\pm }\frac{sm{e}^3{v}^2\tau }{2(1-i\omega \tau )}\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}\frac{v^2{k}_x}{b_{\mathbf{k}}^3}{f}^{\prime \prime }({ϵ}_{\mathbf{k}s}).$$

$$\frac{v^2{k}_x}{b_{\mathbf{k}}}{f}^{\prime \prime }({ϵ}_{\mathbf{k}})=\frac{\partial }{\partial k_x}{f}^{\prime }({ϵ}_{\mathbf{k}s})-s\alpha f^{\prime \prime }({ϵ}_{\mathbf{k}s}).$$

$${\chi }_{xxy}^{(s)}$$

$$-\frac{m\alpha e^3{v}^2\tau }{2(1-i\omega \tau )}\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}\frac{1}{b_{\mathbf{k}}^2}{f}^{\prime \prime }({ϵ}_{\mathbf{k}s}).$$

$${\chi }_{xxy}^{(s)}$$

$$-\frac{m\alpha e^3{v}^2\tau }{2(1-i\omega \tau )}\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}\left\{\frac{f^{\prime }({ϵ}_{\mathbf{k}s})}{b_{\mathbf{k}}^3}+\frac{f^{\prime \prime }({ϵ}_{\mathbf{k}s})}{b_{\mathbf{k}}^2}\right\}.$$

$${\chi }_{xxy}^{(s)}=\frac{e^3\tau }{(1-i\omega \tau )}\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}\left(\frac{sm{v}^2}{b_{\mathbf{k}}^3}\right)\frac{\partial f({ϵ}_{\mathbf{k}s})}{\partial k_x}\equiv \frac{e^3\tau {\mathcal{D}}_x^{(s)}}{1-i\omega \tau },$$

$${\mathcal{D}}_x^{(s)}=\frac{3\alpha m|\mu |{\tilde{\gamma }}^2{v}^3\left(v^2-{\alpha }^2\right)}{4\pi ({\mu }^2{v}^2+m^2{\alpha }^2{)}^{5/2}}$$

$$[{\chi }_{xxy}{]}_{\omega \tau \ll 1}\approx \left(\frac{3e^3}{2\pi }\right)\frac{m\alpha \tau }{{\mu }^2}.$$

Full text

I Introduction

Studies of the Hall effect have long held an important place in both experimental and theoretical condensed matter physics. In the conventional Hall effect, the transverse conductivity $\sigma_{xy}$ describes the system’s response to an external electric field and arises from the orbital motion of conduction electrons coupled to an external magnetic field [1]. In contrast, the anomalous Hall effect originates either from the Zeeman coupling of electrons to an external magnetic field or, in the absence of a field, from a nonzero net magnetization. In clean systems, the finite contribution to the anomalous Hall effect is determined by the Berry curvature of the Bloch wavefunctions, which modifies the electron velocity [2, 3, 4, 5, 6, 7, 8, 9]. In both cases, however, it is evident that within linear response theory a nonzero Hall conductivity requires the breaking of time-reversal symmetry.

The next natural step in the study of the Hall effect has been to generalize the theory to nonlinear response [10, 11, 12, 13, 14, 15, 16]. In particular, it has been shown that a nonlinear Hall-like current can arise in systems lacking inversion symmetry, regardless of whether time-reversal symmetry is preserved [17, 18, 19]. This contribution to the conductivity is geometric in origin and is commonly referred to as the Berry curvature dipole (BCD) contribution [18, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. Most recently, it has also been theoretically proposed that the nonlinear Hall effect can be induced by superconducting fluctuations [31].

During the past few years several papers appeared reporting on the experimental observation of the nonlinear Hall effect in time-reversal invariant electronic systems [32, 33, 34, 35, 36, 37, 38]. Furthermore, giant nonlinear Hall effect has also been observed in an $f$-electron system Ce3Bi4Pd3 [39]. What makes this material special is the presence of strong electronic correlations, which, thanks to the partially filled $f$-electronic orbitals on Ce ions, drive the formation of the local magnetic moments. It has been argued that this material belongs to the class of Weyl-Kondo semimetals [39, 40]. Although it was suggested that BCD may be responsible for the magnitude of the nonlinear Hall response, the onset of the effect below $T^{*}\approx 3$K still remains not fully understood [41, 42].

The intrinsic BCD contribution to the nonlinear Hall response necessarily requires that the corresponding crystal is described by one of the following point groups $\{O,T,C_{1},C_{1v},C_{2v},C_{n},D_{n}\}$ with $n=1,2,3,4,6$  [18]. All of these point groups have a polar axis. For example, Ce3Bi4Pd3 is described by the $T_{d}$ point group which has $C_{3}$ and $C_{2}$ as its subgroups and therefore allows for nonzero Berry curvature dipole. More generally, in two-dimensional systems lacking inversion symmetry, the nonlinear Hall response can be tuned by an external in-plane magnetic field, as illustrated in Fig.1 [43, 44, 45]. Depending on its orientation, the field may either enhance or suppress the intrinsic contribution. The objective of the present work is to provide a detailed microscopic investigation of this effect.

In what follows, we consider a rather simple microscopic model of a two-dimensional disordered electronic system with the nontrivial geometric band structure which gives rise to both nonzero Berry curvature and BCD. Specifically, we will use the model of two tilted massive Dirac cones. To account for the effect of the in-plane magnetic field, we assume that the mass has a weak quadratic dispersion.

Using the theory of nonequilibrium systems, we derive the quantum kinetic equation and compute the second-order response of the current density to an applied electric field in the presence of an in-plane magnetic field. Our approach enables the analysis of both the $ac$- and $dc$-limits. We further identify the conditions under which the magnetic-field contribution can become of the same order as the geometric Berry curvature dipole contribution.

Our paper is organized as follows. In the next section, we introduce the model and provide the details of the derivation of the quantum kinetic equation for the Wigner distribution function. In Section III, we present the calculation of the nonlinear contribution to the current density. In the first part of Section III, we analyze the nonlinear current in the absence of a magnetic field and re-derive the Berry curvature dipole contribution, while in the second part we examine the contribution arising from the external magnetic field. Section IV is devoted to a discussion of our results and conclusions. Finally, Appendices A–E contain the technical details of the transport calculations and model derivation. Throughout this work, we use the units with $\hbar=c=1$.

II Model and Formalism

In this Section, we describe the formalism and provide the list of main equations which are used to calculate the second harmonic of the current density.

II.1 Tilted massive Dirac fermion

We begin by considering the following two-band model Hamiltonian in two spatial dimensions [46, 31]:

[TABLE]

In crystalline topological insulators Pauli matrices $\hat{\sigma}_{i}$ act in a space of underlying Kramers doublets and $m$ is determined by the magnitude of the ferroelectric distortion [47, 48]. Here we focus on the strained transition metal dichalcogenides (TMDs) realization of the Hamiltonian (1) (see App. E for details). In these systems the Pauli matrices act in the space of orbitals comprising the conduction and valence bands. Based on DFT calculations the conduction band at $\pm\bm{K}$ are formed by the $d_{z^{2}}$ atomic transition metal orbitals. The valence band at $\pm\bm{K}$ has a $d_{x^{2}-y^{2}}\pm id_{xy}$ orbital character [49, 50]. Both valence and conduction band wave functions are even in the basal plane mirror operation $(m_{z})$.

The model (1) ignored the actual spin. In principle one has two copies of the Hamiltonian (1) per spin. This Hamiltonian represents a minimal effective model for the surface states in crystalline topological insulator SnTe and monolayer WTe2 as it describes electrons with dispersion determined by the two tilted massive Dirac cones located at momenta $\pm K$ of the two-dimensional Brillouin zone, as shown in Fig. 2 [51, 50, 52].

Furthermore, we assume the scalar potential disorder $U(\mathbf{r})$ with the white-noise Gaussian statistics $\left\langle U(\mathbf{r})U(\mathbf{r}^{\prime})\right\rangle_{\textrm{dis}}={\delta(\mathbf{r}-\mathbf{r}^{\prime})}/({2\pi\nu_{F}\tau})$, where $\tau^{-1}$ is the disorder scattering rate and $\nu_{F}$ is the single-particle density of states per spin at the Fermi level. In the first part of this work, we will be interested in transport properties in the absence of a magnetic field, and for this reason it has not been included into Eq. (1). Lastly, parameter $\alpha$ in Eq. (1) accounts for the tilt of the Dirac cones.

We find it convenient to work in the eigenbasis of Hamiltonian (1), $\hat{H}_{sK}(\mathbf{k})=\sum_{\eta=1}^{2}\varepsilon_{\mathbf{k}s}^{(\eta)}|u_{s\eta}(\mathbf{k})\rangle\langle u_{s\eta}(\mathbf{k})|$, where $s=\pm 1$ is the valley index and the eigenenergies are given by

[TABLE]

The corresponding eigenvectors are

[TABLE]

Here we have introduced the shorthand notation $b_{\mathbf{k}}=\sqrt{m^{2}+v^{2}k^{2}}$ and $\tan\varphi_{\mathbf{k}}=k_{y}/k_{x}$.

Next, we introduce the velocity operator which is defined according to $\hat{v}_{a}^{(s)}(\mathbf{k})=\partial_{k_{a}}\hat{H}_{s}(\mathbf{k})$. Denoting $\mathbf{n}=\mathbf{k}/k$, the matrix elements of the velocity operator $\hat{v}_{x}^{(s)}$ in the basis (3) are given by

[TABLE]

Similarly, for the matrix elements of the velocity operator $\hat{v}_{y}^{(s)}$ we have

[TABLE]

It is well known that the band structure described by Eq. (1) possesses nontrivial geometric phase. Introducing the Berry connection ${\vec{\mathbf{\cal A}}^{(\eta)}}(\mathbf{k})=i\langle\eta\mathbf{k}|{\vec{\nabla}}_{\mathbf{k}}|\eta\mathbf{k}\rangle$ and the Berry curvature $\Omega_{a}^{(\eta)}(\mathbf{k})=\epsilon^{abc}\partial_{b}{\cal A}_{c}^{(\eta)}(\mathbf{k})$, we find that $\Omega_{x}^{(\eta)}=\Omega_{y}^{(\eta)}=0$ and

[TABLE]

As one may have expected, the Berry curvature is independent of the tilt parameter $\alpha$ and is proportional to the energy gap in the Dirac spectrum and the valley index. The latter turns out to be crucial for tje nonlinear Hall effect.

II.2 Quantum kinetic equation

The key quantity in our subsequent analysis is the Wigner distribution function (WDF) which is defined according to

[TABLE]

Here $\psi_{\alpha}^{\dagger}(\mathbf{r},t)$ and $\psi_{\alpha}(\mathbf{r},t)$ are the usual creation and annihilation fermionic fields for a given valley and averaging is performed over the ground state of the Hamiltonian (1). Spatial homogeneity of WDF is restored after averaging over disorder. In the absence of external fields, we find for the WDF computed in the basis of Eq. (3):

[TABLE]

where $f(\epsilon)$ is the Fermi-Dirac distribution function. In the presence of a spatially uniform electric field, the WDF satisfies the following quantum kinetic equation [53, 54]

[TABLE]

Here $\{\hat{f},\hat{g}\}_{\pm}=\hat{f}\hat{g}\pm\hat{g}\hat{f}$ and

[TABLE]

In the kinetic equation (9), we have also neglected the collision integral beyond $\tau$-approximation, which is technically equivalent to neglecting vertex corrections due to disorder scattering. In the limit when electric field is weak, kinetic equation (9) can be solved by perturbation theory. For example, for the monochromatic electric field

[TABLE]

the first order correction to the WDF is

[TABLE]

where $a=(1,2)$, $\overline{1}=2$, $\overline{2}=1$, $z_{\omega}=-i\omega+\tau^{-1}$ and $Z_{\mathbf{k}\omega}^{(a)}=[z_{\omega}+2ib_{\mathbf{k}}\cos(\pi a)]^{-1}$.

II.3 Linear response

Given the Wigner distribution function, the current density can be computed as

[TABLE]

where $\hat{v}_{x}$ and $\hat{v}_{y}$ are the $x$- and $y$-components of the velocity operator with the matrix elements given by Eqs. (4,5). One can check the correctness of this definition by computing the current density in the linear response approximation. Without loss of generality we can ignore the tilt by setting $\alpha=0$. Specifically, in the limit $\omega\ll\mu$ we have reproduced the Drude formula $\sigma_{xx}=e^{2}{D}\nu_{F}/(1-i\omega\tau)$ and also linear intrinsic contribution to the anomalous Hall effect $\sigma_{xy}=se^{2}/4\pi$ (see Appendix A for details).

III Nonlinear response

III.1 Berry curvature dipole contribution

Before analyzing the role of the in-plane magnetic field, we first demonstrate the microscopic derivation of the geometric contribution to the second harmonic using the method of quantum kinetic equation. This contribution is determined by the Berry curvature dipole, which in two dimensions is defined as [18]

[TABLE]

with $i=x,y$. We find that it is proportional to the tilt parameter $\alpha$, Eq. (1).

Given our choice of the model, the Berry curvature dipole will be pointing along the $x$-axis (see Eq. (24) below). We find it more convenient to consider the situation when electric field lies in the $xy$-plane and compute the current along the $x$-axis. In this case the nonlinear Hall contribution to the current is defined by

[TABLE]

An interesting property of the quantum kinetic equation (9) is that the second order corrections to the WDF are formally given by the same expressions as the first order ones, Eq. (12), in which we replace $\hat{w}_{\mathbf{k}\epsilon}^{(0)}\to\hat{w}_{\mathbf{k}\epsilon}^{(1)}$. Straightforward analysis shows that the contribution to $\chi_{xxy}$ from the diagonal matrix elements of the velocity operator and WDF given by $[\hat{v}_{x}]_{11}[\hat{w}_{\mathbf{k}\epsilon}^{(2)}]_{11}+[\hat{v}_{x}]_{22}[\hat{w}_{\mathbf{k}\epsilon}^{(2)}]_{22}\propto s\alpha k_{x}^{2}$ which gives zero upon the summation over valleys, $s=\pm$. Below we will see that in order to get nonzero Hall response one will need to include an in-plane magnetic field. We therefore conclude that nonlinear correction to the current will necessarily originate from the scattering processes which involve the off-diagonal matrix elements of the WDF.

We thus analyze the remaining contribution $[\hat{v}_{x}]_{12}[\hat{w}_{\mathbf{k}\epsilon}^{(2)}]_{21}+[\hat{v}_{x}]_{21}[\hat{w}_{\mathbf{k}\epsilon}^{(2)}]_{12}$. With the help of expressions for $[\hat{w}_{\mathbf{k}\epsilon}^{(2)}]_{a\overline{a}}$ (see Appendix C for details), we find that the combination of velocity matrix elements $[\hat{v}_{x}]_{12}[\hat{v}_{x}]_{21}[\hat{v}_{y}]_{aa}\propto k_{y}$ gives zero upon integration over $k_{y}$. Therefore, the only nonzero contribution will arise from the following combinations: $[\hat{v}_{y}]_{12}[\hat{v}_{x}]_{21}[\hat{v}_{x}]_{aa}$ and $[\hat{v}_{x}]_{12}[\hat{v}_{y}]_{21}[\hat{v}_{x}]_{aa}$. In particular, with the help of expressions (4,5) we find

[TABLE]

Again keeping in mind the subsequent integrations over momentum components, one realizes that the terms $\propto k_{x}k_{y}$ do not contribute to the second harmonic since the resulting expression under the integral will depend only on $k_{x}$ and $k_{x}^{2}+k_{y}^{2}$. As a result, we derive the following expression for the response function $\chi_{xxy}$:

[TABLE]

where function ${\cal W}({\mathbf{k}\epsilon};\omega)$ is defined according to

[TABLE]

Equation (17) is derived in the limit $\alpha\ll v$ and ignores contributions of order $O(\alpha^{2})$.

In Fig. 3, we present the frequency dependence of the real part of the response function, $\chi_{xxy}(\omega)=\chi_{xxy}^{\prime}(\omega)+i\chi_{xxy}^{\prime\prime}(\omega)$. Notably, function $\chi_{xxy}^{\prime}(\omega)$ exhibits a peak at frequency $\omega\approx\mu$. This peak describes the resonant transitions between the valence and conduction band. The location of the peak is deduced from the fact that the denominator in the expression under the integral (17) approaches zero as $\tau\to\infty$. This resonant effect is not captured by the semiclassical Boltzmann equation approach [18] (see also [55] and references therein).

Low frequency limit.

We consider the limit of low frequencies first. Specifically, we take $\omega\ll m$ and for weak disorder $m\tau\gg 1$ we have $z_{2\omega}\ll 2b_{\mathbf{k}}$, so that the expression under the integral (17) simplifies to

[TABLE]

Here we assumed for definiteness that the chemical potential lies in the upper band, $\mu>0$. Expanding the distribution functions in powers of $\omega$ yields

[TABLE]

Next, we express function $\chi_{xxy}$ explicitly in terms of the integral over the Berry curvature. To do that we replace the expression under the integral using the following identity

[TABLE]

At this point it will be convenient to separate the contributions from different valleys, so we write $\chi_{xxy}=\sum\limits_{s=\pm}\chi_{xxy}^{(s)}$ with

[TABLE]

Integrating by parts in the first term and using Eq. (21) again we find:

[TABLE]

Let us first consider the first term in Eq. (23). Given the expressions for the Berry curvature, Eq. (6), and the Berry curvature dipole, Eq. (14), one immediately recognizes that the first integral in Eq. (23) is of geometric origin:

[TABLE]

where we calculated for the Berry curvature dipole per valley

[TABLE]

and we defined the parameter $\tilde{\gamma}=\sqrt{\mu^{2}v^{2}/(v^{2}-\alpha^{2})-m^{2}}$. This contribution does not depend on the valley index $s=\pm 1$, so that ${\cal D}_{x}^{(+)}={\cal D}_{x}^{(-)}={\cal D}_{x}/2$, where ${\cal D}_{x}$ is the total Berry curvature dipole. After subsequent summation over valleys, we find $\chi_{xxy}=e^{3}\tau{\cal D}_{x}/(1-i\omega\tau)$. For details of the evaluation of the momentum integral in Eq. (24) we refer the reader to Appendix B.

As for the remaining integral in Eq. (23), one notices that it is proportional to $\alpha$ and after subsequent rearrangement of the terms it will give the contribution $O(\alpha^{2})$ which we ignore. In principle, it can be shown that this contribution will exactly cancel the corresponding contribution from the similar term in Eq. (17), which we have ignored as well. This cancellation is not accidental and is guaranteed by the time-reversal symmetry of the system.

Equation (25) represents the only nonzero component of the Berry curvature dipole. Indeed, as expected on purely symmetry grounds, we find that the second harmonic response is linearly proportional to the tilt parameter $\alpha$ which breaks the mirror symmetry and gives rise to nonzero Berry curvature dipole contribution ${\vec{\cal D}}={\cal D}_{x}{\vec{e}}_{x}$. Finally, in the dc limit , $\omega\tau\ll 1$, and assuming also $v\gg\alpha$, $\mu\gg m$, we find for the response function after summation over the two valleys

[TABLE]

The calculation presented above can of course be repeated assuming the electric field is aligned along the $x$-axis, while the current is flowing along the $y$-axis, i.e., only the $y$-components of the velocity operator in Eq. (13) needs to be considered. As a result, the general expression for the nonlinear current can be written as

[TABLE]

Absolutely analogous expression can be derived for the dc second-order current, with $\left({\vec{\cal E}}\times{\vec{e}}_{z}\right)$ replaced by $\left({\vec{\cal E^{*}}}\times{\vec{e}}_{z}\right)$. We note that a factor of 2 difference with Ref. [18] originates from a slightly different definition of electric field and current.

High frequency limit.

We now turn our attention to the analysis of $\chi_{xxy}(\omega)$, Eq. (17), in the limit of high frequencies, $\omega\sim m$ and $\omega\tau\gg 1$ while keeping $\omega\ll\mu$. The latter condition still allows us to expand the arguments of the distribution functions up to the second order in powers of $\omega/\mu$. The remaining part of the calculation is very similar to that leading us to Eq. (23). As a result, we find that at high frequencies

[TABLE]

which is purely dissipative.

As the frequency increases, the second-order response exhibits a sharp resonance around $\omega\approx\mu$, as shown in Fig. 3. At even higher frequencies, $\omega\gg\mu$, we find that $[\chi_{xxy}]_{\omega\tau\gg 1}\propto i/\omega^{3}$.

III.2 Effect of in-plane magnetic field

As we have discussed so far, Eq. (27) represents the only nonzero second-order Hall response without an external magnetic field ${\mathbf{B}}$. In order to elucidate the contribution from an in-plane magnetic field, we assume that it couples to spin $\hat{\bm{\sigma}}$ through the Zeeman coupling. We will discuss the following model Hamiltonian

[TABLE]

which represents the extension of the model of the TMDs, Eq. (1), derived in Appendix E based on the symmetry considerations applied to a TMD monolayer. The mass dispersion in Eq. (29), $\beta k^{2}\hat{\sigma}_{z}/2$, which is naturally present in realistic band structures, will turn out to be crucial for our analysis. On a substrate there is a finite admixture of the $m_{z}$ odd orbitals $d_{xz}$ and $d_{yz}$ that is necessary for a finite $\gamma$. In principle, we could have also considered nonlinear-in-momentum corrections to the Dirac dispersion, however, while amenable to a numerical analysis, producing analytically tractable calculations for this case turns out to be quite challenging.

Expressions for the basis functions can be easily generalized by noting that in the absence of the tilt and mass dispersion ($\alpha=\beta=0$) the effect of magnetic field results in the shift of momentum:

[TABLE]

Primarily for computational convenience, we consider the setup in which the direction of the electric field points along the $x$-axis, ${\vec{\cal E}}={\cal E}{\vec{e}}_{x}$, and compute the second harmonic of the current along the $y$-axis $j_{y}^{(2\omega)}$ (see Fig. 1). Given that the Berry curvature dipole is pointing along the $x$-axis, this will also allow us to make a direct comparison with the BCD contribution to the current, Eq. (27).

The weakly dispersive mass term $\propto\beta k^{2}$ modifies the velocity operators according to

[TABLE]

and

[TABLE]

Here we introduced $q=\sqrt{(k_{x}+\gamma B_{y}/v)^{2}+(k_{y}+\gamma B_{x}/v)^{2}}$, $m_{\mathbf{k}}=m+\beta k^{2}/2$, and $\epsilon_{\mathbf{q}}=(v^{2}q^{2}+m_{\mathbf{k}}^{2})^{1/2}$. In what follows, we will be focusing solely on contributions from the magnetic field and for this reason we will set $\alpha=0$.

Second order corrections.

We proceed by looking at the second-order corrections in powers of electric field to the $y$-component of the current. As in the derivation of the Berry curvature dipole contribution to the nonlinear current, one should distinguish between the contributions from the diagonal parts of the WDF $[\hat{w}^{(2)}]_{aa}$ and the off-diagonal ones $[\hat{w}^{(2)}]_{\overline{a}a}$. We show in Appendix C that the off-diagonal contribution has additional small factor $1/(\mu\tau)^{2}\ll 1$ compared to the diagonal one, see Eq. (91).

Thus we focus on the contribution from the diagonal components of the WDF. We have:

[TABLE]

In this expression, we have used the fact that $[\hat{v}_{y}^{(s)}]_{11}+[\hat{v}_{y}^{(s)}]_{22}=0$, see Eq. (32). In the limit of small magnetic fields, we find to linear order in $\bf B$

[TABLE]

The dimensionless function $J(\omega)$ is given by Eq. (75) from Appendix C, and its frequency dependence is shown in Fig. 4. For frequencies $\omega\leq\mu-m$, the momentum integral for $J(\omega)$ can be evaluated analytically, resulting in Eq. (76).

At small frequencies, $\omega\ll\mu$, we obtain a compact expression:

[TABLE]

To elucidate the physical origin of this result, we reproduce it (up to a factor of 2) using semiclassical Boltzmann approach in Appendix D. We show that it is determined solely by the dispersion relation near the Dirac nodes and is not related to the geometry of the band structure. Indeed, in the absence of tilt, the Berry curvature of Hamiltonian (29) has the opposite sign in two valleys even at finite magnetic field, hence the Berry curvature dipole contribution vanishes. However, this cancellation is accidental and is not enforced by any symmetry. We expect that the higher-order couplings to Zeeman field allowed by symmetry will result in finite geometrical contribution, and we leave this problem for future study.

Finally, we find that the longitudinal component of the current equals $j_{x}^{(2\omega)}({\mathbf{B}})=(3B_{y}/B_{x})j_{y}^{(2\omega)}({\mathbf{B}})$. This relation constitutes one of the experimental predictions of our theory.

IV Discussion and conclusions

In addition to the intrinsic magnetic field dependent contribution calculated in this work, realistic materials also exhibit extrinsic mechanisms – such as side-jump and scew-scattering processes – that can contribute comparably to the nonlinear Hall response and are often independent of the disorder scattering rate [56, 57]. Although we did not include disorder-assisted effects in our analysis, we expect the predicted dependence of the nonlinear Hall effect on the magnetic field direction to remain robust and observable at sufficiently high fields. A detailed investigation of these extrinsic contributions is left for future work.

To estimate the magnitude of the effect studied in this work, we can use the following dimensionless parameter

[TABLE]

that determines the ratio between the magnetic field and the Berry curvature dipole contributions to the nonlinear Hall response in the low-frequency limit $\omega\tau\ll 1$. In a disordered system ($m\tau\sim 1$) this parameter is quite small. However, in a relatively clean electronic system the smallness of the energy scale $\beta k_{F}^{2}$ can be partially compensated by the small value of $\alpha k_{F}$. For example, in WSe2 and WTe2 monolayers we can estimate $\mu\tau\approx 100$, $\beta k_{F}^{2}\sim 0.05\mu\approx 80$ meV and $\alpha k_{F}\approx 5$ meV [55], which yields $\eta\approx 10(\gamma B/m)$. For magnetic fields in the range $B\sim 1\div 10$ T, we find $\eta\sim 0.1\div 1$, which should be experimentally observable. The situation for Ce3Bi4Pd3 seems even more favorable, as the electron’s velocity is significantly reduced due to the strong hybridization between the conduction and $f$-electrons of the cerium ions.

Finally, we highlight the advantages of the quantum kinetic equation approach employed in this work. Unlike the semiclassical Boltzmann formalism, this method provides a unified framework that remains valid across a wide range of light frequencies, magnetic field strengths, and chemical potentials. It also naturally incorporates various forms of disorder through the collision integral, enabling a more realistic modeling of experimental conditions. Crucially, the quantum kinetic equation captures resonant interband transitions between the valence and conduction bands – phenomena beyond the reach of Boltzmann-based treatments [54]. Looking ahead, we plan to extend this method to explore orbital effects of strong magnetic fields in topological semimetals and flat-band systems [58, 59, 60], as well as to investigate how generic disorder influences the quantization of the photogalvanic effect [61, 62, 63].

In conclusion, we have explored the nonlinear Hall effect in two-dimensional Dirac materials under the influence of an external in-plane magnetic field. Using the quantum kinetic equation, we identified two distinct contributions to the second harmonic response. The first, of geometric origin, arises intrinsically from the Berry curvature dipole; the second is magnetic-field-induced and increases linearly with the field strength. Notably, the latter contribution exhibits a strong directional dependence, enabling the geometric response to be either enhanced or suppressed by adjusting the orientation of the applied field. Together, these findings suggest a practical and tunable tool for controlling nonlinear Hall phenomena in Dirac systems, offering clear experimental pathways for future investigation.

V Acknowledgments

We thank Leonid Golub for bringing several references to our attention and enlightening discussions. This work was financially supported by the National Science Foundation grant DMR-2400484 (M.D.) M. D. has performed this work in part at Aspen Center for Physics, which is supported by the National Science Foundation grant PHY-2210452. M. K. acknowledges the support of the grant NSF-BSF DMR-2023693. V. K. was supported by the Pittsburgh Quantum Institute Community Collaboration Award.

Appendix A Drude conductivity and anomalous Hall effect

It is useful to briefly discuss transport at the level of the linear response approximation. We start with the calculation of the longitudinal conductivity $\sigma_{xx}=j_{x}/{\cal E}_{x}$. Using the definitions (12) and (13) in the main text, after a short calculation we find

[TABLE]

and $b_{\mathbf{k}}=\sqrt{m^{2}+(vk)^{2}}$. We note that this expression includes both contributions from diagonal as well as off-diagonal elements of $\hat{v}_{\alpha}$. In deriving this formula we assumed the limit of $T\to 0$ and also that the chemical potential is in the conduction band. In the low-frequency limit $\omega\ll vk_{F}$ we recover the Drude formula

[TABLE]

Here ${D}=v_{F}^{2}\tau/2$ is the diffusion coefficient in two dimensions, $v_{F}=(\partial b_{\mathbf{k}}/\partial k)_{k=k_{F}}=v^{2}k_{F}/\mu$ is the Fermi velocity, $\nu_{F}=\mu/(\pi v^{2})$ is the total density of states at the Fermi level and we assumed $\mu\tau\gg 1$.

Let us now evaluate linear Hall conductivity $\sigma_{xy}=j_{y}/{\cal E}_{x}$. Again taking the limit $T\to 0$ we find

[TABLE]

where $\mu$ is the chemical potential. Note that only off-diagonal components of the matrices $\hat{v}_{x}$ and $\hat{v}_{y}$ contribute to $\sigma_{xy}$. The prefactor of $2$ in front appears as a result of the partial cancellation of factor of $4$ which comes in combination with $b_{\mathbf{k}}$ in the numerator and factor of $2$ in the denominator of $[\hat{w}^{(1)}]_{\textrm{a}\overline{\textrm{a}}}$, Eq. (12). We emphasize that integral in Eq. (39) can be rewritten as the integral over the Berry curvature, after performing the integration by parts. In the simplest case $\alpha=0$, the integral can be easily evaluated with the following result

[TABLE]

We consider the static limit $\omega\to 0$ in the absence of disorder, $\tau\to\infty$, but keeping $\omega\tau\ll 1$. Then, setting $\mu=m+0$ in Eq. (40), we readily recover the expression for the intrinsic contribution to the anomalous Hall conductivity for a single valley:

[TABLE]

The subsequent summation over the valley index leads to the cancellation of the corresponding contributions, as guaranteed by the time-reversal symmetry. The second-order Hall effect, which is the main focus of our work, survives even in time-reversal-symmetric systems.

Appendix B Berry curvature dipole contribution

In this Appendix, we provide the details of the calculation of the integral in Eq. (24). We consider the following integral:

[TABLE]

In the low-temperature limit $T\to 0$, the Fermi distribution function can be approximated by the step function

[TABLE]

Thus, the integration in Eq. (42) is performed over the surface

[TABLE]

We will integrate over $k_{x}$ and $k_{y}$ separately. We start with integrating over $k_{y}$ first by introducing a new variable

[TABLE]

It follows that

[TABLE]

where $m_{x}=\sqrt{m^{2}+v^{2}k_{x}^{2}}$. We split Eq. (46) into two integrals:

[TABLE]

Note that $J_{\textrm{B.d.}}^{(+)}+J_{\textrm{B.d.}}^{(-)}=0$ when we set $\alpha=0$. In order to evaluate the first integral, we change the integration variable according to

[TABLE]

Given Eq. (48), it also follows that

[TABLE]

Furthermore, we have

[TABLE]

where

[TABLE]

Keeping in mind that $\alpha<v$, we note that

[TABLE]

It will be convenient to use the shorthand notation

[TABLE]

After we substitute these expressions into the integral and make a change of variables $x=q_{x}\sqrt{v^{2}-\alpha^{2}}$, we have

[TABLE]

with $a=v^{2}u\mu/\alpha^{2}\gamma$. We note that

[TABLE]

which follows from the definition (47). Changing under the integral $q\to-q$, we find

[TABLE]

This expression can be further simplified by using the following identities:

[TABLE]

It then follows that

[TABLE]

We emphasize that this expression should be viewed as the leading-order contribution to the second harmonic at small $\alpha$. Taking the limit $\alpha\ll v$ and $m\ll\mu$, we find

[TABLE]

where we took into account that $[\gamma]_{m\ll\mu}\approx\mu$.

Appendix C Corrections to the Wigner distribution function and nonlinear current in the presence of magnetic field

In the presence of magnetic field, the right hand side of the kinetic equation (9) acquires an additional term -\beta\left\{\mathbf{k}\hat{\sigma}_{z},\tilde{\mbox{\nabla}}\hat{w}_{\mathbf{k}\epsilon}(t)\right\}_{+}/2. Because of this term, the corresponding expressions for the first-order corrections to the Wigner distribution function can be compactly expressed in terms of the matrix elements of the velocity operators as follows:

[TABLE]

Here the matrix elements of the velocity operators are given by Eqs. (31)-(32) of the main text.

The second-order corrections to the WDF can be found similarly to Eq. (59). For the correction to the first diagonal element of the WDF, we have

[TABLE]

Here, $z_{2\omega}=(1-2i\omega\tau)/\tau$. Similarly, for the off-diagonal elements we find

[TABLE]

We can employ these formulas to evaluate the Hall current. Specifically, let us compute the contributions which contain the diagonal elements of the velocity operators and consistent with our discussion in the main text we take the electric field to point along the $x$-axis, ${\vec{\cal E}}={\cal E}{\vec{e}}_{x}$.

We start with the contribution to the $y$-component of the current from the diagonal parts of the WDF:

[TABLE]

Using (59), (60) and (61) we have

[TABLE]

where for convenience we introduced

[TABLE]

Similarly, we find

[TABLE]

Given that matrix $\hat{v}_{y}$ is traceless, Eq. (32), it follows

[TABLE]

Inserting these expressions into (62) and assuming that chemical potential lies in the conduction band $({\mathrm{a}}=2)$ and assuming that the frequenies do not exceed $\mu-m$ we then have

[TABLE]

In order to keep the subsequent expressions compact we will ignore the terms which are proportional to the $y$-component of the magnetic field since, as one can easily check, the current will be proportional to $B_{x}$. We have ($\mathbf{q}=\mathbf{k}+{\mathbf{h}}$):

[TABLE]

where $u=v^{2}/\beta$, $h_{y}=(\gamma/v)B_{x}$, $\epsilon_{\mathbf{k}}=\sqrt{(vk)^{2}+m_{\mathbf{k}}^{2}}$ and used the expansion

[TABLE]

We now have to insert (68) into expression for the current. An analytic formula for the current can be obtained by changing the integration over momentum to the integration over energy $\xi=\sqrt{v^{2}k^{2}+m_{\mathbf{k}}^{2}}$. In particular, we have

[TABLE]

with $k^{2}=2(m_{\xi}-m)/\beta$.

The final result for the current can be expressed as a sum of three separate contributions. The first one originates from the expansion of the arguments of the distribution functions in powers of ${\mathbf{h}}$:

[TABLE]

For our discussion below it proves convenient to represent the delta functions as the derivatives of the Heaviside step functions and perform the integration by parts. It obtains:

[TABLE]

where $\Omega=\mathrm{max}(m,\mu-\omega)$.

The second integral is defined by the second term in (68):

[TABLE]

Lastly, the third integral will be given by

[TABLE]

Comparing (73) and (74) with (72), we see that the sum of the first two terms in (72) equals to $-J_{b}-J_{c}$. Introducing $J=J_{a}+J_{b}+J_{c}$ we have

[TABLE]

This integral can be computed analytically. Assuming that external frequency satisfies $\omega\leq\mu-m$ so that $\Omega=\mu-\omega$ we have

[TABLE]

where we introduced general notation $\overline{a}=a/u$, $r_{\mu}=\sqrt{\mu^{2}+2mu+u^{2}}/u$ and $f_{\mu}=\sinh^{-1}(\mu/\sqrt{u^{2}+2mu})$, and we remind that the parameter $u$ has been defined as $u=v^{2}/\beta$. Expression (76) can be further simplified in the limit of low frequencies $\omega\ll\mu$. As a result, for the current in this limit we find :

[TABLE]

In the limit $\beta\ll\min\{v^{2}/m,v/k_{F}\}$, Eq. (77) can be further simplified to

[TABLE]

Let us now discuss the contribution from the off-diagonal elements of the WDF. We will use the same approximations which lead us to (78). By definition we have

[TABLE]

We would like to remind the reader that off-diagonal parts of the WDF are the ones which gave rise to the Berry dipole contribution to the current. So, let us consider the combination under the integral in more details. Using (59) along with (31) and (32) we find

[TABLE]

It is clear that the terms which are linear in $s$ will give zero upon the summation over the valley index since we took $\alpha=0$. As a result the expression under the integral (79) simplifies to

[TABLE]

On the last step we took into account the definition of $[\hat{v}_{x}]_{22}$, used simple power counting to check that upon carrying out the momentum integration in (79) the result will be proportional to $h_{y}$ and we set $q_{x}=k_{x}$. Thus, expression (79) takes the following form

[TABLE]

and we introduced the shorthand notation

[TABLE]

Note that current from the off-diagonal components of the WDF (82) enters with the overall negative sign compared to (67).

It remains to compute the linear-in-magnetic field correction to $[j_{y}^{(2\omega)}]_{\mathrm{off-diag.}}$. We use the following expansion

[TABLE]

and insert it into into (82). In order to evaluate the integrals it will be useful to consider three different contributions. The first one originates from the second term in (84). After performing the integral over the direction of $\mathbf{k}$ and $\epsilon$ we have

[TABLE]

Here $u=v^{2}/\beta$, we assumed for simplicity $\mu-\omega\geq m$, made a change of variables from $k$ to $\xi=\sqrt{v^{2}k^{2}+m_{\mathbf{k}}^{2}}$ and used $\langle k_{x}^{2}k_{y}^{2}\rangle_{\hat{\mathbf{k}}}=k^{4}/8$. Similarly, the second integral is defined as follows

[TABLE]

Finally, the third integral that we need to consider appears as a result of expanding the distribution functions in powers of $h_{y}$:

[TABLE]

Next we write $\delta(\mu-\xi)=(-d/d\xi)\vartheta(\mu-\xi)$ and integrate by parts. This yields

[TABLE]

The first term on the right hand side equals $-I_{1}$. Furthermore, in the second term we replace $u(u+2m)=(m_{\xi}+u)^{2}-\xi^{2}$. It then follows:

[TABLE]

Taking into account that realistically $\mu\tau\gg 1$ we can approximate

[TABLE]

Since the resulting expression still remains quite cumbersome, we expand it in powers of $\beta$. This yields the following expression for $[j_{y}^{(2\omega)}]_{\mathrm{off-diag.}}$:

[TABLE]

We emphasize that this result is smaller than the diagonal contribution, Eq. (78), by a large factor $(\mu\tau)^{2}\gg 1$, and hence should be neglected to the leading order.

Appendix D Semiclassical Boltzmann approach

It is instructive to discuss how Eq. (78) can be obtained from the semiclassical Boltzmann theory. With the electric field along the $x$-axis, ${\vec{\cal E}}={\cal E}{\vec{e}}_{x}$, the corresponding expression for the second harmonic of the current density is given by [18]

[TABLE]

Here, $f_{0}(\epsilon)$ is the Fermi-Dirac distribution function and $\Omega_{z}(\mathbf{k})$ is the Berry curvature. Neglecting the tilt, $\alpha=0$, and assuming the chemical potential is in the conduction band, the dispersion of Hamiltonian (29) is given by $\epsilon_{\mathbf{q}}=\sqrt{v^{2}q^{2}+m_{\mathbf{k}}^{2}}$, with $m_{\mathbf{k}}=m+\beta k^{2}/2$ and $q$ given by Eq. (30). Note that the dispersion is identical for both valleys in the absence of tilt.

The first term in Eq. (92) describes the Berry curvature dipole contribution. It is straightforward to show that the Berry curvature of Hamiltonian (29) is strictly opposite in different valleys, even in presence of a finite magnetic field $\mathbf{B}$. Combined with the symmetric dispersion at $\alpha=0$, it implies that the Berry curvature dipole vanishes in this case. This cancellation is accidental and is not dictated by any symmetry. We expect the corresponding contribution to be nonzero once other symmetry-allowed terms with higher powers of $\mathbf{B}$ are included in Eq. (29).

The second term in Eq. (92) is determined solely by the band dispersion and vanishes in time-reversal-symmetric systems. Since finite Zeeman field in Eq. (29) breaks time reversal, this term produces a nonzero linear-in-$B$ contribution to the second harmonic. The contributions from the two valleys in this case merely add up resulting in an extra factor of 2.

The evaluation of the integral is most easily done in the zero-temperature limit, such that $f^{\prime}_{0}(\epsilon)=-\delta(\epsilon-\mu)$. After straightforward manipulations with the derivatives of the delta-functions, we find to the linear order in $\mathbf{B}$ (and quadratic in $\cal E$):

[TABLE]

This expression reproduces Eq. (35), up to a difference in factor of 2 which we cannot explain. We again stress that this result is derived in the absence of tilt, $\alpha=0$. We anticipate that a nonzero $\alpha$ will change the leading power of $\beta$ in the final expression.

To make the connection between the Boltzmann and quantum kinetic approaches, we introduce

[TABLE]

Integrating by parts with respect to $k_{x}$ and then with respect to $k_{y}$ one can show that ${\cal J}_{yxx}={\cal J}_{xxy}$ where

[TABLE]

Next, we rewrite the derivative of the distribution function in Eq. (94) into the following form:

[TABLE]

The first term in this expression can be rewritten as

[TABLE]

where $[\hat{v}_{x}^{(s)}]_{22}$, for a given valley $s$, has been defined in Eq. (31). On the other hand, the second term leads to

[TABLE]

After inserting Eq. (96) into Eq. (94) and using Eq. (98), we obtain

[TABLE]

which leads to

[TABLE]

where integration now is over the vicinity of the corresponding valley instead of the whole Brillouin zone. Inserting this result back to Eq. (92), we recover the corresponding expression for the $y$-component of the current in Eq. (67) obtained within the quantum kinetic approach, up to a factor of $2$. The reason for this discrepancy remains unclear.

Appendix E Derivation of the model Hamiltonian

In this appendix we derive the Hamiltonian used in the present paper to compute the BCD and in-plane field contributions using the $\bm{k}\cdot\bm{p}$ approximation [64, 65]. We start with the monolayer TMD ignoring the substrate. In this case the point group symmetry of the lattice is $D_{3h}$ which we represent in the following form, $D_{3h}=C_{3v}\times\{E,m_{z}\}$ where $E$ is the identity, $m_{x,y,z}$ represent the mirror in the plane that is orthogonal to $\hat{x}$, $\hat{y}$ and $\hat{z}$, respectively. And the $C_{3v}$ is the symmetry of the ammonia molecule. It is convenient to view the $C_{3v}$ group as generated by the three fold rotation, $C_{3}$ and the $m_{x}$ mirror. We also ignore the spin for the moment. The group of $\bm{K}$ is $C_{3}$ such that the Bloch states are the eigenstates of the three-fold rotation operation.

The transition metal $d$-orbitals are split into two groups. The orbitals $\{d_{z^{2}},d_{x^{2}-y^{2}},d_{xy}\}$ or even in $m_{z}$ and $\{d_{xz},d_{yz}\}$ are odd. Hence following the DFT results we focus on the first set of three orbitals even in $m_{z}$. At $\bm{K}$ the conduction band electrons has a $d_{z^{2}}$ character, and the valence band has a character $d_{x^{2}-y^{2}}-id_{xy}$ with the respective eigenvalues $1$ and $\exp(i4\pi/3)$ as $C_{3}$ operation eigenvalues. Here $r=\exp(i2\pi/3)$. Following we write the effective Hamiltonian operating in the space formed by the conduction and valence bands, $\Psi^{\pm\bm{K}}_{\bm{k}}=(\Phi^{c}_{\pm\bm{K}+\bm{k}},\Phi^{v}_{\pm\bm{K}+\bm{k}})^{tr}$, where $tr$ stands for the transposition and

[TABLE]

where the summation is over the locations of the $N$ transition metal ions, and the shorthand notation $d_{\pm}(\bm{r})=d_{x^{2}-y^{2}}(\bm{r})\pm id_{xy}(\bm{r})$ has been introduced. The basis at $\pm\bm{K}$ is not unique, yet our choice, (E) is convenient because of the simple relationship between the two basis sets

[TABLE]

The effective two by two Hamiltonian at $\pm\bm{K}$ is by construction,

[TABLE]

where $\mathcal{H}$ is the full Bloch Hamiltonian. We are employing the generic $\bm{k}\cdot\bm{p}$ method of getting $H_{\pm{K}}(\bm{k})$ which amounts to the series expansion of Eq. (103) around $\pm\bm{K}$ constrained by the symmetry of the problem.

We start with the constrains imposed by the group of $\bm{K}$ which contains just three elements and isomorphic $C_{3}$ group. We therefore focus on one, say $\bm{K}$ valley. And we only have to make sure that the Hamiltonian is invariant under the three fold rotation, $C_{3}$. It is convenient to split the Hamiltonian (103) into the diaogonal and off-diagonal parts, $H^{d}$ and $H^{o}$.

The diagonal part contains the zero strain part, $H^{d_{0}}$ and the contribution of the strain, $H^{str}$. At zero strain, it is insensitive to the orbital content, and therefore to the second order in momentum reads,

[TABLE]

In Eq. (104) we only keep the second term as the first one is inconsequential. To add the effect of the strain we classify the strain components by the group of $\bm{K}$ ($C_{3}$),

[TABLE]

Therefore, the invariants allowed by the group of $\bm{K}$ are real and imaginary parts of $k_{+}(u_{xx}-u_{yy}+i2u_{xy})$, and the preliminary Hamiltonian,

[TABLE]

As before, the $m_{x}$ symmetry gives

[TABLE]

and the time reversal gives,

[TABLE]

The consistency of the last two equations enforces, $B=0$, and in summary,

[TABLE]

Since we are interested in the physical conditions with $m_{x}$ preserved, the strain we consider is such that $u_{xy}=0$, therefore we get,

[TABLE]

Let us note that the tensor $B_{i}B_{j}-(1/2)\bm{B}^{2}$ transforms in the same way as the components of the traceless strain tensor considered above. Therefore, on par with Eq. (109) we also have the contribution of the second order in magnetic field,

[TABLE]

which has a meaning of the two-fold modulation of the Dirac cone tilt as the in-plane magnetic field rotates. Such two fold variation is normally associated with the planar Hall effect [66, 67, 68]. The field induced dependent tilt of the Dirac cones has been previously identified as the source of the linear planar Hall effect on a surface of the topological insulator [69]. In these works, however the planar Hall has been induced by the terms in the Hamiltonian of the form, $B_{x}k_{y}-B_{y}k_{x}$. In our model such terms are only allowed if the system is placed on a substrate, while the quadratic terms of Eq. (111) exists regardless of substrate and give the non-linear planar Hall effect.

For the off-diagonal part notice that

[TABLE]

acquires the phase under the $C_{3}$ rotation, $C_{3}\sigma_{+}=\langle C_{3}\Phi_{\bm{K}}^{c}|\mathcal{H}|C_{3}\Phi^{v}_{\bm{K}}\rangle=e^{i4\pi/3}\sigma_{+}$ and similarly, $C_{3}\sigma_{-}=e^{-i4\pi/3}\sigma_{-}$. The immediate consequence of these relations is that the Hamiltonian is diagonal strictly at $\bm{K}$. To fix the Hamiltonian to linear order in $\bm{k}$ note that the combinations, $k_{\pm}=k_{x}\pm ik_{y}$ transform as $C_{3}k_{\pm}=e^{\mp i2\pi/3}k_{\pm}$. In result, the most general invariant to this order allowed by the group of $\bm{K}$ takes the form,

[TABLE]

where $v$ is a complex constant. This, however is not a finite form of the Hamiltonian as it is further constrained by the $m_{x}$ and $\mathcal{T}$ symmetries, and this we explore next.

The $m_{x}$ symmetry tells us that all the Pauli matrices such as Eq. (112) are the same for the two valies, $\bm{k}=\bm{K}$ and $\bm{k}=-\bm{K}$. Away from these high symmetry points, based on Table 1 and (113), $m_{x}$ implies

[TABLE]

Similarly the time reversal implies,

[TABLE]

The consistency between (114) and (115) requires $v$ is a real number. In summary, we have

[TABLE]

Finally we derive the coupling to the in-plane magnetic field, $\bm{B}$. Such term is requires breaking of the basal mirror symmetry, $m_{z}$ which is achieved by placing the whole system on a substrate. As we do not consider the actual spin the coupling to the magnetic field we have in pir problem is via the orbital degrees of freedom in the form of the $\bm{B}\cdot\bm{L}$ coupling. To have $\bm{B}$ in the effective Hamiltonian a finite admixture of the $m_{z}$-odd orbitals to the $m_{z}$-even orbitals $\{d_{z^{2}},d_{x^{2}-y^{2}},d_{xy}\}$. This can only happen if $m_{z}$ symmetry is broken e.g. by a substrate. The combination $B_{x}\pm iB_{y}$ transforms as $k_{x}\pm ik_{y}$. Therefore, the allowed invariants are real and imaginary parts of $B_{+}\sigma_{-}$. The preliminary allowed Hamiltonian reads $H^{B}_{\bm{K}}=\delta(B_{x}\sigma_{x}+B_{y}\sigma_{y})+\gamma(B_{x}\sigma_{y}-B_{y}\sigma_{x})$. The $m_{x}$ symmetry gives

[TABLE]

The $\mathcal{T}$ yields

[TABLE]

The consistency of these equations enforces $\delta=0$, and we have

[TABLE]

We now summarize the model Hamiltonian that contains Eqs. (110), (116), (119), and (104),

[TABLE]

We note that upon the $\pi/2$ rotation, $k_{x}\rightarrow-k_{y}$, $k_{y}\rightarrow k_{x}$, and redefining Eq. (120) becomes at $B=0$,

[TABLE]

This Hamiltonian becomes identical to the one formulated in [18] upon the redefinition $\alpha\rightarrow-\alpha$, at $v_{x}=v_{y}=v$. This explains clearly why in our setting the mirror symmetry preserved by strain is $m_{x}$ while in the setting describing (121) the preserved mirror is $m_{y}$. For this reason our BCD is along $\hat{x}$ and in the (121) is along $\hat{y}$. Correspondingly, in our setting the Hall current from BCD is along $\hat{y}$ for $\bm{E}\parallel\hat{x}$ and in the (121) setting along $\hat{x}$ for $\bm{E}\parallel\hat{y}$.