Scattering theory of the chiral magnetic effect in a Weyl semimetal: Interplay of bulk Weyl cones and surface Fermi arcs

We formulate a linear response theory of the chiral magnetic effect in a finite Weyl semimetal, expressing the electrical current density $j$ induced by a slowly oscillating magnetic field $B$ or chiral chemical potential $\mu$ in terms of the scattering matrix of Weyl fermions at the Fermi level. Surface conduction can be neglected in the infinite-system limit for $\delta j/\delta \mu$, but not for $\delta j/\delta B$: The chirally circulating surface Fermi arcs give a comparable contribution to the bulk Weyl cones no matter how large the system is, because their smaller number is compensated by an increased flux sensitivity. The Fermi arc contribution to $\mu^{-1}\delta j/\delta B$ has the universal value $(e/h)^2$, protected by chirality against impurity scattering --- unlike the bulk contribution of opposite sign.


I. INTRODUCTION
The conduction electrons in a Weyl semimetal have an unusual velocity distribution in the Brillouin zone [1]. The conical band structure (Weyl cone) has a chirality that generates a net current at the Fermi level in the presence of a magnetic field [2]. The Weyl cones come in pairs of opposite chirality, so that the total current vanishes in equilibrium [3][4][5], but a nonzero current I parallel to the field B remains if the cones are offset by an energy µ -slowly oscillating to prevent equilibration [6][7][8][9][10][11]. This is the chiral magnetic effect (CME) from particle physics [12][13][14][15], see Refs. [16][17][18] for recent reviews in the condensed matter setting. In an infinite system the current density has the universal form [19][20][21] j 0 = −(e/h) 2 independent of material parameters. This amounts to a conductance of e 2 /h in the lowest (zeroth) Landau level, multiplied by the degeneracy equal to the enclosed flux in units of the flux quantum. The recent condensed-matter realizations of Weyl semimetals [22][23][24][25] have boosted the search for the chiral magnetic effect [26][27][28][29][30][31][32][33][34]. Future experimental developments may well include nanostructured materials, to minimize effects of disorder. In a finite system, the zeroth Landau level in the bulk hybridizes with the Fermi arcs connecting the two Weyl cones along the surface [35,36]. Previous studies [37,38] have pointed to the importance of boundaries for the chiral magnetic effect -a sign reversal of the current density as one moves from the bulk towards a boundary ensures that zero current flows in response to a static perturbation. Here we wish to study how this interplay of surface and bulk states impacts on the chiral magnetic effect in response to a low-frequency dynamical perturbation. For that purpose we seek a linear response theory that does not assume translational invariance in an infinite system. A scattering formulationà la Landauer seems most appropriate for such a mesoscopic system.
The Landauer approach to electrical conduction considers the current driven between two spatially separated electron reservoirs by a chemical potential difference, and expresses this in linear response by a sum over transmission probabilities at the Fermi level [39][40][41]. The chiral magnetic effect is driven by a nonequilibrium population of the Weyl cones, so in reciprocal space (Brillouin zone) rather than in real space -we will show how to modify the Landauer formula accordingly.
We first apply our scattering formula to a current driven by a slowly oscillating offset µ of the Weyl cones (a so-called "chiral" or "axial" chemical potential [15]), and recover Eq. (1) in the infinite-system limit. We then turn to the more practical scenario of a current driven by a slowly oscillating magnetic field B. We find that the surface Fermi arcs give a contribution to the total induced current equal to minus twice the bulk contribution in the infinite-system limit. That the surface Fermi arc contribution does not vanish relative to the bulk contribution is unexpected and not captured by previous calculations of the chiral magnetic effect.
The outline of the paper is as follows. In the next section II we derive the scattering formula for the chiral magnetic effect, in a general setting. In Sec. IV we apply it to the model Hamiltonian of a Weyl semimetal from Ref. 3, summarized in Sec. III. We evaluate the induced current in response to variations in µ and B, both numerically for a finite system and analytically in the limit of an infinite system size. Finite-size corrections are considered in some detail in Sec. V. We conclude in Sec. VI with a summary and a discussion of the robustness of the results against disorder scattering.

