Quantum Kinetic Theory of the Chiral Anomaly

3 downloads 0 Views 569KB Size Report
Jun 5, 2017 - Berry curvature over the full 2D Brillouin zone divided by 2π, and must be an ...... We will show in Sec. VIC that DB () describes the gen-.
Quantum Kinetic Theory of the Chiral Anomaly Akihiko Sekine,1 Dimitrie Culcer,2 and Allan H. MacDonald1

arXiv:1706.01200v1 [cond-mat.mes-hall] 5 Jun 2017

2

1 Department of Physics, The University of Texas at Austin, Austin, Texas 78712, USA School of Physics and Australian Research Council Centre of Excellence in Low-Energy Electronics Technologies, UNSW Node, The University of New South Wales, Sydney 2052, Australia (Dated: June 6, 2017)

We present a general quantum kinetic theory of low-field magnetotransport in weakly disordered crystals that accounts fully for the interplay between electric-field induced interband coherence, Bloch-state scattering and an external magnetic field. The quantum kinetic equation we derive for the Bloch-state density matrix naturally incorporates the momentum-space Berry phase effects whose influence on Bloch-state wavepacket dynamics is normally incorporated into transport theory in an ad hoc manner. For example, the Berry phase correction to the momentum-space density of states in the presence of an external magnetic field implied by semiclassical wavepacket dynamics appears as an intrinsic density-matrix response to a magnetic field. Our theory provides a simple and general procedure for expanding the linear response of the Bloch-state density matrix to an electric field in powers of magnetic field. As an illustration, we apply our theory to the magnetotransport of Weyl semimetals. We show that the chiral anomaly (positive magnetoconductivity quadratic in magnetic field) associated with the presence of separate Fermi surface pockets surrounding distinct Weyl points only partially survives disorder scattering even in the limit of extremely weak intervalley scattering.

I.

INTRODUCTION

The electrical transport properties of weakly disordered quantum degenerate crystalline conductors are normally controlled by changes in the occupation probabilities of states close to the Fermi surface in response to electric fields or temperature gradients. Transport properties then depend on the shapes of all Fermi surface valleys, and on the distributions of group velocities and the details of the processes that scatter electrons between Bloch states on those valleys. It has recently become more widely appreciated however that intrinsic effects, independent of disorder and related more to Bloch state wave functions than energies, are sometimes important. In special cases these effects are nearly completely characterized by momentum-space Berry curvatures [1], which can be non-zero only in systems that have broken timereversal symmetry or broken inversion symmetry or both. One important example of a transport effect that is often dominated by Berry curvature physics is the anomalous Hall effect [2], i.e.,the Hall effect in the absence of a magnetic field in a crystal in which time-reversal symmetry is broken by magnetic order. This paper is motivated by the need for a practical theory suitable for application to real crystals that can provide a general description of magnetotranport in systems in which Bloch state wavefunction properties play an important role. Large momentum-space Berry curvatures are often associated with nontrivial band topology. For example the anomalous Hall conductivity of two-dimensional (2D) insulators [3, 4] is equal to the quantum unit of conductance, e2 /h, times the sum of the Chern numbers of all occupied bands. The Chern number of a band is a topological index that is simply the integral of its Berry curvature over the full 2D Brillouin zone divided

by 2π, and must be an integer. Indeed, recent interest in the topological classification of crystalline matter started from theories of the quantum Hall effect, and was later generalized to time-reversal invariant topological insulators [5–10]. The classification of topological phases has now been extended to 3D gapless metallic systems like Weyl semimetals [11–13] and Dirac semimetals [14– 17]. Although these systems cannot be characterized by bulk topological invariants, non-degenerate band touching points (Weyl points) are topologically stable and can be regarded as monopoles sources of phase flux in momentum space [18, 19]. Dirac points, which can be viewed as a superposition of two degenerate Weyl points, are not generically stable, but can be stabilized by crystalline symmetries [17]. The experimental identification of Dirac semimetals [20–23] and Weyl semimetals [24–27] has motivated a growing effort, with both theoretical and experimental components, aimed at identifying and exploring novel phenomena in these materials. Three dimensional topological semimetals host a variety of transport and magnetotransport properties that are dependent on Bloch state wavefunction response, and therefore require the type of transport theory discussed in this paper for a proper theoretical description. Topologically nontrivial electronic band structures can lead to an approximate realization of the chiral anomaly in condensed matter physics. In quantum field theory the chiral anomaly, which occurs only in systems with odd space dimensions, refers to the violation of axial current conservation, ∂µ J5µ 6= 0. In the case of U (1) electromagnetic gauge fields, the chiral anomaly is described by a so-called θ term in the electromagnetic action [28]: Sθ =

e2 4π 2 ~c

Z

dt d3 r θ(r, t)E · B,

(1)

