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

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–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–11]. This is the chiral magnetic effect (CME) from particle physics [12–15], see [16–18] for recent reviews in the condensed matter setting. In an infinite system the current density has the universal form [19, 20] ³

$$\begin{eqnarray}&&{j}_{0}=-{(e/h)}^{2}\mu B,\end{eqnarray} \tag{ 1 }$$

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 [21–24] have boosted the search for the CME [25–33]. 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 [34, 35]. Previous studies [36, 37] have pointed to the importance of boundaries for the CME—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 CME 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 [38–40]. The CME 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 equation (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 CME.

The outline of the paper is as follows. In the next section 2 we derive the scattering formula for the CME, in a general setting. In section 4 we apply it to the model Hamiltonian of a Weyl semimetal from [3], summarized in section 3. 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 section 5. We conclude in section 6 with a summary and a discussion of the robustness of the results against disorder scattering.

For a scattering theory of the CME 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 index $n=1,2,\ldots$ and dependent on the wave vector k along the lead. At a given energy ε (measured relative to the equilibrium Fermi level ${E}_{{\rm{F}}}$), each incident mode has wave vector ${k}_{n}(\varepsilon )$ and carries the same current e/h per unit energy interval⁴ .

The scattering matrix $S(\varepsilon )$ 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\times 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 ${ \mathcal P }=\displaystyle \left(\genfrac{}{}{0em}{}{1\;0}{0\;0}\right)$, where each sub-block is an N × N matrix. The current is driven by a set of nonequilibrium occupation numbers $\delta {f}_{n}(\varepsilon )$, with $n=1,2,\ldots N$ for the left lead and $n=N+1,N+2,\ldots 2N$ for the right lead. We collect these numbers in a $2N\times 2N$ diagonal matrix $\delta { \mathcal F }(\varepsilon )$. The net current in the left lead is then given by the difference of incoming and outgoing currents

$$\begin{eqnarray}&&I=\displaystyle \frac{e}{h}\int {\rm{d}}\varepsilon \;{\rm{Tr}}[{ \mathcal P }\delta { \mathcal F }(\varepsilon )-{ \mathcal P }S(\varepsilon )\delta { \mathcal F }(\varepsilon ){S}^{\dagger }(\varepsilon )].\end{eqnarray} \tag{ 2 }$$

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}+\delta X)$. To first order in the perturbation $\delta X$ the energy shift at constant k is

$$\begin{eqnarray}&&{E}_{n}(k| {X}_{0}+\delta X)-{E}_{n}(k| {X}_{0})=\delta X\underset{X\to {X}_{0}}{\mathrm{lim}}\displaystyle \frac{\partial }{\partial X}{E}_{n}(k| X).\end{eqnarray} \tag{ 3 }$$

The corresponding deviation of the occupation number from the equilibrium Fermi function ${f}_{{\rm{eq}}}{(\varepsilon )=(1+{e}^{\varepsilon /{k}_{{\rm{B}}}T})}^{-1}$ is

$$\begin{eqnarray}\delta {f}_{n}(\varepsilon ) & = & {f}_{{\rm{eq}}}({E}_{n}({k}_{n}(\varepsilon )| {X}_{0}))-{f}_{{\rm{eq}}}(\varepsilon )\\ & = & -\delta {{Xf}}_{{\rm{eq}}}^{\prime }(\varepsilon )\underset{X\to {X}_{0}}{\mathrm{lim}}\displaystyle \frac{\partial }{\partial X}{E}_{n}({k}_{n}(\varepsilon )| X),\end{eqnarray} \tag{ 4 }$$

where we have used that ${E}_{n}({k}_{n}(\varepsilon )| {X}_{0}+\delta X)\equiv \varepsilon$.

At zero temperature the derivative ${f}_{{\rm{eq}}}^{\prime }(\varepsilon )\to -\delta (\varepsilon )$, 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