II. SCATTERING FORMULA
For a scattering theory of the chiral magnetic effect we consider a disordered mesoscopic system attached to ideal leads. Such an "electron wave guide" has propagating modes with band structure E n (k), labeled by a mode arXiv:1512.02144v3 [cond-mat.mes-hall] 30 Mar 2016 index n = 1, 2, . . . and dependent on the wave vector k along the lead. At a given energy ε (measured relative to the equilibrium Fermi level E F ), each incident mode has wave vector k n (ε) and carries the same current e/h per unit energy interval [42].
The scattering matrix S(ε) relates amplitudes of incident and outgoing modes. We take a two-terminal geometry (the multi-terminal generalization is straightforward), with N modes each in the left and right lead -so S is a 2N × 2N unitary matrix. The current I through the system can be calculated in the left lead, by current conservation it must be the same through each cross section.
The projection matrix onto the left lead is P = 1 0 0 0 , where each sub-block is an N × N matrix. The current is driven by a set of non-equilibrium occupation numbers δf n (ε), with n = 1, 2, . . . N for the left lead and n = N + 1, N + 2, . . . 2N for the right lead. We collect these numbers in a 2N × 2N diagonal matrix δF(ε). The net current in the left lead is then given by the difference of incoming and outgoing currents, We consider the linear response to a slowly varying parameter X that adiabatically perturbs the system away from its equilibrium state at X = X 0 . We assume that the wave vector k along the lead (say, in the z-direction) is not changed by the perturbation. This requires that the perturbation should neither break the translational invariance along z, nor involve a time-dependent vector potential component A z .
The band structure evolves from E n (k|X 0 ) to E n (k|X 0 + δX). To first order in the perturbation δX the energy shift at constant k is (3) The corresponding deviation of the occupation number from the equilibrium Fermi function f eq (ε) = (1 + e ε/kBT ) −1 is where we have used that E n (k n (ε)|X 0 + δX) ≡ ε. At zero temperature the derivative f eq (ε) → −δ(ε), so the expression (2) for the current contains only Fermi level scattering amplitudes. We may write it in a more explicit form in terms of the transmission probabilities The sign of the derivative ∂E n /∂k distinguishes the rightmoving modes n = 1, 2, . . . N from the left-moving modes n = N + 1, N + 2, . . . 2N . The Landauer conductance formula [39][40][41] is obtained from Eq. (6) if we identify δX = V with the voltage difference between the left and right lead, and then set χ n = 1 for n = 1, 2, . . . N and χ n = 0 for n = N + 1, N + 2, . . . 2N . The chiral magnetic effect is driven by a non-equilibrium population in momentum space, rather than in real space, so modes from both leads contribute -hence the need to sum over 2N rather than N modes.

III. MODEL HAMILTONIAN OF A WEYL SEMIMETAL
A simple model of a Weyl semimetal is given by the four-band Hamiltonian [3] The Pauli matrices σ j and τ j (j = x, y, z) act, respectively, on the spin and orbital degree of freedom. The momentum k varies over the Brillouin zone −π < k j < π of a simple cubic lattice (lattice constant a ≡ 1). The material is layered in the x-y plane, with nearest-neighbor hopping energies t (within the layer) and t z (along the zaxis). The primed terms t , t z indicate hopping with spinorbit coupling. Inversion symmetry, τ x H(−k)τ x = H(k), is broken by strain ∝ γ, while time-reversal symmetry, σ y H * (−k)σ y = H(k), is intrinsically broken by a magnetization ∝ β. Additionally, we may apply a magnetic field in the +z-direction, by substituting k y → k y − eBx/ . The field strength is characterized by the magnetic length l B = /eB. We confine the layers to a W x × W y lattice in the x-y plane, infinite in the z-direction so k z remains a good quantum number. The tight-binding Hamiltonian in this wire geometry is diagonalized with the help of the kwant toolbox [43], see Fig. 1. In zero magnetic field (panels a,c) there are two Weyl cones, gapped by the finite system size. The conical points (Weyl points) are separated along k z by approximately β/t z and they are separated in energy by approximately γ. The precise energy separation µ that governs the chiral magnetic effect was determined from the bandstructure in an infinite system, for our parameter values it differs from γ by a few percent.
As long as µ, M 0 β the Weyl cones remain distinct in an energy interval around E = 0. The Fermi velocity of the massless Weyl fermions is v F = t / in the plane of the layers and v F,z = t z / perpendicular to the layers. Surface states connect the Weyl cones across the Brillouin zone, forming the so-called Fermi arc. The arc states are chiral, spiraling along the wire with a velocity v arc,z = (µ/β)v F,z , as illustrated in Fig. 2. In a magnetic field (panels b,d in Fig. 1) Landau levels develop. The Weyl cones are pushed away from E = 0, but the zeroth Landau level closes the gap. Just like the Fermi arc, the zeroth Landau level propagates along the wire, in opposite direction for the two Weyl cones. Fermi arc spiraling along the wire (cross-sectional area A and perimeter P). Its flux sensitivity is set by the orbital magnetic moment evFA/P, while the number of surface modes at the Fermi level scales ∝ P, so their total contribution to the magnetic response is ∝ A -of the same order as the bulk contribution.