2 where θ(r, t) is a coupling coefficient, and E and B are external electric and magnetic fields. Three-dimensional materials whose electromagnetic effective action includes a finite θ term are referred to as magnetoelectric materials. magnetoelectricity is possibly only in materials with broken time-reversal and and inversion symmetry. A 3D topological insulator acts like a magnetoelectric material with a value of θ that is close to π [29–33] when its timereversal symmetry is weakly broken by ferromagnetic order, and the sign of the magnetization component perpendicular to its surface does not change sign spatially, Some previous works have argued that the electromagnetic effective action of a Weyl semimetal with two Weyl nodes includes a θ term with θ(r, t) = 2(b · r − b0 t) [34– 39], where b is the vector connecting the two Weyl nodes in momentum space, and b0 is a difference between local chemical potentials established in different regions of momentum space. More generally 3D insulators whose electromagnetic Lagrangians contain a θ(r, t) term related to a low-energy degree of freedom with important quantum or thermal fluctuations, for example the N´eel vector order parameter [40] of an antiferromagnet, are referred to as axionic insulators. The θ term implies a charge current ˙ t)B] j(r, t) = δSθ /δA = (e2 /4π 2 ~c)[∇θ(r, t) × E + θ(r, [28]. Here, A is the electromagnetic field’s vector potential. The electric-field induced current is an anomalous Hall effect, since the current is perpendicular to the electric field. The magnetic-field induced current is the so-called chiral magnetic effect [34, 35, 41]. A related dynamical realization of the chiral magnetic effect in axionic insulators was recently proposed theoretically [49, 50]. In 3D topological insulators with θ = π, the presence of the θ term implies that P = δSθ /δE = (e2 /4π~c)B and M = δSθ /δB = (e2 /4π~c)E [29], where P and M are the electric polarization and magnetization, respectively. This is the topological magnetoelectric effect. Our theory is able to realistically address the partial realization of these types of effects in real materials. The chiral magnetic effect in Weyl semimetals refers to charge current generation by a magnetic field when a chemical potential difference is established between different Weyl nodes. Although the possibility of an equilibrium chiral magnetic effect has been raised theoretically [42–48], it is generally agreed that this possibility can be ruled out in crystalline solids [42]. A chemical potential difference can however be established by the combined effect of electric and magnetic fields, as explained below, and when present is responsible via the chiral magnetic effect for a negative magnetoresistance (positive magnetoconductance) that is quadratic in magnetic field [51–54]. A positive magnetoconductance stands in striking contrast to the familiar negative magnetoconductance due to the Lorentz force. Remarkably, negative magnetoresistance has recently been observed in the low magnetic field regime in the Dirac semimetals Na3 Bi [55], Cd3 As2 [56, 57], and ZrTe5 [58], and in the Weyl semimetals TaAs [59] and TaP [60]. The peculiar positive magnetoconductance behavior occurs only for par-

allel electric and magnetic fields, which suggests that it is related to the E · B contribution to the electromagnetic action, i.e.,to the θ term. This very specific effect is a particularly attractive target of the magnetotransport theory developed here. The anomalous positive magnetoconductance found in Weyl semimetals was anticipated theoretically by observing that equations for the semiclassical wavepacket dynamics of Bloch electrons in a crystal predict that electrons can be pumped between Fermi surface pockets when E · B 6= 0. Under free semiclassical dynamics, the number of electrons in a given valley changes at the rate [35, 51] ∂N e2 =Q 2 2 E·B ∂t 4π ~ c

(2)

where Q=

Z

d3 k ∂f0 (εm k) m vk · Ω m k 2π~ ∂εm k

(3)

is the chirality of the valley. In Eq. (3) Ωm k is the Berry curvature of band m (that possesses the Fermi surface), vkm is the Bloch state group velocity, and f0 (εm k ) is the Fermi-Dirac distribution function whose derivative with respect to energy εm k provides a δ-function at the Fermi energy. It is easy to show that the chirality of each valley must be an integer equal to the sum of the chiralities [18] of all the Weyl (band crossing) points enclosed by its Fermi surface, and that the sum of chirality over all Fermi surface valleys vanishes. When Bloch state scattering between valleys is much weaker than scattering within valleys, a circumstance that is common in semimetals, the anomalous E · B pumping leads to differences between the chemical potentials of different valleys and to a corresponding anomalous contribution to the current. Berry curvature corrections to semiclassical electron dynamics capture the modification of Bloch state wavefunctions by an electric field, i.e., electric-field induced mixing of the Bloch states that diagonalize the Hamiltonian in the absence of fields. In this paper, we develop a general quantum kinetic theory of low-field magnetotransport in weakly disordered crystals that accounts fully for electric-field induced interband coherence and for the interplay between Bloch state scattering and the presence of an external magnetic field. We take the effect of magnetic fields into account using a semiclassical approximation that we expect to be accurate when the magnetic field is weak enough that Landau quantization can be neglected. Our main result is a quantum kinetic equation which includes driving terms associated with both external electric and magnetic fields. Our goal is to derive a general transport formalism that will be valuable for understanding the interesting magnetotransport anomalies that appear in any conductor that has band crossings at energies close to the Fermi energy. As an illustration we apply it to the magnetoconductivity of Weyl semimetals. We find that the chiral magnetotransport anomaly is observable only when intervalley

3 scattering at the Fermi energy is very weak compared to intravalley scattering, and that it survives only partially even in this limit. Specifically we obtain the following general expression for the rate of pumping into a valley when the electric and magnetic fields are both along the zˆ-direction:   m m d3 k ∂f0 (εm m k) Ωm vk,x Ωk,x + vk,y k,y . 2π~ ∂εm k (4) m The term in Eq. (3) proportional to vk,z is absent in Eq. (4), because free semiclassical dynamics are disturbed by intravalley scattering of electrons, i.e., those electrons are swept across the valley m Fermi surface by the electric field. This smaller rate of pumping results in a smaller value of the positive quadratic magnetoconductivity induced by the chiral anomaly than that obtained by invoking semiclassical wavepacket dynamics and neglecting scattering. This paper is organized as follows. In Sec. II we derive a general quantum kinetic equation that accounts fully for momentum-dependent Bloch state wavefunctions and associated interband coherence response in the presence of external electric and magnetic fields. In Sec. III we describe a general scheme to apply our theory to magnetotransport phenomena, by performing a systematic low-field expansion. In Sec. IV we explain in detail when the intervalley scattering rate emerges as an important parameter in transport properties of many-valley electronic systems, and when it does not. As an application of our formalism, we explain in Sec. V how our formalism quite generally captures the pumping between valleys that occurs in systems with Fermi surface pockets that enclose Weyl points with a net chirality, and why the rate of pumping differs from the expression obtained when only intrinsic Bloch state dynamics, and not scattering, is considered. Then in Sec. VI apply our quantum kinetic equation to the standard top-model of a two-valley Weyl semimetal, showing that our formalism fully captures the anomalous Hall effect and the chiral magnetic effect in Weyl semimetals. In Sec. VII we report on the derivation of an explicit expression for positive magnetoconductance quadratic in magnetic field for a simplified model of a semimetal with two Weyl points. Finally in Sec. VIII we discuss our results and comment on some other potentially interesting applications of our transport theory. ∂N e2 Ez Bz = ∂t 4π 2 ~2 c

Z

we refer to as the eigenstate basis: H0 |m, ki = εm k |m, ki,

(5)

where H0 is the crystal Hamiltonian, εm k is an eigenvalue of H0 , k is a momentum in the crystal’s Brillouin-zone, and m is a band index. The eigenstates of a crystal Hamiltonian are Bloch states, products of plane-waves with wave-vector k and periodic functions that are eigenfunctions of the k-dependent k · p Hamiltonian. A Bloch state does not, of course, have to be an eigenstate of the crystal Hamiltonian. Our formalism assumes that a good approximation to the Hamiltonian matrix of the perfect crystal is known in a representation of k-independent periodic functions. In the following we assume that the Hamiltonian is known in a representation of Wannier functions, but the formalism could be applied with little change if we used a representation of k · p eigenstates at a particular Brillouin-zone point of interest. Given this starting point, we write X X αα′ ′ ′ H0 = εm HLL ′ |α, Lihα , L | k |m, kihm, k| = k,m

LL′ αα′

(6)



αα ′ ′ with HLL ′ = hα, L|H0 |α , L i. Here |α, Li is a Wannier function associated with orbital α and real-space lattice vector L. For a particular orbital the Bloch and Wannier functions are related by

1 X −ik·L |α, Li = √ e |α, ki, N k 1 X +ik·L |α, ki = √ e |α, Li, N k

(7)

where N is the number of Bravais lattice sites in the crystal. |α, ki is a Bloch state with wavefunction hr|α, ki = eik·r uα k (r),

(8)

where r is position and uα k (r) is the cell-periodic part of the Bloch wave function. Our transport theory is formulated in terms of the single-particle density matrix. The density-matrix operator can be expressed in either the eigenstate or Wannier representations using X X ′ ′ ′ ′ ′ ′ ρ= ρmm ραα kk′ |m, kihm , k | = LL′ |α, Lihα , L |, LL′ αα′

kk′ mm′

(9)



II.

QUANTUM KINETIC EQUATION

In this section, we derive a quantum kinetic equation for Bloch electrons in the presence of electric and magnetic fields. To this end, we start with some general considerations related to Bloch Hamiltonians and basis function choices. Throughout this paper, we work in the basis of the disorder-free Hamiltonian eigenstates, which



′ ′ αα ′ ′ where ρmm kk′ = hm, k|ρ|m , k i and ρLL′ = hα, L|ρ|α , L i. Since we focus on circumstances in which translational symmetry is not broken, the non-equilibrium expectation value of the density matrix operator ρ will always be diagonal in Bloch state wavevector. The eigenstate basis and the Wannier basis are related to each other through the unitary transformation related to P band eigenvectors m in the α-representation: |m, ki = α zα |α, ki with m zα = hα, k|m, ki.

4 A.

Introducing disorder

The main result of this paper is a generic quantum kinetic equation that accounts for disorder and electric and magnetic fields. Throughout this paper we take disorder into account within the Born approximation, implicitly assuming therefore that disorder is weak. We follow a procedure developed in an earlier paper [61], generalizing it to allow for magnetic fields. To establish some notation we first consider the case without external electric and magnetic fields. We start with the quantum Liouville equation i ∂ρ + [H, ρ] = 0, ∂t ~

(10)

where ρ and H are the density-matrix operator and the Hamiltonian of the system, respectively. In the absence of external fields, the total Hamiltonian of the system is H = H0 + U where U is the disorder potential. We decompose the density matrix ρ into two parts, writing ρ = hρi + g, where hρi is the density matrix averaged over disorder configurations, and g is the deviation from this average. The quantum Liouville equation can then be decomposed into coupled equations for hρi and g: ∂hρi i i + [H0 , hρi] + h[U, g]i = 0, ∂t ~ ~

(11)

i i i i ∂g + [H0 , g] + [U, g] − h[U, g]i = − [U, hρi]. (12) ∂t ~ ~ ~ ~ In the Born approximation we can ignore the last two terms on the left hand side of Eq. (12). The equation for g, Eq. (12), can then be solved straightforwardly. By substituting the solution for g into Eq. (11), we obtain ∂hρi i + [H0 , hρi] + K(hρi) = 0, ∂t ~

(13)

where the scattering term K(hρi) is given by [61] 1 K(hρi) = 2 ~

Z

0



dt′

Dh iE ′ ′ U, [e−iH0 t /~ U eiH0 t /~ , hρ(t)i] ,

(14)

where hρ(t)i = e−iH0 t/~ hρieiH0 t/~ . In the following, we do not exhibit the explicit time dependence of hρ(t)i in order to simplify the notation. We separate the density matrix into band-diagonal and band-off-diagonal parts because these two types of components behave quite differently using the notation hρi = hni + hSi, where hni is diagonal in band index and hSi is off-diagonal. The scattering kernel can be separated into four parts which map hni and hSi to band-diagonal and band off-diagonal contributions to ∂hρ(t)i/∂t. This separation is discussed at length later in connection with magnetotransport. Focusing on the elastic-scattering case

and using the notation of Ref. [61] we find that the banddiagonal to band-diagonal part of K(hρi) from hni is [I(hni)]mm = k

′ 2π X mm′ m′ m m m m′ hU ′ Uk′ k i(nk − nm k′ )δ(εk − εk′ ), ~ ′ ′ kk

mk

(15)

with m and m′ being band indices. This is exactly Fermi’s golden rule. Similarly, the band-diagonal to band-off-diagonal part of K(hρi) from hni is ′′ ′ π X mm′ m′ m′′ h m [J(hni)]mm = hUkk′ Uk′ k i (nk − nm k k′ ) ~ ′ ′ mk i ′′ m m′ m′ m′′ m′ ×δ(εk − εk′ ) + (nm − n )δ(ε − ε ) , (16) ′ ′ k k k k

where m 6= m′′ . The full expression for K(hρi), including the contributions the off-diagonal contributions to ρ (hSi) can be found in Ref. [61]. B.

Introducing electric and magnetic fields

Next, we take the effects of magnetic fields into account using a semiclassical approximation that we expect to be accurate when the the weak magnetic fields condition ωc τ ≪ 1 is satisfied and Landau quantization can to be neglected. Here, ωc is the cyclotron frequency and τ is the transport relaxation time discussed at greater length below. In order to derive the kinetic equation for systems in magnetic fields, we apply the Wigner transformation to the quantum Liouville equation (10). We consider a generic single-particle D-dimensional Bloch Hamiltonian H0 (p) with momentum operator p = −i~∇. In the presence of a vector potential A(r, t), minimal coupling results in p → p + eA and adopt the notation that e > 0 is the magnitude of the electron charge. The Wigner distribution function in the presence of a vector potential A is defined by [62] Z hρimn (r) = dD R e−(i/~)P ·R hm, r+ |ρ|n, r− i, (17) p where P = p − eA, r± P = r ± R/2, and |n, ri = P ik·r |n, ki is the Fourier |m, kihm, k|n, ri = ke m,k transform of |n, ki. Note that the sign in front of eA in P is different from the one in minimal coupling [62]. We project the quantum Liouville equation (10) onto the {|m, ri} space. Let us consider the term hm, r+ |[H0 , ρ]|n, r− i Z X = d3 r′ [hm, r+ |H0 |m′ , r ′ ihm′ , r ′ |ρ|n, r− i m′

−hm, r+ |ρ|m′ , r ′ ihm′ , r ′ |H0 |n, r− i] .

(18)

Since the Hamiltonian is diagonal in real space, we have hm, r1 |H0 |m′ , r2 i = (hm′ , r2 |H0 |m, r1 i)∗ = P −iq·r1 iq ′ ·r2 hm, q|H0 (p1 + eA1 )|m′ , q ′ i with p1 = e q,q ′ e

5 −i~∂r1 . We expand p± and  A± up to linear order in R as p± = −i~ 21 ∇ ± ∇R and A± = A ± 12 (R · ∇) A, where ∇R ≡ ∂/∂R. Finally we obtain hm, q|H0 (p+ + A+ )|m′ , q ′ i 1 = hm, q|H0 (q)|m′ , q ′ i − i~ hm, q|∇q H0 (q) · ∇|m′ , q ′ i 2 1 + e hm, q|∇q H0 (q) · [(R · ∇) A] |m′ , q ′ i (19) 2 with q = −i~∇R +eA = p+eA. In reaching Eq. (19) we have observed that it follows from the definition of the Wigner distribution function (17) that R = i~∇p and p = −i~∇R . We perform the Wigner transformation on the quantum Liouville equation (10) as   Z i ∂ρ + [H0 , ρ] |n, r− i = 0. dD R e−(i/~)P ·R hm, r+ | ∂t ~ (20) Note that the scattering term K(hρi) is not changed after the Wigner transformation. We use the following identities e−(i/~)P ·R ∇ = {∇ − (ie/~)[∇(A · R)]} e−(i/~)P ·R ,

e−(i/~)P ·R ∂/∂t = [∂/∂t + (ie/~)E · R] e−(i/~)P ·R , (21)

where E = −∂A/∂t is the electric field. Note that we can obtain the same electric-field dependent term as in Eq. (21) by taking the electric field into account using a scalar potential φ(r) = −E · r, i.e., by adding HE = −eφ(r) = eE · r to the Hamiltonian as H0 → H0 + HE in Eq. (18) [61]. By multiplying Eq. (19) by e−(i/~)P ·R from the left we see that the second term in the right hand side of Eq. (19) generates an additional term proportional to ∇(A · R). Finally, the vector potential can be replaced by the magnetic field using the identity ∇(A · R) − (R · ∇) A = R × B,

(22)

where B = ∇ × A is the magnetic field. Notice that the matrix element of ∇q H0 (q) in Eq. (19) can be written as hm, q|∇q H0 (q)|m′ , q ′ i = δ(q − q ′ )× h i mm′ n′ n′ m′ mn′ n′ m′ ∇q H0q + hum |∇ u iH − H hu |∇ u i , q q q q 0q 0q q q

(23)



mm where H0q = hm, q|H0 (q)|m′ , qi = δmm′ εm q and m |uq i is the periodic part of the Bloch function. The magnetic-field dependent term in Eq. (19) becomes P − 12 e n′ ,q ′′ hm, q|∇q H0 (q)|n′ , q ′′ i · hn′ , q ′′ |R|m′ , q ′ i × B. Here, the term hn′ , q ′′ |R|m′ , q ′ i (with R = i~∇q ) acts on hm′ , r ′ |ρ|n, r− i in Eq. (18), which finally results in ′

n hn′ , q ′′ |R|m′ , q ′ ihρim = i~ δ(q ′′ − q ′ )δ(q ′ − q)× q h i ′ ′ ′ m′ n n′ m′ m′ n ∇q hρinq n + hunq |∇q um ihρi − hρi hu |∇ u i , q q q q q q

(24)

where the third term on the right hand side comes from the Hermiticity of the equation and consistency with the single-band limit. Combining Eq. (13) with the analysis above, we arrive at the final form for the quantum kinetic equation in the presence of disorder, an electric field E, and a magnetic field B:   ∂hρi i 1 DH0 + [H0 , hρi] + · ∇hρi + K(hρi) ∂t ~ 2~ Dk = DE (hρi) + DB (hρi), (25) where k = q/~ is the crystal wave vector, and the angle brackets in hρi imply an impurity average. Here and below {a · b} ≡ a · b + b · a (with a and b being vectors) denotes a symmetrized operator product. In Eq. (25) DE (hρi) and DB (hρi) are the electric and magnetic driving terms: DE (hρi) =

eE Dhρi · ~ Dk

(26)

and DB (hρi) =

e 2~2



DH0 ×B Dk



·

Dhρi Dk



.

(27)

Eqs. (26) and (27) have been simplified by introducing the covariant derivative notation defined by Dhρi = ∇k hρi − i[Rk , hρi], (28) Dk P n where Rk = a=x,y,z Rak ea with [Rak ]mn = ihum k |∂ka uk i being the generalized Berry connection. From another point of view, the covariant derivative is simply the derivative evaluated in the k-independent Wannier state representation. Note that the covariant k-derivative acting on the density matrix accounts for changes in the density matrix due to the fact that its band-eigenstate representation elements are k-dependent and also due to the fact that the band-eigenstates themselves are k-dependent. The latter contributions capture the momentum-space Berry curvature contributions to semi-classical wavepacket dynamics, but now in a formalism that accounts consistently for the scattering terms that are necessary to establish the transport steady state. The covariant derivative involving H0 in Eq. (27) is defined in the same way, simply replacing hρi by H0 in Eq. (28). For our formulation of magnetotransport, it will be important that DB is linear in the density matrix. As we see below the k-dependence of the eigenstates in this derivative play the essential role in capturing the chiral anomaly. The covariant derivatives reduce to simple derivatives in a spin-independent single-band system, for example in a parabolic band system with H0 (k) = ~2 k2 /2m. In the same limit, Eq. (25) reduces to the usual semiclassical Boltzmann equation:   ∂ e + vk · ∇ − (E + vk × B) · ∇k fk = −I(fk ), ∂t ~ (29)

6 where we have defined the velocity vk = (1/~)∇k H0 (k) and I(fk ) is the single-band version of Eq. (15). The quantum kinetic equation (25) we have derived can therefore be understood as a generalization of the simple Boltzmann equation (29) in which the velocity and distribution function scalars are replaced by matrices, the simple derivatives ∇k are replaced by covariant derivatives D/Dk, and scalar products are replaced by symmetrized matrix products 21 { · }. Equation (25) is the principal result of this paper. III.

MAGNETOTRANSPORT THEORY

In this section, we describe a general scheme to apply the quantum kinetic equation (25) to magnetotransport phenomena in systems in which the momentum-space Berry connection plays an important role. We start by obtaining the expression for the density matrix in linear response to an electric field in the absence of a magnetic field. Then we incorporate the effect of a magnetic field by performing a systematic low-magnetic-field expansion. Throughout this paper, we assume that the system is uniform, i.e., that ∇hρi = 0. A.

Expression for hρi without a magnetic field

We first briefly summarize the B = 0 limit of our transport theory. In linear response, we write the electron density matrix hρi as hρi = hρ0 i + hρE i, where hρ0 i is the equilibrium density matrix and hρE i is the correction to hρ0 i in linear order in an electric field E. With this notation we need to solve the kinetic equation in the form ∂hρE i i + [H0 , hρE i] + K(hρE i) = DE (hρ0 i). ∂t ~

(30)

We divide the electron density matrix response hρE i into a diagonal part hnE i and an off-diagonal part hSE i, writing hρE i = hnE i + hSE i. Note that the equilibrium density matrix hρ0 i is diagonal in the band index. When only band-diagonal to band-diagonal terms are included in the scattering kernel, it is easy to solve for the steady-state value of hnE i. The kinetic equation (30) in this limit is [I(hnE i)]mm = [DE (hρ0 i)]mm = eE · vkm k k

∂f0 (εm k) , (31) ∂εm k

mm where vkm = (1/~)∇k εm = f0 (εm k and hρ0 i k ) is the Fermi-Dirac distribution function. The equation for nE is therefore a familiar linear integral equation which is discussed at greater length below and yields hnE im k = m m m eτtrk E · vkm ∂f0 (εm k )/∂εk , where τtrk is the transport lifetime which is often nearly constant across the Fermi surface. Next we consider the solution for the off-diagonal part of the density matrix hSE i, which is independent of weak

disorder. From Eq. (30) the kinetic equation for hSE i is given by i ∂hSE i ′ + [H0 , hSE i] = DE (hρ0 i) − J(hnE i), ∂t ~

(32)

′ where DE (hρ0 i) is the off-diagonal part of the intrinsic driving term: ′

[DE (hρ0 i)]mm =i k

 ′ eE m′ · Rmm f0 (εm k k ) − f0 (εk ) . (33) ~

The off-diagonal part of the driving term is responsible for the Berry phase contribution to the Hall conductivity of systems with broken time-reversal symmetry in the absence of a magnetic field. As we will see, It also plays an essential role in the anomalous magnetoconductivity response. The solution to this equation is [61] Z ∞ ′ ′ hSE i = dt′ e−iH0 t /~ [DE (hρ0 i) − J(hnE i)]eiH0 t /~ , 0

(34)

where we have not explicitly exhibited the time dependences of hρ0 (t − t′ )i and hnE (t − t′ )i. It can be further expanded′ in the eigenstate basis by inserting an infinitesimal e−ηt and taking the limit η → 0 to obtain ′

′ hSE imm k



[DE (hρ0 i)]mm − [J(hnE i)]mm k k = −i~ . ′ m εk − εm k

(35)

Here we have written only the principal value part and omitted δ-function terms which can be important when bands touch, giving rise for example to the Zitterbewegung contribution to the minimum conductivity [61] in graphene. We note that, as shown in Ref. [61], the contribution from J(hnE i) corresponds to the vertex correction in the ladder-diagram approxmation of perturbation theory. B.

Expression for hρi with a magnetic field

Equation (25) can be used to derive a low-magneticfield expansion for the linear response of the densitymatrix, and hence of any single-particle observable, to an electric field E. To explain this we write the steadystate uniform system limit of Eq. (25) in the form (L − DB )hρi = DE (hρ0 i),

(36)

with L ≡ P + K. Here, P hρi ≡

i [H0 , hρi] ~

(37)

is the precession term that accounts for the timedependence of the density matrix in the absence of fields and disorder, and K is the disorder scattering kernel. Note that we have linearized Eq. (36) with respect to the electric field by replacing the density matrix on

7 which DE acts by the equilibrium density matrix hρ0 i. Note also that the matrix P is purely diagonal in both wavevector and in density-matrix element at a given wavevector, and that it is non-zero only for off-diagonal density-matrix elements. We view the three terms P , K, and DB as matrices that act on vectors formed by all eigenstate-representation density-matrix components at all wavevectors, both band diagonal and band offdiagonal, ordered in any convenient way. Given this notation the linear response of the density matrix to an electric field is given in the presence of a magnetic field by hρi = (1 − L−1 DB )−1 L−1 DE (hρ0 i) X = (L−1 DB )N L−1 DE (hρ0 i).

(38)

N ≥0

Note that the N -th term in this expansion is proportional to B N where B is the magnetic field strength, and that the N = 0 term is given by Eqs. (31) and (35). At each order in B, contributions to hρi can quite generally be organized by their order in an expansion in powers of scattering strength λ by letting K → λK and identifying terms with a particular power of λ. The various low-field expansion terms are generated by repeated action of DB and L−1 . In Sec. IV we will take a close look at the properties of L−1 , pointing out that even in systems with weak disorder, the terms that have the highest power of λ−1 do not necessarily dominate. We are now in a position to obtain the magnetoconductivity of a system, i.e., the electrical conductivity in the presence of a magnetic field. With the use of Eq. (38) the total conductivity σµν is obtained from the P N σ definition σµν = Tr[(−e)vµ hρi]/Eν = N ≥0 µν (B ) N (µ, ν = x, y, z). Here, σµν (N ≥ 1) is the magnetoconductivity contribution proportional to B N , which is given by   N σµν = Tr (−e)vµ (L−1 DB )N L−1 DE (hρ0 i) /Eν . (39) In Eq. (39), vµ is the µ component of the velocity operator matrix in the eigenstate representation (i.e., vµ = DH0 /Dkµ ), and Tr indicates both a matrix trace and an integration with respect to k over the Brillouin zone. Note that in Eq. (39) we have allowed for an arbitrary angle between the directions of the electric and magnetic fields. In Sec. VII we will apply Eq. (39) to the positive magnetoconductivity proportional to B 2 , which arises as a consequence of the chiral anomaly, in Dirac and Weyl semimetals. C.

Roles of the magnetic driving term

In closing this section we look more closely at the magnetic driving term (27) by examining linear response to magnetic field in the absence of an electric field. We write the electron density matrix as hρi = hρ0 i + hρB i, where

hρ0 i is the B = 0 density matrix, and hρB i is the correction to hρ0 i at linear order in the magnetic field B. The kinetic equation takes the form ∂hρB i i + [H0 , hρB i] + K(hρB i) = DB (hρ0 i). ∂t ~

(40)

For concreteness, we set B = (0, 0, Bz ). In this case the magnetic driving term is written as     DH0 Dhρ0 i eBz DH0 Dhρ0 i − , , , DB (hρ0 i) = 2~2 Dky Dkx Dkx Dky (41) where { , } indicates a matrix anticommutator. We will show in Sec. VI C that DB (hρ0 i) describes the generation of valley-dependent electrical currents by a magnetic field, i.e.,the so-called chiral magnetic effect. In equilibrium, currents cancel when summed over valleys. We separate the electron density matrix response hρB i into its diagonal hξB i and off-diagonal hSB i parts: hρB i = hξB i + hSB i. By writing the P Hamiltonian and density m matrix of a system as H0 = m,k εk |m, kihm, k| and P m hρ0 i = m,k f0 (εm k )|m, kihm, k|, where εk is an energy m eigenvalue and f0 (εk ) is the Fermi-Dirac distribution function, it can be shown quite generally that only bandoff-diagonal components in DB (hρ0 i) are non-zero, and hence that for weak disorder that K(hρB i) = 0. See Appendix A for a detailed discussion. This means that unlike an electric field [Eq. (31)], a magnetic field does not generate a dissipative current when applied to a conductor, as we know. However both fields do alter the band off-diagonal part of the density matrix. We find that ′

′ hSB imm k

[DB (hρ0 i)]mm k , = −i~ m′ εm k − εk

(42)

where m 6= m′ . See Eq. (35) for the corresponding electric field equation. As mentioned above, [DB (hρ0 i)]mm = 0. However, it k can also be shown quite generally that the band-diagonal components of P −1 DB (hρ0 i) are not zero, i.e. that [P −1 DB (hρ0 i)]mm 6= 0. See Appendix A for a detailed k derivation. We find that in the weak disorder (W τtr ≫ ~ where W is the scale of the typical energetic separation between Bloch bands) limit hξB imm ≡ [P −1 DB (hρ0 i)]mm = k k

e m f0 (εm k ) B · Ωk , (43) ~

m abc where Ωm ih∂kb um k,a = ǫ k |∂kc uk i is the Berry curvature of band m. Note that the magnetic-field induced change in the diagonal density matrix hξB i is quite distinct in character from the changes induced by an electric field hnE i (31) because hξB i is intrinsic (independent of disorder effects) and not related to dissipation. Eq. (43) captures the same physics as, and can be viewed as an alternate derivation of, the Berry phase correction to the

8 density of states implied by semiclassical wavepacket dynamics, (2π)−D [1 + (e/~)B · Ωm ] [63]. In our quantum kinetic formalism, however, it is the the electron density matrix that is modified by a magnetic field, while the density of states remains unchanged. For example, from Eq. (43) we can immediately obtain the Strˇeda formula for the quantum Hall effect with a magnetic field along the z direction [64]: Z e2 X dD k ∂ m Tr[hρB i] = − f0 (εm σxy = −e k )Ωk,z , ∂Bz ~ m (2π)D (44) Finally we note that the same result is obtained when D on any band diagonal density-matrix, hni = B P acts m The calculation details do not m,k gk |m, kihm, k|. appeal to special properties of the equilibrium densitymatrix. In the general case we find that hξB imm = [P −1 DB (hni)]mm = k k

e m g B · Ωm k. ~ k

(45)

This property will play an important role in evaluating contributions to the magnetoconductivity (39). IV.

FERMI SURFACE RESPONSE IN MULTI-VALLEY SYSTEMS

Anomalous transport is often related to band crossings at energies that are close to the Fermi energy and give rise to large Berry curvatures. These are often required by symmetry to have replicas at different Brillouin-zone points, for example at points related by time-reversal operations. For that reason, anomalous transport properties of the type that the present formalism is intended to describe often occur in systems with more than one Fermi surface pocket (In what follows, a pocket refers to a closed Fermi surface segment.) We will follow the common practice in the literature on Weyl and Dirac semimetals by borrowing terminology from semiconductor physics and referring to these clearly separated regions of momentum space that contribute importantly to the physical properties of interest as valleys. It is common that scattering between valleys related to each other by time-reversal, or by some approximate crystal symmetry, is often weaker than scattering within valleys, giving rise to long intervalley relaxation times. In this section, we explain when these long relaxation times are physically relevant. A.

Properties of the operator L−1

The scattering term K(hρi) [Eq. (14)] may be viewed as a linear operator that acts on the density matrix and is the sum of four terms categorized by the following separation: processes that map diagonal hρi to diagonal ∂hρi/∂t (defined as K dd), processes that map diagonal hρi to offdiagonal ∂hρi/∂t (K od ), processes that map off-diagonal

hρi to diagonal ∂hρi/∂t (K do ) and processes that map offdiagonal hρi to off-diagonal ∂hρi/∂t (K oo ). For example, with this notation we can rewrite I(hni) and J(hni) as X ′ [I(hni)]mm = K dd [mk; m′ k′ ]nm k k′ , m′ ,k′

′ [J(hni)]mm k

=

X

K od[mm′ k; m′′ k′ ]nm k′

′′

(46)

m′′ ,k′

with m 6= m′ . Explicit expressions for K dd and K od can be read off from Eqs. (15) and (16). Those for K do and K oo can be found in Ref. [61]. When block-diagonalized to separate band-diagonal and band-off-diagonal components of the density matrix, with this notation the Liouville operator  dd  K K do L=P +K = . (47) K od P + K oo . The disorder (K) contributions to L are proportional to the disorder strength parameter λ, whereas P is intrinsic and independent of disorder strength. Our expressions for magnetoconductivity are formulated in terms of repeated action of the operator  −1 dd  ((L ) (L−1 )do L−1 = . (48) (L−1 )od (L−1 )oo . The leading terms in the disorder strength expansion of L are ∼ λ−1 . As discussed in Ref. [61] it follows that up to the order of λ0   (K dd )−1 −(K dd )−1 K do P −1 L−1 = . (49) −P −1 K od (K dd )−1 P −1 Note that the dd block in L−1 is proportional to disorder strength λ−1 and diverges in the limit of weak disorder, while all other blocks are proportional to λ0 . We see from Eqs. (35) and (42) that the matrix P −1 is diagonal in the density-matrix-component wavevector space. When it acts on a driving term D, it simply multiplies it by a numerical factor: ′

P −1 D = −i~

Dkmm ′ . m εk − εm k

(50)

A similar expression applies when P −1 acts on any Liouville-space density-matrix vector with non-zero offdagonal components. We conclude from these considerations that the leading order term in the disorder strength expansion of the contribution to the response of the diagonal part of the density matrix hnN i at N -th order in magnetic field strength B is obtained by setting L−1 in Eq. (38) to  dd −1  (K ) 0 −1 L → . (51) 0 0 Then it follows that d hnN i = (K dd )−1 DB (hρN −1 i),

(52)

9 d where DB is the diagonal part of the magnetic driving term DB , and hρN −1 i = (L−1 DB )N −1 L−1 DE (hρ0 i) is the density matrix at (N − 1)-th order in magnetic field. On the other hand, the leading order term in the disorder strength expansion of the contribution to the response of the off-diagonal part of the density matrix hSN i is obtained by setting   0 0 L−1 → . (53) −P −1 K od (K dd )−1 P −1

Then it follows that o hSN i =P −1 DB (hρN −1 i)

d (hρN −1 i) − P −1 K od (K dd )−1 DB   −1 o od =P DB (hρN −1 i) − K hnN i ,

(54)

o where DB is the off-diagonal part of the magnetic driving term DB . Notice from Eq. (46) that K od hnN i = J(hnN i). Therefore, we see that Eq. (54) is similar to o the B = 0 equation hSE i = P −1 [DE (hρ0 i) − K od hnE i] [see Eq. (35)]. Here, recall that, as shown in Ref. [61], the contribution from J(hnE i) = K od hnE i corresponds to the vertex correction in the ladder-diagram approximation of perturbation theory. The same argument is applied to the case of magnetotransport. Thus, we conclude that the contribution from K od hnN i to hSN i corresponds to a vertex correction in the ladder-diagram approximation of perturbative magnetotransport calculation. Finally, we note that hSN i is accompanied by a diagonal density matrix o hξN i = P −1 DB (hnN −1 i).

(55)

This is the Berry phase correction to the diagonal part of the density matrix, Eq. (45). hnN i will be dominant in the case where both hnN i and hξN i are present, since (K dd )−1 is of the order of λ−1 while P −1 is of the order of λ0 . However, as will be shown in Sec. VII A, hξN i can make a large contribution to the magnetoconductivity in systems with large momentum-space Berry curvatures.

B.

Scattering times and the total number of electrons in a given valley

We have delayed to this point a discussion of some important properties of K dd and (K dd )−1 that are related to particle number conservation and play an important role in anomalous magnetotransport. The components of these matrices are labelled by band and wavevector labels. Comparing with Eq. (15) we see that  X 2π mn 2 n K dd [mk; m′ k′ ] = δmm′ δkk′ h|Ukp | iδ(εm k − εp ) ~ n,p  mm′ 2 m m′ − h|Ukk′ | iδ(εk − εk′ ) . (56)

P It follows that m′ k′ K dd [mk; m′ k′ ] = 0, i.e.,that the response vector corresponding to total particle number is an eigenvector of K dd with eigenvalue 0. Scattering does not relax total particle number. Since the leading diagd onal response hnE i is given as hnE i = (K dd )−1 DE with d DE the diagonal part of DE [see Eq. (31)], it might seem that the response is divergent. However it is easy to see that the total particle number component of the driving d term, i.e.,the driving term DE summed over all band and wavevector labels is also zero. In fact a stronger statement is valid, namely that the sum of the driving term vanishes when summed over all wavevectors associated with each individual Fermi surface pocket surrounding each relevant valley: X eE ∂f0 (εm ) ∂Nm X k = · = 0, [DE (hρ0 i)]mm = k ∂t ~ ∂k k k (57) where Nm is the total particle number in band m of a given valley. The electric field moves electrons through momentum space; it does not create or destroy particles. In the matrix equations we have been using to express response to electric fields, the total particle number component of the band-diagonal response has implicitly been projected out. The eigenvalues of (K dd )−1 have units of time and characterize the various time scales on which contributions to non-equilibrium populations on the Fermi surface relax. For example, the time scale for relaxation of the non-equilibrium distribution induced by an electric field in a single-valley system, the transport lifetime τtr , is given approximately by P mm 2 m m 1 2π k,k′ h|Ukk′ | i(1 − cos θkk′ )δ(µ − εk )δ(µ − εk′ ) P = m m τtr ~ k δ(µ − εk ) (58)

where θkk′ is the angle between the Bloch state group velocities in valley m at k and k′ (i.e., cos θkk′ = k · k′ /|k||k′ |). In multi-valley systems with weak intervalley scattering it is clear that deviations from equilibrium in the total number of electrons in a valley will relax especially slowly. As we have explained above, these are not driven directly by an electric field, but we will show later that they can be driven by the combined action of electric and magnetic fields. In order to precisely define intervalley relaxation times explicitly it is necessary to separate the Bloch state Hilbert space into valleys in some precise way. The best way to do this depends on the system being studied. For the sake of definiteness, below we use band labels m to distinguish valleys. When K dd is expressed in the representation in which it is diagonal in the absence of intervalley scattering, and intervalley scattering is parametrically weaker than intravalley scattering, its smallest eigenvalues (longest relaxation times) are associated with an M × M block with entries X 1 1 dd − mm′ (59) Kmm ′ = δmm′ mm′′ τ τ inter m′′ inter

10 where (m, m′ ) = 1, . . . , M and 1 ′

mm τinter

2π = ~

P

k,k′





mm 2 m m h|Ukk ′ | iδ(µ − εk )δ(µ − εk′ ) P . (60) m k δ(µ − εk )

Here M is the number of Fermi surface pockets, and for simplicity we have assumed that all pockets have the same Fermi-level density of states because they are related by symmetry. Note that this total valley population block of the scattering kernel still has one zero eigenvalue, corresponding to total population summed over all valleys. In the M = 2 case the non-zero eigenvalue has value 1,2 2/τinter , which is the time scale for relaxing differences in population between the two valleys. We are now in a position to explain when the magnetoresponse of semimetals with atomically smooth disorder is anomalous, and when it is not, by making the following observations: i) When the Fermi surface of a conductor consists of a few small valleys that are well separated in momentum space and disorder is smooth, scattering between valleys is parametrically weaker because it requires large momentum transfers; ii) When intervalley scattering is negligible, the diagonal-component response function (K dd )−1 is divergent not only in the total particle number channel, but also in the valleyprojected particle number channel; iii) because the electric driving term vanishes separately for each valley, very weak intervalley scattering does not lead to anomalous response at B = 0; but iv) in Weyl and Dirac semimetals DB L−1 DE (ρ0 ) does drive the valley-projected total particle number. The total particle number in a given valley is not conserved when the magnetic driving term summed over a given valley is not zero when it acts on the electric-field disturbed density matrix: X ∂N mm = [DB (hρE i)]k 6= 0. ∂t

