Valley-Dependent Magnetoresistance in Two-Dimensional Semiconductors
VValley-Dependent Magnetoresistance in Two-Dimensional Semiconductors
Akihiko Sekine ∗ and Allan H. MacDonald Department of Physics, The University of Texas at Austin, Austin, Texas 78712, USA (Dated: October 17, 2018)We show theoretically that two-dimensional direct-gap semiconductors with a valley degree of freedom,including monolayer transition-metal dichalcogenides and gapped bilayer graphene, have a longitudinalmagnetoconductivity contribution that is odd in valley and odd in the magnetic field applied perpendicular tothe system. Using a quantum kinetic theory we show how this valley-dependent magnetoconductivity arisesfrom the interplay between the momentum-space Berry curvature of Bloch electrons, the presence of a magneticfield, and disorder scattering. We discuss how the e ff ect can be measured experimentally and used as a detectorof valley polarization. Introduction.—
Studies of magnetotransport in metals havea long standing in condensed matter physics. From the view-point of technology the discoveries of giant magnetoresis-tance [1, 2] and tunnel magnetoresistance [3–5] have led todrastic improvements in the performance of magnetic infor-mation storage devices. More generally magnetoresistancestudies can play an important role in characterizing the elec-tronic structure of solids. For example, Shubnikov–de Haasresistance oscillations are routinely used to measure Fermisurfaces. More recently the existence of three-dimensional(3D) Dirac and Weyl semimetals, which have topologicallynontrivial band structures, has been confirmed experimentally[6–11] by measuring a remarkable and characteristic negativelongitudinal magnetoresistance property associated with thechiral anomaly [12–15].This Rapid Communication addresses magnetotransport in2D semiconductors with more than one valley. Valley has re-cently attracted greater attention as an observable degree offreedom of electrons in solids [16–18], in part because ofthe emergence of monolayer transition-metal dichalcogenides(TMDs) and gapped bilayer graphene, both 2D semiconduc-tors in which valence and conduction band extrema occurat the K and K (cid:48) time-reversal partner Brillouin-zone cornerpoints. When intervalley scattering by disorder or phononsis weak, valley remains an approximate quantum numbereven beyond the Bloch band approximation. Weak valleyrelaxation combined with valley-dependent contributions tothe conductivity tensor can lead to observable e ff ects analo-gous to those produced by spin accumulations in conductorswith weak spin-orbit scattering. To date, attention has fo-cused mainly on the valley-dependent anomalous Hall e ff ect[19, 20], which occurs in the absence of a magnetic field and isrelated to the broken time-reversal symmetry of the Hamilto-nian’s projection to a single valley, and to momentum-spaceBerry phase e ff ects. Given the negative magnetoresistancein 3D Dirac and Weyl semimetals, which also involves val-leys related by time-reversal, longitudinal magnetotransporte ff ects should be expected in 2D multi-valley systems. Weapproach this issue theoretically using a massive Dirac modelfor 2D multi-valley semiconductors and a recently developed ∗ Present address: Center for Emergent Matter Science, RIKEN, Wako,Saitama 351-0198, Japan. [email protected] quantum kinetic theory [15, 21]. We find that the longitudinalmagnetoconductivity has a contribution that is odd in valleyand odd in perpendicular magnetic field. Our theoretical pre-dictions can be tested by observing a change from quadraticto linear magnetoresistance in systems in which a finite valleypolarization is induced by optical pumping or valley injection,as schematically illustrated in Fig. 1.
Magnetotransport theory.—
The transport theory we em-ploy is valid in the low magnetic field regime where Landauquantization can be neglected and enables us to systematicallycompute the conductivity tensor in the presence of disorderin arbitrary spatial dimensions. It is based on a quantum ki-netic equation that accounts for disorder, and for electric E and magnetic B fields [15, 22]: ∂ (cid:104) ρ (cid:105) ∂ t + i (cid:126) [ H , (cid:104) ρ (cid:105) ] + K ( (cid:104) ρ (cid:105) ) = D E ( (cid:104) ρ (cid:105) ) + D B ( (cid:104) ρ (cid:105) ) , (1)where (cid:104) ρ (cid:105) is the impurity-averaged Bloch-electron densitymatrix, H is the unperturbed Bloch Hamiltonian, K ( (cid:104) ρ (cid:105) ) isa disorder contribution discussed further below, and D E ( (cid:104) ρ (cid:105) )and D B ( (cid:104) ρ (cid:105) ) are the electric and magnetic driving terms: D E ( (cid:104) ρ (cid:105) ) = e E (cid:126) · D (cid:104) ρ (cid:105) D k , D B ( (cid:104) ρ (cid:105) ) = e (cid:126) (cid:40)(cid:32) D H D k × B (cid:33) · D (cid:104) ρ (cid:105) D k (cid:41) . (2) FIG. 1. (a) Schematic illustration of valley polarization due to achemical potential di ff erence δµ between two valleys in systems withweak intervalley scattering. (b) Schematic of the magnetic-field de-pendence of the low-magnetic-field magnetoresistance. The mag-netoresistance is predicted to have a linear dependence on B z when δµ (cid:44)
0, and a quadratic dependence on B z when δµ = a r X i v : . [ c ond - m a t . m e s - h a ll ] M a y Here, e > k is the crystal wave-vector, { a · b } = a · b + b · a (with a and b being vectors) is a symmetrized operator prod-uct, and we have introduced a covariant derivative notationfor the wave-vector dependence of the matrices X ( = (cid:104) ρ (cid:105) , H )expressed in an eigenstate representation: DXD k = ∂ X ∂ k − i [ R k , X ] , (3)where R k = (cid:80) α = x , y , z R k ,α e α , and [ R k ,α ] mn = i (cid:104) u m k | ∂ k α u n k (cid:105) is ageneralized Berry connection of Bloch electrons.The steady-state linear response of the density matrix toan electric field can be expressed as a formal expansion inpowers of the magnetic field strength B [15, 23]: (cid:104) ρ (cid:105) = (1 −L − D B ) − L − D E ( (cid:104) ρ (cid:105) + (cid:104) Ξ B (cid:105) ) ≡ (cid:104) ρ E (cid:105) + (cid:80) N ≥ (cid:104) ρ B , N (cid:105) , where (cid:104) ρ E (cid:105) = L − D E ( (cid:104) ρ (cid:105) ), (cid:104) ρ B , N (cid:105) = ( L − D B ) N L − D E ( (cid:104) ρ (cid:105) ) + ( L − D B ) N − L − D E ( (cid:104) Ξ B (cid:105) ), and we have defined the Liouvil-lian operator L ≡ P + K with P (cid:104) ρ (cid:105) ≡ ( i / (cid:126) ) [ H , (cid:104) ρ (cid:105) ]. In thisexpansion (cid:104) ρ (cid:105) is the Fermi-Dirac equilibrium density matrixin the absence of both fields, and (cid:104) Ξ B (cid:105) is the equilibrium den-sity matrix in the absence of an electric field. (cid:104) Ξ B (cid:105) accountsfor the Berry phase correction to the density of states impliedby semiclassical wave-packet dynamics [24].Throughout this Rapid Communication, we work in theeigenstate basis for the various contributions to the steady-state density matrix, and decompose (cid:104) ρ B , N (cid:105) into its band-diagonal part (cid:104) n B , N (cid:105) + (cid:104) ξ B , N (cid:105) , and its band-o ff -diagonal part (cid:104) S B , N (cid:105) . We adopt a relaxation time approximation for the dis-order scattering that influences the diagonal part of (cid:104) ρ B , N (cid:105) : (cid:104) n B , N (cid:105) mm k = τ m [ D B ( (cid:104) ρ B , N − (cid:105) )] mm k , (cid:104) ξ B , N (cid:105) mm k = e (cid:126) B · Ω m k (cid:104) n B , N − (cid:105) mm k , (4)where N ≥ (cid:104) ρ B , (cid:105) = (cid:104) ρ E (cid:105) , (cid:104) n B , (cid:105) = (cid:104) n E (cid:105) , and τ m and Ω m k are respectively the scattering time and the Berry curvaturevector for band m . We also have (cid:104) Ξ B (cid:105) mm k = ( e / (cid:126) ) B · Ω m k (cid:104) ρ (cid:105) mm k .In Eq. (4), (cid:104) n B , N (cid:105) is the extrinsic (Lorentz force) contribution,while (cid:104) ξ B , N (cid:105) is the intrinsic (Berry phase) contribution. Theband o ff -diagonal part is given by [15] (cid:104) S B , N (cid:105) mm (cid:48) k = (cid:126) i [ D B ( (cid:104) ρ B , N − (cid:105) )] mm (cid:48) k − [ J ( (cid:104) n B , N (cid:105) )] mm (cid:48) k ε m k − ε m (cid:48) k , (5)where m (cid:44) m (cid:48) and ε m k is the energy eigenvalue of band m . InEq. (5) the term proportional to D B ( (cid:104) ρ B , N − (cid:105) ) is purely a band-structure property expressed in terms of the Berry connec-tion, whereas the term proportional to J ( (cid:104) n B , N (cid:105) ) is a disorder-dependent Fermi-surface response corresponding to a vertexcorrection in the ladder-diagram approximation [15, 21]. Theexplicit form of J ( (cid:104) n B , N (cid:105) ) will be given later. Massive Dirac model.—
We consider 2D semiconductorswith broken inversion symmetry, like monolayer TMDs, thathave two low-energy valleys related by time-reversal. Thelow-energy e ff ective Hamiltonians in these systems normallyhave the massive Dirac form [16, 19, 25] H τ z ( k ) = v F ( τ z k x σ x + k y σ y ) + m σ z . (6)(As we discuss briefly later, gated bilayer graphene is an ex-ception.) In Eq. (6) τ z = ± v F is the Fermi velocity, 2 m is the band gap, and the Paulimatrices σ i act in the space of the retained conduction andvalence bands. The eigenvalues of the Hamiltonian (6) are ± ε k = ± (cid:113) v F ( k x + k y ) + m with the eigenfunctions | u ± k ( τ z ) (cid:105) .From Eq. (3), we see that the wavevector dependence of theeigenfunctions | u ± k ( τ z ) (cid:105) plays an important role in our trans-port theory. In the eigenstate basis the Berry connection vector[ R τ z k ,α ] mn = i (cid:104) u m k ( τ z ) | ∂ k α u n k ( τ z ) (cid:105) with m , n = ± has the explicitform R τ z k , x = k τ z sin θ − ˜ σ z m k ε k τ z sin θ − ˜ σ y v F m ε k cos θ − ˜ σ x v F ε k τ z sin θ, R τ z k , y = − k τ z cos θ + ˜ σ z m k ε k τ z cos θ − ˜ σ y v F m ε k sin θ + ˜ σ x v F ε k τ z cos θ, (7)where e ± i θ = ( k x ± ik y ) / k , k = (cid:113) k x + k y , and ˜ σ α is a Paulimatrix that acts in the eigenstate basis. Also, the Berrycurvature takes the form [ Ω τ z k , z ] ± = i (cid:104) ∂ k x u ± k ( τ z ) | ∂ k y u ± k ( τ z ) (cid:105) − i (cid:104) ∂ k y u ± k ( τ z ) | ∂ k x u ± k ( τ z ) (cid:105) = ∓ τ z v F m / (2 ε k ). Valley-dependent longitudinal magnetoconductivity.—
Weapply our magnetotransport theory to the 2D systems de-scribed by Eq. (6) with a static magnetic field B = (0 , , B z )applied perpendicular to the system. For the moment we ne-glect the vertex correction, i.e., the contribution proportionalto J in Eq. (5). Without loss of generality we may consider anelectron-doped case with positive chemical potential µ . Ourgoal is to compute the xx -component of the magnetoconduc-tivity tensor using σ ( N ) µν ( B z ) = Tr (cid:104) ( − e ) v µ (cid:104) ρ B , N (cid:105) (cid:105) / E ν . (8)In the eigenstate basis the velocity operator reads v x = v F ( ˜ σ z v F k ε k cos θ + ˜ σ y τ z sin θ − ˜ σ x m ε k cos θ ).We first evaluate the magnetoconductivity contributionsproportional to odd powers of B z for µ > m in the conductionband. The linear magnetoconductivity σ (1) xx ( B z ) is determinedby the density matrix (cid:104) ρ B , (cid:105) = L − D B ( (cid:104) ρ E (cid:105) ) + L − D E ( (cid:104) Ξ B (cid:105) ).We find that [26] σ (1) xx ( B z ) = τ z e B z E x (cid:90) [ d k ] v F m ε k ∂∂ k x + v F mk x ε k (cid:104) n E (cid:105) ++ k ≡ τ z e (cid:126) B z v F C ( µ, m ) τ tr , (9)where [ d k ] ≡ d k (2 π ) , (cid:104) n E (cid:105) ++ k = τ tr ( eE x / (cid:126) ) ∂ f ( ε k ) /∂ k x , f ( ε k )is the Fermi-Dirac distribution function, τ tr is the intravalleyscattering time, and C ( µ, m ) < (cid:104) n E (cid:105) arise respectivelyfrom the o ff -diagonal intrinsic contribution (cid:104) S B , (cid:105) and the di-agonal intrinsic contribution (cid:104) ξ B , (cid:105) . In Fig. 2 we show the µ and m dependences of σ (1) xx ( B z ). It should be noted here that, FIG. 2. (a) Chemical potential µ dependence of σ (1) xx ( B z ) [Eq. (9)]and σ (1)ver xx ( B z ) [Eq. (14)] for m = . σ (1) xx ( B z ) and σ (1)ver xx ( B z )are proportional to m /µ in the case of µ (cid:29) m and thus approachzero in the limit µ (cid:29) m . (b) Band gap m dependence of σ (1) xx ( B z ) for µ/ m = .
1. In both (a) and (b), we set B z = . v F = · Å, τ tr = . T = as seen from Fig. 2(b), σ (1) xx ( B z ) is significantly enhanced in asmaller gap system. Similarly, we find the cubic magnetocon-ductivity obtained from (cid:104) ρ B , (cid:105) [26] σ (3) xx ( B z ) = τ z e B z E x (cid:90) [ d k ] v F m ε k ∂∂ k x + v F mk x ε k (cid:104) n B , (cid:105) ++ k ≡ τ z e (cid:126) B z v F C ( µ, m ) τ , (10)where (cid:104) n B , (cid:105) ++ k = ( eB z τ tr ) ( ∂ε k ∂ k y ∂∂ k x − ∂ε k ∂ k x ∂∂ k y ) (cid:104) n E (cid:105) ++ k , and C ( µ, m ) >
0. In the case of µ (cid:29) m , we find that C ( µ, m ) ∝ m /µ . For the material parameters used inFig. 2(a), | σ (3) xx ( B z ) /σ (1) xx ( B z ) | ∼ − ( B z [T]) . The two termsin square brackets of Eq. (10) acting on the extrinsic response (cid:104) n B , (cid:105) arise respectively from the o ff -diagonal intrinsic con-tribution (cid:104) S B , (cid:105) and the diagonal intrinsic contribution (cid:104) ξ B , (cid:105) .There are no valley-independent contributions to the linearand cubic magnetoconductivities, as required by time-reversalsymmetry. Higher-order odd-power terms have the generalform σ ( N ) xx ( B z ) = τ z e N + (cid:126) B Nz v NF C N ( µ, m ) τ N tr , (11)where N = , , · · · is an odd integer and C N ( µ, m ) has thedimension of [Energy] − N .Next, we consider the magnetoconductivity contributionsproportional to even powers of B z , which cannot be valley de-pendent due to time-reversal symmetry [27]. We find that thequadratic magnetoconductivity obtained from (cid:104) ρ B , (cid:105) is domi-nated by the Lorentz-force contribution σ (2) xx ( B z ) ≈ − e B z τ E x (cid:90) [ d k ] v F k x ε k (cid:18) ∂ε k ∂ k y ∂∂ k x − ∂ε k ∂ k x ∂∂ k y (cid:19) (cid:104) n E (cid:105) ++ k = − σ (0) xx ( ω c τ tr ) , (12)where σ (0) xx = ( − e / E x ) (cid:82) [ d k ]( v F k x /ε k ) (cid:104) n E (cid:105) ++ k is the Drudeconductivity and ω c = eB z v F /µ is the cyclotron frequency.Intrinsic contributions to the quadratic magnetoconductivity are not zero, but they are suppressed by ∼ / ( µτ tr ) (cid:28) Vertex corrections.—
From Eq. (5) the vertex correctioncontribution to the density matrix linear in B z is given by (cid:104) S (cid:48) B , (cid:105) mm (cid:48)(cid:48) k ≡ i (cid:126) [ J ( (cid:104) n B , (cid:105) )] mm (cid:48)(cid:48) k / ( ε m k − ε m (cid:48)(cid:48) k ), where [21][ J ( (cid:104) n (cid:105) )] mm (cid:48)(cid:48) k = π (cid:126) (cid:88) m (cid:48) k (cid:48) (cid:104) U mm (cid:48) kk (cid:48) U m (cid:48) m (cid:48)(cid:48) k (cid:48) k (cid:105) (cid:104) ( n m k − n m (cid:48) k (cid:48) ) δ ( ε m k − ε m (cid:48) k (cid:48) ) + ( n m (cid:48)(cid:48) k − n m (cid:48) k (cid:48) ) δ ( ε m (cid:48)(cid:48) k − ε m (cid:48) k (cid:48) ) (cid:105) . (13)Here, m (cid:44) m (cid:48)(cid:48) and (cid:104) n (cid:105) = diag[ n m k ] is an arbitrary band-diagonal density matrix. We assume short-range disorderof the form U ( r ) = U (cid:80) i δ ( r − r i ) with (cid:104) U ( r ) U ( r (cid:48) ) (cid:105) = n imp U δ ( r − r (cid:48) ), where n imp is the impurity density. Aftera lengthy calculation [26], we find the vertex correction to thelinear magnetoconductivity σ (1)ver xx ( B z ) = Tr (cid:104) ( − e ) v x (cid:104) S (cid:48) B , (cid:105) (cid:105) / E x ≡ τ z e (cid:126) B z v F C ver1 ( µ, m ) τ tr , (14)where C ver1 ( µ, m ) <
0. This means that the vertex correc-tion enhances the valley-dependent linear magnetoconductiv-ity [see Fig. 2(a)]. This contrasts with its well-known influ-ence on the the spin Hall [28] and anomalous Hall [29] con-ductivities in certain Rashba models, i.e., the suppression ofthese conductivities by the vertex correction. Here, we notethat usual golden-rule intraband scattering rates proportionalto (cid:104) U ++ kk (cid:48) U ++ k (cid:48) k (cid:105) are not valley dependent. Nonzero contribu-tions to Eq. (14) require interband scattering matrix elementslike (cid:104) U ++ kk (cid:48) U + − k (cid:48) k (cid:105) , which are valley dependent. The vertex cor-rection is valley dependent because it is due to interband co-herence induced by the magnetic field. Discussion.—
A longitudinal total magnetoconductivityproportional to odd powers of magnetic field can occur onlyin systems with broken time-reversal symmetry [30, 31]. Inmonolayer TMDs and gated bilayer graphene, spatial in-version symmetry is broken but time-reversal symmetry isretained. The valley-dependent magnetoconductivity con-tributes to transport only in the presence of a finite valley po-larization, for example one due to a chemical potential dif-ference between the two valleys, that explicitly breaks time-reversal symmetry. Valley polarization in TMDs can be re-alized by applying circularly polarized light [16, 32–34] togenerate an excess population of carriers in one valley. Whenintravalley scattering is much stronger than intervalley scatter-ing, equilibration will occur within valleys to establish valley-dependent chemical potentials. This approach has been usedpreviously to measure the valley Hall e ff ect [20], and is themost direct way to measure the valley-dependent magneto-conductivity derived in this Rapid Communication. In dis-cussing the results of such a measurement below, we assumethat the contribution to transport from photo-generated holesis negligible.Including terms up to order of B z and allowing for a chem-ical potential di ff erence between valleys, the total magneto- FIG. 3. (a) Valley polarization δµ dependence of the magnetore-sistance ∆ ρ xx ( B z , δµ ) [Eq. (16)] in the case of | σ (0) xx /σ (0) xy | (cid:29) m = . µ = µ = .
82 eV for the δµ = µ = .
83 eVand µ = .
82 eV for the δµ (cid:44) ∆ ρ xx ( B z , δµ ) in the caseof | σ (0) xx /σ (0) xy | (cid:28) m =
30 meV, corresponding to low-carrier-density gated bilayer graphene. We used µ = µ =
33 meV for the δµ = µ =
36 meV and µ =
33 meV for the δµ (cid:44) ff ect is much stronger in a smallergap system. In both (a) and (b), we set v F = · Å, τ tr = . T = conductivity of a two-valley system reads σ Bxx ( µ , µ ) = e (cid:126) B z v F [ C tot1 ( µ , m ) − C tot1 ( µ , m )] τ tr − σ (0) xx ( µ )( ω c τ tr ) − σ (0) xx ( µ )( ω c τ tr ) , (15)where µ i ( i = ,
2) is the chemical potential of valley i , C tot1 ( µ i , m ) = C ( µ i , m ) + C ver1 ( µ i , m ), and ω ci = eB z v F /µ i .In the low-field limit σ Bxx ∝ B z when δµ = µ − µ (cid:44) σ Bxx ∝ B z when δµ =
0. The resistivity is defined by ρ xx ( B z ) = σ xx / ( σ xx + σ xy ) with σ µν ≡ σ (0) µν + σ B µν . It followsthat the low-field magnetoresistance ∆ ρ xx ≡ ρ xx ( B z ) − ρ xx (0) ρ xx (0) ≈ ∓ σ Bxx ( µ , µ ) σ (0) xx ( µ ) + σ (0) xx ( µ ) . (16)Here, the − ( + ) sign applies in the | σ (0) xx /σ (0) xy | (cid:29) | σ (0) xx /σ (0) xy | (cid:28) σ (0) xx is not valley dependent. Thus the change in magnetic-field dependence from B z to B z when illuminated by circularlypolarized light, illustrated schematically in Fig. 1, should be readily observable [36]. Interestingly, we can always makethe δµ (cid:44) δµ = ff ects discussed in this Rapid Com-munication are much stronger, for a given Fermi velocity,in Dirac models with a smaller gap, and we expect themto be much more easily observed experimentally in bilayergraphene systems than in monolayer TMDs. Bilayer grapheneis described by a generalized Dirac model with chirality J = J =
1, and has quadratic dispersion in the absenceof a gap [37]. For a given gap the size of the magnetoresis-tance e ff ect in bilayer graphene will exceed the J = µ isonly slightly larger than the gap m ), the current partitioningbetween valleys corresponding to equal electro-chemical po-tential gradients changes across interfaces at which the carrierdensity changes. If intervalley scattering is weak, valley accu-mulation will persist within a valley relaxation length of anysuch interface, and should be detectable via Kerr microscopy[38]. Summary.—
To summarize, we have theoretically demon-strated the existence of a valley-dependent longitudinal mag-netoconductivity in 2D semiconductors with a valley de-gree of freedom. The e ff ect arises from the interplay be-tween the momentum-space Berry curvature of Bloch elec-trons, the presence of a magnetic field, and disorder scatter-ing. Our prediction can be verified by measuring the influ-ence of circularly-polarized light illumination on magnetore-sistance in the low-field limit. Because the magnetoresistanceis proportional to valley polarization it can be used as a valleydetector. We predict that these magnetoresistance e ff ects willbe significantly enhanced in bilayer graphene samples withsmall gaps and small carrier densities.We thank H. Chen and D. Culcer for valuable discus-sions. This work was supported by the Department of En-ergy, O ffi ce of Basic Energy Sciences under Contract No. DE-FG02-ER45958 and by the Welch foundation under Grant No.TBF1473. A.S. is supported by the JSPS Overseas ResearchFellowship. [1] M. N. Baibich, J. M. Broto, A. Fert, F. Nguyen Van Dau, F.Petro ff , P. Etienne, G. Creuzet, A. Friederich, and J. Chazelas,Phys. Rev. Lett. , 2472 (1988).[2] G. Binasch, P. Gr¨unberg, F. Saurenbach, and W. Zinn, Phys.Rev. B , 4828(R) (1989).[3] M. Julliere, Phys. Lett. , 225 (1975).[4] T. Miyazaki and N. Tezuka, J. Magn. Magn. Mater. , L231(1995).[5] J. S. Moodera, L. R. Kinder, T. M. Wong, and R. Meservey, Phys. Rev. Lett. , 3273 (1995).[6] J. Xiong, S. K. Kushwaha, T. Liang, J. W. Krizan, M.Hirschberger, W. Wang, R. J. Cava, and N. P. Ong, Science ,413 (2015).[7] C.-Z. Li, L.-X. Wang, H. Liu, J. Wang, Z.-M. Liao, and D.-P.Yu, Nat. Commun. , 10137 (2015).[8] X. Huang, L. Zhao, Y. Long, P. Wang, D. Chen, Z. Yang, H.Liang, M. Xue, H. Weng, Z. Fang, X. Dai, and G. Chen, Phys.Rev. X , 031023 (2015). [9] H. Li, H. He, H.-Z. Lu, H. Zhang, H. Liu, R. Ma, Z. Fan, S.-Q.Shen, and J. Wang, Nat. Commun. , 10301 (2016).[10] Q. Li, D. E. Kharzeev, C. Zhang, Y. Huang, I. Pletikosi, A. V.Fedorov, R. D. Zhong, J. A. Schneeloch, G. D. Gu, and T. Valla,Nat. Phys. , 550 (2016).[11] F. Arnold, C. Shekhar, S.-C. Wu, Y. Sun, R. D. dos Reis, N. Ku-mar, M. Naumann, M. O. Ajeesh, M. Schmidt, A. G. Grushin,J. H. Bardarson, M. Baenitz, D. Sokolov, H. Borrmann, M.Nicklas, C. Felser, E. Hassinger, and B. Yan, Nat. Commun. , 11615 (2016).[12] D. T. Son and B. Z. Spivak, Phys. Rev. B , 104412 (2013).[13] A. A. Burkov, Phys. Rev. Lett. , 247203 (2014).[14] B. Z. Spivak and A. V. Andreev, Phys. Rev. B , 085107(2016).[15] A. Sekine, D. Culcer, and A. H. MacDonald, Phys. Rev. B ,235134 (2017).[16] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Phys. Rev.Lett. , 196802 (2012).[17] X. Xu, W. Yao, D. Xiao, and T. F. Heinz, Nat. Phys. , 343(2014).[18] K. F. Mak and J. Shan, Nat. Photonics , 216 (2016).[19] D. Xiao, W. Yao, and Q. Niu, Phys. Rev. Lett. , 236809(2007).[20] K. F. Mak, K. L. McGill, J. Park, and P. L. McEuen, Science , 1489 (2014).[21] D. Culcer, A. Sekine, and A. H. MacDonald, Phys. Rev. B ,035106 (2017).[22] The quantum kinetic equation (1) reduces to the usual semiclas-sical Boltzmann equation in spin-independent single-band sys-tems, since the covariant derivatives reduce to simple deriva-tives in such systems. In this work we treat disorder in theBorn approximation, although we can in principle go beyondthis in our formalism. In the approximation we employ, the re-sults obtained by our quantum kinetic formalism are equivalentto those obtained by the Kubo formalism with vertex correc-tions included in the ladder-diagram approximation and bothband diagonal and band o ff -diagonal terms taken into account.If only band diagonal terms were retained, our approach would be equivalent to the standard semiclassical Boltzmann theory.One of the advantages of our formalism is that Berry phase anddisorder contributions can be incorporated in a unified fashionat each order of electric and magnetic field strengths.[23] See Supplemental Material for the detailed derivation.[24] D. Xiao, J. Shi, and Q. Niu, Phys. Rev. Lett. , 137204 (2005).[25] C.-C. Liu, W. Feng, and Y. Yao, Phys. Rev. Lett. , 076802(2011).[26] See Supplemental Material for a detailed calculation of σ (1) xx ( B z ), σ (2) xx ( B z ), and σ (3) xx ( B z ).[27] Under time reversal valleys change place and the field directionis reversed. The response must be odd in time and odd in fieldor even in time and even in field.[28] J. Inoue, G. E.W. Bauer, and L.W. Molenkamp, Phys. Rev. B , 041303 (2004).[29] J. Inoue, T. Kato, Y. Ishikawa, H. Itoh, G. E. W. Bauer, and L.W. Molenkamp, Phys. Rev. Lett. , 046604 (2006).[30] L. Onsager, Phys. Rev. , 405 (1931).[31] H. Chen, Y. Gao, D. Xiao, A. H. MacDonald, and Q. Niu,arXiv:1511.02557.[32] T. Cao, G. Wang, W. Han, H. Ye, C. Zhu, J. Shi, Q. Niu, P. Tan,E. Wang, B. Liu, and J. Feng, Nat. Commun. , 887 (2012).[33] H. Zeng, J. Dai, W. Yao, D. Xiao, and X. Cui, Nat. Nanotechnol. , 490 (2012).[34] K. F. Mak, K. He, J. Shan, and T. F. Heinz, Nat. Nanotechnol. , 494 (2012).[35] Note that the anomalous Hall conductivity σ (0) xy takes the max-imum value τ z e / h when the system is insulating, i.e., when µ i < m in our model. See Supplemental Material for a detailedderivation of ∆ ρ xx .[36] Note that Fig. 1(b) is an enlarged view of the absolute valuesof the magnetoresistance ∆ ρ xx at the magnetic field scale B z ∼ .
01 T.[37] H. Min and A. H. MacDonald, Phys. Rev. B , 155416 (2008).[38] J. Lee, K. F. Mak, and J. Shan, Nat. Nanotechnol. , 421(2016). Supplemental Material
I. STEADY-STATE LINEAR RESPONSE OF THE DENSITY MATRIX TO AN ELECTRIC FIELD
Let us consider a density matrix in the presence of electric and magnetic fields. We write the electron density matrix as (cid:104) ρ (cid:105) = (cid:104) ρ (cid:105) + (cid:104) ρ (cid:105) F , where (cid:104) ρ (cid:105) is the Fermi-Dirac equilibrium density matrix in the absence of fields, and (cid:104) ρ (cid:105) F is the field-induceddensity matrix. Then we can rewrite the steady-state uniform system limit of Eq. (1) in the form( L − D E − D B ) (cid:104) ρ (cid:105) F = ( D E + D B ) (cid:104) ρ (cid:105) , (S1)where we have defined the Liouvillian operator L ≡ P + K with P (cid:104) ρ (cid:105) ≡ ( i / (cid:126) ) [ H , (cid:104) ρ (cid:105) ] and have used the fact that L(cid:104) ρ (cid:105) =
0. Itfollows that (cid:104) ρ (cid:105) F = (cid:2) − L − ( D E + D B ) (cid:3) − L − ( D E + D B ) (cid:104) ρ (cid:105) = (cid:88) N ≥ (cid:2) L − ( D E + D B ) (cid:3) N L − ( D E + D B ) (cid:104) ρ (cid:105) . (S2)Equation (S2) can be used to derive a low-magnetic-field expansion for the linear response of the density matrix, and hence ofany single-particle observable, to an electric field E . From Eq. (S2) we have (cid:104) ρ (cid:105) F = (cid:88) N , N (cid:48) ≥ ( L − D B ) N L − D E ( L − D B ) N (cid:48) (cid:104) ρ (cid:105) . (S3)Here, the N = N (cid:48) = E (which is linear in E ): (cid:104) ρ E (cid:105) = L − D E ( (cid:104) ρ (cid:105) ).It is convenient to define a density matrix induced solely by the magnetic field B as (cid:104) ρ B (cid:105) ≡ (cid:80) N ≥ ( L − D B ) N (cid:104) ρ (cid:105) . Note that thisexpression is the generalized solution for ( L − D B ) (cid:104) ρ B (cid:105) = D B ( (cid:104) ρ (cid:105) ). Since we are assuming that the magnetic field is very weak,we retain only the term linear in B and neglect the higher-order terms, i.e., we may set (cid:104) ρ B (cid:105) = (cid:104) Ξ B (cid:105) , where (cid:104) Ξ B (cid:105) accounts for theBerry phase correction to the density of states implied by semiclassical wave-packet dynamics. This means that the correctionto the Fermi-Dirac distribution function due to B is given by the Berry phase correction (cid:104) Ξ B (cid:105) . Finally, we obtain (cid:104) ρ (cid:105) F = (cid:104) ρ E (cid:105) + (cid:88) N ≥ (cid:104) ρ B , N (cid:105) , (S4)where (cid:104) ρ B , N (cid:105) = ( L − D B ) N L − D E ( (cid:104) ρ (cid:105) ) + ( L − D B ) N − L − D E ( (cid:104) Ξ B (cid:105) ). II. GENERAL EXPRESSION FOR THE MAGNETIC DRIVING TERM
Let us consider a generic two-band system with the energy eigenvalues ε ± k = ± ε k and a magnetic field applied along the z direction such that B = (0 , , B z ). In what follows we work in the eigenstate basis where matrices are written in the basis of (cid:34) ++ + −− + −− (cid:35) . The magnetic driving term obtained from a diagonal density matrix (cid:104) n (cid:105) = (cid:34) n + k n − k (cid:35) is given by D B ( (cid:104) n (cid:105) ) = eB z (cid:34)(cid:40) D H Dk y , D (cid:104) n (cid:105) Dk x (cid:41) − (cid:40) D H Dk x , D (cid:104) n (cid:105) Dk y (cid:41)(cid:35) = eB z ∂ε k ∂ k y ∂ n + k ∂ k x − ∂ε k ∂ k x ∂ n + k ∂ k y − i ε k ( R + − x ∂∂ k y − R + − y ∂∂ k x )( n + k + n − k ) i ε k ( R − + x ∂∂ k y − R − + y ∂∂ k x )( n + k + n − k ) − (cid:18) ∂ε k ∂ k y ∂ n − k ∂ k x − ∂ε k ∂ k x ∂ n − k ∂ k y (cid:19) . (S5)On the other hand, the magnetic driving term obtained from an o ff -diagonal density matrix (cid:104) S (cid:105) = (cid:34) a k b k (cid:35) is given by D B ( (cid:104) S (cid:105) ) = eB z (cid:34)(cid:40) D H Dk y , D (cid:104) S (cid:105) Dk x (cid:41) − (cid:40) D H Dk x , D (cid:104) S (cid:105) Dk y (cid:41)(cid:35) = eB z (cid:34) A k − B k A k − B k (cid:35) (S6)with A k = − i ∂ε k ∂ k y ( R + − x b k − R − + x a k ) + i R + − y ε k (cid:34) i ( R ++ x − R −− x ) b k + ∂ b k ∂ k x (cid:35) − i R − + y ε k (cid:34) − i ( R ++ x − R −− x ) a k + ∂ a k ∂ k x (cid:35) , B k = − i ∂ε k ∂ k x ( R + − y b k − R − + y a k ) + i R + − x ε k (cid:34) i ( R ++ y − R −− y ) b k + ∂ b k ∂ k y (cid:35) − i R − + x ε k (cid:34) − i ( R ++ y − R −− y ) a k + ∂ a k ∂ k y (cid:35) . (S7)Especially in the case of (cid:104) S (cid:105) = c k σ y (i.e., a k = − ic k and b k = ic k ), we have A k = ( R + − x + R − + x ) ∂ε k ∂ k y c k − i ( R + − y − R − + y )( R ++ x − R −− x ) ε k c k − ( R + − y + R − + y ) ε k ∂ c k ∂ k x , B k = ( R + − y + R − + y ) ∂ε k ∂ k x c k − i ( R + − x − R − + x )( R ++ y − R −− y ) ε k c k − ( R + − x + R − + x ) ε k ∂ c k ∂ k y . (S8)Similarly in the case of (cid:104) S (cid:105) = c k σ x (i.e., a k = c k and b k = c k ), we have A k = − i ( R + − x − R − + x ) ∂ε k ∂ k y c k − ( R + − y + R − + y )( R ++ x − R −− x ) ε k c k + i ( R + − y − R − + y ) ε k ∂ c k ∂ k x , B k = − i ( R + − y − R − + y ) ∂ε k ∂ k x c k − ( R + − x + R − + x )( R ++ y − R −− y ) ε k c k + i ( R + − x − R − + x ) ε k ∂ c k ∂ k y . (S9) III. THEORETICAL MODEL
Let us consider the two-component massive Dirac fermion model whose Hamiltonian is given by H τ z ( k ) = v F ( τ z k x σ x + k y σ y ) + m σ z (S10)with τ z = ± | u ± k ( τ z ) (cid:105) = √ (cid:113) ± m ε k ± τ z e i τ z θ (cid:113) ∓ m ε k , (S11)where ε ± k = ± ε k = ± (cid:113) v F ( k x + k y ) + m are the eigenvalues and e ± i θ = ( k x ± ik y ) / k with k = (cid:113) k x + k y . The Berry connection[ R τ z k ,α ] mn = i (cid:104) u m k ( τ z ) | ∂∂ k α u n k ( τ z ) (cid:105) with m , n = ± has the explicit form R τ z k , x = k τ z sin θ − ˜ σ z m k ε k τ z sin θ − ˜ σ y v F m ε k cos θ − ˜ σ x v F ε k τ z sin θ, R τ z k , y = − k τ z cos θ + ˜ σ z m k ε k τ z cos θ − ˜ σ y v F m ε k sin θ + ˜ σ x v F ε k τ z cos θ. (S12)Also, the Berry curvature takes the form [ Ω τ z k , z ] ± = i (cid:104) ∂ k x u ± k ( τ z ) | ∂ k y u ± k ( τ z ) (cid:105) − i (cid:104) ∂ k y u ± k ( τ z ) | ∂ k x u ± k ( τ z ) (cid:105) = ∓ τ z v F m / (2 ε k ). The velocityoperator in the eigenstate representation is obtained as v x ( τ z ) = (cid:104) u m k ( τ z ) | ∂ H τ z ∂ k x | u n k ( τ z ) (cid:105) = v F (cid:32) ˜ σ z v F k ε k cos θ + ˜ σ y τ z sin θ − ˜ σ x m ε k cos θ (cid:33) . (S13) IV. LINEAR MAGNETOCONDUCTIVITY
Let us consider a case where an electric field is applied along the x direction and a magnetic field is applied along the z direction such that E = ( E x , ,
0) and B = (0 , , B z ). The longitudinal linear magnetoconductivity is given by σ (1) xx ( B z ) = Tr (cid:104) ( − e ) v x (cid:110) L − D B ( (cid:104) ρ E (cid:105) ) + L − D E ( (cid:104) Ξ B (cid:105) ) (cid:105)(cid:111) / E x , (S14)where (cid:104) ρ E (cid:105) = L − D E ( (cid:104) ρ (cid:105) ) is the density matrix linear in the electric field, and (cid:104) Ξ B (cid:105) mm k = ( e / (cid:126) ) f ( ε m k ) B · Ω m k is the density matrixlinear in the magnetic field, which corresponds to the Berry phase correction to the density of states in semiclassical wave-packetdynamics. Here, (cid:104) ρ (cid:105) = diag[ f ( ε + k ) , f ( ε − k )] with µ > (cid:104) ρ E (cid:105) , which is given by (cid:104) n E (cid:105) = τ tr eE x ∂ f ( ε + k ) ∂ k x
00 0 , (S15)where τ tr is the intravalley scattering time. Next we compute the magnetic driving term obtained from (cid:104) n E (cid:105) . From Eq. (S5) wehave D B ( (cid:104) n E (cid:105) ) = eB z ∂ε k ∂ k y ∂ n + E k ∂ k x − ∂ε k ∂ k x ∂ n + E k ∂ k y − i ε k (cid:18) R + − x ∂ n + E k ∂ k y − R + − y ∂ n + E k ∂ k x (cid:19) i ε k (cid:18) R − + x ∂ n + E k ∂ k y − R − + y ∂ n + E k ∂ k x (cid:19) . (S16)Then the corresponding diagonal and o ff -diagonal parts of the density matrix (cid:104) ρ B (cid:105) ( = (cid:104) S B (cid:105) + (cid:104) n B (cid:105) + (cid:104) ξ B (cid:105) ) are obtained as (cid:104) S B (cid:105) = L − [ D B ( (cid:104) n E (cid:105) )] = − i (cid:126) [ D B ( (cid:104) n E (cid:105) )] mm (cid:48) k ε m k − ε m (cid:48) k = eB z (cid:40) ∂ n + E k ∂ k x (cid:34) R + − k , y R − + k , y (cid:35) − ∂ n + E k ∂ k y (cid:34) R + − k , x R − + k , x (cid:35)(cid:41) = eB z ˜ σ y − ∂ n + E k ∂ k x v F m ε k sin θ + ∂ n + E k ∂ k y v F m ε k cos θ + ˜ σ x (cid:32) ∂ n + E k ∂ k x τ z v F ε k cos θ + ∂ n + E k ∂ k y τ z v F ε k sin θ (cid:33) , (S17) FIG. S1. Schematic illustration of longitudinal quadratic magnetoconductivity contributions. (cid:104) n (cid:105) and (cid:104) ξ (cid:105) indicates band-diagonal densitymatrix components, and (cid:104) S (cid:105) indicates band-o ff -diagonal density matrix components (see the main text). Note that, for a generic two-bandmodel, the magnetic driving terms acting on purely o ff -diagonal density matrices are purely diagonal and proportional to the identity matrix[see Eq. (S6)]. and (cid:104) n B (cid:105) = L − [ D B ( (cid:104) n E (cid:105) )] = τ tr [ D B ( (cid:104) n E (cid:105) )] mm k = τ tr eB z ∂ε k ∂ k y ∂ n + E k ∂ k x − ∂ε k ∂ k x ∂ n + E k ∂ k y
00 0 . (S18)Note that (cid:104) S B (cid:105) is accompanied by the Berry curvature contribution term that is purely diagonal: (cid:104) ξ B (cid:105) ++ k = eB z Ω + k , z n + E k = − τ z eB z v F m ε k n + E k , (cid:104) ξ B (cid:105) −− k = . (S19)The evaluation of L − D E ( (cid:104) Ξ B (cid:105) ) is quite similar to that of (cid:104) n E (cid:105) , which gives L − D E ( (cid:104) Ξ B (cid:105) ) = (cid:104) ξ B (cid:105) . Note that the o ff -diagonal partof L − D E ( (cid:104) Ξ B (cid:105) ) does not contribute to the linear magnetoconductivity.Finally, the total longitudinal linear magnetoconductivity is calculated to be σ (1) xx ( B z ) = Tr (cid:2) ( − e ) v x ( (cid:104) S B (cid:105) + (cid:104) ξ B (cid:105) ) (cid:3) / E x = τ z e B z E x (cid:90) d k (2 π ) v F m ε k ∂∂ k x + v F mk x ε k n + E k ≡ τ z e (cid:126) B z v F C ( µ, m ) τ tr , (S20)where C ( µ, m ) <
0. Note that the contribution from (cid:104) n B (cid:105) vanishes because it is an odd function of k y , and that the contributionfrom (cid:104) n (cid:48) B (cid:105) in Fig. S1 also vanishes because (cid:104) n (cid:48) B (cid:105) is proportional to the identity matrix. In the case of µ (cid:29) m , we find that C ( µ, m ) ∝ m /µ . Using this relation, we can check that the dimension of σ (1) xx ( B z ) is indeed the dimension of two-dimensionalelectrical conductivity: τ z e (cid:126) B z v F m µ τ tr = (A · s) · s V · sm m s
1J s = A s VJ = AV , (S21)where we have used J = V · A · s and (cid:126) = J · s. V. QUADRATIC MAGNETOCONDUCTIVITY
Let us calculate the quadratic magnetoconductivity, which is given by σ (2) xx ( B z ) = Tr (cid:110) ( − e ) v x (cid:104) ( L − D B ) L − D E ( (cid:104) ρ (cid:105) ) + L − D B L − D E ( (cid:104) Ξ B (cid:105) ) (cid:105)(cid:111) / E x . (S22)First we consider the contribution from the Lorentz force. This contribution is purely extrinsic, i.e., comes from the diagonalpart of D B . From Eq. (S5) the diagonal part of the density matrix obtained from (cid:104) n B (cid:105) reads (cid:104) n B (cid:105) = L − [ D B ( (cid:104) n B (cid:105) )] = τ tr eB z ∂ε k ∂ k y ∂ n + B k ∂ k x − ∂ε k ∂ k x ∂ n + B k ∂ k y
00 0 . (S23)Then we immediately get σ (2)LF xx ( B z ) = − e B z τ E x (cid:90) d k (2 π ) v F k x ε k (cid:18) ∂ε k ∂ k y ∂∂ k x − ∂ε k ∂ k x ∂∂ k y (cid:19) n + E k ≡ − σ (0) xx ( ω c τ tr ) , (S24)where σ (0) xx = − e / E x (cid:82) d k (2 π ) ( v F k x /ε k ) n + E k is the Drude conductivity and ω c = eB z v F /µ is the cyclotron frequency.There also exist intrinsic contributions to the quadratic magnetoconductivity (see Fig. S1). After a calculation we find that theintrinsic contributions take the form σ (2)in xx ( B z ) = e B z v F C ( µ, m ) τ tr , (S25)which is not dependent on τ z . We numerically find that C ( µ, m ) ∝ m /µ . Then we see that σ (2)in xx ( B z ) is quite small compared to σ (2)LF xx ( B z ): σ (2)in xx ( B z ) σ (2)LF xx ( B z ) ∼ e B z v F ( m /µ ) τ tr e B z v F (1 /µ ) τ ∼ µτ tr ) m µ (cid:28) , (S26)where we have used the fact that the condition ( µτ tr ) (cid:29) σ (2) xx ( B z ) ≈ σ (2)LF xx ( B z ) = − σ (0) xx ( ω c τ tr ) . (S27)Note that there are no τ z -dependent contributions to σ (2) xx ( B z ), as required by time-reversal symmetry. VI. CUBIC MAGNETOCONDUCTIVITY
Let us calculate the cubic magnetoconductivity, which is given by σ (3) xx ( B z ) = Tr (cid:110) ( − e ) v x (cid:104) ( L − D B ) L − D E ( (cid:104) ρ (cid:105) ) + ( L − D B ) L − D E ( (cid:104) Ξ B (cid:105) ) (cid:105)(cid:111) / E x . (S28)The contributions to σ (3) xx ( B z ) are obtained from the magnetic driving terms acting on the density matrices of the order of B z thatare shown in Fig. S1.The o ff -diagonal density matrix (cid:104) S B (cid:105) obtained from the o ff -diagonal part of D B ( (cid:104) n B (cid:105) ) is given by [see Eq. (S17) for a similarcalculation] (cid:104) S B (cid:105) = L − D B ( (cid:104) n B (cid:105) ) = eB z ˜ σ y − ∂ n + B k ∂ k x v F m ε k sin θ + ∂ n + B k ∂ k y v F m ε k cos θ + ˜ σ x (cid:32) ∂ n + B k ∂ k x τ z v F ε k cos θ + ∂ n + B k ∂ k y τ z v F ε k sin θ (cid:33) , (S29)which is accompanied by the Berry curvature contribution term (cid:104) ξ B (cid:105) ++ k = eB z Ω + z (cid:104) n B (cid:105) ++ k = − τ z eB z v F m ε k n + B k , (cid:104) ξ B (cid:105) −− k = , (S30)where n + B k is given by Eq. (S23). Then one of the contributions to the cubic magnetoconductivity is calculated to be σ (3) xx ( B z ) = Tr (cid:2) ( − e ) v x ( (cid:104) S B (cid:105) + (cid:104) ξ B (cid:105) ) (cid:3) / E x = τ z eB z (cid:90) ∞−∞ d k (2 π ) mv F ε k ∂∂ k x + v F mk x ε k n + B k ≡ τ z e (cid:126) B z v F D ( µ, m ) τ , (S31)0where D ( µ, m ) >
0. In the case of µ (cid:29) m , we numerically find that D ( µ, m ) ∝ m /µ .There is another contribution to the cubic magnetoconductivity. The diagonal density matrix (cid:104) n (cid:48) B (cid:105) obtained from (cid:104) S B (cid:105) = L − D B ( (cid:104) n E (cid:105) ) = c ˜ σ y + c ˜ σ x [Eq. (S17)] is calculated from Eq. (S6) to be (cid:104) n (cid:48) B (cid:105) = L − [ D B ( (cid:104) S B (cid:105) )] = τ tr eB z G k , (S32)where G k = − τ z v F k ε k c k − τ z v F m ε k k c k − τ z v F cos θ ∂ c k ∂ k x − τ z v F sin θ ∂ c k ∂ k y − v F m ε k sin θ ∂ c k ∂ k x + v F m ε k cos θ ∂ c k ∂ k y , (S33)with c k = eB z − ∂ n + E k ∂ k x v F m ε k sin θ + ∂ n + E k ∂ k y v F m ε k cos θ , (S34) c k = eB z (cid:32) ∂ n + E k ∂ k x τ z v F ε k cos θ + ∂ n + E k ∂ k y τ z v F ε k sin θ (cid:33) . (S35)Here, note that all the terms in G k are dependent on τ z . Then, the diagonal density matrix (cid:104) n (cid:48) B (cid:105) obtained from the diagonal partof D B ( (cid:104) n (cid:48) B (cid:105) ) is written as (cid:104) n (cid:48) B (cid:105) = L − [ D B ( (cid:104) n (cid:48) B (cid:105) )] = ( eB z ) τ tr (cid:32) ∂ε k ∂ k y ∂ G k ∂ k x − ∂ε k ∂ k x ∂ G k ∂ k y (cid:33) ˜ σ z , (S36)from which the contribution to the cubic magnetoconductivity is calculated to be σ (3) xx ( B z ) = Tr (cid:104) ( − e ) v x (cid:104) n (cid:48) B (cid:105) (cid:105) / E x = − e B z τ tr (cid:90) ∞−∞ d k (2 π ) v F k x ε k (cid:32) ∂ε k ∂ k y ∂ G k ∂ k x − ∂ε k ∂ k x ∂ G k ∂ k y (cid:33) = τ z e (cid:126) B z v F D ( µ, m ) τ . (S37)The other two contributions come from the diagonal parts of the density matrices at each power of B z acting on (cid:104) n E (cid:105) : oneof the three actions of D B is the Berry curvature contribution eB z Ω + z , and the other two actions of D B are the Lorentz forcecontribution τ tr eB z [ ∂ε k ∂ k y ∂∂ k x − ∂ε k ∂ k x ∂∂ k y ]. After a calculation we find that the other contribution to the cubic magnetoconductivity is σ (3) xx ( B z ) = τ z e (cid:126) B z v F D ( µ, m ) τ . (S38)Finally, from Eqs. (S31), (S37), and (S38) we obtain the total longitudinal cubic magnetoconductivity σ (3) xx ( B z ) = τ z e (cid:126) B z v F C ( µ, m ) τ , (S39)where C ( µ, m ) = (8 / D ( µ, m ) <
0. Note that there are no τ z -independent contributions to σ (3) xx ( B z ), as required by time-reversalsymmetry. VII. DISORDER EFFECT IN THE BORN APPROXIMATION
Let us consider a short-range (on-site) disorder of the form U ( r ) = U (cid:80) i δ ( r − r i ), and assume that the correlationfunction satisfies (cid:104) U ( r ) U ( r (cid:48) ) (cid:105) = n imp U δ ( r − r (cid:48) ) with n imp the impurity density. The eigenvectors of the Hamiltonian H τ z ( k ) = v F ( τ z k x σ x + k y σ y ) + m σ z with eigenvalues ε ± k = ± ε k = ± (cid:113) v F ( k x + k y ) + m are given by | u ± k ( τ z ) (cid:105) = √ (cid:113) ± m ε k ± τ z e i τ z θ (cid:113) ∓ m ε k , (S40)1where e ± i θ = ( k x ± ik y ) / k with k = (cid:113) k x + k y . From these eigenstates, we immediately get U ++ kk (cid:48) = U (cid:104) u + k ( τ z ) | u + k (cid:48) ( τ z ) (cid:105) = U (cid:114)(cid:18) + m ε (cid:19) (cid:18) + m (cid:48) ε (cid:48) (cid:19) + e i τ z γ (cid:114)(cid:18) − m ε (cid:19) (cid:18) − m (cid:48) ε (cid:48) (cid:19) , U + − kk (cid:48) = U (cid:104) u + k ( τ z ) | u − k (cid:48) ( τ z ) (cid:105) = U (cid:114)(cid:18) + m ε (cid:19) (cid:18) − m (cid:48) ε (cid:48) (cid:19) − e i τ z γ (cid:114)(cid:18) − m ε (cid:19) (cid:18) + m (cid:48) ε (cid:48) (cid:19) , U − + kk (cid:48) = U (cid:104) u − k ( τ z ) | u + k (cid:48) ( τ z ) (cid:105) = U (cid:114)(cid:18) − m ε (cid:19) (cid:18) + m (cid:48) ε (cid:48) (cid:19) − e i τ z γ (cid:114)(cid:18) + m ε (cid:19) (cid:18) − m (cid:48) ε (cid:48) (cid:19) , U −− kk (cid:48) = U (cid:104) u − k ( τ z ) | u − k (cid:48) ( τ z ) (cid:105) = U (cid:114)(cid:18) − m ε (cid:19) (cid:18) − m (cid:48) ε (cid:48) (cid:19) + e i τ z γ (cid:114)(cid:18) + m ε (cid:19) (cid:18) + m (cid:48) ε (cid:48) (cid:19) , (S41)where γ = θ (cid:48) − θ . Then we obtain (cid:104) U ++ kk (cid:48) U + − k (cid:48) k (cid:105) = n imp U (cid:114)(cid:18) + m ε (cid:19) (cid:18) + m (cid:48) ε (cid:48) (cid:19) + e i τ z γ (cid:114)(cid:18) − m ε (cid:19) (cid:18) − m (cid:48) ε (cid:48) (cid:19) (cid:114)(cid:18) + m (cid:48) ε (cid:48) (cid:19) (cid:18) − m ε (cid:19) − e − i τ z γ (cid:114)(cid:18) − m (cid:48) ε (cid:48) (cid:19) (cid:18) + m ε (cid:19) = n imp U (cid:34) v F k ε m (cid:48) ε (cid:48) + (cid:18) i τ z sin γ − m ε cos γ (cid:19) v F k (cid:48) ε (cid:48) (cid:35) , (cid:104) U + − kk (cid:48) U −− k (cid:48) k (cid:105) = n imp U (cid:114)(cid:18) + m ε (cid:19) (cid:18) − m (cid:48) ε (cid:48) (cid:19) − e i τ z γ (cid:114)(cid:18) − m ε (cid:19) (cid:18) + m (cid:48) ε (cid:48) (cid:19) (cid:114)(cid:18) − m (cid:48) ε (cid:48) (cid:19) (cid:18) − m ε (cid:19) + e − i τ z γ (cid:114)(cid:18) + m (cid:48) ε (cid:48) (cid:19) (cid:18) + m ε (cid:19) = n imp U (cid:34) − v F k ε m (cid:48) ε (cid:48) + (cid:18) − i τ z sin γ + m ε cos γ (cid:19) v F k (cid:48) ε (cid:48) (cid:35) = −(cid:104) U ++ kk (cid:48) U + − k (cid:48) k (cid:105) , (cid:104) U −− kk (cid:48) U − + k (cid:48) k (cid:105) = (cid:104) U ++ kk (cid:48) U + − k (cid:48) k (cid:105) ( m → − m , m (cid:48) → − m (cid:48) ) = − (cid:104) (cid:104) U ++ kk (cid:48) U + − k (cid:48) k (cid:105) (cid:105) ∗ , (cid:104) U − + kk (cid:48) U ++ k (cid:48) k (cid:105) = (cid:104) U + − kk (cid:48) U −− k (cid:48) k (cid:105) ( m → − m , m (cid:48) → − m (cid:48) ) = (cid:104) (cid:104) U ++ kk (cid:48) U + − k (cid:48) k (cid:105) (cid:105) ∗ , (S42)From the definition of J ( (cid:104) n (cid:105) ) + − k with (cid:104) n (cid:105) = diag[ n + k , n − k ], we have[ J ( (cid:104) n (cid:105) )] + − k = π (cid:88) k (cid:48) (cid:104) U ++ kk (cid:48) U + − k (cid:48) k (cid:105) (cid:104) ( n + k − n + k (cid:48) ) δ ( ε + k − ε + k (cid:48) ) + ( n − k − n + k (cid:48) ) δ ( ε − k − ε + k (cid:48) ) (cid:105) + π (cid:88) k (cid:48) (cid:104) U + − kk (cid:48) U −− k (cid:48) k (cid:105) (cid:104) ( n + k − n − k (cid:48) ) δ ( ε + k − ε − k (cid:48) ) + ( n − k − n − k (cid:48) ) δ ( ε − k − ε − k (cid:48) ) (cid:105) = π (cid:88) k (cid:48) (cid:104) U ++ kk (cid:48) U + − k (cid:48) k (cid:105) (cid:104) ( n + k − n + k (cid:48) ) δ ( ε + k − ε + k (cid:48) ) − ( n − k − n − k (cid:48) ) δ ( ε − k − ε − k (cid:48) ) (cid:105) , (S43)where we have used δ ( ε − k − ε + k (cid:48) ) = δ ( ε + k − ε − k (cid:48) ) =
0. Similarly, we get[ J ( (cid:104) n (cid:105) )] − + k = π (cid:88) k (cid:48) (cid:104) U − + kk (cid:48) U ++ k (cid:48) k (cid:105) (cid:104) ( n − k − n + k (cid:48) ) δ ( ε − k − ε + k (cid:48) ) + ( n + k − n + k (cid:48) ) δ ( ε + k − ε + k (cid:48) ) (cid:105) + π (cid:88) k (cid:48) (cid:104) U −− kk (cid:48) U − + k (cid:48) k (cid:105) (cid:104) ( n − k − n − k (cid:48) ) δ ( ε − k − ε − k (cid:48) ) + ( n + k − n − k (cid:48) ) δ ( ε + k − ε − k (cid:48) ) (cid:105) = π (cid:88) k (cid:48) (cid:104) U ++ kk (cid:48) U + − k (cid:48) k (cid:105) ∗ (cid:104) ( n + k − n + k (cid:48) ) δ ( ε + k − ε + k (cid:48) ) − ( n − k − n − k (cid:48) ) δ ( ε − k − ε − k (cid:48) ) (cid:105) = (cid:2) J ( (cid:104) n (cid:105) ) + − (cid:3) ∗ . (S44) A. Correction to the linear magnetoconductivity due to disorder scattering
Let us calculate the contribution from J ( (cid:104) n (cid:105) ) to the linear magnetoconductivity, which corresponds to the vertex correction inthe ladder-diagram approximation. The correction to (cid:104) S B (cid:105) [Eq. (S17)] is given by (cid:104) S (cid:48) B (cid:105) nn (cid:48) k = i [ J ( (cid:104) n B (cid:105) )] nn (cid:48) k ε n k − ε n (cid:48) k = i J + − ε + k − ε − k − i ( J + − ) ∗ ε + k − ε − k = − ˜ σ y Re (cid:110) [ J ( (cid:104) n B (cid:105) )] + − k (cid:111) ε k − ˜ σ x Im (cid:110) [ J ( (cid:104) n B (cid:105) ) + − k (cid:111) ε k , (S45)2where Re (cid:110) [ J ( (cid:104) n B (cid:105) )] + − k (cid:111) = π n imp U (cid:88) k (cid:48) (cid:34) v F k ε m (cid:48) ε (cid:48) − m ε v F k (cid:48) ε (cid:48) (cos θ (cid:48) cos θ + sin θ sin θ (cid:48) ) (cid:35) ( n + B k − n + B k (cid:48) ) δ ( ε + k − ε + k (cid:48) ) = π n imp U v F k ε n + B k (cid:88) k (cid:48) m (cid:48) ε (cid:48) δ ( ε + k − ε + k (cid:48) ) + π n imp U m ε sin θ (cid:88) k (cid:48) v F k (cid:48) ε (cid:48) sin θ (cid:48) n + B k (cid:48) δ ( ε + k − ε + k (cid:48) ) , (S46)and Im (cid:110) [ J ( (cid:104) n B (cid:105) )] + − k (cid:111) = τ z π n imp U (cid:88) k (cid:48) (sin θ (cid:48) cos θ − cos θ (cid:48) sin θ ) v F k (cid:48) ε (cid:48) ( n + B k − n + B k (cid:48) ) δ ( ε + k − ε + k (cid:48) ) = − τ z π n imp U θ (cid:88) k (cid:48) sin θ (cid:48) v F k (cid:48) ε (cid:48) n + B k (cid:48) δ ( ε + k − ε + k (cid:48) ) , (S47)with n + B k = τ tr eB z ∂ε k ∂ k y ∂ n + E k ∂ k x − ∂ε k ∂ k x ∂ n + E k ∂ k y = τ e E x B z (cid:32) ∂ε k ∂ k y ∂∂ k x − ∂ε k ∂ k x ∂∂ k y (cid:33) ∂ f ( ε + k ) ∂ k x , (S48)and n + E k = τ tr eE x ∂ f ( ε + k ) ∂ k x [see Eq. (S18)]. Here, we have used the fact that n + B k is an odd function of k y , and sin γ = sin θ (cid:48) cos θ − cos θ (cid:48) sin θ and cos γ = cos θ (cid:48) cos θ + sin θ sin θ (cid:48) . Then we have v x (cid:104) S (cid:48) B (cid:105) = v F (cid:32) ˜ σ z v F k ε k cos θ + ˜ σ y τ z sin θ − ˜ σ x m ε k cos θ (cid:33) − ˜ σ y Re (cid:110) [ J ( (cid:104) n B (cid:105) )] + − k (cid:111) ε k − ˜ σ x Im (cid:110) [ J ( (cid:104) n B (cid:105) ) + − k (cid:111) ε k = τ z v F π n imp U − m ε k (cid:88) k (cid:48) sin θ (cid:48) v F k (cid:48) ε (cid:48) n + B k (cid:48) δ ( ε + k − ε + k (cid:48) ) − v F k ε sin θ n + B k (cid:88) k (cid:48) m (cid:48) ε (cid:48) δ ( ε + k − ε + k (cid:48) ) + (traceless terms) . (S49)Finally, we obtain the valley-dependent vertex correction due to disorder scattering to the linear magnetoconductivity σ (1)ver xx ( B z ) = Tr (cid:2) ( − e ) v x (cid:104) S (cid:48) B (cid:105) (cid:3) / E x = (2 e ) τ z v F π n imp U (cid:90) ∞−∞ d k (2 π ) (cid:90) ∞−∞ d k (cid:48) (2 π ) δ ( ε + k − ε + k (cid:48) ) m ε k sin θ (cid:48) v F k (cid:48) ε (cid:48) n + B k (cid:48) + v F k ε sin θ n + B k m (cid:48) ε (cid:48) ≡ τ z e (cid:126) B z v F C dis1 ( µ, m ) τ tr , (S50)where C dis1 ( µ, m ) < τ tr = (4 (cid:126) v F / n imp U µ )(1 + m /µ ) − is the intravalley scattering time (see the next section for a detailedcalculation). Note that n + B k ∝ τ , while the coe ffi cient in front of the integral is proportional to τ − . B. Calculation of the intravalley scattering time
Let us calculate the intravalley scattering time (or equivalently transport lifetime) in the two-band massive Dirac model wehave studied. The intravalley scattering time of band m can be calculated from the following well-known equation:1 τ m tr = π (cid:126) (cid:90) d k (cid:48) (2 π ) (cid:104)| U mm kk (cid:48) | (cid:105) (1 − cos θ kk (cid:48) ) δ ( µ − ε m k (cid:48) ) , (S51)where cos θ kk (cid:48) = k · k (cid:48) / | k || k (cid:48) | . Because (cid:104)| U ++ kk (cid:48) | (cid:105) = (cid:104)| U −− kk (cid:48) | (cid:105) , we may drop the band index and consider the case of µ >
0. FromEq. (S41) we get (cid:104)| U ++ kk (cid:48) | (cid:105) = n imp U (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:114)(cid:18) + m ε (cid:19) (cid:18) + m (cid:48) ε (cid:48) (cid:19) + (cos γ + i τ z sin γ ) (cid:114)(cid:18) − m ε (cid:19) (cid:18) − m (cid:48) ε (cid:48) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = n imp U + m ε k ε k (cid:48) + cos γ v F kk (cid:48) ε k ε k (cid:48) . (S52)3We define the Fermi wave number k = k (cid:48) = k F ≡ (cid:112) µ − m / v F ( ε k = ε k (cid:48) = µ ) and θ kk (cid:48) ≡ φ . Also, without loss of generalitywe can set k = ( k F ,
0) (i.e., θ = γ = cos φ . Finally, we obtain1 τ tr = π (cid:126) n imp U π (cid:90) ∞ k (cid:48) dk (cid:48) (cid:90) π d φ + m ε k ε k (cid:48) + cos γ v F kk (cid:48) ε k ε k (cid:48) δ ( µ − ε m k (cid:48) ) , = n imp U µ (cid:126) v F (cid:32) + m µ (cid:33) . (S53) VIII. APPROXIMATE FORMS OF THE LONGITUDINAL MAGNETORESISTANCE ∆ ρ xx Let us derive approximate forms of the longitudinal magnetoresistance ∆ ρ xx ( B z ) ≡ [ ρ xx ( B z ) − ρ xx (0)] /ρ xx (0), where ρ xx ( B z ) = σ xx / ( σ xx + σ xy ) with σ µν ≡ σ (0) µν + σ B µν , in the low-field limit with the condition σ (0) xx (cid:29) σ (0) xy or σ (0) xx (cid:28) σ (0) xy for the massive Diracmodel (S10). In this model, the anomalous Hall conductivity σ (0) xy takes the maximum value τ z e / h when the Fermi level lies inthe gap.First, let us consider the case of σ (0) xx (cid:29) σ (0) xy , i.e., the case of high-carrier-desity semiconductors ( µ > m ). In this case, we have ρ xx ( B z ) = σ (0) xx + σ Bxx ( σ (0) xx + σ Bxx ) + ( σ (0) xy + σ Bxy ) ≈ σ (0) xx + σ Bxx ( σ (0) xx ) − σ (0) xx σ Bxx + σ (0) xy σ Bxy ( σ (0) xx ) ≈ σ (0) xx + σ Bxx ( σ (0) xx ) − σ (0) xx σ Bxx + σ (0) xy σ Bxy ( σ (0) xx ) ≈ σ (0) xx − σ Bxx ( σ (0) xx ) , (S54)where we have used that ( σ B µν /σ (0) xx ) (cid:28) σ (0) xy /σ (0) xx )( σ Bxy /σ (0) xx ) (cid:28)
1. We also have ρ xx (0) = σ (0) xx ( σ (0) xx ) + ( σ (0) xy ) ≈ /σ (0) xx . Then thelongitudinal magnetoresistance ∆ ρ xx ( B z ) in the low-field limit is obtained as ∆ ρ xx ( B z ) ≈ − σ Bxx σ (0) xx . (S55)Next, let us consider the case of σ (0) xx (cid:28) σ (0) xy , i.e., the case of low-carrier-density semiconductors (with small m ) or the case of µ ≈ m . In this case, we have ρ xx ( B z ) = σ (0) xx + σ Bxx ( σ (0) xx + σ Bxx ) + ( σ (0) xy + σ Bxy ) ≈ σ (0) xx + σ Bxx ( σ (0) xy ) − σ (0) xx σ Bxx + σ (0) xy σ Bxy ( σ (0) xy ) ≈ σ (0) xx + σ Bxx ( σ (0) xy ) − σ (0) xx σ (0) xx σ Bxx + σ (0) xy σ Bxy ( σ (0) xy ) ≈ σ (0) xx + σ Bxx ( σ (0) xy ) , (S56)where we have used that ( σ B µν /σ (0) xy ) (cid:28) σ (0) xx /σ (0) xy )( σ Bxy /σ (0) xy ) (cid:28)
1. We also have ρ xx (0) = σ (0) xx ( σ (0) xx ) + ( σ (0) xy ) ≈ σ (0) xx / ( σ (0) xy ) . Thenthe longitudinal magnetoresistance ∆ ρ xx ( B z ) in the low-field limit is obtained as ∆ ρ xx ( B z ) ≈ + σ Bxx σ (0) xx ..