IV. INDUCED CURRENT IN LINEAR RESPONSE
A. Numerical results from the scattering formula We have calculated the current density δj flowing along the wire in response to a slowly varying µ or B. In the former case we fix B at l B = 25 and increase µ ≡ X from X 0 ≡ 0 to δX ≡ δµ, in the latter case we fix µ = 0.196 t and increase B ≡ X from X 0 ≡ 0 to δX ≡ δB. We obtain the CME coefficients in linear response, from the scattering formula (6), with T n ≡ 1 (no disorder, so unit transmission for all modes). The Fermi level is set at E F = 0. Results are shown in Fig. 3. We see that the numerical data points [44] lie close to the dashed lines given by The CME coefficient J µ agrees with the expected value from Eq. (1), while the CME coefficient J B has the opposite sign and is smaller by a factor of two. Inspection of the contributions from individual modes, plotted in Fig.  4, indicates that surface states are behind the different response, as we now explain in some detail.
B. Why surface Fermi arcs contribute to the magnetic response in the infinite-system limit Consider the propagating modes through a wire of diameter W . The number of surface modes scales ∝ W ,  while the number of bulk modes scales ∝ W 2 , so one might surmise that surface contributions to the current density I/W 2 can be neglected in the limit W → ∞. This is correct for J µ -but not for J B , because each surface mode individually contributes an amount ∝ W , so the total surface contribution scales ∝ W 2 , just like the bulk contribution.
To make this argument more precise, we consider the effective Hamiltonian of the surface Fermi arcs, with p s the component of the momentum along the perimeter of the wire in the x-y plane, of length P enclosing a flux Φ = BA in an area A. The energy spacing of the surface states at given momentum p z along the wire is δE = hv F /P, so an energy separation µ of the Weyl cones pushes N arc = µ/δE = µP/hv F surface modes through the Fermi level. Each contributes to the induced current, which is just its orbital magnetic moment. The total surface contribution takes on the universal value J B,arc = (e/h)N arc χ n /µA = (e/h) 2 .
The red dotted line in Fig. 4b confirms this reasoning.

C. Bulk Weyl cone contribution to the magnetic response
The numerical data in Fig. 3 indicates that the bulk Weyl cones contribute to the CME coefficient induced by a magnetic field, for a total J B,bulk + J B,arc = 1 2 (e/h) 2 . We have not found a simple intuitive argument for Eq. (14), but we do have an explicit analytical calculation, see App. A.
The difference between J µ and J B goes against the original expectation [6] that the low-frequency response to small variations in µ at fixed B should be the same as to small variations in B at fixed µ. That there is no such reciprocity was found recently in two studies [10,11] of currents induced by an oscillating magnetic field in an infinite isotropic system. Their bulk response has a different numerical coefficient than our Eq. (14) (1/3 instead of 1/2), possibly because of the intrinsic anisotropy of a wire geometry.

D. Interplay of surface Fermi arcs with bulk
Landau levels So far we have considered the magnetic response in the zero-magnetic field limit, when the bulk contribution arises from Weyl cones. We can also ask for the current density δj in response to a slow variation δX ≡ δB around some nonzero X 0 ≡ B 0 , all at fixed µ. As shown in Fig. 5, the magnetic response is the same whether we vary B around zero or nonzero B 0 . This is remarkable, because the bulk states are entirely different -Weyl cones versus Landau levels, compare the band structures in Figs. 1c and 1d. The individual modes also contribute FIG. 5: Same as Fig. 3, but for a larger range of widths Wx at fixed Wy = 5000 (periodic boundary conditions in the ydirection). The data for JB = µ −1 δj(B)/δB is shown at γ = 0.1 t ⇒ µ = 0.098 t in the limit B → 0 and at a large magnetic field in the Landau level regime (lB = 50). The data for Jµ = B −1 δj(µ)/δµ is shown at lB = 50 in the limit µ → 0 and for a large µ = 0.098 t.
FIG. 6: Same as Fig. 4b, but now at a large magnetic field corresponding to the bandstructure in Fig. 1d. The contribution from the surface Fermi arcs, close to A/P, goes to zero when they hybridize with the zeroth Landau level (for which ∂En/∂B = 0). very differently to J B , compare Figs. 4b and 6, and yet the net contribution is still close to 1 2 × (e/h) 2 . We have not succeeded in an analytical derivation of this numerical result.