(61)

m,k

Here, hρE i = L−1 DE (hρ0 i) represents a density matrix linear in electric field. As a consequence of the total particle number non-conservation in a given valley, the intervalley scattering time τinter appears as an eigenvalue of (K dd )−1 . So far we have focused on systems where valleys are not degenerate, i.e., are separated in momentum space, such as Weyl semimetals. Similar considerations apply to systems where valleys are degenerate, such as Dirac semimetals. 3D Dirac semimetals such as Cd3 As2 and Na3 Bi have two Dirac points which are protected by crystalline symmetry. The effective Hamiltonian around a Dirac point in such Dirac semimetals can be written in the block-diagonal form [17]     HAA HAB H+ (q) 0 , (62) HDirac (q) = = 0 H− (q) HBA HBB where A and B denote two states after a unitary transformation, and H± (q) is a 2 × 2 Weyl Hamiltonian with

chirality ±1. As is seen from Eq. (62), two Weyl Hamiltonians that belong to different states A and B are degenerate around the Dirac point. In other words, two different valleys are degenerate. In general, atomically smooth disorder does not allow scattering processes between two different states A and B. Namely, intervalley scattering is much weaker than intravalley scattering, i.e., we have τinter ≫ τintra . Thus we can apply the same argument as in the case of non-degenerate valleys to the case of degenerate valleys: the total particle number in a given valley is not conserved if the magnetic driving term summed over a given valley is not zero [see Eq. (61)]. We will demonstrate in Sec. VII that the particle number non-conservation in a given valley characterized by Eq. (61) occurs in the process of calculating the magnetoconductivity proportional to B 2 in Dirac and Weyl semimetals, which is identified as a consequence of the chiral anomaly.

V. CHIRAL ANOMALY IN THREE-DIMENSIONAL METALS

In Weyl semimetals the term “chiral anomaly” is sometimes used to refer to the property that the total number of electrons in a given valley is not conserved in the presence of simultaneous electric and magnetic fields, and sometimes to the enhanced magnetoconductivity that this lack of separate particle number conservation produces. In this section we focus on pumping of charge between valleys without referring to a specific microscopic model, which leads to observable effects only when intervalley scattering times are much longer than intravalley scattering times. For long intervalley scattering times, anomalous pumping leads to a difference between the effective chemical potentials of different valleys and thus to a enhanced current via the magnetoelectric effect, as explained in great deal for a Weyl metal toy model in the following sections. a general model with the Hamiltonian H0 = PWe study m ε |m, kihm, k| and the equilibrium density mam,k k P m m trix hρ0 i = m,k f0 (εk )|m, kihm, k|, where εk is an energy eigenvalue and f0 (εm k ) is the Fermi-Dirac distribution function. For concreteness, we consider the case of E = (0, 0, Ez ) and B = (0, 0, Bz ). We consider the quantity [Eq. (61)]   ∂N = Tr DB L−1 DE (hρ0 i) ∂t

(63)

evaluated for a particular Fermi surface pocket associated with a particular valley. When L−1 acts on ∂N/∂t it yields the contribution to the density matrix that is linear in both electric and magnetic fields. It follows that ∂N/∂t is the rate of change of the density matrix that must be compensated by scattering and free evolution in order to yield a steady-state density matrix. As we have explained earlier, hρE i = L−1 DE (hρ0 i), the

