SciPost Phys. 11, 046 (2021)
Chiral anomaly trapped in Weyl metals:
Nonequilibrium valley polarization at zero magnetic field
Pablo M. Perez-Piskunow1? , Nicandro Bovenzi2† ,
Anton R. Akhmerov3‡ and Maxim Breitkreiz4◦
1 Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST,
Campus UAB, Bellaterra, 08193 Barcelona, Spain
2 Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands
3 Kavli Institute of Nanoscience, Delft University of Technology,
P.O. Box 4056, 2600 GA Delft, The Netherlands
4 Dahlem Center for Complex Quantum Systems and Fachbereich Physik,
Freie Universität Berlin, 14195 Berlin, Germany
Abstract
In Weyl semimetals, the application of parallel electric and magnetic fields leads to valley
polarization—an occupation disbalance of valleys of opposite chirality—a direct conse-
quence of the chiral anomaly. In this work, we present numerical tools to explore such
nonequilibrium effects in spatially confined three-dimensional systems with a variable
disorder potential, giving exact solutions to leading order in the disorder potential and
the applied electric field. Application to a Weyl-metal slab shows that valley polarization
also occurs without an external magnetic field as an effect of chiral anomaly “trapping”:
Spatial confinement produces chiral bulk states, which enable the valley polarization
in a similar way as the chiral states induced by a magnetic field. Despite its finite-size
origin, the valley polarization can persist up to macroscopic length scales if the disorder
potential is sufficiently long ranged, so that direct inter-valley scattering is suppressed
and the relaxation then goes via the Fermi-arc surface states.
Copyright Perez Piskunow et al. Received 21-04-2021
This work is licensed under the Creative Commons Accepted 16-08-2021
Check for
Attribution 4.0 International License. Published 31-08-2021 updates
Published by the SciPost Foundation. doi:10.21468/SciPostPhys.11.2.046
Contents
1 Introduction 1
2 Quantum transport approach 3
2.1 General quantum transport approach 3
2.2 Application to a slab model 5
3 Weyl-semimetal slab model 6
4 Scattering 7
4.1 Disorder potential 7
4.2 Scattering lengths 9
1
SciPost Phys. 11, 046 (2021)
5 Numerical results 10
5.1 Conductivity 10
5.2 Valley polarization 11
6 Discussion 13
6.1 Impurity-range dependence 13
6.2 Valley polarization 14
6.3 Conductivity in case of a large number of bulk states 15
6.4 Conductivity in case of a small number of bulk states 16
7 Conclusion 17
References 19
1 Introduction
The most famous effect associated with Weyl Fermions is the chiral anomaly [1, 2]—magnetic-
field induced chiral states moving parallel or antiparallel to the field, depending on the chirality
of the Weyl Fermion. In Weyl semimetals [3, 4] the two chiralities occur pairwise, ensuring an
equal number of forward- and backward-propagating states, and the chiralities are connected
by Fermi-arc surface states. Existing Weyl-semimetal materials typically have a small but finite
Fermi momentum k F measured from the Weyl node and a much larger momentum-space
separation ∆k of valleys that host the opposite chiralities.
The valley degree of freedom [5] plays a central role in the transport behavior of Weyl
semimetals. Parallel electric and magnetic fields produce a difference in the non-equilibrium
occupation of valleys [3, 6, 7]. A direct consequence of this valley polarization is an enhanced
conductivity parallel to the magnetic field due to the polarization-enhanced occupation disbal-
ance of countermoving chiral states [8–13]. Experimental observations, although obscured
by the competing current-jetting effect [14], support the general feature of a chiral-anomaly
enhancement of the conductivity [15–17]. Other manifestations of the valley degree of freedom
are found in nonlocal transport measurements [6,18] and in the photogalvanic response [19,20].
Crucial in understanding “valleytronic” transport is to explore the effect of disorder and
the finite size of the crystal, which are two unaviodable properties of real materials. Disorder
plays a subtle role if the Fermi level lies at the Weyl nodes, where it may or may not destroy the
ideal semimetal phase by inducing a finite density of states [21–24] or, for finite inter-valley
scattering, drive the system into an insulating phase [25]. At finite chemical potentials, well-
separated Weyl nodes, and a weak disorder potential the Weyl-semimetal phase has proven to
be robust, allowing for a perturbative treatment of disorder, which will be employed in this
work.
Finite-size effects in meso- and macroscopic Weyl semimetals (crystal dimensions much
larger than the lattice constant), have also been explored with the focus on the role of topological
Fermi-arc surface states [26–36]. Peculiarities are rooted in the specifics of the momentum-space
structure of Fermi arcs connecting valleys of opposite chirality and their unidirectional motion
at a single surface, see Fig. 1. Separation of countermoving Fermi arcs to opposite surfaces
explains their relevance at large systems sizes, most prominently in the intrinsic anomalous
Hall effect [34, 37–39]. The relevance of finite-size effects for the valley degree of freedom, on
the other hand, is much less obvious, since the valleys consist of extended bulk states, lacking
2
SciPost Phys. 11, 046 (2021)
Fermi arc
(a) upper surface z = W/2
(b) chiral bulk states
ky 2kF
AAAB7HicbVBNS8NAEJ3Ur1q/qh69LBbBU0mqoMeiIB4rmLbQhrLZbtqlm03YnQil9Dd48aCIV3+QN/+N2zYHbX0w8Hhvhpl5YSqFQdf9dgpr6xubW8Xt0s7u3v5B+fCoaZJMM+6zRCa6HVLDpVDcR4GSt1PNaRxK3gpHtzO/9cS1EYl6xHHKg5gOlIgEo2glv0ZGvbteueJW3TnIKvFyUoEcjV75q9tPWBZzhUxSYzqem2IwoRoFk3xa6maGp5SN6IB3LFU05iaYzI+dkjOr9EmUaFsKyVz9PTGhsTHjOLSdMcWhWfZm4n9eJ8PoOpgIlWbIFVssijJJMCGzz0lfaM5Qji2hTAt7K2FDqilDm0/JhuAtv7xKmrWqd1GtPVxW6jd5HEU4gVM4Bw+uoA730AAfGAh4hld4c5Tz4rw7H4vWgpPPHMMfOJ8/4EaOEg==
kx
2AAAB7nicbVDLSgNBEOz1GeMr6tHLYBA8hd0o6DHoxWME84BkCbOTTjJk9sFMrxCWfIQXD4p49Xu8+TdOkj1oYkFDUdVNd1eQKGnIdb+dtfWNza3twk5xd2//4LB0dNw0caoFNkSsYt0OuEElI2yQJIXtRCMPA4WtYHw381tPqI2Mo0eaJOiHfBjJgRScrNSqsm6AxHulsltx52CrxMtJGXLUe6Wvbj8WaYgRCcWN6XhuQn7GNUmhcFrspgYTLsZ8iB1LIx6i8bP5uVN2bpU+G8TaVkRsrv6eyHhozCQMbGfIaWSWvZn4n9dJaXDjZzJKUsJILBYNUsUoZrPfWV9qFKQmlnChpb2ViRHXXJBNqGhD8JZfXiXNasW7rFQfrsq12zyOApzCGVyAB9dQg3uoQwMEjOEZXuHNSZwX5935WLSuOfnMCfyB8/kDjSGPDg==
ky
kx
z = -W/2
Fermi arc
lower surface
Figure 1: (a) Schematic picture of the considered Weyl-slab model in a mixed
momentum-/real-space illustration. The plot shows the slab spectrum as a func-
tion of in-plane momenta. Top (blue) and bottom (red) surface states are indicated
with their position and propagation direction in real space. Chiral bulk states connect-
ing surface states of both surfaces are depicted in green. (b) Top-view on equi-energy
contours of the slab spectrum of (a).
the spatial separation.
In this work, we show that in a disordered Weyl semimetal slab valley polarization can
be induced without external magnetic fields as a finite-size effect at mesoscopic slab widths,
possibly extending to even larger sizes. Crucial turn out to be confinement-induced chiral bulk
states [33, 36]: At zero magnetic field and a finite Fermi momentum k F there is a residual
density of chiral bulk states, which must remain to reconnect the two Fermi-arc surface states
as shown in Fig. 1. The density of chiral bulk states of a single valley, relative top
the density
of magnetic-field induced chiral states, is k F /l B−2 W , where l B = ħ
p
h/eB ≈ 26 nm 1/B[T] is
the magnetic length, and W the width of the slab. Taking an experimentally realistic value
−1
of k F = 0.01 Å , the density of anomalous chiral states is larger than that of magnetic-field
induced ones for W ® 100 nm /B[T]. At a mesoscopic width W ∼ 100 nm, the effect of
anomalous chiral states is thus comparable to that of the field-induced chiral states at B ∼ 1 T,
which relevance is commonly accepted and experimentally supported [15–17].
We find that the confinement-induced valley polarization and the presence of surface states
can lead to conductivity enhancements by several orders of magnitude, compared to that of
the infinite system. This conductivity enhancement is suppressed with increasing width W as
1/W , simply due to the decreasing density of confinement-induced chiral states. The valley
polarization, however, turns out to remain unsupressed up to widths set by the probability of
direct inter-valley scattering,
which in case of Gaussian-type disorder potetials is exponentially
enhanced ∼ exp (ξ∆k)2 .
To reveal this effect, we develop a two-part numerical approach combining the full quantum
mechanical calculations of a multiband slab dispersion and wavefunctions, with a numerical
solution for the non-equilibrium corrections to the density matrix [40]. The resulting non-
equilibrium density matrix is exact to leading order in the disorder potential and the applied
external electric field.
The paper is organized as follows. In Section II we start with a derivation of transport
3
SciPost Phys. 11, 046 (2021)
equations for multiband systems with a large number of bands and discuss its validity regime.
In Section III we introduce the model of the Weyl slab, the impurity potential, and calculate
the scattering rate. In Section IV we present the transport results obtained from solving the
transport equations numerically, which we discuss in Section V by comparing with simplified
analytic calculations. We conclude in Section VI.
2 Quantum transport approach
In the first part of this section we recapitulate the general transport formalism in the presence
of weak disorder following Kohn and Luttinger [41]. This is necessary to identify the validity
regime of this formalism when applied to a spatially confined system, which we do in the
second part.
2.1 General quantum transport approach
We separate a general single-particle Hamiltonian into the free-particle part H0 , the additional
weak scattering potential V , and a time-dependent electric-field term eE · r est with the position
operator r and an adiabatic time dependence est with s → 0+ ,
H = H0 + V + est eE · r . (1)
The scattering is due to a random configuration of impurities with a vanishing impurity-averaged
potential ⟪V ⟫ = 0.
We make the ansatz for the full density matrix
ρ = p + g est , (2)
where p is the equilibrium density matrix
e−β(H0 +V )
p= (3)
Tr e−β(H0 +V )
and g the non-equilibrium correction. Inserting into the von Neumann equation for the density
matrix,
i∂ t ρ = [H, ρ] , (4)
and expanding to first order in E we obtain
i s g = [eE · r , p] + [H0 + V, g] . (5)
The following analysis consists in expanding (5) in powers of V . We write (5) in terms of
its matrix elements in the basis of H0 eigenstates |κ〉, where κ combines the quantum numbers.
Off-diagonal and diagonal elements read, respectively,
X
(Eκ − Eκ0 − is) gκκ0 = (gκ − gκ0 )Vκκ0 + Cκκ0 + (gκκ00 Vκ00 κ0 − Vκκ00 gκ00 κ0 ) , (6)
κ00 6=κ0 ,κ
X
−is gκ = Cκ + (gκκ0 Vκ0 κ − Vκκ0 gκ0 κ ) , (7)
κ0 6=κ
where the notation is Aκκ0 = 〈κ|A|κ0 〉, Aκ = Aκκ , Eκ = 〈κ|H0 |κ〉. The field-dependent term
Cκκ0 = eE · [r , p]κκ0 (8)
4
SciPost Phys. 11, 046 (2021)
expands in powers of V starting with the zeroth order,
(0)
Cκκ0 = eE · rκκ0 [n F (Eκ ) − n F (Eκ0 )] , (9)
where n F (E) is the Fermi distribution. From (6) and (7) we see that the leading order of the
off-diagonals of g are of order V −1 , while the diagonals are of order V −2 . To leading order, the
latter two terms in (6) can thus be neglected, leading to
(gκ − gκ0 )Vκκ0
gκκ0 = . (10)
Eκ − Eκ0 − is
Inserting into (7), taking the adiabatic limit s → 0+ , and applying disorder averaging ⟪. . .⟫ we
obtain
δ(Eκ − Eκ0 ) ⟪|Vκκ0 |2 ⟫ (gκ0 − gκ ) .
X
Cκ(0) = i 2π (11)
κ0 0 6=κ
If the electric field points in a direction in which the system is infinite, let it be x and y, the
eigenstates can be chosen as momentum k = (k x , k y ) eigenstates. The field term (9) becomes
Cκ(0) = i eE · vκ n0F (Eκ ) , (12)
where vκ = ∂κ Eκ is the velocity.
Making the ansatz
gκ = −eE · Λκ n0F (Eκ ) (13)
(11) simplifies to
n0F (Eκ )δ(Eκ − Eκ0 ) ⟪|Vκκ0 |2 ⟫ (Λκ − Λκ0 ) ,
X
n0F (Eκ )vκ = 2π (14)
κ0 6=κ
known as Boltzmann equation (BE), to be solved with respect to the vector-valued state-resolved
transport mean free paths Λκ , which we will refer to as transport length in short. The average
magnitude of the transport length scales with the strength of the impurity potential as V −2 .
Note that since the summation operator acting on Λκ in (14) has an eigenvalue zero for a
κ independent vector, the solution is generally determined up to a constant
P 0
κ n (Eκ )Λκ
c = P F0 . (15)
κ n F (Eκ )
Particle conservation however requires κ gκ = 0, which fixes the constant to c = 0.
P
The current-density expectation value reads
1X
j = −e v g , (16)
V κ κ κ
where V is the system volume. The conductivity tensor σ, defined as j = σE, becomes,
using (13),
1X 0
σ = e2 n (E )v ⊗ Λκ . (17)
V κ F κ κ
Note that the BE (14) is exact in the weak-disorder limit, giving a conductivity that scales with
the squared inverse strength of the disorder potential. Leading corrections, which will not be
considered here, are of zeroth order in the impurity potential, they include, e.g., the anomalous
Hall effect.
5
SciPost Phys. 11, 046 (2021)
2.2 Application to a slab model
We now discuss the validity regime of (14) when applied to a slab model. We consider a system
that is infinite in two spatial directions, x and y (as specified above), and confined in direction
z to −W /2 < z < W /2. The slab energy eigenspace κ = (k, b), where k is the in-plane
momentum, has the particularity that the number of bands (band index b) is potentially very
large, scaling with the width W of the system. Since the BE that we have just derived relies
only on considering the leading order in the scattering potential V , it can still be applied to
the slab, provided that V can be taken to be arbitrary small. For the slab, a problem arises if
we want to consider such a large width W that the effect of boundaries becomes smaller than
that of the impurity scattering, which can invalidate the expansion in powers of V . We now
examine when exactly the width becomes “too large” in that sense.
The large width W enters our above formalism through the position matrix elements in
Cκκ0 , see Eq. (8). Let us thus repeat the above steps without neglecting higher V orders in
C since they might still be large due to W . In this case the BE obtains an extra term on the
right-hand side (rhs), so that the extended BE reads
X X
Cκ = δ(Eκ − Eκ0 )|Vκκ0 |2 (gκ0 − gκ ) + δ(Eκ − Eκ0 )(Cκκ0 Vκ0 κ − Vκκ0 Cκ0 κ ) . (18)
κ0 6=κ κ0 6=κ
Now expanding C in powers of V the rhs term with C (0) vanishes upon impurity averaging
since the mean potential due to impurities is zero. We thus consider the next order term,
X n (E 00 ) − n F (Eκ0 ) n (E ) − n F (Eκ00 )
(1)
Cκκ0 = eE · rκκ00 Vκ00 κ0 F κ −rκ00 κ0 Vκκ00 F κ .
κ00
Eκ00 − Eκ0 Eκ − Eκ00
While it also vanishes upon averaging on the left-hand side (lhs) of (18), the new term on the
rhs becomes
XX n F (Eκ ) − n F (Eκ00 )
eE · δ(Eκ − Eκ0 ) [Vκ00 κ0 Vκ0 κ rκκ00 + Vκκ0 Vκκ00 rκ00 κ0 ] .
κ0 6=κ κ00
Eκ − Eκ00
In this sum there are terms that are proportional to |Vκκ0 |2 , which certainly do not vanish upon
impurity averaging. Compared to the first term on the rhs of (18) with the general ansatz (13),
the new term is generally smaller if the position matrix elements are smaller than the typical
values of the transport length, let us denote them by Λ̄. This correction can thus be neglected
only if the width is much smaller,
W Λ̄ . (19)
Higher order terms due to the expansion of C are of order W times a higher power of V and
thus give even smaller corrections.
Summarizing section III, Eq. (19) characterizes the validity regime of the BE (14) if applied
to the slab. In words, one is allowed to consider impurity scattering as a weak perturbation to
a free propagation in the slab treated as a two-dimensional multiband system as long as the
mean free path is much larger than the width.
3 Weyl-semimetal slab model
We consider a minimal lattice model of a Weyl semimetal [42],
H0 (κ) = tσ x sin kz + tσ y sin k y + mκ σz , (20)
6
SciPost Phys. 11, 046 (2021)
0.040
0.035
0.030 ky
0.025 13 contours
| |2 kx
0.020
0.015
0.010
0.005
0.000
100 50 0 50 100
z [a]
Figure 2: Position-z-resolved probability density |ψ|2 = |〈z|ψ〉|2 of Fermi-level states.
Colors indicate surface states at the lower boundary (red) and upper boundary (blue),
the confinement-induced chiral bulk states (green), and normal bulk states (gray).
The inset shows the position of the states in momentum space. Other parameters are
E F = 0.3 t, W = 201, µ b = 0.3 t.
where mκ = t(2 + cos β − cos k x − cos k y − cos kz ) and σi are pseudospin Pauli matrices
(corresponding to an arbitrary degree of freedom), t is the hopping amplitude, and the lattice
constant is set to unity. The two Weyl nodes are placed at k = ±β, where β = β x̂ corresponds
to a time-reversal breaking magnetization. We consider a “good” Weyl semimetal with a cone
separation β ∼ 1.
The Hamiltonian of the slab is given by the lattice Hamiltonian (20) but for a finite number
W of sites in the z direction. Transformation into the site basis in the z direction replaces
cos kz → (δi, j+1 + δi, j−1 )/2, where i = 0, 1, . . . , (W − 1) is the site number, corresponding to
the discrete position in z,
W −1
z≡i− , (21)
2
in units of the lattice constant which is set to one. We furthermore add a boundary potential µ b
at the surface layers of the slab, which main effect is to bend the Fermi-arc surface states. We
label the eigenstates by κ = (k, b) where k = (k x , k y ) are the continuous in-plane momenta
and b denotes the 2W modes at each value of k.
The eigenstates |ψκ 〉 and eigenenergies Eκ of the slab are obtained from exact diagonaliza-
tion of the Hamiltonian at a fixed in-plane momentum k using standard methods of numerical
diagonalization [40]. For our transport considerations we need to take into account all Fermi-
level states, which are continuous contours in the space of the in-plane momentum k. We
numerically [40] determine the contours by means of the marching squares algorithm [43],
whereby the contours are discretized. The precision level of the discretization is improved until
full convergence of the results. Figure 2 illustrates typical results of numerical diagonalization
giving the Fermi-level contours (inset, see also Fig. 1) and the wavefunction probability density
|ψ|2 ≡ |〈z|ψ〉|2 .
Most bulk states form closed contours located at one of the valleys. Additionally there is
the special contour that connects the valleys by wrapping around them. In between it consists
7
SciPost Phys. 11, 046 (2021)
Figure 3: Plot of the function M (z, z 0 , ξ, W ) (a) and dz 0 M (z, z 0 , ξ, W ) (b) at ξ = 3
R
and W = 100. The width of the diagonal peak in (a) and the region of reduced weight
in (b) are both set by ξ.
of surface states but at places where the contour touches the bulk contours, the states are
delocalized — we call them chiral bulk states, since they are unidirectional, moving parallel or
antiparallel to the intrinsic magnetization (here x direction) depending on the valley. While
the number of bulk contours ∼ k F W /π increases with the width, there is always only a single
contour that contains the surface and chiral bulk states.
4 Scattering
4.1 Disorder potential
We model the disorder by static Gaussian potentials,
X 2
/2ξ2
V= Uα φ(r − rα ) , φ(r ) = e−|r | , (22)
α
where the sum runs over the Gaussian’s with a characteristic width ξ, random and uncorrelated
potential magnitudes Uα ∈ [−δ, δ], and random positions rα .
The disorder potential enters the BE (14) in the form of the scattering rate between two
energy eigenstates Q(κ, κ0 ) = 2πδ(Eκ − E F )q(κ, κ0 ), where
q(κ, κ0 ) = ⟪|〈ψκ |V |ψκ0 〉|2 ⟫ (23)
and ⟪. . .⟫ denotes disorder average.
In the slab model, as compared to a translation-invariant system, the scattering rate is not
a simple Gaussian as a function of the momentum difference. Inserting the impurity potential
and averaging over the disorder configurations within the slab we obtain
(2πξ2 )3 δ2 ni −ξ2 (k−k 0 )2 /2 X
q(κ, κ0 ) = e M (z, z 0 , ξ, W )ψ†κ0 (z)ψκ (z)ψ†κ (z 0 )ψκ0 (z 0 ) , (24)
3L x L y z,z 0
where ψκ (z) = 〈z|ψκ 〉 and ni is the impurity concentration. A detailed derivation can be found
8
SciPost Phys. 11, 046 (2021)
in the Appendix. The function M (z, z 0 , ξ, W ) is given by
0 2 2
e−(z−z ) /4ξ W + z + z0 W − z − z0
0
M (z, z , ξ, W ) = p erf + erf ,
4 πξ 2ξ 2ξ
p Rx 2
where the error function is defined as erf (x) = (2/ π) 0 e−t dt. As illustrated in Fig. 3,
M (z, z 0 , ξ, W ) is mainly the (z − z 0 ) dependent part of the Gaussian impurity potential, which
magnitude however reduces by approximately a factor of 2 at the slab surfaces, where the
possible impurity positions obviously fill only half of the space. Another effect of the finite size
is that the z dependence of the wavefunctions is not plane-wave like, hence the sum in (24)
does not reduce to a Fourier transformation. In particular, note that the wavefunction factor
in (24) strongly suppresses the scattering rate between surface states of opposite surfaces when
ξ and the penetration depth of the surface states are both much smaller than W due to a
vanishing overlap between the surface states.
In the limit ξ 1 one obtains M (z, z 0 , ξ, W ) = δ(z − z 0 ), in which case
ξ 1 (2πξ2 )3 δ2 ni −ξ2 (k−k 0 )2 X
q(κ, κ0 ) −−→ e |ψκ0 (z)|2 |ψκ (z)|2 .
3L x L y z
Despite the rather complex form, the impurity scattering is fully determined by two param-
eters — the real-space impurity width ξ and the overall impurity strength (set by δ2 ni ), the
latter will be in the following quantified by the average mean free path l, defined below. The
impurity width ξ essentially sets the momentum-space range of most scattering processes to
|k − k 0 | ® ξ−1 , hence for a large value of ξ scattering between states that are far apart from
each other in the in-plane momentum k is exponentially suppressed ∼ exp[−ξ2 (k − k 0 )2 ].
4.2 Scattering lengths
In the presence of surface and bulk states it is interesting to quantify averaged scattering
rates and scattering probabilities between different types of states, which will be helpful to
understand the numerical transport results.
To quantify the overall strength of impurity scattering, we define the averaged mean free
path * −1 +
X
l = |vκ | δ(Eκ − E F )q(κ, κ0 ) , (25)
κ0
where the Fermi-surface average is given by
1X X
〈. . . 〉 = δ(Eκ − E F )(. . . ) ; N= δ(Eκ − E F ) , (26)
N κ κ
and N is the number of states at the Fermi level. Note that for a nearly constant |vκ | that we
have, the mean free path is inversely proportional to the total scattering probability.
To quantify the scattering probability between different types of states i ∈ [b, s], where b
(s) denotes bulk (surface) states, we define the scattering length
* !−1 +
j
X
0
li j = |vκ | δ(Eκ − E F )q(κ, κ ) , (27)
κ0 i
where the sum runs overP thePtype j of states and the averaging is analogous to (26) but only
i
over type i of states ( κ → κ ).
9
SciPost Phys. 11, 046 (2021)
We now want to determine the dependence of the scattering probability on the width W
of the slab in two limiting cases (i) the number of bulk states being much larger than that of
surface states, N Ns and (ii) the number of states being dominated by surface states, N ∼ Ns .
In regime (ii) the number of bulk states Nb = N − Ns , is given by the number of chiral bulk
states. Thus the scaling of Nb with the slab width reads
¨
W (i)
Nb ∼ (28)
1 (ii).
The scaling of the scattering rate (24) is governed by the z dependence of the wavefunctions.
A normalized surface state with a penetration depth λ ∼ β −1 ∼ 1 and a normalized bulk state
are of the form
v
t2 p
ψs (z) = e−z/λ , ψ b (z) = 1/ W , (29)
λ
respectively. Consequently, for W λ (which we always consider), we can estimate the scaling
of the scattering rate between the different types of states as
1 same surface
q(κ, κ ) ∼ W1 bulk-surface, bulk-bulk
0
(30)
0 opposite surfaces.
From this, (28), and (27) the width dependence of the scattering probabilities summarizes to
l b b /l lss /l lsb /l l bs /l
Nb Nb
(i) 1 Ns W ∼1 1 Ns ∼ W (31)
Ns Ns
(ii) W Nb 1 W Nb W
Most importantly, in the regime (i) the scattering probability from bulk to surface (∝ l −1 bs
) is
a factor W smaller than other scattering probabilities, due to the ratio of the number of bulk
states to surface states, which is large and increases with W . In the regime (ii) instead, the
number of bulk states (consisting only of the chiral bulk states) does not depend on W . The
peculiarity of this regime is that surface states scatter most probably within the surface states,
which is due to the larger overlap of surface wavefunctions.
5 Numerical results
We calculate the nonequilibrium occupation function (13) at zero temperature by numerically
solving Eq. (14) with respect to the transport length on the basis of the numerical solution of
the discretized slab spectrum discussed in Section III [40]. The nonequilibrium occupation
function (or, equivalently, the transport length) determines the conductivity, given in Eq. (17),
and the valley polarization, defined in (5.2) below.
5.1 Conductivity
We consider the conductivity in units of the standard Drude estimate given by the mean free
path l, the density of states at the Fermi level n = N /V , and the Fermi velocity v ≈ t/ħ
h,
e2 n l v
σ0 = . (32)
3
10
SciPost Phys. 11, 046 (2021)
(a) W = 201 (b) xx / 0 yy / 0
1
10 30
20
ky [ /a]
0 13 contours 2.6
5 2.3
10
1
2 1 0 1 2 0 2 4 0 2 4
kx [ /a] [a] [a]
EF = 0.3
(c) W = 21 (d) xx / 0 yy / 0
1
20 100
ky [ /a]
0
2.2 50 2.1
1 0
2 1 0 1 2 0 2 4 0 2 4
kx [ /a] [a] [a]
(e) W = 201 (f) xx / 0 yy / 0
1
7.5 15000
0.1
EF = 0.03
ky [ /a]
0 5.0 10000
2.7
2.5 5000
1
2 1 0 1 2 0 2 4 0 2 4
kx [ /a] [a] [a]
Figure 4: Impurity-range dependence of the conductivity. (a),(c),(e): Discretized
Fermi-level contours for the three parameter choices (E F = 0.3 t, W = 201),
(E F = 0.3 t, W = 21), and (E F = 0.03 t, W = 201). The number of finite-size
induced states (surface states and chiral bulk states), Ns + Nc , relative to the total
number of states N is (a) (Ns + Nc )/N = 0.14, (c) (Ns + Nc )/N = 0.68, and (e)
(Ns + Nc )/N = 1. Color indicates surface states at the upper (blue) and lower (red)
surfaces, bulk states (gray), and the chiral bulk states (green), cf. Fig. 2. (b), (d), (f):
Conductivity as a function of the impurity width ξ. The green curves show power-law
fits, the exponent is indicated in the plot. Other parameters are β = 1.5, µ b = 0.3 t.
This is the result one would expect to find for a system with point-like impurities (ξ → 0) and
only bulk states arranged in form of a spherical Fermi surface.
The dependence of the slab conductivity on the width of the impurity potential ξ is sum-
marized in Fig. 4. For ξ > 1 the conductivity is well fitted by a power-law dependence on ξ,
with an exponent between 2 and 3. An exception is found for σ y y at E F t, which shows a
weak ξ dependence at a strongly enhanced conductivity at all ξ. In total, the magnitude of the
conductivity, especially in the direction of motion of Fermi arcs, may be enhanced by several
orders of magnitude, either due to a wide impurity range, or if the Fermi energy is close to the
Weyl points.
To gain further insight, in Fig. 5 we consider the width dependence of the conductivity.
Figure 5(a) shows that in the case of a large number of bulk states and point-like impurities,
the conductivity is nearly independent of W and is close to σ0 — in this regime the slab
thus resembles a conventional metal. At a large ξ, however, the conductivity enhancement
decreases antiproportional to the width, indicating that the conductivity enhancement at large
E F is related to the presence of the confinement-induced surface and chiral bulk states, which
number, relative to the total number of states N , is antiproportional to W . At small E F , however,
when there is only one Fermi-level contour mainly consisting of surface states, the total number
of states is nearly independent of W . In this case, the large conductivity in the direction of
motion of surface states, σ y y , linearly increases with the width W .
11
SciPost Phys. 11, 046 (2021)
= 0.1 =5
(a) xx / 0 yy / 0 (b) xx / 0 yy / 0
2.8 14
W 1 W 1
EF = 0.3 t
1.26
2.6 40
13 contours 12
ky 2.4
1.25 10 30
kx
100 W [a] 200 100 W [a] 200 100 W [a] 200 100 W [a] 200
(c) xx / 0 yy / 0 (d) xx / 0 yy / 0
10000 9.5 4000
W
EF = 0.03 t
1.50 7500 9.0
ky 2000
1.25 W 5000 W W
kx 8.5
100 W [a] 200 100 W [a] 200 100 W [a] 200 100 W [a] 200
Figure 5: Width dependence of the conductivity at (a) E F = 0.3 t, ξ = 0.1, (b)
E F = 0.3 t, ξ = 5, (c) E F = 0.03 t, ξ = 0.1, and E F = 0.03 t, ξ = 5. Other parameters
are β = 1.5, µ b = 0.3 t. The dashed lines in (b) indicate W −1 dependence and in (c)
and (d) linear dependence.
xx / 0 yy / 0
21
225
20
200
19
175
18
150
17
125
16
15 W [a] 35 15 W [a] 35
Figure 6: Width dependence of the conductivity at ξ = 5 and E F = 0.15 t, µ b = 0.3 t.
The transition from a single contour to three contours is indicated by the dotted lines
and the contour-plots below the W axis.
To explore more carefully the transition from an increasing to a decreasing W dependence,
in Fig. 6 we plot the conductivity at smaller width, when the first bulk contours appear. This
plot shows that the strong W enhancement of the surface conductivity requires the normal bulk
contours to vanish, which happens if W k F /π ® 1. As soon as at least one normal bulk contour
appears, the conductivity jumps to a significantly lower value and becomes decreasing in W .
12
SciPost Phys. 11, 046 (2021)
150 150
100 100
50 2.5
50
0 0
0 5 100 200
[a] W [a]
Figure 7: Valley polarization χ, defined in Eq. (5.2). (a) Dependence of χ on the
impurity range ξ at a width W = 201. (b) Width dependence of χ at ξ = 5. Other
parameters are µ b = 0.3 t, β = 1.5, and E F = 0.3 t.
5.2 Valley polarization
Besides the conductivity it is interesting to explore the non-equilibrium occupation difference
P± the valley separation, E = x̂ E.
of the valleys, which occurs when the electric field points along
The average occupation of the valley at k x = ±β is given by κ gκ /(Nb /2), where the sum runs
P±
over all bulk states at the valley ±, and Nb /2 = κ is the number of those states. We quantify
the valley polarization by the difference of the valley occupations relative to the standard
occupation difference of states δ(E F − Eκ )eEl due to the mean free motion in the electric field,
P+ P−
κ gκ κ gκ 〈Λκx 〉+ − 〈Λκx 〉−
χ ≡ P+ − P− = , (33)
κ δ(E F − Eκ )eEl κ δ(E F − Eκ )eEl l
where inPthe second line wePused Eq. (13) and defined the average over valley bulk states
± ±
〈. . . 〉± = κ δ(E F − Eκ ) . . . / κ δ(E F − Eκ ).
A representative result for the impurity range and width dependence of the valley polariza-
tion is shown in Fig. 7. The valley polarization shows the power-law ξ dependence, similar to
the conductivity. However, unlike for the conductivity, there is no significant width dependence
at large ξ (here already at ξ ¦ 1) in the numerically accessible width range. This is very
surprising as it seems to imply a presence of valley polarization for sufficiently large impurity
widths in arbitrary large systems, contradicting previous predictions based on infinite-system
calculation [3]. Below we will show that the valley polarization in fact does decay but for
2
width above ¦ l e(ξ∆k) , which becomes exponentially large for ξ∆k ¦ 1.
6 Discussion
The numerical results of the previous section show enhancements of the slab conductivity by
several orders of magnitude (compared to the expectation for a conventional metal (32)) and
a substantial valley polarization in a wide region of the parameter space. The characteristic
dependencies on ξ and W allow to identify the main mechanisms of these effects, which we
now systematically discuss.
13
SciPost Phys. 11, 046 (2021)
6.1 Impurity-range dependence
When the impurity range ξ increases, the scattering rate (24) between two countermovers
separated by ∆k becomes exponentially suppressed, ∼ exp[−(ξ∆k)2 ]. The transport length
does not inherit the exponential enhancement though, since relaxation happens via multiple
small-angle scattering processes. To illustrate this, we consider a toy model of a closed chain of
N states labeled i ∈ [1, N ] with arbitrary velocities v i . The BE (14) is of the form (summation
over repeated indices assumed)
v i = Mi j Λ j , (34)
where Mi j = δi j k qik − qi j is given by the scattering rates between states, qi j . We assume
P
scattering only between the nearest neighbors, with the rate qnn , and direct scattering between
countermovers with the rate qd , in which case the matrix M becomes
Mi j = (2qnn + 2qd )δi j − qnn (δi j+1 + δi, j−1 ) , (35)
where we used that the direct-scattering part of qi j cancels when multiplied with Λ. For the
nearest-neighbor part Mnn of M there is a left pseudo-inverse,
1 |i − j| (|i − j| − N )
Pi j = , (36)
qnn 2N
so that P Mnn Λ = Λ (note Λi = 0 due to particle conservation). With its help, the full
P
i
solution of (34) becomes
Λ = (1 + 2qd P)−1 P v . (37)
We first ignore direct scattering. From the form of P, it is clear that the solution can de-
pend on N to maximally the power N 2 . In particular, for a circular velocity arrangement
v = {− cos[2π(i − 1)/N ], sin[2π(i − 1)/N ]} the vector mean free path for N 1 assumes
the value
1 (ξ̃∆k)2
Λ= v, (38)
qnn 4
where we have written the number of nearest neighbors as N = π ξ̃∆k, in terms of the spacing
between nearest neighbors ξ̃−1 and the distance between countermovers ∆k.
Considering the full result with direct scattering in Eq. (37), we see that the nearest-neighbor
scattering contributions to Λ dominate as long as
(ξ̃∆k) p 1
, (39)
qnn qd
where the exponent p depends on the velocity arrangement. In the opposite limit we instead
obtain
1
Λ= v. (40)
2qd
Transferring this insight to the Weyl slab model, the distance between the nearest neighbors
ξ̃−1 corresponds to ξ−1 . The mean free path corresponds to the inverse diagonal of M times
velocity, l ∼ v/(qd + qnn ). Interpolating between the two regimes of dominant direct scattering
and nearest-neighbor scattering the typical values of the transport length may be well estimated
as
Λ̄ = l [1 + (ξ∆k) p ] . (41)
This explains the general power-law conductivity enhancement with ξ in Fig. 4, except for σ y y
in Fig. 4(f), which we discuss separately.
14
SciPost Phys. 11, 046 (2021)
6.2 Valley polarization
The averaged occupation of the two valleys ± can be expressed as
P±
gκ
µ± = κ , (42)
Nb /2
P±
where κ δ(E F − Eκ ) = Nb /2. The occupation difference is related to the valley polarization χ
defined in Eq. (5.2),
µ+ − µ− = χ eE x l . (43)
A single valley has a total velocity in the x direction — the unbalanced velocity of the chiral
states vNc /2, where Nc is the total number of chiral states (Nc /2 in each valley). An electric
field in the x direction thus pumps charge between the valleys with the rate veE x Nc , which in
a steady state must be counterbalanced by scattering. We can write down a simple balance
equation as a condition for a steady valley occupation,
dµ± ∂ µ±
eNb = 0 = ±veE x Nc + eNb , (44)
dt ∂ t scat
where we used Nb Nc . The time-change of µ+ − µ− due to scattering is proportional to the
occupation difference itself and the scattering probability, which we quantify by the scattering
length l c , hence
∂ (µ+ − µ− ) µ+ − µ− χ eEl
= −2 = −2 . (45)
∂t scat l c /v l c /v
Together with (44) we obtain
Nc l c
χ= . (46)
Nb l
For a point-like disorder potential (ξ → 0) the ratio l c /l goes to one and the valley polarization
is small. For larger ξ, direct scattering between the valleys becomes strongly suppressed and
the relaxation of µ+ − µ− must go via surface states. Thereby the relaxation along the arcs
and the scattering from surface to bulk is much faster than the scattering from bulk to surface,
see (31). The scattering length l c is thus set by the bulk-surface scattering length l bs , which is
proportional to the ratio Nb /Ns so that
Nc
χ∝ . (47)
Ns
This explains the surprising result that the valley polarization does not depend on the width,
as seen numerically in Fig. 6. This is surprising since the origin of the valley polarization are
the chiral bulk states, which number Nc is a factor ∝ W smaller than the total number of bulk
states Nb . The explanation is that the valley relaxation also becomes suppressed ∝ W −1 since
the probability to scatter into Fermi arcs decreases with an increasing number of bulk states.
For a much larger width, when the probability of relaxation via Fermi arcs becomes smaller
than the probability of relaxation via direct inter-valley scattering, the valley polarization
will ultimately go to zero like W −1 . In case of a Gaussian potential, the amplitude of direct
inter-valley cattering is however exponentially suppressed ∼ exp[−(ξ∆k)2 ] so that the width
independence can easily extend to arbitrary macroscopic sizes for realistic values of the impurity
width and the cone separation ξ∆k ¦ 1.
Regarding the strong enhancement of valley polarization with ξ, one is tempted to under-
stand it as a consequence of an increasing number of scattering events needed for a relaxation
along the Fermi arcs. However, the contribution of such a process to χ would enter in the form
lss (ξ∆) p /l ∝ 1/W , with a clear 1/W dependence, which we do not observe. It rather must be
15
SciPost Phys. 11, 046 (2021)
the suppression of the bulk-surface scattering probability with ξ, which is plausible in view of
the typical effect of an increasing ξ to increase the relaxation time. We note, however, that
the simplified analytical calculation of Section 6.1 does not apply since in this case the valleys
correspond to two-dimensional pools of states, while the model of Section 6.1 only considers
one-dimensional chains.
6.3 Conductivity in case of a large number of bulk states
For bulk states (excluding chiral bulk states) the small separation of countermovers in momen-
tum space ∼ k F makes their transport length close to the mean free path l and not significantly
enhanced with ξ. We thus approximate the current contribution of bulk states as
jn = σ0 E. (48)
For an electric field in the y direction there is additionally the contribution of surface states,
which transport length is mainly set by surface-surface scattering, enhanced by ξ according
to (41),
js = e2 ns v lss (ξ∆k) p y E y , (49)
where ns = Ns /V is the density of surface states.
The current contribution of chiral bulk states is negligible compared to (48), except if the
valley polarization becomes large, in which case
jc = ev x̂ nc (µ+ − µ− )/2 = e2 v x̂ nc χ l E/2 , (50)
where nc = Nc /V is the density of chiral bulk states. Adding all the current contributions, we
obtain
σ(i)
xx 3 Nc
=1+ χ, (51)
σ0 2 Nb
σ(i)
yy Ns lss (ξ∆k) p 1
=1+3 ∼ 1 + (ξ∆k) p , (52)
σ0 Nb l W
where we again used (31). This rough estimate is in qualitative agreement with the numerical
results. Note that Nb increases proportional to W , which explains the W dependence of σ x x
and σ y y in Fig. 5(b). We now understand that the enhancement of σ x x and σ y y are mainly
due to chiral bulk states and surface states, respectively.
In conventional metals the conductivity is width independent as long as the width is much
larger than the mean free path; when the width becomes smaller, the conductivity tends to
decrease due to additional scattering at boundaries. Our work shows that in Weyl semimetals
away from charge neutrality the opposite trend of increasing conductivity with shrinking the
width may occur due to chiral bulk states and surface states. A qualitatively similar width
dependence may occur also in the regime l W (we consider the opposite limit l W ),
which has been considered in Ref. [34]. Observations of enhanced conductivity for reduced
widths have been reported in Ref. [35], where the Weyl semimetal nanobelts, according to
estimates of the mean free path, are presumably in the regime of our work, l W , or in the
crossover regime l ∼ W .
6.4 Conductivity in case of a small number of bulk states
We now come to the case (ii) when the number of states is dominated by surface states, while in
the bulk only chiral bulk states are present. The conductivity in the y direction can be written
in the form
σ(ii)
yy Ns Λs
=1+ , (53)
σ0 N l
16
SciPost Phys. 11, 046 (2021)
where Ns ≈ N is the number and Λs the transport length of surface states. Scattering within
surface states at the same surface does not lead to relaxation of motion in the y direction
since the average velocity v y of those states is not zero. Countermoving surface states, on
the other hand, have no overlap with each other, direct scattering between them is blocked.
The relaxation of Fermi arc states must thus go via the small number of chiral bulk states,
so that Λs is set by the surface-bulk scattering probability, Λs ∼ lsb . In Section 4.2 we found
ls b /l ∼ W Ns /Nc —Eq. (31)—which leads to
σ(ii)
yy Ns
∼W . (54)
σ0 Nc
Both factors are of order 100 for parameters in Fig. 4(f), which explains the large magnitude.
Also the W dependence in Fig. 5 is consistent since both Ns and Nc are W independent in this
case.
For large ξ, surface-bulk scattering becomes limited to small regions at the nodes and the
full relaxation must involve nearest-neighbor scattering along the Fermi arc. According to
Section 6.1, the latter should elongate the full transport length by an additional lss (ξ∆k) p . For
the considered parameters, ls b is much larger than this additional part since (ξ∆k) p W Ns /Nc ,
which explains the weak ξ dependence in Fig. 4(f).
The transition from (ii) to (i) spoils the strong enhancement of σ y y in two ways: First,
the ratio Ns /Nc changes to Ns /(Nc + Nb ) and thus becomes smaller and second, according to
Eq. (31), ls b /l ∼ W Ns /Nb → 1 is no longer width dependent, which in Fig. 6 explains the jump
and the change of slope.
The conductivity in the x direction is also governed by the dominant number of surface
states. Relaxation however happens via scattering within the same surface, since vx averages
to zero at each surface separately. Since lss /l ∼ 1, the conductivity in the x direction is not
significantly enhanced.
7 Conclusion
In conclusion, we have studied linear-response properties of a finite Weyl semimetal slab (width
W ) in the presence of long-ranged disorder (disorder potential width ξ). Our work highlights
the remarkable property of Weyl semimetals to realize valleys of opposite chirality that are
well separated in momentum space (∆k) and continuously connected only via surface states.
For a Fermi energy that is not exactly at the Weyl nodes, the surface states occur together
with confinement-induced chiral bulk states. In the presence of an electric field parallel to
cone separation they allow to violate chiral charge conservation even without an external
magnetic field. This peculiarity stabilizes an anomalous valley polarization at zero magnetic
field. If the potential width is substantially larger than the inverse separation of valleys, the
valley polarization persists up to very large slab width. This is explained by the fact that direct
inter-valley scattering is strongly suppressed and the relaxation must go via Fermi arcs, which
is however also increasingly ineffective owing to their vanishing density with an increasing
width. The resulting width independence of the confintenment-induced valley polarization
persists up to a width, for which relaxation via direct inter-valley scattering becomes more
effective than relaxation via Fermi arcs. For Gaussian-type disorder potentials this maximum
width is exponentially enhanced by exp[(ξ∆k)2 ] and can thus easily reach macroscopic length
scales at realistic values of cone separation ∆k and inpurity-potential widths ξ.
The valley polarization and Fermi-arc surface states lead to a conductivity enhancement
which increases with an increasing width of the disorder potential and a decreasing width of
the slab (σ ∝ 1/W ). Moreover, if the Fermi energy is reduced towards charge neutrality such
17
SciPost Phys. 11, 046 (2021)
that normal bulk states vanish completely, the conductivity in the direction of motion of surface
states becomes strongly enhanced because relaxation of surface states can only go via bulk
states which number becomes strongly reduced.
Methodologically our work performs first steps in the application of the weak-disorder
transport formalism to a multilayer system with a large (¦ 100) number of layers and con-
sequently a similarly large number of bands in the in-plane Brillouin zone. The numerical
code [40] is designed to be easily applicable to an arbitrary lattice model and can thus be used
to explore in detail the confinement-induced valley polarization in various Weyl-metal models.
In this work, we introduced the formalism by considering a minimal two-Weyl-cone model. We
find that the qualitative aspects of the valley polarization are robust to lattice details such as
boundary potentials (which give the Fermi arcs a finite curvature) or velocity anisotropy of
the Weyl cones. Our analytical discussion shows that the valley polarization depends on the
mere presence of chiral bulk states and surface states (which is topological) and the ratio of
the inverse separation of Weyl cones vs. the width of the scattering potential, which explains
the robustness of this effect.
An interesting application of the introduced tools is to consider lattice models of existing
Weyl semimetals. For the case of several pairs of Weyl nodes that are sufficiently separated in
momentum space, such as in the TaAs material family, we expect valley polarization to occur in
each pair which cone sepration aligns with the electric field, similarly to the two-cone case. The
reason is that due to the large pair separation, scattering between pairs should be negligible
compared to the intervalley scattering within a single pair, making each pair independent and
thus reduce the problem to the two-cone case.
General limitations of the introduced numerical tools are the restriction to slab width
being smaller than the mean free path and the restriction to the leading order in the disorder
potential. Both the fate of confinement-induced effects for larger widths as well as corrections
of higher order in the disorder potential, which are known to start with Berry phase effects [44],
constitute interesting directions to extend this formalism.
Data availability
All the code and data used to produce the reported results is available in Ref. [40].
Author contributions
M.B. formulated the project idea, developed the theory with input from N.B. and A.A., performed
and analyzed numerical calculations, and developed the analytical model. P.M.P.P. developed the
numerical code with input from N.B., A.A., and M.B., and performed and analyzed numerical
calculations. The manuscript was written by M.B. with input from P.M.P.P. and A.A.
Acknowledgments
This research was supported by the European Union Horizon 2020 research and innova-
tion programme under Grant Agreement No. 824140, Grant No. 18688556 of the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation), ERC Starging Grant 638760,
and the Netherlands Organisation for Scientific Research (NWO/OCW), as part of the Frontiers
of Nanoscience program.
18
SciPost Phys. 11, 046 (2021)
Supplementary Material
Derivation of scattering amplitudes
We consider Gaussian-type static impurity potentials,
X 2
/2ξ2
V= Uα φ(r − rα ), φ(r ) = e−|r | , (S1)
α
where the sum runs over impurities with a characteristic width ξ, random and uncorrelated
potential magnitudes Uα ∈ [−δ, δ], and random positions rα .
For our transport consideration we consider the scattering rate Q(κ, κ0 ) = 2πδ(Eκ − E F ) ×
q(κ, κ0 ) between energy eigenstates |ψκ 〉 and |ψκ0 〉, which we calculate using Fermi’s Golden
Rule,
q(κ, κ0 ) = ⟪|〈ψκ |V |ψκ0 〉|2 ⟫ , (S2)
where the disorder average is defined as
YZ δ Z
d Uα d rα
⟪. . .⟫ = (. . . ) . (S3)
α −δ 2δ V
We write the normalized wavefunctions as
1
〈r |ψκ 〉 = Æ e ik·ρ ψκ (z) , (S4)
Lx L y
where ρ = (x, y) is the in-plane position and L x L y the in-plane volume and ψκ (z) is the
normalized eigenvector of numerical diagonalization of the lattice model, z denoting the
discrete sites in the z direction.
The expectation value of the impurity then calculates to
2πξ2 −ξ2 (k−k 0 )2 /2 X X
〈ψκ0 |V |ψκ 〉 = e Uα φ(z − zα ) ψ†κ0 (z)ψκ (z) . (S5)
Lx L y α z
Inserting into (S2) and using (S3), we obtain
4π2 ξ4
Z
q(κ, κ ) = 2 2 e−ξ (k−k ) ⟪ Uα Uβ φ(z − zα )φ(z 0 − zβ ) ψ†κ0 (z)ψκ (z)ψ†κ (z 0 )ψκ0 (z 0 )⟫
2 0 2 X dz
α
X
0
Lx L y α,β z,z 0
W
4π ξ δ −ξ2 (k−k 0 )2
2 4 2
Z
X X
† † 0 0 dzα
= 2 2
e ψκ0 (z)ψκ (z)ψκ (z )ψκ0 (z ) φ(z − zα )φ(z 0 − zα )
3L x L y α z,z 0
W
4π ξ δ ni −ξ2 (k−k 0 )2 X †
2 4 2
Z
= e ψκ0 (z)ψκ (z)ψκ (z )ψκ0 (z ) dzi φ(z − zi )φ(z 0 − zi )
† 0 0
3L x L y z,z 0
(2πξ2 )3 δ2 ni −ξ2 (k−k 0 )2 X †
= e ψκ0 (z)ψκ (z)ψ†κ (z 0 )ψκ0 (z 0 )M (z, z 0 , ξ, W ) ,
3L x L y z,z 0
(S6)
where ni = α /L x L y W is the impurity concentration, furthermore we used
P
δ2
⟪Uα Uβ . . . ⟫ = δαβ ⟪. . .⟫ , (S7)
3
and defined the function
Z
0 1
M (z, z , ξ, W ) = dzi φ(z − zi )φ(z 0 − zi ) (S8)
2πξ2
0 2 2
e−(z−z ) /4ξ W + z + z0 W − z − z0
= p erf + erf , (S9)
4 πξ 2ξ 2ξ
19
SciPost Phys. 11, 046 (2021)
p Rx 2
plotted in Fig. 3; the error function is defined as erf (x) = (2/ π) 0 e−t dt. Note that
Z
dz 0 M (z, z 0 , ξ, W ) (S10)
is a weak function of z, equal to 1 in the middle of the slab, and going down to ≈ 0.4 at the
edges in the range ξ.
References
[1] S. L. Adler, Axial-vector vertex in spinor electrodynamics, Phys. Rev. 177, 2426 (1969),
doi:10.1103/PhysRev.177.2426.
[2] J. S. Bell and R. Jackiw, A PCAC puzzle: π0 → γγ in the σ-model, Nuov Cim A 60, 47
(1969), doi:10.1007/BF02823296.
[3] N. P. Armitage, E. J. Mele and A. Vishwanath, Weyl and Dirac semimet-
als in three-dimensional solids, Rev. Mod. Phys. 90, 015001 (2018),
doi:10.1103/RevModPhys.90.015001.
[4] B. Yan and C. Felser, Topological materials: Weyl semimetals, Annu. Rev. Condens. Matter
Phys. 8, 337 (2017), doi:10.1146/annurev-conmatphys-031016-025458.
[5] A. Rycerz, J. Tworzydło and C. W. J. Beenakker, Valley filter and valley valve in graphene,
Nat. Phys. 3, 172 (2007), doi:10.1038/nphys547.
[6] S. A. Parameswaran, T. Grover, D. A. Abanin, D. A. Pesin and A. Vishwanath, Probing the
chiral anomaly with nonlocal transport in three-dimensional topological semimetals, Phys.
Rev. X 4, 031035 (2014), doi:10.1103/PhysRevX.4.031035.
[7] G. Bednik, K. S. Tikhonov and S. V. Syzranov, Magnetotransport and intern-
odal tunnelling in Weyl semimetals, Phys. Rev. Research 2, 023124 (2020),
doi:10.1103/physrevresearch.2.023124.
[8] H. B. Nielsen and M. Ninomiya, The Adler-Bell-Jackiw anomaly and Weyl fermions in a
crystal, Phys. Lett. B 130, 389 (1983), doi:10.1016/0370-2693(83)91529-0.
[9] K. Fukushima, D. E. Kharzeev and H. J. Warringa, Chiral magnetic effect, Phys. Rev. D 78,
074033 (2008), doi:10.1103/PhysRevD.78.074033.
[10] D. T. Son and N. Yamamoto, Kinetic theory with Berry curvature from quantum field theories,
Phys. Rev. D 87, 085016 (2013), doi:10.1103/PhysRevD.87.085016.
[11] A. A. Burkov, Chiral anomaly and diffusive magnetotransport in Weyl metals, Phys. Rev.
Lett. 113, 247203 (2014), doi:10.1103/PhysRevLett.113.247203.
[12] B. Z. Spivak and A. V. Andreev, Magnetotransport phenomena related to the chiral anomaly
in Weyl semimetals, Phys. Rev. B 93, 085107 (2016), doi:10.1103/PhysRevB.93.085107.
[13] J. Behrends and J. H. Bardarson, Strongly angle-dependent magnetoresistance in
Weyl semimetals with long-range disorder, Phys. Rev. B 96, 060201 (2017),
doi:10.1103/PhysRevB.96.060201.
20
SciPost Phys. 11, 046 (2021)
[14] R. D. dos Reis, M. O. Ajeesh, N. Kumar, F. Arnold, C. Shekhar, M. Naumann, M. Schmidt, M.
Nicklas and E. Hassinger, On the search for the chiral anomaly in Weyl semimetals: the neg-
ative longitudinal magnetoresistance, New J. Phys. 18, 085006 (2016), doi:10.1088/1367-
2630/18/8/085006.
[15] J. Xiong, S. K. Kushwaha, T. Liang, J. W. Krizan, M. Hirschberger, W. Wang, R. J. Cava
and N. P. Ong, Evidence for the chiral anomaly in the Dirac semimetal Na3 Bi, Science 350,
413 (2015), doi:10.1126/science.aac6089.
[16] M. Hirschberger, S. Kushwaha, Z. Wang, Q. Gibson, S. Liang, C. A. Belvin, B. A. Bernevig,
R. J. Cava and N. P. Ong, The chiral anomaly and thermopower of Weyl fermions in the
half-Heusler GdPtBi, Nat. Mater. 15, 1161 (2016), doi:10.1038/nmat4684.
[17] S. Liang, J. Lin, S. Kushwaha, J. Xing, N. Ni, R. J. Cava and N. P. Ong, Experimental tests
of the chiral anomaly magnetoresistance in the Dirac-Weyl semimetals Na3 Bi and GdPtBi,
Phys. Rev. X 8, 031002 (2018), doi:10.1103/PhysRevX.8.031002.
[18] C. Zhang et al., Room-temperature chiral charge pumping in Dirac semimetals, Nat. Commun.
8, 13741 (2017), doi:10.1038/ncomms13741.
[19] Q. Ma et al., Direct optical detection of Weyl fermion chirality in a topological semimetal,
Nat. Phys. 13, 842 (2017), doi:10.1038/nphys4146.
[20] C.-K. Chan, N. H. Lindner, G. Refael and P. A. Lee, Photocurrents in Weyl semimetals, Phys.
Rev. B 95, 041104 (2017), doi:10.1103/PhysRevB.95.041104.
[21] B. Sbierski, G. Pohl, E. J. Bergholtz and P. W. Brouwer, Quantum transport of dis-
ordered Weyl semimetals at the nodal point, Phys. Rev. Lett. 113, 026602 (2014),
doi:10.1103/PhysRevLett.113.026602.
[22] M. Buchhold, S. Diehl and A. Altland, Vanishing density of states in weakly disordered Weyl
semimetals, Phys. Rev. Lett. 121, 215301 (2018), doi:10.1103/PhysRevLett.121.215301.
[23] S. V. Syzranov and L. Radzihovsky, High-dimensional disorder-driven phenomena in Weyl
semimetals, semiconductors, and related systems, Annu. Rev. Condens. Matter Phys. 9, 35
(2018), doi:10.1146/annurev-conmatphys-033117-054037.
[24] J. H. Wilson, D. A. Huse, S. Das Sarma and J. H. Pixley, Avoided quantum criticality in
exact numerical simulations of a single disordered Weyl cone, Phys. Rev. B 102, 100201
(2020), doi:10.1103/PhysRevB.102.100201.
[25] C.-Z. Chen, J. Song, H. Jiang, Q.-f. Sun, Z. Wang and X. C. Xie, Disorder and
metal-insulator transitions in Weyl semimetals, Phys. Rev. Lett. 115, 246603 (2015),
doi:10.1103/PhysRevLett.115.246603.
[26] E. V. Gorbar, V. A. Miransky, I. A. Shovkovy and P. O. Sukhachov, Origin of dis-
sipative Fermi arc transport in Weyl semimetals, Phys. Rev. B 93, 235127 (2016),
doi:10.1103/PhysRevB.93.235127.
[27] A. Igarashi and M. Koshino, Magnetotransport in Weyl semimetal nanowires, Phys. Rev. B
95, 195306 (2017), doi:10.1103/PhysRevB.95.195306.
[28] G. Resta, S.-T. Pi, X. Wan and S. Y. Savrasov, High surface conductivity of Fermi-arc electrons
in Weyl semimetals, Phys. Rev. B 97, 085142 (2018), doi:10.1103/PhysRevB.97.085142.
21
SciPost Phys. 11, 046 (2021)
[29] N. Bovenzi, M. Breitkreiz, T. E. O’Brien, J. Tworzydło and C. W. J. Beenakker, Twisted Fermi
surface of a thin-film Weyl semimetal, New J. Phys. 20, 023023 (2018), doi:10.1088/1367-
2630/aaaa90.
[30] E. V. Gorbar, V. A. Miransky, I. A. Shovkovy and P. O. Sukhachov, Hydrodynamic electron
flow in a Weyl semimetal slab: Role of Chern-Simons terms, Phys. Rev. B 97, 205119 (2018),
doi:10.1103/PhysRevB.97.205119.
[31] P. O. Sukhachov, M. V. Rakov, O. M. Teslyk and E. V. Gorbar, Fermi arcs and dc trans-
port in nanowires of Dirac and Weyl semimetals, Ann. Phys. 532, 1900449 (2020),
doi:10.1002/andp.201900449.
[32] V. Kaladzhyan and J. H. Bardarson, Quantized Fermi arc mediated trans-
port in Weyl semimetal nanowires, Phys. Rev. B 100, 085424 (2019),
doi:10.1103/PhysRevB.100.085424.
[33] J. Behrends, R. Ilan and J. H. Bardarson, Anomalous conductance scal-
ing in strained Weyl semimetals, Phys. Rev. Research 1, 032028 (2019),
doi:10.1103/PhysRevResearch.1.032028.
[34] M. Breitkreiz and P. W. Brouwer, Large Contribution of Fermi Arcs to the
conductivity of topological metals, Phys. Rev. Lett. 123, 066804 (2019),
doi:10.1103/PhysRevLett.123.066804.
[35] C. Zhang et al., Ultrahigh conductivity in Weyl semimetal NbAs nanobelts, Nat. Mater. 18,
482 (2019), doi:10.1038/s41563-019-0320-9.
[36] M. Breitkreiz, Parabolic Hall effect due to copropagating surface modes, Phys. Rev. Research
2, 012071 (2020), doi:10.1103/PhysRevResearch.2.012071.
[37] A. A. Burkov, Anomalous Hall effect in Weyl metals, Phys. Rev. Lett. 113, 187202 (2014),
doi:10.1103/PhysRevLett.113.187202.
[38] T. Suzuki, R. Chisnell, A. Devarakonda, Y.-T. Liu, W. Feng, D. Xiao, J. W. Lynn and J. G.
Checkelsky, Large anomalous Hall effect in a half-Heusler antiferromagnet, Nat. Phys. 12,
1119 (2016), doi:10.1038/nphys3831.
[39] H. Li and G. Ouyang, Nonlocal transport in superconducting heterostructures based on Weyl
semimetals, Phys. Rev. B 100, 085410 (2019), doi:10.1103/PhysRevB.100.085410.
[40] P. M. Perez-Piskunow, N. Bovenzi, A. R. Akhmerov and M. Breitkreiz, Code
and data associated with the paper “Chiral Anomaly Trapped in Weyl Metals:
Nonequilibrium Valley Polarization at Zero Magnetic Field”, Zenodo v1.0.1, (2021),
doi:https://doi.org/10.5281/zenodo.4704325.
[41] W. Kohn and J. M. Luttinger, Quantum theory of electrical transport phenomena, Phys. Rev.
108, 590 (1957), doi:10.1103/PhysRev.108.590.
[42] K.-Y. Yang, Y.-M. Lu and Y. Ran, Quantum Hall effects in a Weyl semimetal:
Possible application in pyrochlore iridates, Phys. Rev. B 84, 075129 (2011),
doi:10.1103/PhysRevB.84.075129.
[43] W. E. Lorensen and H. E. Cline, Marching cubes: A high resolution 3D surface construction
algorithm, SIGGRAPH Comput. Graph. 21, 163 (1987), doi:10.1145/37402.37422.
[44] D. Xiao, M.-C. Chang and Q. Niu, Berry phase effects on electronic properties, Rev. Mod.
Phys. 82, 1959 (2010), doi:10.1103/RevModPhys.82.1959.
22