V. FINITE-SIZE EFFECTS
We have seen that the surface Fermi arcs modify the magnetic response δj/δB even in the limit that the size of the system tends to infinity. The response δj/δµ to an energy displacement µ of the Weyl cones is unaffected by the surface states in the infinite-system limit, given by Eq. (1) in that limit. There are however finite-size effects from the surface state contributions, which we consider in this section.
As shown in Fig. 7, finite-size effects on J µ = B −1 δj/δµ are sensitive to whether or not the Fermi level E F is symmetrically arranged between the Weyl points (E F = 0 in Fig. 1). The earlier plots (Figs. 3 and 5) were for E F = 0, when finite-size effects are small. A variation of E F away from the symmetry point has little effect on the magnetic response J B = µ −1 δj/δB, provided |E F | |µ|. In contrast, the Fermi level displacement introduces a substantial size-dependence in J µ .
Inspection of the band structure in Fig. 1b shows that the degeneracy of the zeroth Landau level increases with increasing E F , because surface modes are converted into bulk modes, at a rate given by the inverse of the level spacing δE = hv F /P (cf. Sec. IV B). Each bulk mode contributes −(e/h)δµ/A to the induced current density δj, so we expect a finite-size correction to J µ = B −1 δj/δµ equal to −(e/h)(BA) −1 (E F /δE), hence As seen in Fig. 7, the slope of the E F dependence of J µ is accurately described by Eq. (15).