11 linear-response density matrix in the absence of a magnetic field, contains both band-diagonal hnE i and band off-diagonal hSE i contributions, as hρE i = hnE i + hSE i. Then we have Tr [DB L−1 DE (hρ0 i)] = Tr [DB (hnE i)] + Tr [DB (hSE i)], since DB is linear in the density matrix. mm First we evaluate Pthe diagonal element [DB (hnE i)]k )]|m, kihm, k|. In this with hnE i = eEz m,k [∂kz f0 (εm k case we can use the result of [DB (hρ0 i)]mm by replacing k f0 by eEz ∂kz f0 . Namely, from Eq. (A6) we get [DB (hnE i)]mm k  m  ∂εk ∂ ∂εm ∂f0 (εm 2 k ∂ k) = e Ez Bz − . ∂ky ∂kx ∂kx ∂ky ∂kz

(64)

It is obvious that this is an odd function of kx , ky , and kz , which means that Tr [DB (hnE i)] = 0. Then we see that the non-conservation of the total electron number in a given valley induced by the chiral anomaly is associated with the action of the DB operator on the off-diagonal density matrix hSE i: Z ∂N d3 k X = [DB (hSE i)]mm , k 3 ∂t CA FS (2π) m

(65)

where FS represents the integration on the Fermi surface of the valley. After a calculation (see Appendix B for a detailed calculation), we find a general expression that represents the Fermi surface contribution X m

mm

[DB (hSE i)]k 2

= e Ez Bz

X  ∂f0 (εm ) k

m

∂kx

Ωm k,x

 ∂f0 (εm k) m + Ωk,y , ∂ky

(66)

m abc where Ωm ih∂kb um k,a = ǫ k |∂kc uk i is the Berry curvature of band m. Then Eq. (65) is rewritten as

Z 3   m m ∂N d k ∂f0 (εm e2 Ez Bz m k) vk,x Ωk,x + vk,y Ωm = k,y , ∂t CA 4π 2 2π ∂εm k (67) where vkm is the Bloch state group velocity, and we have assumed that only the band m possesses the Fermi surm face, i.e., ∂f0 (εnk )/∂εnk = δmn ∂f0 (εm k )/∂εk , which can in general apply to multi-valley systems. Equation (67) is in sharp contrast to the result obtained in free semiclassical dynamics, which shows the number of electrons in a given valley changes at the rate [35, 51] e2 Ez Bz ∂N = ∂t 4π 2

Z

d3 k ∂f0 (εm k) m vk · Ω m k. 2π ∂εm k

(68)

m The term proportional to vk,z in Eq. (68) is absent in Eq. (67). This result can be interpreted as follows: free semiclassical dynamics are disturbed by intravalley scattering of electrons, i.e., those electrons are swept across the valley m Fermi surface by the electric field.

VI.

CHIRAL ANOMALY IN WEYL SEMIMETALS

So far we have formulated a general theory of lowfield magnetotransport that accounts for the wavevectordependence of Bloch states in multiband systems. In this section, we illustrate the power of our theory by showing that it correctly captures the physics of the chiral anomaly in which large momentum-space Berry curvatures play essential roles. We show that the electric driving term (26) and the magnetic driving term (27) properly describe electron transport phenomena that arise from topologically nontrivial electronic band structures: the anomalous Hall effect and chiral magnetic effect, respectively, in Weyl semimetals. We set ~ = 1 in the rest of this paper.

A.

Theoretical model

We start with the 4 × 4 continuum toy model Hamiltonian for two-node Weyl semimetals with broken timereversal symmetry [39, 42, 65, 66] H0 (k) = vF τz k · σ + ∆τx + bσz ,

(69)

where the Pauli matrices τi and σi are Weyl-node and spin degrees of freedom, respectively, vF is the Dirac velocity, and ∆ is the mass of 3D Dirac fermions. The term bσz can be regarded as a magnetic interaction of the z component of spin such as s-d coupling and Zeeman coupling. Therefore, it can be said that the above Hamiltonian represents a 3D (topological or normal) insulator doped with magnetic impurities. As will be shown just below, a Weyl semimetal is realized when |b/∆| > 1. The time-reversal T and spatial inversion (parity) P operators are given by T = −iσy C with C the complex conjugation operator and P = τx . Here, we have the following relations: T −1 = iσy C, P −1 = τx , T 2 = −1, and P 2 = 1. Then it is easily seen that time-reversal symmetry of the system is broken such that T H0 (k)T −1 6= H0 (−k), but inversion symmetry of the system is preserved such that PH0 (k)P −1 = H0 (−k). Performing a canonical transformation such that σx,y → τz σx,y and τx,y → σz τx,y , Eq. (69) can be rewrit˜ 0 (k) = vF (kx σx + ky σy ) + (b + vF τz kz + ∆τx )σz . ten as H In the representation of eigenstates of b + vF τz kz + ∆τx , Eq. (69) is block-diagonal with two 2 × 2 Hamiltonians given by [39] H± (k) = vF (kx σx + ky σy ) + m± (kz )σz

(70)

p with m± (kz ) = b ± vF2 kz2 + ∆2 . The two √ Weyl nodes are located at W± = (0, 0, ±k0 ) with k0 = b2 − ∆2 /vF . Near the Weyl nodes, we have m− (qz ) ≈ ∓vF2 (k0 /b)qz with momentum q = k − W± (q 2 ≪ 1). The eigenvectors of Ht (k) (t = ±) with eigenvalues ε± tk = ±εtk =

Energy

12

FIG. 1. Schematic figure of a two-node Weyl semimetal. Q± = ±1 is the chirality of a Weyl node. 2b and 2µ5 are the the momentum-space distance and the chemical potential difference between the Weyl nodes, respectively. The dashed line indicates the Fermi level of the system.

±

q vF2 (kx2 + ky2 ) + m2t are given by

 q  mt (kz ) 1 ± 1  εtk , q |u± tk i = √ 2 ±eiθ 1 ∓ mt (kz ) εtk

where eiθ = (kx + iky )/k⊥ with k⊥ =

(71)

q kx2 + ky2 . We see

that H− (k) describes a subsystem with two Weyl nodes, while H+ (k) has a fully gapped spectrum. For convenience, we omit the subscripts t = ± at intermediate steps of the calculations below and restore them in the final results. The generalized Berry connection in the eigenstate repn resentation is given by [Rk,α ]mn = ihum k |∂kα uk i with α = x, y, z and m, n = ±. The individual components are given explicitly by 1 m vF m 1 sin θ − σ ˜z sin θ − σ ˜y 2 cos θ 2k⊥ 2k⊥ εk 2εk vF −σ ˜x sin θ, 2εk 1 m vF m 1 cos θ + σ ˜z cos θ − σ ˜y 2 sin θ =− 2k⊥ 2k⊥ εk 2εk vF +σ ˜x cos θ, 2εk vF k⊥ ∂m =˜ σy , (72) 2ε2k ∂kz

Rk,x =

Rk,y

Rk,z

where σ ˜α are the Pauli matrices in the eigenstate basis. Also, the individual components of the Berry curvature, ± abc Ω± ih∂kb u± k,a = ǫ k |∂kc uk i, are given by Ω± k,x = ∓

∂m vF2 ky vF2 m ∂m vF3 kx ± ± , Ω = ∓ , Ω = ∓ . k,y k,z ∂kz 2ε2k ∂kz 2ε3k 2ε3k (73)

Here, we briefly review the derivation of the θ term from the microscopic four-band model (69). In order to

describe a more generic Weyl semimetal, we add the term µ5 τz to the Hamiltonian, which generates a chemical potential difference 2µ5 between the two Weyl nodes, and replace bσz by b · σ. Note that this term breaks inversion symmetry. Figure 1 shows a schematic diagram of a two-node Weyl semimetal. We set ∆ = 0 for simplicity. The action of the system in the presence of an external electromagnetic potential Aµ = (A0 , −A) is given by Z S = dt d3 r ψ † {i(∂t − ieA0 ) − [H0 (k + eA) − µ5 τz ]} ψ Z ¯ µ (∂µ − ieAµ − ibµ γ 5 )ψ, = dt d3 r ψiγ (74) where e > 0, ψ is a four-component spinor, γ¯ = ψ † γ 0 , γ 0 = τx , γ j = τx τz σj = −iτy σj , γ 5 = iγ 0 γ 1 γ 2 γ 3 = τz , and bµ = (µ5 , −b). After applying the Fujikawa method [67], i.e., applying an infinitesimal gauge transformation 5 such that ψ → e−idφθ(r,t)γ /2 ψ with θ(r, t) = −2xµ bµ = 2(b · r − µ5 t) and φ ∈ [0, 1] for infinite times, the action of the system becomes [34, 39] Z ¯ µ (∂µ − ieAµ )]ψ + Sθ , S = dt d3 r ψ[iγ (75) where Sθ is the θ term [Eq. (1)] with θ(r, t) = 2(b · r − µ5 t). Note that the first term in Eq. (75) represents the (trivial) action of massless Dirac fermions. In general, regardless of whether the system is gapless or gapped, the four-current density j ν can be obtained from the variation of the θ term with respect to the four-potential Aν as j ν = δSθ /δAν = −(e2 /4π 2 )ǫµνρλ [∂µ θ(r, t)]∂ρ Aλ [28]. The induced current density in the presence of external electric and magnetic fields is given by j(r, t) =

i e2 h ˙ t)B , ∇θ(r, t) × E + θ(r, 4π 2

(76)

e2 (b × E − µ5 B) , 2π 2

(77)

where θ˙ = ∂θ(r, t)/∂t, E = −∇A0 − ∂A/∂t, and B = ∇ × A. In the present case of Weyl semimetals with θ(r, t) = 2(b · r − µ5 t), we obtain in the ground state a static current of the form j=

where the electric-field induced and magnetic-field induced terms are the anomalous Hall effect and chiral magnetic effect, respectively [34–39, 41, 42, 49]. The anomalous Hall effect in Weyl semimetals can be understood straightforwardly. It is well known that the Hall conductivity of 2D massive Dirac fermions (such as those on the surface of 3D topological insulators) of the ± (kz ) = sgn[m± (kz )]e2 /2h, when form (70) is given by σxy the chemical potential lies in the energy gaps. This (half)quantized value holds valid even in the presence of disorder [31, 68]. We see that the 3D Hall conductivity is nonzero only in the region −k0 ≤ kz ≤ k0 , which gives Rk 3D ± σxy = −k0 0 (kz /2π)σxy (kz ) = k0 e2 /πh. This value is

13 exactly the same as the first term in Eq. (77). On the other hand, the chiral magnetic effect in Weyl semimetals looks like a peculiar phenomenon. The chiral magnetic effect indicates a direct current generation along a static magnetic field (without electric fields) in the presence of a chemical potential difference between the two nodes, as is seen from the second term in Eq. (77). If the static chiral magnetic effect exists in real materials, there will be substantial possible applications. Discussions on the existence of the static chiral magnetic effect in Weyl semimetals have continued theoretically [42–48]. However, the existence of the static chiral magnetic effect would be ruled out in crystalline solids (i.e., lattice systems) as discussed in Ref. [42], which is consistent with our understanding that static magnetic fields do not generate equilibrium currents. Lastly, we refer to these two phenomena from the symmetry viewpoint. The term b · σ breaks time-reversal symmetry but preserves inversion symmetry: T σi T −1 = −σi and Pσi P −1 = σi (i = x, y, z). This indicates that the occurrence of the anomalous Hall effect in Weyl semimetals requires the breaking of time-reversal symmetry, which is consistent with the generic requirement for the realization of the anomalous Hall effect. On the other hand, the term µ5 τz breaks inversion symmetry but preserves time-reversal symmetry: T τz T −1 = τz and Pτz P −1 = −τz . This indicates that the occurrence of the chiral magnetic effect in Weyl semimetals requires the breaking of inversion symmetry. Note that both the anomalous and chiral magnetic effects are possible (at least theoretically) in systems with broken time-reversal and inversion symmetries.

B.

±

(εk −µ)/T where f0 (ε± + 1] is the Fermi-Dirac disk ) = 1/[e tribution function. The velocity operator is given by vx = ∂Ht /∂kx = vF σx . We transform to the eigenstate n basis, in which [vx ]mn = vF hum k |σx |uk i (m, n = ±) is given by   vF k⊥ m ˜z vx = vF σ cos θ + σ ˜y sin θ − σ ˜x cos θ . (80) εk εk I The anomalous Hall conductivity σxy is calculated from I the definition σxy = Tr [(−e)vx hSE i] /Ey . Combining the above ingredients, we obtain X Z ∞ d3 k mt   − I σxy = 2e2 f0 (ε+ tk ) − f0 (εtk ) 3 3 (2π) 4εtk t=± −∞ # " Z X mt e2 ∞ d3 k m− + f0 (ε−k ) − = 2 −∞ (2π)3 ε3−k ε3 t=± tk Z Z e2 δ d3 q kb0 qz − kb0 qz e2 X ∞ mt ≈ − dkz 2 −δ (2π)3 ε3−q 4π 2 t=± 0 |mt |

e2 k0 , (81) 2π 2 where we have used the fact that m− (kz ) ≈ ∓vF2 (k0 /b)qz around the Weyl nodes with q = k − W± (q 2 ≪ 1). Then we see that the contribution from ε+ tk always vanishes as long as the chemical potential µ lies sufficiently close to the Weyl nodes. Also, it should be noted that the quantity Ωtk,z = mt /2ε3tk is the z component of the Berry curvature of the two-band Hamiltonian (70). The anomalous Hall conductivity (81) is exactly the same as the value obtained from a field-theoretical approach, the first term in Eq. (77). Next, we need to evaluate the extrinsic contribution II σxy , which originates from the Fermi surface and the presence of disorder. As described in Sec. II, we treat disorder within the Born approximation. We consider a short-range P (on-site) disorder potential of the form U (r) = U0 i δ(r − ri ), and assume that the correlation function satisfies hU (r)U (r ′ )i = nimp U02 δ(r − r ′ ) with nimp the impurity density. Then we obtain =−

Anomalous Hall effect in Weyl metals

Let us see that the electric driving term (26) properly describes the anomalous Hall effect in Weyl metals. As described in Sec. III, the diagonal part of the electric driving term DE results in usual Drude conductivity. Hence, we need the explicit matrix expression for the off-diagonal part of the electron density hSE i to obtain the anomalous Hall conductivity. As is seen from Eq. (35), there are contributions from the electric driving term DE (hρ0 i) and the anomalous driving term J(hnE i) to the off-diagonal part of the density matrix hSE i. Let us denote the resultant total anomalous Hall conductivity of the system as I II σxy = σxy + σxy ,

I First, we evaluate the intrinsic contribution σxy , which originates from the Fermi sea of the electronic band. Using the expression for DE (hρ0 i) (33) and the Berry connection (72), the off-diagonal part of the electron density matrix hSE i is given by  −  f0 (ε+ 0 R+− k ) − f0 (εk ) k,y hSE i = eEy R−+ 0 2εk k,y   + − f0 (εk ) − f0 (εk ) m σ ˜ cos θ − σ ˜ sin θ , v = eEy x y F 4ε2k εk (79)

(78)

I II where σxy is the Hall conductivity from DE (hρ0 i) and σxy is that from J(hnE i). In the following, we consider the case where an electric field is applied along the y direction as E = Ey ey , and set µ5 = 0. Also, we add a positive small chemical potential µ which lies sufficiently close to the Weyl nodes. Note that the sign of µ is not essential, since the Hamiltonian (70) has particle-hole symmetry.

++ +− hUkk ′ U k′ k i

=

   ′  nimp U02 k⊥ m(kz′ ) m(kz ) k⊥ , + i sin γ − cos γ 2 ε k ε k′ εk ε k′ (82)

14 ′

′ where γ = θ′ − θ with eiθ = (kx′ + iky′ )/k⊥ . We also have +− −− ++ +− −− −+ the relations hUkk′ Uk′ k i = −hUkk′ Uk′ k i, hUkk ′ U k′ k i = ++ +− ∗ −+ ++ ++ +− ∗ −hUkk′ Uk′ k i , and hUkk′ Uk′ k i = hUkk′ Uk′ k i . From the expression for J(hnE i) [Eq. (16)] with hnE i = − diag[n+ Ek , nEk ], we get X  + ++ +− + + [J(hnE i)]+− =π hUkk ′ Uk′ k i nEk δ(εk − εk′ ) k

where Qν and µν are the chirality and chiral chemical potential of a Weyl node, respectively. In two-node Weyl semimetals, we have Q± = ±1 and 2µ5 = µ+ − µ− (see Fig. 1). As is obvious, the chiral magnetic effect does not occur when there is no chemical potential difference between the nodes as µ+ = µ− .

+ 2 where n+ Ek = eτ+ Ey ∂f0 (εk )/∂ky with τ+ ∝ 1/nimp U0 + the scattering rime for the band εk . Here, we have − − + used δ(ε+ k − εk′ ) = δ(εk − εk′ ) = 0 and the fact that − − nEk = nEk′ = 0, because we have considered the case of µ > 0. Similarly, we obtain [J(hnE i)]−+ = k ∗ {[J(hnE i)]+− } . Then the off-diagonal density matrix k resulting from J(hnE i) is given by   Im [J(hnE i)]+− Re [J(hnE i)]+− k k ′ −σ ˜x , hSE i = −˜ σy 2εk 2εk (84)

First, it is informative to check the absence of the chiral magnetic effect in our formulation when µ+ = µ− . Since there is no electric fields, we start from the di− agonal density matrix hρ0 i = diag[f0 (ε+ k ), f0 (εk )] with ± (εk −µ)/T f0 (ε± + 1]. Without loss of generalk ) = 1/[e ity we can consider the case of B = (0, 0, Bz ) and can set µ > 0. Explicit matrix expression for DB (hρ0 i) can be found in Appendix D. Using Eqs. (42) and (D9), the linear-response off-diagonal electron density matrix induced by the magnetic field is obtained as   −  ∂[f0 (ε+ 0 R+− e k ) + f0 (εk )] k,y hSB i = Bz R−+ 0 2 ∂kx k,y   − ∂[f0 (ε+ 0 R+− k ) + f0 (εk )] k,x − . (87) −+ Rk,x 0 ∂ky

k′

 + + −n+ Ek′ δ(εk − εk′ ) ,

where Re and Im represent the real and imaginary part, respectively. The calculation of the conductivity from II ′ the definition σxy = Tr [(−e)vx hSE i] /Ey is somewhat long. See Appendix C for a detailed calculation. Finally, II we find that σxy = 0 when the chemical potential µ lies sufficiently close to the Weyl nodes. In the end, the total anomalous Hall conductivity of the Weyl metal reads I II σxy = σxy + σxy =−

e2 k0 , 2π 2

(85)

which means that the anomalous Hall conductivity of Weyl metals is purely intrinsic, i.e., is determined by the Berry curvature of the filled band, as long as the chemical potential µ lies sufficiently close to the Weyl nodes. This plateau-like behavior of the anomalous Hall conductivity is consistent with a calculation by Burkov [52]. As mentioned in Sec. III A the contribution from J(hni), i.e., II σxy , corresponds to the vertex correction in the ladder approximation. Therefore it can be said that the vertex correction to the anomalous Hall conductivity of Weyl metals described by the Hamiltonian (69) is absent, as long as the chemical potential µ lies sufficiently close to the Weyl nodes. C.

Chiral magnetic effect in Weyl metals

Let us see that the magnetic driving term (27) properly describes the chiral magnetic effect in Weyl metals. The chiral magnetic effect, the second term in Eq. (77), is written as 2

j=

e B 4π 2

X ν

Q ν µν ,

1.

(83)

(86)

The case of µ+ = µ−

By substituting the Berry connections (72) into this equation, we get a concrete expression for hSB i. Also, the intrinsic diagonal density matrix induced by the magnetic field (43) is given by   + f0 (ε+ 0 k )Ωk,z . (88) hξB i = eBz − 0 f0 (ε− k )Ωk,z The velocity operator is written in the eigenstate basis as   ∂m m vF k⊥ vz = (89) σ ˜z + σ ˜x . ∂kz εk εk We are in a position to calculate the electric current. Since the chemical potential µ lies sufficiently close to the Weyl nodes, we do not need to take into account the contribution from the t = + band. Finally, we find that the electric current along the magnetic field B vanishes as expected: jz =Tr[(−e)vz hSB i] + Tr[(−e)vz hξB i] Z ∂f0 (ε+ d3 k vF4 kz 1 X e2 −k ) k = − Bz a 2 2 (2π)3 m− − b ε−k a=x,y ∂ka Z  d3 k vF4 kz m2−  − 2 + e Bz f0 (ε+ −k ) + f0 (ε−k ) 4 3 (2π) m− − b 2ε−k =0,

(90)

where we have used the fact that the integrand is an odd function of kz . This shows that the chiral magnetic effect does not occur when there is no chemical potential difference between the Weyl nodes, i.e., when µ+ = µ− .

15 to be

The case of µ+ 6= µ−

2.

Next, let us consider the case where there is a chemical potential difference between the two nodes as 2µ5 = µ+ − µ− . As is shown above, the expression for the chiral magnetic effect (86) is derived analytically from the four-band model (69) with nonzero µ5 τz term. However, it is difficult to obtain analytically the eigenvalues and eigenstates in the presence of the µ5 τz term, unlike in the case of the bσz term where we can obtain the eigenstates analytically as Eq. (71). To determine whether the magnetic driving term (27) can describe the chiral magnetic effect, we consider two-band Hamiltonians around each Weyl node separately and sum each contribution to the electric current. The Hamiltonian around each Weyl node reads ˜ ν (q) = vF (qx σx + qy σy + Qν qz σz ), H

(91)

q qx2 + qy2 + qz2 . which gives the eigenvalues ε± q = ±vF Here, ν = ± is the node index and Q± = ±1 is the chirality of each Weyl node. The Fermi-Dirac distribution function around each node is given by f0ν (ε± q) = ±

1/[e(εq −˜µν )/T +1] with µ ˜ ν = µ+µν , where µ is a uniform chemical potential and µν is a chiral chemical potential such that µ+ 6= µ− . Without loss of generality, we can set µ ˜ν > 0. In linear response, the diagonal part hξB,ν i and the off-diagonal part hSB,ν i of the density matrix induced by the magnetic field are given respectively by ν ± Eqs. (88) and (87), with f0 (ε± q ) replaced by f0 (εq ). The velocity operator is given by vz,ν = Qν vF



 vF q⊥ Qν vF qz σ ˜z + σ ˜x . εq εq

(92)

In the present case, we may approximate the total electric current by the sum of the contributions from each node. Then the total electric current along the magnetic P S ξ + jz,ν ), where field B is calculated as jz = ν=± (jz,ν X jz,ν ≡ Tr[(−e)vz,ν hXB,ν i] (X = S, ξ). The contribution from hSB,ν i is calculated to be ∂f0ν (ε+ d3 q vF3 X q) q a 3 2 ∂qa δΛν (2π) εq a=x,y Z Z π δqν e2 Bz vF = Q dq q sin3 θδ(q − µ ˜ν /vF ), dθ ν 8π 2 0 0 2 e2 = Bz Q ν µ ˜ν , (93) 3 4π 2

S jz,ν =−

e2 Bz Q ν 2

Z

whereqδΛν is a small momentum space around each node, + q = qx2 + qy2 + qz2 , and we have used ∂f0ν (ε+ q )/∂εq =

−δ(εq − µ ˜ν ). The contribution from hξB,ν i is calculated

ξ jz,ν

 d3 q vF5 qz2  ν + f0 (εq ) + f0ν (ε− q) 3 4 (2π) 2εq Z π "Z µ˜ν Z Λ # 2 e Bz = Qν dθ dq + dq sin θ cos2 θ 8π 2 0 0 0 2

= e Bz Q ν

=

Z

1 e2 Bz Qν (˜ µν + Λ), 3 4π 2

(94)

where Λ is the cutoff of the Fermi sea of Weyl cones. Combining Eqs. (93) and (94), we arrive at the final expression for the chiral magnetic effect: X e2 Bz Q ν µν , (95) 2 4π ν P P ˜ν = where we have used ν=± Qν µν and ν=± Qν µ P Q Λ = 0. This is in full agreement with the result ν=± ν by a field-theoretical approach (86). Here, we comment on an important observation in this study. The portion 2/3 comes from the Fermi surface contribution, as seen from Eq. (93). The portion 1/3 comes from the Fermi sea contribution, as see from Eq. (94). This observation is in contrast to a derivation by semiclassical wavepacket dynamics [51, 54] in which the chiral magnetic effect comes entirely from the Fermi surface contribution. jz =

VII.

NEGATIVE MAGNETORESISTANCE IN WEYL AND DIRAC METALS

We have shown that the electric driving term DE (26) and magnetic driving term DB (27) properly describe the chiral-anomaly-induced phenomena in Weyl metals, the anomalous Hall effect and chiral magnetic effect, respectively. Another phenomenon manifested by the chiral anomaly is a negative magnetoresistance (or equivalently positive magnetoconductance) quadratic in B for parallel electric and magnetic fields in Weyl and Dirac metals [51–54]. Such an unusual negative magnetoresistance has recently been experimentally observed in Dirac semimetals Na3 Bi [55], Cd3 As2 [56, 57], and ZrTe5 [58], and in the Weyl semimetals TaAs [59] and TaP [60]. Here, note that the usual magnetoresistance due to Lorentz force is positive. As denoted in the introduction, R the chiral anomaly is described by the θ term Sθ = dt d3 r (e2 /4π 2 )θ(r, t)E · B. As we can see from this expression for the θ term, the E · B contribution becomes largest in the case of parallel electric and magnetic fields. The positive quadratic magnetoconductivity arising from the chiral anomaly first derived by Son and Spivak reads [51] σzz (Bz2 ) =

e2 (eBz )2 vF3 τ, 4π 2 µ2

(96)

where µ is the chiral potential and τ is the intervalley scattering time. Since their derivation is based on semiclassical wavepacket dynamics, it is difficult to take into

16 confirmed that only the magnetoconductances originating from the chiral anomaly and the usual Lorentz force survive. We divide the calculation of the matrix (L−1 DB )2 L−1 DE (hρ0 i) in Eq. (97) into three steps as follows: ′



n n hSE i = L−1 DE (hρ0 i) = [DE (hρ0 i)]nn k /i(εk − εk ) ∝ eEz ,

hnEB i = L−1 DB (hSE i) = τ [DB (hSE i)]nn k ∝ e2 Ez Bz τ,





hSEB 2 i = L−1 DB (hnEB i) = [DB (hnEB i)]knn /i(εnk − εnk ) FIG. 2. Schematic flow to obtain the longitudinal quadratic magnetoconductivity σzz (Bz2 ) (97) for E k B in Weyl metals described by the Hamiltonian (70). hni and hSi indicate a diagonal density matrix and an off-diagonal density matrix, respectively. τ and τintra are the intervalley and intravalley scattering times, respectively. σzz (Bz2 ) is given by the sum of the CA contribution from the chiral anomaly σzz (Bz2 ) and that from LF 2 2 CA LF the Lorentz force σzz (Bz ) as σzz (Bz ) = σzz (Bz2 ) + σzz (Bz2 ). ′ ′ The contributions from hnEB2 i and hSEB2 i to the conductivity vanish because they are odd functions of kx and ky . The contribution from hn′′EB2 i vanishes because the matrix vz hn′′EB2 i is traceless.

account the detailed band structure of the system. Also, the microscopic origin of the appearance of the intervalley scattering time in Eq. (96) is not clear. In this section, we apply the magnetotransport theory we have developed so far to the positive quadratic magnetoconductivity, starting from a microscopic continuum model of Weyl semimetals (69). It should be mentioned that the positive quadratic magnetoconductivity has been experimentally observed in the low magnetic field regime. Thus, our semiclassical treatment of magnetic fields can be justified. In our theory, the formal expression for the µν-component of the quadratic magnetoconductivity is written from Eq. (39) in the form   σµν (B 2 ) = Tr (−e)vµ (L−1 DB )2 L−1 DE (hρ0 i) /Eν . (97) Here, note that the angle between the electric and magnetic fields is arbitrary in this formalism. Our calculation of Eq. (97) is a little tricky. Our calculation is inspired by (i) that the expression of σzz (Bz2 ) (96) has a linear dependence on τ , and (ii) that this expression is derived via the chiral magnetic effect [51– 54, 58]. We use the fact that a magnetic driving term obtained from an off-diagonal density matrix is purely diagonal, and that a magnetic driving term obtained from a diagonal density matrix has the diagonal and off-diagonal components that lead to the chiral magnetic effect as we have seen in Sec. VI C (see Appendix D for explicit expressions for the magnetic driving term). Indeed, as the result of an explicit calculation shown in Fig. 2, we have

∝ e3 Ez Bz2 τ,

(98)

− ± where hρ0 i = diag[f0 (ε+ k ), f0 (εk )] with f0 (εk ) = (±εk −µ)/T 1/[e + 1] being the Fermi-Dirac distribution function, hni and hSi indicate diagonal and off-diagonal matrices, respectively, and τ is the intervalley scattering time. See Eqs. (52) and (54) for the formal expressions for hni and hSi in a magnetic field. As we shall show below, the appearance of the intervalley scattering time is due to a special property such that the magnetic driving term DB (hSE i) has a nonzero value when integrated over a given valley (Weyl cone). In the following, we consider the low-temperature case where T ≪ µ, with T and µ > 0 being the temperature and chemical potential of the system, respectively.

A.

~ kB ~ Quadratic magnetoconductivity for E 1.

Contribution from the chiral anomaly

Let us consider the case of parallel electric and magnetic fields as E = (0, 0, Ez ) and B = (0, 0, Bz ). We start by obtaining the off-diagonal part of the density matrix induced solely by the electric field, hSE i in Eq. (98). As we have seen in Sec. VI B, this off-diagonal density matrix results in an intrinsic effect, i.e., the anomalous Hall effect. Using the expression for the electric driving term DE (hρ0 i) (33) and the Berry connection (72), we get − hSE i = σ ˜y eEz [f0 (ε+ k ) − f0 (εk )]

vF k⊥ ∂m . 4ε3k ∂kz

(99)

At this stage we recall from Eq. (35) that there exists an extrinsic contribution to the off-diagonal density matrix hSE i from the anomalous driving term J(hnE i). Here, − ± ± hnE i = diag[n+ Ek , nEk ] with nEk = eτintra Ez ∂f0 (εk )/∂kz (τintra is the intravalley scattering time). Let us consider the same setup as in the evaluation of the anomalous Hall conductivity in Sec. VI B. Namely, we consider the case of a short-range (on-site) disorder potential of the form P U (r) = U0 i δ(r − ri ), and assume that the correlation function satisfies hU (r)U (r ′ )i = nimp U02 δ(r − r ′ ) with nimp the impurity density. In this case, we can show

17 explicitly from Eqs. (C2) and (C3) that   Re [J(hnE i)]+− = Im [J(hnE i)]+− = 0, k k

where (100)

and hence [J(hnE i)]+− = [J(hnE i)]−+ = 0, as long as k k the chemical potential µ lies sufficiently close to the Weyl nodes. Hence, the expression for hSE i, Eq. (99), remains valid even after the anomalous driving term is taken into account in the case of sufficiently small µ. Second, we compute the diagonal density matrix hnEB i proportional to Ez Bz in Eq. (98). This hnEB i is the most important ingredient in our derivation of the magnetoconductivity induced by the chiral anomaly, as is understood from its form similar to E · B. As described in Appendix D, the magnetic driving term obtained from an off-diagonal matrix is purely diagonal. Using the expressions for the magnetic driving term DB (hSE i) (D15) and (D16) and the Berry connection (72), we have DB (hSE i) = e2 Ez Bz Fk 1

(101)

with Fk = −

∂ck ∂ck vF m2 vF3 k⊥ ck − vF cos θ − vF sin θ , c − k ε2k ε2k k⊥ ∂kx ∂ky (102)

where 1 is the 2 × 2 identity matrix and ck = [f0 (ε+ k) − 3 f0 (ε− )](v k /4ε )∂m/∂k . Here, we show that this F ⊥ z k k DB (hSE i) has a special property. Around a Weyl node with momentum q = k − W± (q 2 ≪ 1), m− (kz ) is approximated as m− (qz ) ≈ ∓vF2 (k0 /b)qz . Then we find that Fq is an even function of qx , qy , and qz . Accordingly, we find that the integral of DB (hSE i) over the Fermi surface of a given valley has a nonzero R d3 q P mm value: FS (2π) 6= 0. As discussed in 3 m [DB (hSE i)]q Sec. IV, this is a consequence of the total particle number non-conservation in a given Weyl cone. Namely we obtain the rate of pumping Z ∂N e2 Ez Bz d3 q = 2 Fq . (103) ∂t 4π 2 FS 2π In the case of isotropic Weyl fermions [i.e., m− (qz ) ≡ vF qz ], Eq. (102) can be simplified to be  X 1   ∂ − + 2 + iso f0 (ε+ + 3q (Ωq,a ) Ωq,a Fq = q ) − f0 (εq ) , 2 ∂qa a=x,y (104)

Ω+ q,a

3

where = −qa /(2q ) is the Berry curvature and q = q qx2 + qy2 + qz2 . Integrating over the Fermi surface we have Z 3 Z d k ∂f0 (ǫm 2 2 d3 q m k) m 2 Fqiso = m vk · Ωk = 3 . (105) 2π 3 2π ∂ǫ FS k It follows from Eq. (52) that the intervalley scattering time τ appears when L−1 acts on DB (hSE i). In the end, the diagonal density matrix hnEB i is obtained as hnEB i = L−1 DB (hSE i) = e2 Ez Bz F˜k τ 1,

(106)

  ∂f0 (εm 1 X ∂f0 (εm k) m k) m Ωk,x + Ωk,y F˜k = 2 m=± ∂kx ∂ky

(107)

is the component which represents the Fermi surface response in Eq. (102). Note that the component which represents the Fermi sea response in Eq. (102) is irrelevant, since we are focusing on the scattering on the Fermi surface. Third, we compute the off-diagonal density matrix hSEB 2 i proportional to Ez Bz2 in Eq. (98). Before doing this, notice from Eq. (D7) that there also are the diagonal components in DB (hnEB i) as [DB (hnEB i)]++ = k ˜k ˜k ∂εk ∂ F ∂εk ∂ F 3 2 −[DB (hnEB i)]−− = e E B τ ( − ), which z z k ∂ky ∂kx ∂kx ∂ky represent the Lorentz force contribution. Since these components are odd functions of kx and ky around a Weyl node, DB (hnEB i) does not drive the total particle number in a given Weyl cone. Hence the diagonal density matrix hn′EB 2 i obtained from these [DB (hnEB i)]nn k contains the intravalley scattering time as hn′EB 2 i = ˜ ˜k ∂εk ∂ F k ∂ Fk e3 Ez Bz2 τ τintra ( ∂ε σz . However, we see ∂ky ∂kx − ∂kx ∂ky )˜ that hn′EB 2 i does not contribute to the magnetoconductivity proportional to Bz2 , since it is an odd function of kx and ky and thus vanishes when integrated over momentum space. From Eqs. (54) and (D9), we obtain the relevant off-diagonal density matrix as "   vF vF m ∂ F˜k 3 2 σ ˜x cos θ − σ ˜y 2 sin θ hSEB 2 i =e Ez Bz τ ∂kx 2εk 2εk  # ∂ F˜k vF vF m + σ ˜x sin θ + σ ˜y 2 cos θ . ∂ky 2εk 2εk (108) Here, let us recall that there exists an extrinsic contribution to the off-diagonal density matrix hSEB 2 i from the anomalous driving term J(hn′EB 2 i), as in the case of hSE i. We can show explicitly by replacing hnE i in Eqs. (C2) and (C3) by hn′EB 2 i that   Re [J(hn′EB 2 i)]+− = Im [J(hn′EB 2 i)]+− = 0, (109) k k

and hence [J(hn′EB 2 i)]+− = [J(hn′EB 2 i)]−+ = 0, as long k k as the chemical potential µ lies sufficiently close to the Weyl nodes. Hence, the expression for hSEB 2 i, Eq. (108), remains valid even after the anomalous driving term is taken into account in the case of sufficiently small µ. Notice from Eqs. (45) and (55) that hSEB 2 i is accompanied by an intrinsic Berry phase contribution to the diagonal part of the density matrix proportional to Ez Bz2 . From Eq. (45) we readily obtain hξEB 2 i = −e3 Ez Bz2 F˜k τ

vF2 m σ ˜z . 2ε3k

(110)

We are now in a position to evaluate the zzcomponent of the magnetoconductivity proportional to

18

FIG. 3. k0 dependence of C(b, ∆, T ) in Eq. (111) with ∆ = 0, T = 0.005, and µ = 0.1. In our model Hamiltonian for a twonode Weyl metal (70), the distance between √ the two Weyl nodes in momentum space is given by 2k0 = 2 b2 − ∆2 /vF . CA Bz2 , σzz (Bz2 ) = Tr{(−e)vz [hSEB 2 i + hξEB 2 i]}/Ez . Since we have considered the scattering of electrons at the Fermi surface, we do not need to take into account the contribution from the t = + band. From Eqs. (89), (108) and (110), an explicit expression for σzz (Bz2 ) at low temperatures such that T ≪ µ is obtained as Z d3 k ∂m− vF2 X ∂ F˜k CA σzz (Bz2 ) = − e4 Bz2 τ ka 2 3 (2π) ∂kz ε−k a=x,y ∂ka Z d3 k ∂m− vF2 m2− ˜ Fk + e4 Bz2 τ (2π)3 ∂kz ε4−k



e2 (eBz )2 vF3 C(b, ∆, T ) τ 4π 2 µ2

CA Finally, we note that the contribution to σzz (Bz2 ) from a single isotropic Weyl cone is obtained by setting C(b, ∆, T ) ≃ 0.205 in Eq. (111) (see Fig. 3), since our Hamiltonian (70) with b = ∆ = 0 (i.e, k0 = 0) possesses a single isotropic Weyl cone. Note that the saturated value of C ≃ 0.41 in the region where k0 is large indicates the contribution from two isotropic Weyl cones. The above value from a single isotropic Weyl cone, C ≃ 0.205, is a little smaller than the value obtained by semiclassical wavepacket dynamics (96) that corresponds to C = 1. This is because the rate of pumping into a valley (105), ∂N/∂t = (2/3) e2 Ez Bz /4π 2 , is smaller than that obtained by semiclassical wavepacket dynamics (68), ∂N/∂t = e2 Ez Bz /4π 2 .

2.

Here, we briefly discuss the case of Dirac metals. To this end let us consider the contribution from a single Dirac cone, i.e., two degenerate Weyl cones. Time reversal symmetry is preserved in Dirac metals, and therefore we have to set b = 0. Also, for the existence of a Dirac point we have to set ∆ = 0. As mentioned above, Eq. (111) with b = ∆ = 0 represents the contribution CA from a single Weyl cone. Hence, the value of σzz (Bz2 ) coming from a single Dirac cone is obtained by setting b = ∆ = 0 and multiplying a factor 2 in Eq. (111): CA σzz (Bz2 )Dirac =

(111)

CA We numerically find that σzz (Bz2 ) is proportional to 2 3 1/µ and vF . This is consistent with the result by Son and Spivak, Eq. (96). We also find that, as in the case of the chiral magnetic effect (95), the contributions from CA hSEB 2 i and hξEB 2 i to σzz (Bz2 ) are respectively 2/3 and 1/3 of the total value (111). C(b, ∆, T ) in Eq. (111) is a dimensionless parameter of the order of 1. We show the k0 dependence √ of C(b, ∆, T ) with ∆ = 0 in Fig. 3. Here, 2k0 = 2 b2 − ∆2 /vF is the distance between the two Weyl nodes in momentum space. It is found that C is an increasing function of k0 and is saturated in the region where k0 is large. Here, note that we have not considered the k0 dependence of the intervalley scattering time τ . Although it is not easy to take into account such a k0 dependence in τ , it is expected that the value of τ will become larger as the value of k0 becomes larger, i.e., as the distance between the two Weyl nodes becomes longer. This is because the number of intervalley scattering processes will decrease as the value of k0 becomes larger, since intervalley scattering processes with large k0 require large momentum transfers. Thus, based on these viewpoints, we conclude CA that the value of σzz (Bz2 ) in Weyl metals becomes larger as the distance between the two Weyl nodes becomes longer.

The case of Dirac metals

(eBz )2 vF3 e2 C(b = ∆ = 0, T ) τDirac , 2 2π µ2 (112)

where τDirac is the intervalley scattering time between two Weyl cones that are degenerate in momentum space. As discussed in Sec. IV the effective Hamiltonian around a Dirac point in Dirac semimetals which are protected by crystalline symmetry can be written in the blockdiagonal form [17]     HAA HAB H+ (q) 0 , (113) HDirac (q) = = HBA HBB 0 H− (q) where A and B denote two states after a unitary transformation, and H± (q) is a 2 × 2 Weyl Hamiltonian with chirality ±1. In general, atomically smooth disorder does not allow scattering processes between two different states A and B. Therefore, the intervalley scattering time τDirac is much larger than intravalley scattering time.

3.

Contribution from the Lorentz force

Next, we compute the magnetoconductivity proportional to Bz2 due to the Lorentz force. Such a magnetoconductivity does not originate from the intrinsic (i.e., the Berry curvature) contribution but originates solely

19 from the Fermi surface contribution. Thus, only the diagonal components in the electric and magnetic driving terms are relevant to obtaining it. We start by obtaining the diagonal part of the density matrix induced by the electric field. It is easily seen that the integral of DE (hρ0 i) over the Fermi surface of a given Weyl cone is zero. Then, from the discussion in Sec. IV, the diagonal density matrix obtained from DE (hρ0 i) contains the intravalley scattering time as hnE i = τintra DE (hρ0 i), + − where n+ Ek = eEz τintra [∂f0 (εk )/∂kz ] and nEk = 0 (since we consider the case of µ > 0). From Eq. (D7), the diagonal components in DB (hnE i) are given by [DB (hnE i)]++ = eBz k



∂εk ∂n+ ∂εk ∂n+ Ek Ek − ∂ky ∂kx ∂kx ∂ky



(114)

and [DB (hnE i)]−− = 0. Note that Eq. (114) indeed k represents the Lorentz force term as [DB (hnE i)]++ = k e(vk × B) · ∇k n+ with v = ∇ ε . Again we find that k k k Ek the integral of DB (hnE i) over the Fermi surface of a given Weyl cone is zero. Then we see from Eq. (52) the diagonal density matrix obtained from DB (hnE i) contains the intravalley scattering time as hn′EB i = τintra DB (hnE i). The same procedure can be applied to obtain the diagonal density matrix proportional to Ez Bz2 , and we obtain ∂εk ∂ 2 hnEB 2 i = τintra DB (hnEB i) = diag[(eBz )2 τintra ( ∂k − y ∂kx ∂εk ∂ 2 + ∂kx ∂ky ) nEk , 0].

Finally, noting that only the t = − band is relevant, the magnetoconductivity proportional to Bz2 induced by the Lorentz force at low temperatures such that T ≪ µ is calculated to be d3 k ∂m− m− (2π)3 ∂kz ε−k 2  ∂ε−k ∂ ∂f0 (ε−k ) ∂ε−k ∂ − × ∂ky ∂kx ∂kx ∂ky ∂kz

3 LF (Bz2 ) = − e4 Bz2 τintra σzz

Z

2

0 ≈ − σzz (ωc τintra ) ,

When the contribution from the chiral anomaly dominates that from the Lorentz force, we have the total magnetoconductivity in the low-field limit as CA σzz (Bz2 ) ≈ σzz (Bz2 ) > 0.

We have confirmed that the magnetoconductivity linear in B vanishes in our formalism. Namely, the magnetoconductivity that has the lowest B-dependence is the one proportional to B 2 , which is in agreement with experimental results on the magnetoconductivity in the low-field limit in the Dirac semimetal Cd3 As2 [57]. On other hand, note that in this study we have not considered the contribution from the weak antilocalization such √ that ∼ a B with a < 0 in the low-field limit which is observed in the Weyl semimetal TaAs [59].

B.

2

LF σzz τintra 2 σ CA ∼ τ (µτintra ) . zz

DB (hSE i) = e2 Ez Bx Gk 1

2

which takes a negative value. Here, ≈ e (µ /vF )τintra is the Drude conductivity, and ωc = eBz vF2 /µ is the cyclotron frequency. Combining Eqs. (111) and (115), the total quadratic magnetoconductivity in parallel electric and magnetic CA LF fields is written as σzz (Bz2 ) = σzz (Bz2 ) + σzz (Bz2 ). The ratio of these two contributions to the magnetoconductivity is estimated as (116)

In general, the intervalley scattering time τ is much larger than the intravalley scattering time τintra , i.e., τintra /τ ≪ 1, since the intervalley scatterings require large momentum transfers, i.e., the number of intervalley scattering processes is much smaller than that of intravalley scattering processes. On the other hand, in usual high-mobility semiconductors the condition (µτintra )2 ≫ 1 is satisfied.

~ ⊥B ~ Quadratic magnetoconductivity for E

Let us consider the case of perpendicular electric and magnetic fields as E = (0, 0, Ez ) and B = (Bx , 0, 0). Even in this case, the intrinsic contribution to the magnetoconductivity proportional to Bx2 , which resembles that from the chiral anomaly (111), takes nonzero value in our formalism. However, it is not characterized by the intervalley scattering time τ which appeared in the case of parallel fields, but characterized by the intravalley scattering time τintra . We show that as the result the contribution from the Lorentz force dominates the intrinsic contribution, which is consistent with experimental results. We compute the intrinsic contribution to the quadratic magnetoconductivity. The exactly same calculation as in the case of parallel electric and magnetic fields can be applied here. After a calculation, we obtain the magnetic driving term which is proportional to Ez Bx as

(115) 0 σzz

(117)

(118)

with Gk = vF cos θ

∂ ∂kz

 

 − vF k⊥ ∂m f0 (ε+ k ) − f0 (εk ) 4ε3k ∂kz



.

(119)

Here, we show that this DB (hSE i) has the opposite property of the the case of parallel electric and magnetic fields. Around a Weyl node with momentum q = k − W± (q 2 ≪ 1), we find that Gq is an odd function of qx and qz . Accordingly, we find that the integral of DB (hSE i) over the Fermi surface of a given valley is R d3 q P mm zero: FS (2π) = 0. As discussed in 3 m [DB (hSE i)]q Sec. IV, this is a consequence of the total particle number conservation in a given Weyl cone. Namely we obtain Z e2 Ez Bx d3 q ∂N = 2 Gq = 0. (120) ∂t 4π 2 FS 2π

20 Then it follows from Eq. (52) that the intravalley scattering time τintra appears when L−1 acts on DB (hSE i). In the end, the diagonal density matrix hnEB i is obtained as (121) hnEB i = L−1 DB (hSE i) = e2 Ez Bx G˜k τintra 1, + − ∂ where G˜k = 12 Ω− k,x ∂kz [f0 (εk ) − f0 (εk )] is the component which represents the Fermi surface response in Eq. (119). The off-diagonal part hSEB 2 i and diagonal part hξEB 2 i of the density matrix proportional to Ez Bx2 can be obtained in a similar way as Eqs. (108) and (110), respectively. Finally, noting that only the t = − band is relevant, the zz-component of the quadratic magnetoconductivity from the intrinsic contribution at low temperatures such that T ≪ µ is obtained as Z d3 k ∂m− vF2 kx ∂ G˜k In σzz (Bx2 ) =e4 Bx2 τintra (2π)3 ∂kz ε2−k ∂kz Z d3 k ∂m− vF3 kx m− ˜ + e4 Bx2 τintra Gk (2π)3 ∂kz ε4−k



e2 ′ (eBz )2 vF3 C (b, ∆, T ) τintra 4π 2 µ2

(122)

with C ′ (b, ∆, T ) ∼ O(0.1). As mentioned above, this form resembles that from the chiral anomaly (111). However, the appearance of the intravalley scattering time indicates that the total particle number in a given valley is conserved, while the chiral anomaly requires nonconservation of the total particle number in a given valley. Thus, the magnetoconductivity (122) does not originate from the chiral anomaly. The quadratic magnetoconductivity due to the Lorentz force in the case of perpendicular electric and magnetic fields takes the same form as Eq. (115). Therefore, combining Eqs. (115) and (122), the total quadratic magnetoconductivity in perpendicular electric and magnetic fields In LF is written as σzz (Bx2 ) = σzz (Bx2 ) + σzz (Bx2 ). The ratio of these two contributions to the magnetoconductivity is estimated as LF In σ /σ ∼ (µτintra )2 ≫ 1, (123) zz

zz

where we have used that in usual high-mobility semiconductors the condition (µτintra )2 ≫ 1 is satisfied. Thus, we have the total magnetoconductivity in the low-field limit as LF σzz (Bx2 ) ≈ σzz (Bx2 ) < 0,

(124)

which is in agreement with experimental results in Dirac and Weyl semimetals [55–60]. As in the case of parallel electric and magnetic fields, we have confirmed that the magnetoconductivity linear in B vanishes in our formalism. VIII.

DISCUSSIONS

From the quantum kinetic theory viewpoint, we have clarified the microscopic origin of the appearance of the

intervalley scattering time τ in the positive quadratic magnetoconductivity induced by the chiral anomaly for parallel electric and magnetic fields, whereas preceding works based on semiclassical wavepacket dynamics [51, 54] introduced phenomenologically the intervalley scattering time. The appearance of the intervalley scattering time τ in the present study is due to the argument that the integral of the magnetic driving term DB over the Fermi surface of a given valley has a nonzero value (see Sec. IV). Although this argument is obtained under the assumption of atomically smooth disorder and sufficiently weak intervalley scattering compared to intravalley scattering (i.e., τ ≫ τintra ), such an assumption can in general be justified when the presence of magnetic impurities is neglected. In addition to the case of Weyl semimetals, our argument is naturally extended to the case of Dirac semimetals in which two valleys are degenerate. We have newly introduced the magnetic driving term DB (27), which plays a crucial role in this study. The magnetic driving term is understood as the natural matrix generalization of the well-known scalar formula, −e(∇k εk × B) · ∇k fk . The magnetic driving term has a property such that it maps an equilibrium density matrix hρ0 i (i.e., the Fermi-Dirac distribution function) to a purely off-diagonal density matrix hSB i. In other words, DB (hρ0 i) is a purely off-diagonal matrix. On the other hand, the matrix P −1 DB (hρ0 i) has diagonal components hξB i, which give rise to the Berry phase correction to hρ0 i. As is shown in Sec. VI C, these hSB i and hξB i result in the chiral magnetic effect. Here, recall from Sec. VII A that the density matrices hSEB 2 i and hξEB 2 i, which result in the positive quadratic magnetoconductivity in Weyl and Dirac metals, have similar forms to hSB i and hξB i. Namely, although our derivation of the positive quadratic magnetoconductivity was a little complicated, our study have shown that the positive quadratic magnetoconductivity is a consequence of the chiral magnetic effect, which is in agreement with preceding theoretical works [51–54, 58]. This study is also an illustration of the theory developed in Ref. [61] which accounts for the intraband and interband contributions to a physical quantity on an equal footing. Here, the intraband and interband contributions are described respectively by the diagonal and off-diagonal parts of a density matrix. It should be noted that the theory in Ref. [61], which corresponds to the case of B = 0 in the present theory, properly reproduces various physical quantities with the vertex correction in the ladder-diagram approximation taken into account, such as spin density [69], vanishing spin Hall conductivity [69], and vanishing anomalous Hall conductivity [70] in the Rashba models. This means that the contribution from J(hni) to conductivity corresponds to the vertex correction in the ladder-diagram approximation. As discussed below Eq. (54) the same argument is applied to the case of magnetotransport. Therefore, it is understood that the vertex corrections are absent in the

21 anomalous Hall conductivity (85) and positive quadratic magnetoconductivity (111) in Weyl metals, as long as the chemical potential lies sufficiently close to the Weyl points. Furthermore, our study shows that the anomalous Hall effect and chiral magnetic effect are obtained from the off-diagonal parts of density matrices, i.e., the off-diagonal parts of the electric and magnetic driving terms, respectively. This indicates that transport phenomena induced by the chiral anomaly are also important examples of the interband coherence in conductors. Lastly, we briefly discuss possible applications of our theory. We note that it is also possible to calculate other physical observables such as spin density and spin current in our formalism. More generally, all single-particle observables O maintain their crystal periodicity when they respond to an electric field and therefore have expectation values of the form hOi = Tr [Ohρi] ,

(125)

where hρi is a density matrix. Since the eigenstate basis and the Wannier basis can be used interchangeably, our theory can be incorporated into numerical methods such as first-principles calculations. Of course our theory can be applied to lattice models as well as continuum models. It can also be applied to other models such as insulators and nodal-line semimetals. In the case of insulators there is no Fermi surface so that only the offdiagonal components of density matrices are relevant. It will be interesting to study 2D multivalley systems such as graphene and transition metal dichalcogenides. The kinetic equation approach can also be applicable to interacting systems [71–73]. Thus it will be interesting to investigate from the quantum kinetic theory viewpoint the interplay between nontrivial band topology and electron correlations in the presence of electric and magnetic fields. IX.

Berry phase effects of Bloch electrons into transport phenomena, which are often discussed in the context of semiclassical wavepacket dynamics. Our theory is able to simply explain and predict when large momentum-space Berry curvatures yield important corrections to the standard Boltzmann theory of Fermi-surface dominated magnetotransport. As an illustration of our theory, we have applied our theory to transport phenomena in Weyl metals induced by the chiral anomaly: the anomalous Hall effect, chiral magnetic effect, and positive magnetoconductivity quadratic in magnetic field. We have shown that the anomalous Hall effect in Weyl metals is purely intrinsic, which means that the vertex correction in the ladderdiagram approximation is absent. We have also shown, by constructing the linear response theory to a magnetic field in the absence of electric fields, that the magnetic driving term we have newly introduced describes properly the chiral magnetic effect in a continuum model of Weyl metals. We have obtained an explicit expression for the positive quadratic magnetoconductivity that includes the intervalley scattering time in parallel electric and magnetic fields. In the process of obtaining this expression, we have shown that the vertex correction in the ladder-diagram approximation is absent. Our study indicates that the positive quadratic magnetoconductivity is a consequence of the chiral magnetic effect. On the other hand, in the case of perpendicular electric and magnetic fields, we have obtained a negative quadratic magnetoconductivity due to the Lorentz force that includes the intravalley scattering time, which is in agreement with experimental results. We have clarified by a general discussion that the chiral magnetotransport anomaly is observable only when intervalley scattering at the Fermi energy is very weak compared to intravalley scattering, and that it survives only partially even in this limit.

SUMMARY ACKNOWLEDGMENTS

In summary, we have developed a general quantum kinetic theory of low-field magnetotransport in weakly disordered multiband systems. By applying the Wigner transformation to the quantum Liouville equation, we have derived the quantum kinetic equation (25), which is the principle result of this study. From this equation we have obtained a generic expression for the magnetoconductivity that is applicable to the cases where the angle between electric and magnetic fields is arbitrary. Our theory naturally incorporates the momentum-space

Work at Austin was supported by the Department of Energy, Office of Basic Energy Sciences under Contract No. DE-FG02-ER45118 and by the Welch foundation under Grant No. TBF1473. A.S. is supported by the JSPS Overseas Research Fellowship. D.C. is supported by the Australian Research Council Centre of Excellence in Future Low-Energy Electronics Technologies (Project No. CE170100039). A.H.M. acknowledges stimulating conversations with Anton Burkov and Boris Spivak.

Appendix A: General properties of the magnetic driving term DB

In this Appendix, we consider general properties of the magnetic driving term (27). We write the Hamiltonian and P P the density matrix of a system as H0 = m εm |mihm| and hρ0 i = m f0m |mihm| with m a band index and f0m a

22 Fermi-Dirac distribution function, where we have omitted the wavevector dependences in these equations to simplify the notation. For concreteness, we consider the case of B = (0, 0, Bz ). Then the driving term (27) is written as        eBz DH0 DH0 Dhρ0 i Dhρ0 i DH0 Dhρ0 i e = − . (A1) , , ×B · DB (hρ0 i) = 2 2~ Dk Dk 2~2 Dky Dkx Dkx Dky First, we consider a diagonal component of DB (hρ0 i), i.e., hm|DB (hρ0 i)|mi. We immediately get X X   DH0 = ∂y εm′ |m′ ihm′ | + εm′ |∂y m′ ihm′ | + |m′ ih∂y m′ | , Dky m′ m′ X X   Dhρ0 i = ∂x f0n′ |n′ ihn′ | + f0n′ |∂x n′ ihn′ | + |n′ ih∂x n′ | , Dkx ′ ′ n

(A2)

n

where ∂a = ∂/∂ka . Then we have X X DH0 Dhρ0 i X = ∂y εm′ ∂x f0m′ |m′ ihm′ | + ∂y εm′ f0n′ |m′ ihm′ |∂x n′ ihn′ | + ∂y εm′ f0m′ |m′ ih∂x m′ | Dky Dkx m′ m′ n′ m′ X X X ′ ′ ′ ′ ′ ′ + εm′ ∂x f0m′ |∂y m ihm | + εm′ f0n′ |∂y m ihm |∂x n ihn | + εm′ f0m′ |∂y m′ ih∂x m′ | m′

+

X

m′ n′

m′ n′









εm′ ∂x f0n′ |m ih∂y m |n ihn | +

m′

X

m′ n′







εm′ f0n′ |m ih∂y m |∂x n ihn′ | +

X

m′ n′

εm′ f0n′ |m′ ih∂y m′ |n′ ih∂x n′ |, (A3)

from which a diagonal component is obtained as hm|

X X DH0 Dhρ0 i |mi =∂y εm ∂x f0m + εm′ f0m hm|∂y m′ ihm′ |∂x mi + εm′ f0m′ hm|∂y m′ ih∂x m′ |mi Dky Dkx ′ ′ m m X + εm f0m h∂y m|∂x mi + εm f0n′ h∂y m|n′ ih∂x n′ |mi,

(A4)

n′

where we have used hm′ |∂a n′ i + h∂a m′ |n′ i = ∂a (δm′ n′ ) = 0. Similarly we have hm|

X X Dhρ0 i DH0 |mi =∂y εm ∂x f0m + f0n′ εm hm|∂x n′ ihn′ |∂y mi + f0n′ εn′ hm|∂x n′ ih∂y n′ |mi Dkx Dky n′ n′ X ′ + f0m εm h∂x m|∂y mi + f0m εm′ h∂x m|m ih∂y m′ |mi.

(A5)

m′

Then we see that hm|DB (hρ0 i)|mi = eBz (∂y εm ∂x f0m − ∂x εm ∂y f0m ) = 0.

(A6)

Next, we consider a diagonal component of P −1 DB (hρ0 i), i.e., hm|P −1 DB (hρ0 i)|mi. Here, the matrix P −1 is defined by P hρi ≡ ~i [H0 , hρi], i.e., P −1 is the precession term that accounts for the time dependence of the density matrix in the absence of fields and disorder. Recall from Eq. (35) that when the matrix P −1 acts on a driving term D it simply multiplies it by a numerical factor: P −1 D = −i~

Dkmn n. εm k − εk

(A7)

P Let us take a closer look at the second term in the third line of Eq. (A3). Noticing the fact that ∂a ( m′ |m′ ihm′ |) = 0, we can replace εm′ by (εm′ − εn′ ) in Eq. (A2). Then we can rewrite it as X DH0 Dhρ0 i = (εm′ − εn′ )f0n′ |m′ ih∂y m′ |∂x n′ ihn′ |, Dky Dkx ′ ′

(A8)

mn

from which we have P −1

X εm′ − εn′ X DH0 Dhρ0 i = −i~ f0n′ |m′ ih∂y m′ |∂x n′ ihn′ | = −i~ f0n′ |m′ ih∂y m′ |∂x n′ ihn′ |. Dky Dkx εm′ − εn′ ′ ′ ′ ′ mn

mn

(A9)

23 0 Dhρ0 i Then a diagonal element of this matrix is obtained as hm|P −1 DH Dky Dkx |mi = −i~f0m h∂y m|∂x mi. Similarly we have 0 i DH0 hm|P −1 Dhρ Dkx Dky |mi = +i~f0m h∂x m|∂y mi. In the end, we see that

hm|P −1 DB (hρ0 i)|mi = (e/~)f0m Bz Ωm z ,

(A10)

where Ωm z = ih∂x m|∂y mi − ih∂y m|∂x mi is the z component of the Berry curvature of band m. Note that only the terms of the form of (A8) give rise to nonzero contributions to the diagonal component of P −1 DB (hρ0 i). Appendix B: General expression for Tr [DB (hSE i)] in parallel electric and magnetic fields

In this Appendix, we derive an explicit P expression for the diagonal matrix element [DB (hSE i)]mm in a general system P with H0 = m f0m |mihm| with m a band index and f0m a Fermi-Dirac distribution m εm |mihm| and hρ0 i = function, where we have omitted the wavevector k dependences in these equations to simplify the notation. For concreteness, we consider the case of E = (0, 0, Ez ) and B = (0, 0, Bz ). An off-diagonal component of the electric driving term (26) reads X   hn|DE (hρ0 i)|n′ i = eEz f0m′ hn| |∂z m′ ihm′ | + |m′ ih∂z m′ | |n′ i = eEz (f0n′ − f0n )hn|∂z n′ i (B1) m′

with ∂a = ∂/∂ka , from which we obtain the off-diagonal part of the density matrix induced by the electric field hSE i = (−i)eEz

X f0n′ − f0n |nihn|∂z n′ ihn′ |, ′ ε − ε n n ′

(B2)

nn

where n 6= n′ . On the other hand, the magnetic driving term (27) is written as        e eBz DH0 DH0 DhSE i DhSE i DH0 DhSE i DB (hSE i) = = − . , , ×B · 2 Dk Dk 2 Dky Dkx Dkx Dky

(B3)

Since we are focusing on the Fermi surface response, we consider only the terms proportional to ∂x f0 in DhSE i/Dkx : X ∂x f0n′ − ∂x f0n DhSE i = (−i)eEz |nihn|∂z n′ ihn′ |. ′ Dkx ε − ε n n ′

(B4)

nn

We also have the relevant term X   DH0 = (εm′ − εn′ ) |∂y m′ ihm′ | + |m′ ih∂y m′ | , Dky ′

(B5)

m

P where we have used ∂a ( m′ |m′ ihm′ |) = 0. Note that the terms proportional to ∂y ε in DH0 /Dky do not contribute to the diagonal part of DB (hSE i). Then we obtain hm|

XX   DH0 DhSE i ∂x f0n′ − ∂x f0n |mi = (−i)eEz δn′ m δm′ n hm|∂y nihn|∂z mi + δmm′ h∂y m|nihn|∂z mi (εm′ − εn′ ) Dky Dkx εn − εn′ nn′ m′ X = (−i)eEz (∂x f0m − ∂x f0n )hm|∂y nihn|∂z mi n

= ieEz ∂x f0m h∂y m|∂z mi − ieEz

X n

∂x f0n h∂z n|mihm|∂y ni.

(B6)

Similarly we have hm|

X DhSE i DH0 |mi = −ieEz ∂x f0m h∂z m|∂y mi + ieEz ∂x f0n h∂y n|mihm|∂z ni. Dkx Dky n

Finally we arrive at a general expression for the rate of pumping from the Fermi surface contribution: X  ∂N m = Tr [DB (hSE i)] = e2 Ez Bz ∂x f0m Ωm x + ∂y f0m Ωy , ∂t m,k

abc where Ωm ih∂b m|∂c mi is the Berry curvature of band m. a = ǫ

(B7)

(B8)

24 II Appendix C: Evaluation of the anomalous Hall conductivity σxy

II In this Appendix, we calculate explicitly the anomalous Hall conductivity σxy which results from the extrinsic (i.e., Fermi surface) contribution, and show that it vanishes when the chemical potential µ lies sufficiently close to the Weyl points. We have considered the case where an electric field is applied along the y direction as E = Ey ey . In this case, − the diagonal part of the density matrix induced by the electric field hnE i = diag[n+ Ek , nEk ] is given by + ∂ε+ k⊥ k ∂f0 (εk ) sin θ δ(εk − µ) and n− (C1) = −eτ + Ey vF2 Ek = 0, + ∂k ∂εk εk q = ±εk = ± vF2 (kx2 + ky2 ) + m2 and considered the case of µ > 0. Then from Eq. (83) we

+ n+ Ek = −eτ E ·

where we have used ε± k

obtain the imaginary part of [J(hnE i)]+− as k

  nimp U02 X k′  + + + =π δ(ε+ − ε+ [sin θ′ cos θ − cos θ′ sin θ] ⊥ n+ Im [J(hnE i)]+− ′ ) − nEk′ δ(εk − εk′ ) Ek k k k 2 ε k′ ′ k

X nimp U02 + k′ =π eτ Ey vF2 cos θ sin2 θ′ ⊥ δ(εk′ − µ)δ(εk − εk′ ), 2 ε k′ ′

(C2)

k

′ ′ ′ ′ ′ ′ where we have used that the contribution from n+ Ek vanishes since sin θ = ky /k⊥ and cos θ = kx /k⊥ are odd functions. Similarly, we obtain the real part of [J(hnE i)]+− as k   2 ′ ′   +  nimp U0 X k⊥ m m k⊥ + + + ′ Re [J(hnE i)]+− = π − cos(θ − θ) nEk δ(ε+ − ε+ ′ ) − nEk′ δ(εk − εk′ ) k k k 2 ε k ε k′ ε k ε k′ k′ # " ′ X X m′ k m k nimpU02 + ⊥ sin2 θ′ ⊥ δ(εk′ − µ)δ(εk − εk′ ) . δ(εk − µ) δ(εk − εk′ ) + eτ Ey vF2 sin θ = −π 2 εk ε k′ εk ′ ε k′ ′ k

k

(C3)

II II Finally, the extrinsic contribution to the anomalous Hall conductivity σxy is calculated from the definition σxy = ′ Tr[(−e)vx hSE i]/Ey . Using Eqs. (80) and (84), in the case of small chemical potential µ we obtain  X X 1   mt +− +− II sin θ − 2 Im [J(hnE i)]k cos θ Re [J(hnE i)]k σxy =2e 2εtk 2εtk t=± k  ′ 2 ′ X  m− k⊥ 2 ′ k⊥ 2 m− 2 nimp U0 + 2 = − πe δ(εk′ − µ) + 2 sin θ δ(εk − µ) δ(εk − εk′ ) τ vF sin θ 2 ε2k ε k′ εk ε k′ k,k′ Z δ 3 Z δ 3 ′ 2 d q d q 2 nimp U0 + 4 ≈ − πe δ(εq − εq′ ) τ vF 3 3 2 (2π) (2π) −δ −δ " # k0 k0 k0 ′ ′ qz − kb0 qz′ q⊥ 2 ′ q⊥ 2 b qz − b qz b × δ(εq ′ − µ) + 2 sin θ δ(εq − µ) sin θ ε2q εq ′ εq εq ′

=0,

(C4)

where we have omitted the subscript − in ε−k in the third through last lines. Note that the contributions only from t = − survive since the chemical potential µ lies at ε+ used the fact that m− (kz ) ≈ ∓vF2 (k0 /b)qz −k . Also, we have √ around the Weyl nodes, where q = k − W± = (kx , ky , kz ∓ k0 ) with k0 = b2 − ∆2 /vF and q 2 ≪ 1. Appendix D: Explicit matrix forms of the magnetic driving term DB

In this Appendix, we show the explicit matrix forms of the magnetic driving term (27) obtained from arbitrary diagonal and off-diagonal density matrices. For concreteness, let us consider the case of a two-band model with the two eigenvalues ε± and the density matrix hρi. We define the diagonal and off-diagonal parts of hρi as hni and hSi, respectively. The total density matrix of the system is given by hρi = hni + hSi. In the eigenstate basis, we have  +   +    εk − µ 0 nk 0 0 ak H0 = , hni ≡ , (D1) , hSi ≡ bk 0 0 ε− 0 n− k −µ k

25 where µ is a chemical potential. We immediately get     0 R+− (n− − n+ ) 0 R+− (ε− − ε+ ) α α k k k k , (D2) , [Rα , hni] = [Rα , H0 ] = − + − + −R−+ 0 −R−+ 0 α (nk − nk ) α (ε − εk )   +− k −+ −− Rα bk − Rα ak (R++ α − Rα )ak , (D3) [Rα , hSi] = ++ −− +− −(Rα − Rα )bk −(Rα bk − R−+ α ak ) P where R = α=x,y,z Rα eα is the Berry connection vector. First, let us calculate explicitly the magnetic driving term originating from the diagonal part of the density matrix hni. For concreteness, we consider the case of B = (0, 0, Bz ). Then the driving term (27) is written as        1 DH0 Dhni Dhni DH0 Dhni DH0 1 = eBz − . (D4) , , ×B · DB (hni) = e 2 Dk Dk 2 Dky Dkx Dkx Dky After a straightforward calculation, we obtain  " # ∂ε+ ∂n+ + +− − + k +− − k −iR (ε − ε ) (n − n ) −iR DH0 Dhni  y k k x ∂ky k k ∂kx  = ∂n− ∂ε− + − + −+ − k k Dky Dkx iR (n − n ) iR−+ (ε − ε ) x y k k k k ∂ky ∂kx   + + + ∂ε ∂n− ∂εk ∂nk + − + + + +− −+ − +− − +− − k k + R R (ε − ε )(n − n ) −i R (n − n ) − i R (ε − ε ) y x x y k k k k k k k k ∂k ∂k ∂k ∂k y x , =  ∂ny+ x − ∂ε− ∂ε− − + + + − + −+ − −+ +− − k k ∂nk i ∂kkx R−+ (ε − ε ) + i R (n − n ) + R R (ε − ε )(n − n ) y x y x k k k k k k k k ∂ky ∂ky ∂kx ∂n+ k ∂kx − iR−+ x (nk −

Dhni DH0 = Dkx Dky

"

=



which results in

n+ k)

− −iR+− x (nk ∂n− k ∂kx



n+ k)

#

∂ε+ k ∂ky  − + iR−+ y (εk − εk )

+ ∂ε+ + − + +− − k ∂nk + R−+ y Rx (εk − εk )(nk − nk ) ∂k ∂k y  − x + ∂n ∂εk − + + −+ − i ∂kkx R−+ y (εk − εk ) + i ∂ky Rx (nk − nk )



DH0 Dhni , Dky Dkx





DH0 Dhni , Dkx Dky



where

=



2 ∂kky

=



2 ∂kkx

∂ε+ ∂n+ ++ k ∂kx + gxy  + ∂ + − iRy−+ (ε− k − εk ) ∂kx (nk + nk )

∂ε+ ∂n+ ++ k ∂ky + gyx  + ∂ + − iRx−+ (ε− k − εk ) ∂ky (nk + nk )

−+ + gxy

−+ + gyx

− + −iR+− y (εk − εk ) ∂ε− k ∂ky

 

 ∂ε− ∂n+ − + + +− − k −i ∂kky R+− (n − n ) − i R (ε − ε ) x y k k k k  ∂kx , − ∂ε− + − + +− −+ − k ∂nk + R R (ε − ε )(n − n ) y x k k k k ∂ky ∂kx − + ∂ + − +− −iR+− y (εk − εk ) ∂kx (nk + nk ) + gxy ∂ε− ∂n− k ∂kx

2 ∂kky

−− + gxy

− + ∂ + − +− −iR+− x (εk − εk ) ∂ky (nk + nk ) + gyx ∂ε− ∂n− k ∂ky

2 ∂kkx

(D5)

−− + gyx

(D6)



,



,

(D7)

+ − + + − + ++ −+ − −+ +− − gxy = R+− y Rx (εk − εk )(nk − nk ) + Ry Rx (εk − εk )(nk − nk ),

− + − + + − + −− +− −+ − gxy = Ry−+ R+− x (εk − εk )(nk − nk ) + Ry Rx (εk − εk )(nk − nk ), ∂ − + +− (ε+ + ε− gxy = −iR+− x (nk − nk ) k ), ∂ky k ∂ − + −+ (ε+ + ε− gxy = iR−+ x (nk − nk ) k ). ∂ky k

(D8)

± ± ++ ++ −− −− Here, note that gxy = gyx and gxy = gyx . In the case of n± k = f0 (εk ) where f0 (εk ) is the Fermi-Dirac distribution

function, we have

± ∂ε± k ∂f0 (εk ) ∂ky ∂kx



± ∂ε± k ∂f0 (εk ) = 0. Namely, the magnetic driving term (D4) is purely ∂kx ∂ky ± of εk = ±εk , which applies to Weyl semimetals, we have g +− =

case. Furthermore, in the case such a case the magnetic driving term is simplified to be " 0 DB (hni) = eBz + − ∂ −+ ∂ iεk (R−+ − R x ∂ky y ∂kx )[f0 (εk ) + f0 (εk )]

off-diagonal in this g −+ = 0. Then in

# + − ∂ +− ∂ −iεk (R+− x ∂ky − Ry ∂kx )[f0 (εk ) + f0 (εk )] . 0

(D9)

26 Next, let us calculate explicitly the magnetic driving term originating from the off-diagonal part of the density matrix hSi. For concreteness, we consider the case of B = (0, 0, Bz ). Then the driving term (27) is written as        1 DH0 DhSi DhSi DH0 DhSi 1 DH0 = eBz − . (D10) , , ×B · DB (hSi) = e 2 Dk Dk 2 Dky Dkx Dkx Dky After a straightforward calculation, we obtain    ∂ε+ + +− − k ∂a  −iR (ε − ε ) −iF −i∆a + ∂k DH0 DhSi  y k k  ∂ky x = − ∂b ∂εk − + iF i∆b + ∂k Dky Dkx iR−+ (ε − ε ) x y k k ∂ky   + + ∂εk ∂ε − + + ∂b ∂a +− − (ε − ε )[i∆b + ] [−i∆a + ] + R (ε − ε )F −i ∂kky F − iR+− y y k k k k ∂kx ∂ky ∂kx , = ∂ε− ∂ε− + + ∂b ∂a −+ − k k Ry−+ (ε− − ε )F + [i∆b + ] iR (ε − ε )[−i∆a + ] + i F y k k k k ∂ky ∂kx ∂kx ∂ky   ∂ε+ + +− − k ∂a  (ε − ε ) −iR −iF −i∆a + DhSi DH0 y k k  ∂ky ∂kx  = ∂b ∂ε− + −+ − i∆b + ∂k iF k Dkx Dky iRy (εk − εk ) x ∂ky   + − ∂ε ∂εk − + + ∂a ∂a +− − −i ∂kky F + iR−+ (ε − ε )[−i∆a + ] [−i∆a + ] − R (ε − ε )F y y k k k k ∂kx ∂ky ∂kx , = ∂ε+ ∂ε− + + ∂b ∂b −+ − +− − k −Ry (εk − εk )F + ∂ky [i∆b + ∂kx ] −iRy (εk − εk )[i∆b + ∂kx ] + i ∂kky F

(D11)



(D12)

± −+ where ∆ = R++ − R−− and F = (R+− x x x b − Rx a). In the case of εk = ±εk , which applies to Weyl semimetals, we have         1 DH0 DhSi 1 DH0 DhSi Ak 0 B 0 = , = k (D13) , , 0 Ak 0 Bk 2 Dky Dkx 2 Dkx Dky

with  ∂εk +− ++ +− − R−− (Rx bk − R−+ a ) + iR ε k k i(Rx x )bk + x y ∂ky  ∂εk +− ++ +− − R−− (Ry bk − R−+ a ) + iR ε Bk = −i k k i(Ry y )bk + y x ∂kx

Ak = −i

  ∂bk ++ −− − iR−+ ε k −i(Rx − Rx )ak + y ∂kx   ∂bk ++ − R−− − iR−+ ε k −i(Ry y )ak + x ∂ky

 ∂ak , ∂kx  ∂ak . ∂ky (D14)

Especially in the case of ak = −ick and bk = ick , i.e., hSik = ck σy , we have ∂εk ∂ck ++ +− ck − i(R+− − R−+ − R−− + R−+ , y y )(Rx x )εk ck − (Ry y )εk ∂ky ∂kx ∂εk ∂ck ++ +− Bk = (R+− + R−+ ck − i(R+− − R−+ − R−− + R−+ . y y ) x x )(Ry y )εk ck − (Rx x )εk ∂kx ∂ky

Ak = (R+− + R−+ x x )

Finally, we obtain the magnetic driving term which is purely diagonal,      DH0 DhSi DH0 DhSi 1 F − = eBz k , , DB (hSi) = eBz 0 2 Dky Dkx Dkx Dky

0 Fk



(D15)

(D16)

with Fk = Ak − Bk .

[1] D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys. 82, 1959 (2010).

[2] N. Nagaosa, J. Sinova, S. Onoda, A. H. MacDonald, and N. P. Ong, Rev. Mod. Phys. 82, 1539 (2010).

27 [3] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982). [4] F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988). [5] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005). [6] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005). [7] B. A. Bernevig and S.-C. Zhang, Phys. Rev. Lett. 96, 106802 (2006). [8] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010). [9] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011). [10] Y. Ando, J. Phys. Soc. Jpn. 82, 102001 (2013). [11] G. E. Volovik, The Universe in a Helium Droplet (Clarendon, Oxford, U.K., 2003). [12] S. Murakami, New J. Phys. 9, 356 (2007). [13] X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Phys. Rev. B 83, 205101 (2011). [14] S. M. Young, S. Zaheer, J. C. Y. Teo, C. L. Kane, E. J. Mele, and A. M. Rappe, Phys. Rev. Lett. 108, 140405 (2012). [15] Z. Wang, H. Weng, Q. Wu, X. Dai, and Z. Fang, Phys. Rev. B 88, 125427 (2013). [16] Z. Wang, Y. Sun, X.-Q. Chen, C. Franchini, G. Xu, H. Weng, X. Dai, and Z. Fang, Phys. Rev. B 85, 195320 (2012). [17] B.-J. Yang and N. Nagaosa, Nat. Commun. 5, 4898 (2014). [18] A. H. MacDonald and Qian Niu, Physics World 17, 18 (2004). [19] Z. Fang, N. Nagaosa, K. S. Takahashi, A. Asamitsu, R. Mathieu, T. Ogasawara, H. Yamada, M. Kawasaki, Y. Tokura, and K. Terakura, Science, 302, 92 (2003). [20] Z. K. Liu, B. Zhou, Y. Zhang, Z. J. Wang, H. M. Weng, D. Prabhakaran, S.-K. Mo, Z. X. Shen, Z. Fang, X. Dai, Z. Hussain, and Y. L. Chen, Science 343, 864 (2014). [21] Z. K. Liu, J. Jiang, B. Zhou, Z. J. Wang, Y. Zhang, H. M. Weng, D. Prabhakaran, S.-K. Mo, H. Peng, P. Dudin, T. Kim, M. Hoesch, Z. Fang, X. Dai, Z. X. Shen, D. L. Feng, Z. Hussain, and Y. L. Chen, Nat. Mater. 13, 677 (2014). [22] M. Neupane, S.-Y. Xu, R. Sankar, N. Alidoust, G. Bian, C. Liu, I. Belopolski, T.-R. Chang, H.-T. Jeng, H. Lin, A. Bansil, F. Chou, and M. Z. Hasan, Nat. Commun. 5, 3786 (2014). [23] S. Borisenko, Q. Gibson, D. Evtushinsky, V. Zabolotnyy, B. Bchner, and R. J. Cava, Phys. Rev. Lett. 113, 027603 (2014). [24] B. Q. Lv, N. Xu, H. M. Weng, J. Z. Ma, P. Richard, X. C. Huang, L. X. Zhao, G. F. Chen, C. E. Matt, F. Bisti, V. N. Strocov, J. Mesot, Z. Fang, X. Dai, T. Qian, M. Shi, and H. Ding, Nat. Phys. 11, 724 (2015). [25] L. X. Yang, Z. K. Liu, Y. Sun, H. Peng, H. F. Yang, T. Zhang, B. Zhou, Y. Zhang, Y. F. Guo, M. Rahn, D. Prabhakaran, Z. Hussain, S.-K. Mo, C. Felser, B. Yan, and Y. L. Chen, Nat. Phys. 11, 728 (2015). [26] S.-Y. Xu, I. Belopolski, N. Alidoust, M. Neupane, G. Bian, C. Zhang, R. Sankar, G. Chang, Z. Yuan, C.-C. Lee, S.-M. Huang, H. Zheng, J. Ma, D. S. Sanchez, B. Wang, A. Bansil, F. Chou, P. P. Shibayev, H. Lin, S. Jia, and M. Z. Hasan, Science 349, 613 (2015). [27] L. Lu, Z. Wang, D. Ye, L. Ran, L. Fu, J. D. Joannopoulos, and M. Solja i, Science 349, 622 (2015).