$$\begin{eqnarray}{T}_{n}=\left\{\begin{array}{cc}{\displaystyle \sum }_{m=N+1}^{2N}| {S}_{{mn}}{| }^{2} & {\rm{for}}\ \ 1\leqslant n\leqslant N,\\ {\displaystyle \sum }_{m=1}^{N}| {S}_{{mn}}{| }^{2} & {\rm{for}}\ \ N+1\leqslant n\leqslant 2N,\end{array}\right.\end{eqnarray} \tag{ 5 }$$

evaluated at $\varepsilon =0$. (The two cases correspond to transmission from left to right or from right to left.) Since ${\sum }_{m=1}^{2N}| {S}_{{mn}}{| }^{2}=1$ because of unitarity, we have

$$\begin{eqnarray}I & = & \displaystyle \frac{e}{h}\delta X\displaystyle \sum _{n=1}^{2N}{\chi }_{n}{T}_{n},\\ {\chi }_{n} & = & \underset{k\to {k}_{n}}{\mathrm{lim}}\;\underset{X\to {X}_{0}}{\mathrm{lim}}\displaystyle \frac{\partial {E}_{n}(k| X)}{\partial X}\times {\rm{sign}}\;\displaystyle \frac{\partial {E}_{n}(k| X)}{\partial k}.\end{eqnarray} \tag{ 6 }$$

The sign of the derivative $\partial {E}_{n}/\partial k$ distinguishes the right-moving modes $n=1,2,\ldots N$ from the left-moving modes $n=N+1,N+2,\ldots 2N$.

The Landauer conductance formula [38–40]

$$\begin{eqnarray}&&G=I/V=({e}^{2}/h)\displaystyle \sum _{n=1}^{N}{T}_{n}\end{eqnarray} \tag{ 7 }$$

is obtained from equation (6) if we identify $\delta X=V$ with the voltage difference between the left and right lead, and then set ${\chi }_{n}=1$ for $n=1,2,\ldots N$ and ${\chi }_{n}=0$ for $n=N+1,N+2,\ldots 2N$. The CME is driven by a nonequilibrium 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.

A simple model of a Weyl semimetal is given by the four-band Hamiltonian [3]

$$\begin{eqnarray}H({\boldsymbol{k}}) & = & {t}^{\prime }{\tau }_{z}({\sigma }_{x}\mathrm{sin}{k}_{x}+{\sigma }_{y}\mathrm{sin}{k}_{y})+{t}_{z}^{\prime }{\tau }_{y}\mathrm{sin}{k}_{z}\\ & & +M({\boldsymbol{k}}){\tau }_{x}+\tfrac{1}{2}\gamma {\tau }_{y}{\sigma }_{z}+\tfrac{1}{2}\beta {\sigma }_{z},\\ M({\boldsymbol{k}}) & = & {M}_{0}+t(2-\mathrm{cos}{k}_{x}-\mathrm{cos}{k}_{y})+{t}_{z}(1-\mathrm{cos}{k}_{z}).\end{eqnarray} \tag{ 8 }$$

The Pauli matrices ${\sigma }_{j}$ and ${\tau }_{j}$ ($j=x,y,z$) act, respectively, on the spin and orbital degree of freedom. The momentum ${\boldsymbol{k}}$ varies over the Brillouin zone $-\pi \lt {k}_{j}\lt \pi$ of a simple cubic lattice (lattice constant $a\equiv 1$). The material is layered in the x–y plane, with nearest-neighbor hopping energies t (within the layer) and t_z (along the z-axis). The primed terms ${t}^{\prime },{t}_{z}^{\prime }$ indicate hopping with spin–orbit coupling. Inversion symmetry, ${\tau }_{x}H(-{\boldsymbol{k}}){\tau }_{x}=H({\boldsymbol{k}})$, is broken by strain $\propto \gamma$, while time-reversal symmetry, ${\sigma }_{y}{H}^{*}(-{\boldsymbol{k}}){\sigma }_{y}=H({\boldsymbol{k}})$, is intrinsically broken by a magnetization $\propto \beta$. Additionally, we may apply a magnetic field in the $+z$-direction, by substituting ${k}_{y}\mapsto {k}_{y}-{eBx}/{\hslash }$. The field strength is characterized by the magnetic length ${l}_{B}=\sqrt{{\hslash }/{eB}}$.

We confine the layers to a ${W}_{x}\times {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 [41], see figure 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 $\beta /{t}_{z}^{\prime }$ and they are separated in energy by approximately γ. The precise energy separation μ that governs the CME was determined from the bandstructure in an infinite system, for our parameter values it differs from γ by a few percent.

Zoom In Zoom Out Reset image size

As long as $\mu ,{M}_{0}\ll \beta$ the Weyl cones remain distinct in an energy interval around E = 0. The Fermi velocity of the massless Weyl fermions is ${v}_{{\rm{F}}}={t}^{\prime }/{\hslash }$ in the plane of the layers and ${v}_{{\rm{F}},z}={t}_{z}^{\prime }/{\hslash }$ 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}_{{\rm{arc}},z}=(\mu /\beta ){v}_{{\rm{F}},z}$, as illustrated in figure 2. In a magnetic field (panels (b), (d) in figure 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.

Zoom In Zoom Out Reset image size

4.1. Numerical results from the scattering formula

We have calculated the current density $\delta 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 $\mu \equiv X$ from ${X}_{0}\equiv 0$ to $\delta X\equiv \delta \mu$, in the latter case we fix $\mu =0.196\;t$ and increase $B\equiv X$ from ${X}_{0}\equiv 0$ to $\delta X\equiv \delta B$. We obtain the CME coefficients in linear response

$$\begin{eqnarray}&&{{ \mathcal J }}_{\mu }\equiv {B}^{-1}\delta j/\delta \mu ,\ \ {{ \mathcal J }}_{B}={\mu }^{-1}\delta j/\delta B,\end{eqnarray} \tag{ 9 }$$

from the scattering formula (6), with ${T}_{n}\equiv 1$ (no disorder, so unit transmission for all modes). The Fermi level is set at ${E}_{{\rm{F}}}=0$. Results are shown in figure 3. We see that the numerical data points⁵ lie close to the dashed lines given by

$$\begin{eqnarray}&&{{ \mathcal J }}_{\mu }=-{(e/h)}^{2},\ \ {{ \mathcal J }}_{B}=\tfrac{1}{2}\times {(e/h)}^{2}.\end{eqnarray} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

The CME coefficient ${{ \mathcal J }}_{\mu }$ agrees with the expected value from equation (1), while the CME coefficient ${{ \mathcal J }}_{B}$ has the opposite sign and is smaller by a factor of two. Inspection of the contributions from individual modes, plotted in figure 4, indicates that surface states are behind the different response, as we now explain in some detail.

Zoom In Zoom Out Reset image size

4.2. 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 $\propto W$, while the number of bulk modes scales ∝W², so one might surmise that surface contributions to the current density $I/{W}^{2}$ can be neglected in the limit $W\to \infty$. This is correct for ${{ \mathcal J }}_{\mu }$—but not for ${{ \mathcal J }}_{B}$, because each surface mode individually contributes an amount $\propto W$, so the total surface contribution scales ∝W², just like the bulk contribution.

To make this argument more precise, we consider the effective Hamiltonian of the surface Fermi arcs

$$\begin{eqnarray}&&{H}_{{\rm{arc}}}={v}_{{\rm{F}}}({p}_{s}-e{\rm{\Phi }}/{ \mathcal P })-{v}_{{\rm{arc}},z}{p}_{z},\end{eqnarray} \tag{ 11 }$$

with p_s the component of the momentum along the perimeter of the wire in the x–y plane, of length ${ \mathcal P }$ enclosing a flux ${\rm{\Phi }}=B{ \mathcal A }$ in an area ${ \mathcal A }$. The energy spacing of the surface states at given momentum p_z along the wire is $\delta E={{hv}}_{{\rm{F}}}/{ \mathcal P }$, so an energy separation μ of the Weyl cones pushes ${N}_{{\rm{arc}}}=\mu /\delta E=\mu { \mathcal P }/{{hv}}_{{\rm{F}}}$ surface modes through the Fermi level. Each contributes

$$\begin{eqnarray}&&{\chi }_{n}={\rm{sign}}(\partial {E}_{n}/\partial {p}_{z})\times \partial {E}_{n}/\partial B={{ev}}_{{\rm{F}}}{ \mathcal A }/{ \mathcal P }\end{eqnarray} \tag{ 12 }$$

to the induced current, which is just its orbital magnetic moment. The total surface contribution takes on the universal value

$$\begin{eqnarray}&&{{ \mathcal J }}_{B,{\rm{arc}}}={(e/h){N}_{{\rm{arc}}}{\chi }_{n}/\mu { \mathcal A }=(e/h)}^{2}.\end{eqnarray} \tag{ 13 }$$

The red dotted line in figure 4(b) confirms this reasoning.

4.3. Bulk Weyl cone contribution to the magnetic response

The numerical data in figure 3 indicates that the bulk Weyl cones contribute

$$\begin{eqnarray}&&{{ \mathcal J }}_{B,{\rm{bulk}}}=-\tfrac{1}{2}{(e/h)}^{2}\end{eqnarray} \tag{ 14 }$$

to the CME coefficient induced by a magnetic field, for a total ${{ \mathcal J }}_{B,{\rm{bulk}}}+{{ \mathcal J }}_{B,{\rm{arc}}}=\tfrac{1}{2}{(e/h)}^{2}$. We have not found a simple intuitive argument for equation (14), but we do have an explicit analytical calculation, see appendix.

The difference between ${{ \mathcal J }}_{\mu }$ and ${{ \mathcal 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 equation (14) (1/3 instead of 1/2), possibly because of the intrinsic anisotropy of a wire geometry.

4.4. 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 $\delta j$ in response to a slow variation $\delta X\equiv \delta B$ around some nonzero ${X}_{0}\equiv {B}_{0}$, all at fixed μ. As shown in figure 5, the magnetic response is the same whether we vary B around zero or nonzero B₀. This is remarkable, because the bulk states are entirely different—Weyl cones versus Landau levels, compare the band structures in figures 1(c) and (d). The individual modes also contribute very differently to ${{ \mathcal J }}_{B}$, compare figures 4(b) and 6, and yet the net contribution is still close to $\tfrac{1}{2}\times {(e/h)}^{2}$. We have not succeeded in an analytical derivation of this numerical result.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have seen that the surface Fermi arcs modify the magnetic response $\delta j/\delta B$ even in the limit that the size of the system tends to infinity. The response $\delta j/\delta \mu$ to an energy displacement μ of the Weyl cones is unaffected by the surface states in the infinite-system limit, given by equation (1) in that limit. There are however finite-size effects from the surface state contributions, which we consider in this section.

As shown in figure 7, finite-size effects on ${{ \mathcal J }}_{\mu }={B}^{-1}\delta j/\delta \mu$ are sensitive to whether or not the Fermi level ${E}_{{\rm{F}}}$ is symmetrically arranged between the Weyl points (${E}_{{\rm{F}}}=0$ in figure 1). The earlier plots (figures 3 and 5) were for ${E}_{{\rm{F}}}=0$, when finite-size effects are small. A variation of ${E}_{{\rm{F}}}$ away from the symmetry point has little effect on the magnetic response ${{ \mathcal J }}_{B}={\mu }^{-1}\delta j/\delta B$, provided $| {E}_{{\rm{F}}}| \lesssim | \mu |$. In contrast, the Fermi level displacement introduces a substantial size-dependence in ${{ \mathcal J }}_{\mu }$.

Zoom In Zoom Out Reset image size

Inspection of the band structure in figure 1(b) shows that the degeneracy of the zeroth Landau level increases with increasing ${E}_{{\rm{F}}}$, because surface modes are converted into bulk modes, at a rate given by the inverse of the level spacing $\delta E={{hv}}_{{\rm{F}}}/{ \mathcal P }$ (see section 4.2). Each bulk mode contributes $-(e/h)\delta \mu /{ \mathcal A }$ to the induced current density $\delta j$, so we expect a finite-size correction to ${{ \mathcal J }}_{\mu }={B}^{-1}\delta j/\delta \mu$ equal to $-{(e/h)(B{ \mathcal A })}^{-1}({E}_{{\rm{F}}}/\delta E)$, hence

$$\begin{eqnarray}&&{{ \mathcal J }}_{\mu }=-{(e/h)}^{2}\left(1+\displaystyle \frac{{ \mathcal P }}{{ \mathcal A }}\displaystyle \frac{{E}_{{\rm{F}}}}{{{eBv}}_{{\rm{F}}}}\right).\end{eqnarray} \tag{ 15 }$$

As seen in figure 7, the slope of the ${E}_{{\rm{F}}}$ dependence of ${J}_{\mu }$ is accurately described by equation (15).

Figure 5 summarizes our main finding: it is known [6–11] that the CME 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 CME, that we have developed as an alternative to the established Kubo formulation [42, 43]. 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 equation (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 figure 4(b) 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 ${{ \mathcal J }}_{B,{\rm{bulk}}}=-\tfrac{1}{2}{(e/h)}^{2}$ to the magnetic response, without affecting the arc state contribution ${{ \mathcal J }}_{B,{\rm{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 ${{ \mathcal J }}_{B}$ may be fully suppressed, leaving a B-induced current density equal to $j={(e/h)}^{2}\mu 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 [34, 35]. It has been argued [10, 11] that the CME 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 inversion-symmetry breaking that sets μ is hardly adjustable, preventing a direct measurement of $\delta j/\delta \mu$, while the magnetic field induced current $\delta j/\delta 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 quasi-one-dimensional system long-range impurity scattering may localize the bulk states, without significantly affecting the chiral surface states. One would be searching for an oscillating current $I(\omega )\mathrm{cos}\omega 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(\omega )\propto \omega$.

A final word on nomenclature. The nontopological 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 'CME' for both⁶ .

We have benefited from discussions with I Adagideli, T E O'Brien, V V Cheianov, and B Tarasinski. This research was supported by the Foundation for Fundamental Research on Matter (FOM), the Netherlands Organization for Scientific Research (NWO/OCW), and an ERC Synergy Grant.

We wish to derive the result (14) for the contribution

$$\begin{eqnarray}&&{\mu }^{-1}\delta {j}_{{\rm{bulk}}}/\delta B=-\tfrac{1}{2}{(e/h)}^{2}\end{eqnarray} \tag{ A1 }$$

from the bulk Weyl cones to the magnetic response. We assume that the two Weyl cones are nonoverlapping at the Fermi energy (as they are in figure 1(c)), so we can consider them separately.

A single Weyl cone has Hamiltonian

$$\begin{eqnarray}&&{H}_{{\rm{Weyl}}}={v}_{x}{k}_{x}{\sigma }_{x}+{v}_{y}{k}_{y}{\sigma }_{y}+{v}_{z}{k}_{z}{\sigma }_{z}.\end{eqnarray} \tag{ A2 }$$

(We have set ${\hslash }\equiv 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 $\partial {E}_{n}({k}_{z})/\partial 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

$$\begin{eqnarray}&&{k}_{y}=2\pi n/{W}_{y}\ \ {\rm{or}}\ \ {k}_{y}=2\pi \left(n+\displaystyle \frac{1}{2}\right)/{W}_{y},\\ &&\qquad \qquad n=0,\pm 1,\pm 2,\ldots ,\end{eqnarray} \tag{ A3 }$$

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 $\psi (x,y)$ as an eigenfunction of ${\sigma }_{y}$ at the two surfaces x = 0 and $x={W}_{x}$,

$$\begin{eqnarray}&&\underset{x\to 0}{\mathrm{lim}}\psi (x,y)=f(y)\left(\begin{array}{c}1\\ i\end{array}\right),\ \ \underset{x\to {W}_{x}}{\mathrm{lim}}\psi (x,y)=g(y)\left(\begin{array}{c}1\\ i\end{array}\right),\end{eqnarray} \tag{ A4 }$$

for arbitrary complex functions f (y), g (y). This boundary condition corresponds to confinement by a mass term $\propto {\sigma }_{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 $\psi ={{\rm{e}}}^{{\rm{i}}{k}_{y}y}\displaystyle \left(\genfrac{}{}{0em}{}{1}{i}\right)$, $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}_{{\rm{Weyl}}}\psi =E\psi$ that satisfies the boundary condition (A4) is given by

$$\begin{eqnarray}\psi (x,y) & = & \displaystyle \frac{1}{Z}{{\rm{e}}}^{{\rm{i}}{k}_{y}y}\left[(E-{\rm{i}}{v}_{x}{k}_{x}-{v}_{y}{k}_{y}+{v}_{z}{k}_{z})\left(\begin{array}{c}E+{v}_{z}{k}_{z}\\ {v}_{x}{k}_{x}+{\rm{i}}{v}_{y}{k}_{y}\end{array}\right){{\rm{e}}}^{{\rm{i}}{k}_{x}x}\right.\\ & & \qquad \quad -\left.(E+{\rm{i}}{v}_{x}{k}_{x}-{v}_{y}{k}_{y}+{v}_{z}{k}_{z})\left(\begin{array}{c}E+{v}_{z}{k}_{z}\\ -{v}_{x}{k}_{x}+{\rm{i}}{v}_{y}{k}_{y}\end{array}\right){{\rm{e}}}^{-{\rm{i}}{k}_{x}x}\right],\end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray}&&{k}_{x}=\pi m/{W}_{x},\ \ m=1,2,3,\ldots .\end{eqnarray} \tag{ A6 }$$

The band structure ${E}_{{nm}}({k}_{z})$ is determined by the dispersion relation

$$\begin{eqnarray}&&{E}^{2}={(\pi {{mv}}_{x}/{W}_{x})}^{2}+{(2\pi {{nv}}_{y}/{W}_{y})}^{2}+{({v}_{z}{k}_{z})}^{2}.\end{eqnarray} \tag{ A7 }$$

Normalization $\langle \psi | \psi \rangle =1$ gives

$$\begin{eqnarray}&&{Z}^{2}=8{W}_{x}{W}_{y}E{(E-{v}_{y}{k}_{y})(E+{v}_{z}{k}_{z})}^{2}.\end{eqnarray} \tag{ A8 }$$

The magnetic response is induced by the vector potential ${A}_{y}=B(x+{X}_{0})$, with an offset X₀ that accounts for the flux ${{BW}}_{y}{X}_{0}$ enclosed by the inner surface of the cylinder. The magnetic moment results from

$$\begin{eqnarray}\displaystyle \frac{\partial E}{\partial B} & = & \langle \psi | \partial H/\partial B| \psi \rangle =-{{ev}}_{y}\langle \psi | (x+{X}_{0}){\sigma }_{y}| \psi \rangle \\ & = & -\displaystyle \frac{{{ev}}_{x}{v}_{y}{v}_{z}{k}_{z}}{2E(E-{v}_{y}{k}_{y})}-\displaystyle \frac{{{ev}}_{y}^{2}{k}_{y}}{E}({X}_{0}+W/2).\end{eqnarray} \tag{ A9 }$$

The second term is the magnetic moment of a charge e circulating along the inner surface of the cylinder with velocity $\partial E/\partial {k}_{y}={v}_{y}^{2}{k}_{y}/E$. It drops out when we sum the contributions from $+{k}_{y}$ and $-{k}_{y}$, producing the magnetic moment

$$\begin{eqnarray}&&\displaystyle \sum _{\pm {k}_{y}}\displaystyle \frac{\partial E}{\partial B}=-\displaystyle \frac{{{ev}}_{x}{v}_{y}{v}_{z}{k}_{z}}{{E}^{2}-{v}_{y}^{2}{k}_{y}^{2}}=-\displaystyle \frac{{{ev}}_{x}{v}_{y}{v}_{z}{k}_{z}}{{v}_{x}^{2}{k}_{x}^{2}+{v}_{z}^{2}{k}_{z}^{2}}\end{eqnarray} \tag{ A10 }$$

plotted in figure 8.

Zoom In Zoom Out Reset image size

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 over the magnetic moment of the modes with ${k}_{z}\gt 0$ and ${k}_{z}\lt 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 equations (A3) and (A6), can be replaced by an integration over the k_x–k_y plane

$$\begin{eqnarray}&&{\displaystyle \sum _{{nm}}}^{\prime }\mapsto \displaystyle \frac{{W}_{x}{W}_{y}}{2{\pi }^{2}}{\int }_{0}^{\infty }{\rm{d}}{k}_{x}{\int }_{-\infty }^{\infty }{\rm{d}}{k}_{y}\;\theta ({E}^{2}-{v}_{x}^{2}{k}_{x}^{2}-{v}_{y}^{2}{k}_{y}^{2}),\end{eqnarray} \tag{ A11 }$$

with $\theta (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 equation (A9) is readily evaluated in polar coordinates

$$\begin{eqnarray}{\displaystyle \sum _{{nm}}}^{\prime }{\left.\displaystyle \frac{\partial {E}_{{nm}}}{\partial B}\right|}_{{E}_{{nm}}=E} & = & -({\rm{sign}}\;{k}_{z})\displaystyle \frac{{{eW}}_{x}{W}_{y}}{4{\pi }^{2}E}\\ & & \times {\displaystyle \int }_{-\pi /2}^{\pi /2}{\rm{d}}\phi {\displaystyle \int }_{0}^{| E| }r{\rm{d}}r\;\displaystyle \frac{\sqrt{{E}^{2}-{r}^{2}}}{E-r\mathrm{sin}\phi }\\ & = & -({\rm{sign}}\;{k}_{z})\displaystyle \frac{{{eW}}_{x}{W}_{y}| E| }{8\pi }.\end{eqnarray} \tag{ A12 }$$

The quantity ${\chi }_{{nm}}$ that determines the magnetically induced current $\delta I/\delta B$ according to equation (6) is the magnetic moment $\partial {E}_{{nm}}/\partial B$ times the sign of the velocity $\partial {E}_{{nm}}/\partial {k}_{z}$ in the z-direction. The sign of the velocity in a single Weyl cone (with Weyl point at ${\boldsymbol{k}}=0$, E = 0) equals the product of the sign of k_z and the sign of E_nm, hence

$$\begin{eqnarray}\displaystyle \sum _{{nm}}{\chi }_{{nm}} & = & \displaystyle \sum _{{nm}}\left({\rm{sign}}\;\displaystyle \frac{\partial {E}_{{nm}}}{\partial {k}_{z}}\right){\left.\displaystyle \frac{\partial {E}_{{nm}}}{\partial B}\right|}_{{E}_{{nm}}=E}\\ & = & -\displaystyle \frac{{{eW}}_{x}{W}_{y}E}{4\pi {\hslash }}.\end{eqnarray} \tag{ A13 }$$

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 $\delta j=\delta I/{W}_{x}{W}_{y}$ from a single Weyl cone at energy $E=\mu /2$ is

$$\begin{eqnarray}&&\delta j=\displaystyle \frac{e}{h}\displaystyle \frac{\delta B}{{W}_{x}{W}_{y}}\displaystyle \sum _{{nm}}{\chi }_{{nm}}=-\tfrac{1}{4}{(e/h)}^{2}\mu \delta B.\end{eqnarray} \tag{ A14 }$$

The other Weyl cone contributes the same amount, for a total CME coefficient given by equation (A1).