VI. CONCLUSION AND DISCUSSION OF DISORDER EFFECTS
Fig. 5 summarizes our main finding: It is known [6][7][8][9][10][11] that the chiral magnetic effect in a Weyl semimetal can be driven either by a slowly varying inversion-symmetry breaking µ or by a slowly varying magnetic field B. Contrary to the expectation from an infinite system, we find for a finite system that the induced current in the two cases has opposite sign. The difference is due to the surface Fermi arcs, but it is not a finite-size effect: The surface modes and the bulk modes give comparable contribution to the magnetic response no matter how large the system is, because the smaller number of surface modes is compensated by their stronger B-sensitivity.
This finding results from a scattering formulation of the chiral magnetic effect, that we have developed as an alternative to the established Kubo formulation [45,46]. Similarly to the Landauer formula for electrical conduction, the scattering formula (6) is ideally suited to describe finite and disordered systems, without translational invariance. Here we focused on the surface effects in a finite system, but in closing we briefly consider the disorder effects.
A qualitative prediction can be made without any calculation. In Eq. (6) disorder reduces the contribution from each mode n by its transmission probability T n . Assume that the disorder potential is smooth on the scale of the lattice constant a, so that it predominantly couples nearby modes in the Brillouin zone (with k n 's differing by much less than 1/a). This coupling can only lead to backscattering, reducing T n below unity, if it involves both left-moving and right-moving modes. Inspection of Fig. 4b shows that the surface modes are insensitive to backscattering, because they all move in the same direction along the wire, in contrast to the bulk Weyl cones. Disorder will therefore reduce the Weyl cone contribution J B,bulk = − 1 2 (e/h) 2 to the magnetic response, without affecting the arc state contribution J B,arc = (e/h) 2 . Since these contributions have opposite sign, we predict that disorder will increase the magnetically induced current.
For sufficiently strong disorder the bulk contribution to J B may be fully suppressed, leaving a B-induced current density equal to j = (e/h) 2 µB, carried entirely by the surface Fermi arc. This has the same topological origin as the zeroth Landau level that carries the µ-induced current (1) -both the Fermi arc and the zeroth Landau level connect Weyl cones of opposite chirality [35,36]. It has been argued [10,11] that the chiral magnetic effect produced by an oscillating B is fundamentally different from that produced by an oscillating µ, because the former lacks the topological protection that is the hallmark of the latter. By including surface conduction we can now offer an alternative perspective: Both the µ-response and the B-response are similarly protected by chirality, the difference is that one is a bulk current and the other a surface current.
From an experimental point of view, the inversionsymmetry breaking that sets µ is hardly adjustable, preventing a direct measurement of δj/δµ, while the magnetic field induced current δj/δB seems readily accessible. We note that Landau levels are not required for the B-response, so one can work with a nanowire of width small compared to the magnetic length. In such a quasione-dimensional system long-range impurity scattering may localize the bulk states, without significantly affect-ing the chiral surface states. One would be searching for an oscillating current I(ω) cos ωt along the wire in response to an oscillating parallel magnetic field. The frequency ω should be below µ and above the inelastic relaxation rate of the surface modes. The magnetic response is quasi-dc, showing a plateau in this ω-range that would distinguish it from any electrically induced ac current I(ω) ∝ ω.
A final word on nomenclature. The non-topological bulk contribution to the B-induced current has been termed the "gyrotropic magnetic effect" [11]. Because the B-induced current in the surface Fermi arc originates from the same chiral anomaly as the µ-induced current in the zeroth Landau level, we use the name "chiral magnetic effect" for both [47].
from the bulk Weyl cones to the magnetic response. We assume that the two Weyl cones are non-overlapping at the Fermi energy (as they are in Fig. 1c), so we can consider them separately.
A single Weyl cone has Hamiltonian (We have set ≡ 1 for ease of notation, but we will reinstate it at the end.) For generality, we allow for an anisotropic velocity (v x , v y , v z ). The modes propagating along the cylindrical wire have energy E n (k z ). We seek the magnetic moment ∂E n (k z )/∂B for an infinitesimal magnetic field B in the z-direction (along the axis of the cylinder).
For sufficiently large transverse dimensions W x , W y the boundary conditions should be irrelevant for the bulk response, and we use this freedom to simplify the calculation. To isolate the bulk contribution we prefer a boundary condition that does not bind a surface state.
In the y-direction we impose periodicity, so that k y is a good quantum number. The system is then represented by a hollow cylinder of circumference W y , with an inner and an outer surface at x = 0 and x = W x . We can use periodic or antiperiodic boundary conditions, k y = 2πn/W y or k y = 2π(n + 1 2 )/W y , n = 0, ±1, ±2, . . . , in the large-W y limit it makes no difference.
In the x-direction we choose a zero-current boundary condition. A simple choice is to take the spinor ψ(x, y) as an eigenfunction of σ y at the two surfaces x = 0 and for arbitrary complex functions f (y), g(y). This boundary condition corresponds to confinement by a mass term ∝ σ z of infinite magnitude and opposite sign at the two surfaces. No surface state is produced by mass confinement. For k z = 0 the sign change of the mass term does produce a spurious chiral state ψ = e ikyy 1 i , E = v y k y , which carries no current in the z-direction and can therefore be ignored.
The solution of the eigenvalue equation H Weyl ψ = Eψ that satisfies the boundary condition (A4) is given by The band structure E nm (k z ) is determined by the dispersion relation Normalization ψ|ψ = 1 gives The magnetic response is induced by the vector potential A y = B(x + X 0 ), with an offset X 0 that accounts for the flux BW y X 0 enclosed by the inner surface of the cylinder. The magnetic moment results from The second term is the magnetic moment of a charge e circulating along the inner surface of the cylinder with velocity ∂E/∂k y = v 2 y k y /E. It drops out when we sum the contributions from +k y and −k y , producing the magnetic moment plotted in Fig. 8. We fix the energy E nm = E, adjusting k z accordingly for each n and m. Both +k z and −k z satisfy the dispersion relation (A7). We consider separately the sum Wx, so the discrete modes merge into a continuous curve.
over the magnetic moment of the modes with k z > 0 and k z < 0, so that we can distinguish left-movers from right-movers in the scattering formula (6). In the large-W limit the sum over k x and k y , quantized by Eqs. (A3) and (A6), can be replaced by an integration over the k xk y plane, with θ(s) the unit step function. The prime in the summation indicates that we include only half of the modes, with a given sign of k z . The integral over Eq. (A9) is readily evaluated in polar coordinates, The quantity χ nm that determines the magnetically induced current δI/δB according to Eq. (6) is the magnetic moment ∂E nm /∂B times the sign of the velocity ∂E nm /∂k z in the z-direction. The sign of the velocity in a single Weyl cone (with Weyl point at k = 0, E = 0) equals the product of the sign of k z and the sign of E nm , hence In the last equation we have reinstated the that we had previously set to unity. There is no prime in the summation because now both signs of k z are included.
We conclude that the contribution to the induced current density δj = δI/W x W y from a single Weyl cone at energy E = µ/2 is δj = e h δB W x W y nm χ nm = − 1 4 (e/h) 2 µδB.
The other Weyl cone contributes the same amount, for a total CME coefficient given by Eq. (A1).