[28] F. Wilczek, Phys. Rev. Lett. 58, 1799 (1987). [29] X.-L. Qi, T. L. Hughes, and S.-C. Zhang, Phys. Rev. B 78, 195424 (2008). [30] A. M. Essin, J. E. Moore, and D. Vanderbilt, Phys. Rev. Lett. 102, 146805 (2009). [31] K. Nomura and N. Nagaosa, Phys. Rev. Lett. 106, 166802 (2011). [32] D. Baasanjav, O. A. Tretiakov, and K. Nomura, Phys. Rev. B 90, 045149 (2014). [33] K. Everschor-Sitte, M. Sitte, and A. H. MacDonald, Phys. Rev. B 92, 245118 (2015). [34] A. A. Zyuzin and A. A. Burkov, Phys. Rev. B 86, 115133 (2012). [35] D. T. Son and N. Yamamoto, Phys. Rev. Lett. 109, 181602 (2012). [36] A. Grushin, Phys. Rev. D 86, 045001 (2012). [37] Z. Wang and S.-C. Zhang, Phys. Rev. B 87, 161107 (2013). [38] P. Goswami and S. Tewari, Phys. Rev. B 88, 245107 (2013). [39] A. A. Burkov, J. Phys. Condens. Matter 27, 113201 (2015). [40] R. Li, J. Wang, X.-L. Qi, and S.-C. Zhang, Nat. Phys. 6, 284 (2010). [41] K. Fukushima, D. E. Kharzeev, and H. J. Warringa, Phys. Rev. D 78, 074033 (2008). [42] M. M. Vazifeh and M. Franz, Phys. Rev. Lett. 111, 027201 (2013). [43] M.-C. Chang and M.-F. Yang, Phys. Rev. B 91, 115203 (2015). [44] P. V Buividovich, M. Puhr, and S. N. Valgushev, Phys. Rev. B 92, 205122 (2015). [45] N. Yamamoto, Phys. Rev. D 92, 085011 (2015). [46] J. Ma and D. A. Pesin, Phys. Rev. B 92, 235205 (2015). [47] S. Zhong, J. E. Moore, and I. Souza, Phys. Rev. Lett. 116, 077201 (2016). [48] M. A. Zubkov, Phys. Rev. D 93, 105036 (2016). [49] A. Sekine and K. Nomura, Phys. Rev. Lett. 116, 096401 (2016). [50] A. Sekine and T. Chiba, Phys. Rev. B 93, 220403(R) (2016). [51] D. T. Son and B. Z. Spivak, Phys. Rev. B 88, 104412 (2013). [52] A. A. Burkov, Phys. Rev. Lett. 113, 247203 (2014). [53] A. A. Burkov, Phys. Rev. B 91, 245157 (2015). [54] B. Z. Spivak and A. V. Andreev, Phys. Rev. B 93, 085107 (2016). [55] J. Xiong, S. K. Kushwaha, T. Liang, J. W. Krizan, M. Hirschberger, W. Wang, R. J. Cava, and N. P. Ong, Science 350, 413 (2015). [56] C.-Z. Li, L.-X. Wang, H. Liu, J. Wang, Z.-M. Liao, and D.-P. Yu, Nat. Commun. 6, 10137 (2015). [57] H. Li, H. He, H.-Z. Lu, H. Zhang, H. Liu, R. Ma, Z. Fan, S.-Q. Shen, and J. Wang, Nat. Commun. 7, 10301 (2016). [58] 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. 12, 550 (2016). [59] 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 5, 031023 (2015). [60] F. Arnold, C. Shekhar, S.-C. Wu, Y. Sun, R. D. dos Reis, N. Kumar, M. Naumann, M. O. Ajeesh, M. Schmidt, A. G. Grushin, J. H. Bardarson, M. Baenitz, D. Sokolov, H.

28

[61] [62] [63] [64] [65] [66] [67]

Borrmann, M. Nicklas, C. Felser, E. Hassinger, and B. Yan, Nat. Commun. 7, 11615 (2016). D. Culcer, A. Sekine, and A. H. MacDonald, arXiv:1612.06865. F. T. Vasko and O. E. Raichev, Quantum Kinetic Theory and Applications (Springer, New York, 2005). D. Xiao, J. Shi, and Q. Niu, Phys. Rev. Lett. 95, 137204 (2005). P. Strˇeda, J. Phys. C 15, L717 (1982). A. A. Burkov, M. D. Hook, and L. Balents, Phys. Rev. B 84, 235126 (2011). A. Sekine and K. Nomura, J. Phys. Soc. Jpn. 83, 094710 (2014). K. Fujikawa, Phys. Rev. Lett. 42, 1195 (1979); K. Fu-

jikawa, Phys. Rev. D 21, 2848 (1980). [68] When the chemical potential lies on the conduction band, the situation is quite a bit more complicated. See, for example, D. Culcer and S. Das Sarma, Phys. Rev. B 83, 245441 (2011). [69] J. Inoue, G. E. W. Bauer, and L. W. Molenkamp, Phys. Rev. B 70, 041303 (2004). [70] J. Inoue, T. Kato, Y. Ishikawa, H. Itoh, G. E. W. Bauer, and L. W. Molenkamp, Phys. Rev. Lett. 97, 046604 (2006). [71] D. Culcer, Phys. Rev. B 84, 235411 (2011). [72] J. Wang and D. Culcer, Phys. Rev. B 88, 125140 (2013). [73] H. Liu, W. E. Liu, and D. Culcer, arXiv:1610.05866.