#  Chiral Anomaly Trapped in Weyl Metals: Nonequilibrium Valley Polarization
at Zero Magnetic Field

Pablo M. Perez-Piskunow  Catalan Institute of Nanoscience and Nanotechnology
(ICN2), CSIC and BIST, Campus UAB, Bellaterra, 08193 Barcelona, Spain
Nicandro Bovenzi  Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300
RA Leiden, The Netherlands  Anton R. Akhmerov  Kavli Institute of Nanoscience,
Delft University of Technology, P.O. Box 4056, 2600 GA Delft, The Netherlands
Maxim Breitkreiz<EMAIL_ADDRESS>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 consequence 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.

##  I 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 $\Delta 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 [6, 3, 7]. A direct
consequence of this valley polarization is an enhanced conductivity parallel
to the magnetic field due to the polarization-enhanced occupation disbalance
of countermoving chiral states [8, 9, 10, 11, 12, 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, 16, 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, 22, 23, 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.

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
function 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 connecting surface states of both surfaces are depicted in green.
(b) Top-view on equi-energy contours of the slab spectrum of (a).

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, 27,
28, 29, 30, 31, 32, 33, 34, 35, 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 [37, 38, 39, 34]. 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 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 to the density of magnetic-field induced chiral
states, is $k_{F}/l_{B}^{-2}W$, where $l_{B}=\sqrt{\hbar/eB}\approx
26\,\textrm{nm}\sqrt{1/B[\mathrm{T}]}$ is the magnetic length, and $W$ the
width of the slab. Taking an experimentally realistic value of
$k_{F}=$0.01\text{\,}\mathrm{\SIUnitSymbolAngstrom}^{-1}$$, the density of
anomalous chiral states is larger than that of magnetic-field induced ones for
$W\lesssim 100\,\textrm{nm}$ $/B[\mathrm{T}]$. At a mesoscopic width $W\sim
100\,\textrm{nm}$, the effect of anomalous chiral states is thus comparable to
that of the field-induced chiral states at $B\sim 1\,\mathrm{T}$, which
relevance is commonly accepted and experimentally supported [15, 16, 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 $\sim\exp\big{[}(\xi\Delta
k)^{2}\big{]}$.

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 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.

##  II 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.

###  II.1 General quantum transport approach

We separate a general single-particle Hamiltonian into the free-particle part
$H_{0}$, the additional weak scattering potential $V$, and a time-dependent
electric-field term $e\bm{E}\cdot\bm{r}\,e^{st}$ with the position operator
$\bm{r}$ and an adiabatic time dependence $e^{st}$ with $s\to 0^{+}$,

$H=H_{0}+V+e^{st}e\bm{E}\cdot\bm{r}.$ (1)

The scattering is due to a random configuration of impurities with a vanishing
impurity-averaged potential $\llangle V\rrangle=0$.

We make the ansatz for the full density matrix

$\rho=p+ge^{st},$ (2)

where $p$ is the equilibrium density matrix

$p=\frac{e^{-\beta(H_{0}+V)}}{\mathrm{Tr}\,e^{-\beta(H_{0}+V)}}$ (3)

and $g$ the non-equilibrium correction. Inserting into the von Neumann
equation for the density matrix,

$i\partial_{t}\rho=[H,\rho],$ (4)

and expanding to first order in $E$ we obtain

$i\,s\,g=[e\bm{E}\cdot\bm{r},p]+[H_{0}+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 $H_{0}$ eigenstates
$|{\bm{\kappa}}\rangle$, where ${\bm{\kappa}}$ combines the quantum numbers.
Off-diagonal and diagonal elements read, respectively,

$\displaystyle\begin{split}(E_{\bm{\kappa}}-E_{{\bm{\kappa^{\prime}}}}-is)\,g_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}=(g_{\bm{\kappa}}-g_{\bm{\kappa^{\prime}}})V_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}+C_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}\\\
+\sum_{{\bm{\kappa}}^{\prime\prime}\neq{\bm{\kappa^{\prime}}},{\bm{\kappa}}}\left(g_{{\bm{\kappa}}{\bm{\kappa}}^{\prime\prime}}V_{{\bm{\kappa}}^{\prime\prime}{\bm{\kappa}}^{\prime}}-V_{{\bm{\kappa}}{\bm{\kappa}}^{\prime\prime}}g_{{\bm{\kappa}}^{\prime\prime}{\bm{\kappa}}^{\prime}}\right),\end{split}$
(6) $\displaystyle-
isg_{\bm{\kappa}}=C_{\bm{\kappa}}+\sum_{{\bm{\kappa^{\prime}}}\neq{\bm{\kappa}}}(g_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}V_{{\bm{\kappa^{\prime}}}{\bm{\kappa}}}-V_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}g_{{\bm{\kappa^{\prime}}}{\bm{\kappa}}}),$
(7)

where the notation is
$A_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}=\langle{\bm{\kappa}}|A|{\bm{\kappa^{\prime}}}\rangle$,
$A_{\bm{\kappa}}=A_{{\bm{\kappa}}{\bm{\kappa}}}$,
$E_{\bm{\kappa}}=\langle{\bm{\kappa}}|H_{0}|{\bm{\kappa}}\rangle$. The field-
dependent term

$\displaystyle C_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}$
$\displaystyle=e\bm{E}\cdot\left[\bm{r},p\right]_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}$
(8)

expands in powers of $V$ starting with the zeroth order,

$C^{(0)}_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}=e\bm{E}\cdot\bm{r}_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}\left[n_{F}(E_{\bm{\kappa}})-n_{F}(E_{{\bm{\kappa^{\prime}}}})\right],$
(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_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}=\frac{(g_{\bm{\kappa}}-g_{\bm{\kappa^{\prime}}})V_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}}{E_{\bm{\kappa}}-E_{{\bm{\kappa^{\prime}}}}-is}.$
(10)

Inserting into (7), taking the adiabatic limit $s\to 0^{+}$, and applying
disorder averaging $\llangle\dots\rrangle$ we obtain

$C^{(0)}_{\bm{\kappa}}=i\,2\pi\sum_{{\bm{\kappa^{\prime}}}^{\prime}\neq{\bm{\kappa}}}\delta(E_{\bm{\kappa}}-E_{{\bm{\kappa^{\prime}}}})\left\llangle|V_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}|^{2}\right\rrangle(g_{{\bm{\kappa^{\prime}}}}-g_{{\bm{\kappa}}}).$
(11)

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
$\bm{k}=(k_{x},k_{y})$ eigenstates. The field term (9) becomes

$C^{(0)}_{\bm{\kappa}}=i\,e\bm{E}\cdot\bm{v}_{\bm{\kappa}}\,n_{F}^{\prime}(E_{\bm{\kappa}}),$
(12)

where $\bm{v}_{\bm{\kappa}}=\partial_{\bm{\kappa}}E_{\bm{\kappa}}$ is the
velocity.

Making the ansatz

$g_{\bm{\kappa}}=-e\bm{E}\cdot\bm{\Lambda}_{\bm{\kappa}}\,n_{F}^{\prime}(E_{\bm{\kappa}})$
(13)

(11) simplifies to

$n_{F}^{\prime}(E_{\bm{\kappa}})\bm{v}_{\bm{\kappa}}=2\pi\sum_{{\bm{\kappa^{\prime}}}\neq{\bm{\kappa}}}n_{F}^{\prime}(E_{\bm{\kappa}})\delta(E_{\bm{\kappa}}-E_{{\bm{\kappa^{\prime}}}})\\\
\times\left\llangle|V_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}|^{2}\right\rrangle(\bm{\Lambda}_{{\bm{\kappa}}}-\bm{\Lambda}_{{\bm{\kappa^{\prime}}}}),$
(14)

known as Boltzmann equation (BE), to be solved with respect to the vector-
valued state-resolved transport mean free paths $\bm{\Lambda}_{\bm{\kappa}}$,
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 $\bm{\Lambda}_{\bm{\kappa}}$
in (14) has an eigenvalue zero for a ${\bm{\kappa}}$ independent vector, the
solution is generally determined up to a constant

$\bm{c}=\frac{\sum_{\bm{\kappa}}n_{F}^{\prime}(E_{\bm{\kappa}})\bm{\Lambda}_{\bm{\kappa}}}{\sum_{\bm{\kappa}}n_{F}^{\prime}(E_{\bm{\kappa}})}.$
(15)

Particle conservation however requires $\sum_{\bm{\kappa}}g_{\bm{\kappa}}=0$,
which fixes the constant to $\bm{c}=0$.

The current-density expectation value reads

$\bm{j}=-e\frac{1}{V}\sum_{\bm{\kappa}}\bm{v}_{\bm{\kappa}}\,g_{\bm{\kappa}},$
(16)

where $V$ is the system volume. The conductivity tensor $\sigma$, defined as
$\bm{j}=\sigma\bm{E}$, becomes, using (13),

$\sigma=e^{2}\frac{1}{V}\sum_{\bm{\kappa}}n_{F}^{\prime}(E_{\bm{\kappa}})\bm{v}_{\bm{\kappa}}\otimes\bm{\Lambda}_{{\bm{\kappa}}}.$
(17)

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.

###  II.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 ${\bm{\kappa}}=(\bm{k},b)$, where $\bm{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_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}$, 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

$C_{\bm{\kappa}}=\sum_{{\bm{\kappa}}^{\prime}\neq{\bm{\kappa}}}\delta(E_{\bm{\kappa}}-E_{{\bm{\kappa^{\prime}}}})|V_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}|^{2}(g_{{\bm{\kappa^{\prime}}}}-g_{{\bm{\kappa}}})\\\
+\sum_{{\bm{\kappa^{\prime}}}\neq{\bm{\kappa}}}\delta(E_{\bm{\kappa}}-E_{{\bm{\kappa^{\prime}}}})(C_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}V_{{\bm{\kappa^{\prime}}}{\bm{\kappa}}}-V_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}C_{{\bm{\kappa^{\prime}}}{\bm{\kappa}}}).$
(18)

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,

$C^{(1)}_{{\bm{\kappa}}{\bm{\kappa^{\prime}}}}=e\bm{E}\cdot\sum_{{\bm{\kappa}}^{\prime\prime}}\left[\bm{r}_{{\bm{\kappa}}{\bm{\kappa}}^{\prime\prime}}V_{{\bm{\kappa}}^{\prime\prime}{\bm{\kappa}}^{\prime}}\frac{n_{F}(E_{{\bm{\kappa}}^{\prime\prime}})-n_{F}(E_{{\bm{\kappa}}^{\prime}})}{E_{{\bm{\kappa}}^{\prime\prime}}-E_{{\bm{\kappa}}^{\prime}}}\right.\\\
\left.-\bm{r}_{{\bm{\kappa}}^{\prime\prime}{\bm{\kappa}}^{\prime}}V_{{\bm{\kappa}}{\bm{\kappa}}^{\prime\prime}}\frac{n_{F}(E_{{\bm{\kappa}}})-n_{F}(E_{{\bm{\kappa}}^{\prime\prime}})}{E_{{\bm{\kappa}}}-E_{{\bm{\kappa}}^{\prime\prime}}}\right].$

While it also vanishes upon averaging on the left-hand side (lhs) of (18), the
new term on the rhs becomes

$e\bm{E}\cdot\sum_{{\bm{\kappa}}^{\prime}\neq{\bm{\kappa}}}\sum_{{\bm{\kappa}}^{\prime\prime}}\delta(E_{\bm{\kappa}}-E_{{\bm{\kappa}}^{\prime}})\frac{n_{F}(E_{{\bm{\kappa}}})-n_{F}(E_{{\bm{\kappa}}^{\prime\prime}})}{E_{{\bm{\kappa}}}-E_{{\bm{\kappa}}^{\prime\prime}}}\\\
\times\left[V_{{\bm{\kappa}}^{\prime\prime}{\bm{\kappa}}^{\prime}}V_{{\bm{\kappa}}^{\prime}{\bm{\kappa}}}\bm{r}_{{\bm{\kappa}}{\bm{\kappa}}^{\prime\prime}}+V_{{\bm{\kappa}}{\bm{\kappa}}^{\prime}}V_{{\bm{\kappa}}{\bm{\kappa}}^{\prime\prime}}\bm{r}_{{\bm{\kappa}}^{\prime\prime}{\bm{\kappa}}^{\prime}}\right].$

In this sum there are terms that are proportional to
$|V_{{\bm{\kappa}}{\bm{\kappa}}^{\prime}}|^{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 $\bar{\Lambda}$. This correction can thus be neglected
only if the width is much smaller,

$W\ll\bar{\Lambda}.$ (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.

##  III Weyl-semimetal slab model

Figure 2:  Position-$z$-resolved probability density $|\psi|^{2}=|\langle
z|\psi\rangle|^{2}$ of Fermi-level states. Colors indicate surface states at
the lower boundary (red) and upper boundary (blue), the anomalous 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,\;\mu_{b}=0.3\,t$.

We consider a minimal lattice model of a Weyl semimetal [42],

$\displaystyle H_{0}({\bm{\kappa}})$
$\displaystyle=t\sigma_{x}\sin{k_{z}}+t\sigma_{y}\sin{k_{y}}+m_{\bm{\kappa}}\sigma_{z},$
(20)

where $m_{\bm{\kappa}}=t(2+\cos{\beta}-\cos{k_{x}}-\cos{k_{y}}-\cos{k_{z}})$
and $\sigma_{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 $\bm{k}=\pm\bm{\beta}$, where
$\bm{\beta}=\beta\hat{\bm{x}}$ corresponds to a time-reversal breaking
magnetization. We consider a “good” Weyl semimetal with a cone separation
$\beta\sim 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
k_{z}\to(\delta_{i,j+1}+\delta_{i,j-1})/2$, where $i=0,1,\dots,(W-1)$ is the
site number, corresponding to the discrete position in $z$,

$z\equiv i-\frac{W-1}{2},$ (21)

in units of the lattice constant which is set to one. We furthermore add a
boundary potential $\mu_{b}$ at the surface layers of the slab, which main
effect is to bend the Fermi-arc surface states. We label the eigenstates by
${\bm{\kappa}}=(\bm{k},b)$ where $\bm{k}=(k_{x},k_{y})$ are the continuous in-
plane momenta and $b$ denotes the $2W$ modes at each value of $\bm{k}$.

The eigenstates $|\psi_{\bm{\kappa}}\rangle$ and eigenenergies
$E_{\bm{\kappa}}$ of the slab are obtained from exact diagonalization of the
Hamiltonian at a fixed in-plane momentum $\bm{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 $\bm{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
$|\psi|^{2}\equiv|\langle z|\psi\rangle|^{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 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 $\sim k_{F}W/\pi$ increases
with the width, there is always only a single contour that contains the
surface and chiral bulk states.

##  IV Scattering

###  IV.1 Disorder potential

Figure 3: Plot of the function $M(z,z^{\prime},\xi,W)$ (a) and $\int
dz^{\prime}\,M(z,z^{\prime},\xi,W)$ (b) at $\xi=3$ and $W=100$. The width of
the diagonal peak in (a) and the region of reduced weight in (b) are both set
by $\xi$.

We model the disorder by static Gaussian potentials,

$\displaystyle V$
$\displaystyle=\sum_{\alpha}U_{\alpha}\phi(\bm{r}-\bm{r}_{\alpha}),$
$\displaystyle\phi(\bm{r})$ $\displaystyle=e^{-|\bm{r}|^{2}/2\xi^{2}},$ (22)

where the sum runs over the Gaussian’s with a characteristic width $\xi$,
random and uncorrelated potential magnitudes $U_{\alpha}\in[-\delta,\delta]$,
and random positions $\bm{r}_{\alpha}$.

The disorder potential enters the BE (14) in the form of the scattering rate
between two energy eigenstates
$Q({\bm{\kappa}},{\bm{\kappa^{\prime}}})=2\pi\delta(E_{\bm{\kappa}}-E_{F})q({\bm{\kappa}},{\bm{\kappa^{\prime}}})$,
where

$q({\bm{\kappa}},{\bm{\kappa^{\prime}}})=\left\llangle|\langle\psi_{\bm{\kappa}}|V|\psi_{\bm{\kappa^{\prime}}}\rangle|^{2}\right\rrangle$
(23)

and $\llangle\dots\rrangle$ 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

$q({\bm{\kappa}},{\bm{\kappa^{\prime}}})=\frac{(2\pi\xi^{2})^{3}\delta^{2}n_{i}}{3L_{x}L_{y}}e^{-\xi^{2}(\bm{k}-\bm{k}^{\prime})^{2}/2}\sum_{z,z^{\prime}}M(z,z^{\prime},\xi,W)\\\
\times\psi_{\bm{\kappa^{\prime}}}^{\dagger}(z)\psi_{\bm{\kappa}}^{\phantom{\dagger}}(z)\psi_{\bm{\kappa}}^{\dagger}(z^{\prime})\psi_{\bm{\kappa^{\prime}}}^{\phantom{\dagger}}(z^{\prime})$
(24)

where $\psi_{\bm{\kappa}}^{\phantom{\dagger}}(z)=\langle
z|\psi_{\bm{\kappa}}\rangle$ and $n_{i}$ is the impurity concentration. A
detailed derivation can be found in the Appendix. The function
$M(z,z^{\prime},\xi,W)$ is given by

$M(z,z^{\prime},\xi,W)=\frac{e^{-(z-z^{\prime})^{2}/4\xi^{2}}}{4\sqrt{\pi}\xi}\\\
\times\left[\mathrm{erf}\left(\frac{W+z+z^{\prime}}{2\xi}\right)+\mathrm{erf}\left(\frac{W-z-z^{\prime}}{2\xi}\right)\right],$

where the error function is defined as
$\mathrm{erf}\left(x\right)=(2/\sqrt{\pi})\int_{0}^{x}e^{-t^{2}}\mathrm{d}t$.
As illustrated in Fig. 3, $M(z,z^{\prime},\xi,W)$ is mainly the
$(z-z^{\prime})$ 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 $\xi$ 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 $\xi\ll 1$ one obtains
$M(z,z^{\prime},\xi,W)=\delta(z-z^{\prime})$, in which case

$q({\bm{\kappa}},{\bm{\kappa^{\prime}}})\xlongrightarrow{\xi\ll
1}\frac{(2\pi\xi^{2})^{3}\delta^{2}n_{i}}{3L_{x}L_{y}}e^{-\xi^{2}(\bm{k}-\bm{k}^{\prime})^{2}}\\\
\times\sum_{z}\,|\psi_{\bm{\kappa^{\prime}}}^{\phantom{\dagger}}(z)|^{2}|\psi_{\bm{\kappa}}^{\phantom{\dagger}}(z)|^{2}.$

Despite the rather complex form, the impurity scattering is fully determined
by two parameters — the real-space impurity width $\xi$ and the overall
impurity strength (set by $\delta^{2}n_{i}$), the latter will be in the
following quantified by the average mean free path $l$, defined below. The
impurity width $\xi$ essentially sets the momentum-space range of most
scattering processes to $|\bm{k}-\bm{k}^{\prime}|\lesssim\xi^{-1}$, hence for
a large value of $\xi$ scattering between states that are far apart from each
other in the in-plane momentum $\bm{k}$ is exponentially suppressed
$\sim\exp[-\xi^{2}(\bm{k}-\bm{k}^{\prime})^{2}]$.

###  IV.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

$l=\left\langle\,|{\bm{v}_{\bm{\kappa}}}|\left(\sum_{{\bm{\kappa^{\prime}}}}\delta(E_{\bm{\kappa}}-E_{F})q({\bm{\kappa}},{\bm{\kappa^{\prime}}})\right)^{-1}\,\right\rangle,$
(25)

where the Fermi-surface average is given by

$\displaystyle\langle\dots\rangle$
$\displaystyle=\frac{1}{N}\sum_{\bm{\kappa}}\delta(E_{\bm{\kappa}}-E_{F})(\dots);$
$\displaystyle N$
$\displaystyle=\sum_{\bm{\kappa}}\delta(E_{\bm{\kappa}}-E_{F}),$ (26)

and $N$ is the number of states at the Fermi level. Note that for a nearly
constant $|{\bm{v}_{\bm{\kappa}}}|$ 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\in[b,\;s]$, where $b$ ($s$) denotes bulk (surface) states, we define the
scattering length

$l_{ij}=\left\langle|{\bm{v}_{\bm{\kappa}}}|\left(\sum^{j}_{{\bm{\kappa^{\prime}}}}\delta(E_{\bm{\kappa}}-E_{F})q({\bm{\kappa}},{\bm{\kappa^{\prime}}})\right)^{-1}\right\rangle_{i}$
(27)

where the sum runs over the type $j$ of states and the averaging is analogous
to (26) but only over type $i$ of states
($\sum_{\bm{\kappa}}\to\sum^{i}_{\bm{\kappa}}$).

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\gg N_{s}$ and (ii) the
number of states being dominated by surface states, $N\sim N_{s}$. In regime
(ii) the number of bulk states $N_{b}=N-N_{s}$, is given by the number of
chiral bulk states. Thus the scaling of $N_{b}$ with the slab width reads

$N_{b}\sim\begin{cases}W&\text{(i)}\\\ 1&\text{(ii)}.\end{cases}$ (28)

The scaling of the scattering rate (24) is governed by the $z$ dependence of
the wavefunctions. A normalized surface state with a penetration depth
$\lambda\sim\beta^{-1}\sim 1$ and a normalized bulk state are of the form

$\displaystyle\psi_{s}(z)$
$\displaystyle=\sqrt{\frac{2}{\lambda}}\,e^{-z/\lambda},$
$\displaystyle\psi_{b}(z)$ $\displaystyle=1/\sqrt{W},$ (29)

respectively. Consequently, for $W\gg\lambda$ (which we always consider), we
can estimate the scaling of the scattering rate between the different types of
states as

$q({\bm{\kappa}},{\bm{\kappa}}^{\prime})\sim\begin{cases}1&\text{same
surface}\\\ \frac{1}{W}&\text{bulk-surface, bulk-bulk}\\\ 0&\text{opposite
surfaces}.\end{cases}$ (30)

From this, (28), and (27) the width dependence of the scattering probabilities
summarizes to

| $l_{bb}/l$ | $l_{ss}/l$ | $l_{sb}/l$ | $l_{bs}/l$
---|---|---|---|---
(i) | $1$ | $\tfrac{N_{b}}{N_{s}W}\sim 1$ | $1$ | $\tfrac{N_{b}}{N_{s}}\sim W$
(ii) | $W\tfrac{N_{s}}{N_{b}}$ | $1$ | $W\tfrac{N_{s}}{N_{b}}$ | $W$
(31)

Most importantly, in the regime (i) the scattering probability from bulk to
surface ($\propto 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.

##  V 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 (33) below.

###  V.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\approx t/\hbar$,

$\sigma_{0}=\frac{e^{2}n\,l\,v}{3}.$ (32)

This is the result one would expect to find for a system with point-like
impurities ($\xi\to 0$) and only bulk states arranged in form of a spherical
Fermi surface.

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), $N_{s}+N_{c}$, relative to the total number of
states $N$ is (a) $(N_{s}+N_{c})/N=0.14$, (c) $(N_{s}+N_{c})/N=0.68$, and (e)
$(N_{s}+N_{c})/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
$\xi$. The green curves show power-law fits, the exponent is indicated in the
plot. Other parameters are $\beta=1.5,\;\;\mu_{b}=0.3\,t$.

The dependence of the slab conductivity on the width of the impurity potential
$\xi$ is summarized in Fig. 4. For $\xi>1$ the conductivity is well fitted by
a power-law dependence on $\xi$, with an exponent between 2 and 3. An
exception is found for $\sigma_{yy}$ at $E_{F}\ll t$, which shows a weak $\xi$
dependence at a strongly enhanced conductivity at all $\xi$. 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.

Figure 5:  Width dependence of the conductivity at (a) $E_{F}=0.3\,t$,
$\xi=0.1$, (b) $E_{F}=0.3\,t$, $\xi=5$, (c) $E_{F}=0.03\,t$, $\xi=0.1$, and
$E_{F}=0.03\,t$, $\xi=5$. Other parameters are $\beta=1.5,\;\mu_{b}=0.3\,t$.
The dashed lines in (b) indicate $W^{-1}$ dependence and in (c) and (d) linear
dependence.

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 $\sigma_{0}$ — in this regime the slab thus resembles a
conventional metal. At a large $\xi$, 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,
$\sigma_{yy}$, linearly _increases_ with the width $W$.

Figure 6:  Width dependence of the conductivity at $\xi=5$ and
$E_{F}=0.15\,t$, $\mu_{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 $Wk_{F}/\pi\lesssim 1$. As soon as at least one normal bulk contour
appears, the conductivity jumps to a significantly lower value and becomes
decreasing in $W$.

###  V.2 Valley polarization

Besides the conductivity it is interesting to explore the non-equilibrium
occupation difference of the valleys, which occurs when the electric field
points along the valley separation, $\bm{E}=\hat{\bm{x}}E$. The average
occupation of the valley at $k_{x}=\pm\beta$ is given by
$\sum^{\pm}_{\bm{\kappa}}g_{\bm{\kappa}}/(N_{b}/2)$, where the sum runs over
all bulk states at the valley $\pm$, and $N_{b}/2=\sum^{\pm}_{\bm{\kappa}}$ 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 $\delta(E_{F}-E_{\bm{\kappa}})eEl$ due to the mean free
motion in the electric field,

$\displaystyle\chi$
$\displaystyle\equiv\dfrac{\sum^{+}_{{\bm{\kappa}}}g_{\bm{\kappa}}}{\sum^{+}_{{\bm{\kappa}}}\delta(E_{F}-E_{\bm{\kappa}})eEl}-\dfrac{\sum^{-}_{{\bm{\kappa}}}g_{\bm{\kappa}}}{\sum^{-}_{{\bm{\kappa}}}\delta(E_{F}-E_{\bm{\kappa}})eEl}$
(33)
$\displaystyle=\frac{\langle\Lambda^{x}_{\bm{\kappa}}\rangle_{+}-\langle\Lambda^{x}_{\bm{\kappa}}\rangle_{-}}{l},$

where in the second line we used Eq. (13) and defined the average over valley
bulk states
$\langle\dots\rangle_{\pm}=\sum^{\pm}_{{\bm{\kappa}}}\delta(E_{F}-E_{\bm{\kappa}})\dots/\sum^{\pm}_{{\bm{\kappa}}}\delta(E_{F}-E_{\bm{\kappa}})$.

Figure 7:  Valley polarization $\chi$, defined in Eq. (33). (a) Dependence of
$\chi$ on the impurity range $\xi$ at a width $W=201$. (b) Width dependence of
$\chi$ at $\xi=5$. Other parameters are $\mu_{b}=0.3\,t$, $\beta=1.5$, and
$E_{F}=0.3\,t$.

A representative result for the impurity range and width dependence of the
valley polarization is shown in Fig. 7. The valley polarization shows the
power-law $\xi$ dependence, similar to the conductivity. However, unlike for
the conductivity, there is no significant width dependence at large $\xi$
(here already at $\xi\gtrsim 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
width above $\gtrsim l\ e^{(\xi\Delta k)^{2}}$, which becomes exponentially
large for $\xi\Delta k\gtrsim 1$.

##  VI 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 $\xi$ and
$W$ allow to identify the main mechanisms of these effects, which we now
systematically discuss.

###  VI.1 Impurity-range dependence

When the impurity range $\xi$ increases, the scattering rate (24) between two
countermovers separated by $\Delta k$ becomes exponentially suppressed,
$\sim\exp[-(\xi\Delta 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 ${\cal N}$ states labeled $i\in[1,{\cal N}]$ with arbitrary
velocities $\bm{v}_{i}$. The BE (14) is of the form (summation over repeated
indices assumed)

$\bm{v}_{i}=M_{ij}\bm{\Lambda}_{j},$ (34)

where $M_{ij}=\delta_{ij}\sum_{k}q_{ik}-q_{ij}$ is given by the scattering
rates between states, $q_{ij}$. We assume scattering only between the nearest
neighbors, with the rate $q_{\mathrm{nn}}$, and direct scattering between
countermovers with the rate $q_{\mathrm{d}}$, in which case the matrix $M$
becomes

$M_{ij}=(2q_{\mathrm{nn}}+2q_{\mathrm{d}})\delta_{ij}-q_{\mathrm{nn}}(\delta_{ij+1}+\delta_{i,j-1}),$
(35)

where we used that the direct-scattering part of $q_{ij}$ cancels when
multiplied with $\bm{\Lambda}$. For the nearest-neighbor part
$M_{\mathrm{nn}}$ of $M$ there is a left pseudo-inverse,

$P_{ij}=\frac{1}{q_{\mathrm{nn}}}\frac{|i-j|\left(|i-j|-{\cal
N}\right)}{{2\cal N}},$ (36)

so that $PM_{\mathrm{nn}}\bm{\Lambda}=\bm{\Lambda}$ (note
$\sum_{i}\bm{\Lambda}_{i}=0$ due to particle conservation). With its help, the
full solution of (34) becomes

$\bm{\Lambda}=\left(1+2q_{\mathrm{d}}P\right)^{-1}P\,\bm{v}.$ (37)

We first ignore direct scattering. From the form of $P$, it is clear that the
solution can depend on ${\cal N}$ to maximally the power ${\cal N}^{2}$. In
particular, for a circular velocity arrangement
$\bm{v}=\\{-\cos[2\pi(i-1)/{\cal N}],\;\sin[2\pi(i-1)/{\cal N}]\\}$ the vector
mean free path for ${\cal N}\gg 1$ assumes the value

$\bm{\Lambda}=\frac{1}{q_{\mathrm{nn}}}\frac{(\tilde{\xi}\Delta
k)^{2}}{4}\bm{v},$ (38)

where we have written the number of nearest neighbors as ${\cal
N}=\pi\,\tilde{\xi}\Delta k$, in terms of the spacing between nearest
neighbors $\tilde{\xi}^{-1}$ and the distance between countermovers $\Delta
k$.

Considering the full result with direct scattering in Eq. (37), we see that
the nearest-neighbor scattering contributions to $\bm{\Lambda}$ dominate as
long as

$\frac{(\tilde{\xi}\Delta
k)^{p}}{q_{\mathrm{nn}}}\ll\frac{1}{q_{\mathrm{d}}},$ (39)

where the exponent $p$ depends on the velocity arrangement. In the opposite
limit we instead obtain

$\bm{\Lambda}=\frac{1}{2q_{\mathrm{d}}}\bm{v}.$ (40)

Transferring this insight to the Weyl slab model, the distance between the
nearest neighbors $\tilde{\xi}^{-1}$ corresponds to $\xi^{-1}$. The mean free
path corresponds to the inverse diagonal of $M$ times velocity, $l\sim
v/(q_{\mathrm{d}}+q_{\mathrm{nn}})$. 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

$\bar{\Lambda}=l\,\left[1+(\xi\Delta k)^{p}\right].$ (41)

This explains the general power-law conductivity enhancement with $\xi$ in
Fig. 4, except for $\sigma_{yy}$ in Fig. 4(f), which we discuss separately.

###  VI.2 Valley polarization

The averaged occupation of the two valleys $\pm$ can be expressed as

$\mu_{\pm}=\frac{\sum^{\pm}_{{\bm{\kappa}}}g_{\bm{\kappa}}}{N_{b}/2},$ (42)

where $\sum^{\pm}_{{\bm{\kappa}}}\delta(E_{F}-E_{\bm{\kappa}})=N_{b}/2$. The
occupation difference is related to the valley polarization $\chi$ defined in
Eq. (33),

$\mu_{+}-\mu_{-}=\chi\,eE_{x}l.$ (43)

A single valley has a total velocity in the $x$ direction — the unbalanced
velocity of the chiral states $vN_{c}/2$, where $N_{c}$ is the total number of
chiral states ($N_{c}/2$ in each valley). An electric field in the $x$
direction thus pumps charge between the valleys with the rate $veE_{x}N_{c}$,
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,

$eN_{b}\frac{d\mu_{\pm}}{dt}=0=\pm
veE_{x}N_{c}+eN_{b}\left(\frac{\partial\mu_{\pm}}{\partial
t}\right)_{\text{scat}},$ (44)

where we used $N_{b}\gg N_{c}$. The time-change of $\mu_{+}-\mu_{-}$ due to
scattering is proportional to the occupation difference itself and the
scattering probability, which we quantify by the scattering length $l_{c}$,
hence

$\left(\frac{\partial(\mu_{+}-\mu_{-})}{\partial
t}\right)_{\text{scat}}=-2\frac{\mu_{+}-\mu_{-}}{l_{c}/v}=-2\frac{\chi
eEl}{l_{c}/v}.$ (45)

Together with (44) we obtain

$\chi=\frac{N_{c}}{N_{b}}\,\frac{l_{c}}{l}.$ (46)

For a point-like disorder potential ($\xi\to 0$) the ratio $l_{c}/l$ goes to
one and the valley polarization is small. For larger $\xi$, direct scattering
between the valleys becomes strongly suppressed and the relaxation of
$\mu_{+}-\mu_{-}$ 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 $N_{b}/N_{s}$ so that

$\chi\propto\frac{N_{c}}{N_{s}}.$ (47)

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
$N_{c}$ is a factor $\propto W$ smaller than the total number of bulk states
$N_{b}$. The explanation is that the valley relaxation also becomes suppressed
$\propto 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
$\sim\exp[-(\xi\Delta 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 $\xi\Delta k\gtrsim 1$.

Regarding the strong enhancement of valley polarization with $\xi$, one is
tempted to understand 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 $\chi$ would enter in the form
$l_{ss}(\xi\Delta)^{p}/l\propto 1/W$, with a clear $1/W$ dependence, which we
do not observe. It rather must be the suppression of the bulk-surface
scattering probability with $\xi$, which is plausible in view of the typical
effect of an increasing $\xi$ to increase the relaxation time. We note,
however, that the simplified analytical calculation of Section VI.1 does not
apply since in this case the valleys correspond to two-dimensional pools of
states, while the model of Section VI.1 only considers one-dimensional chains.

###  VI.3 Conductivity in case of a large number of bulk states

For bulk states (excluding chiral bulk states) the small separation of
countermovers in momentum space $\sim k_{F}$ makes their transport length
close to the mean free path $l$ and not significantly enhanced with $\xi$. We
thus approximate the current contribution of bulk states as

$\bm{j}_{n}=\sigma_{0}\bm{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 $\xi$ according to (41),

$\bm{j}_{s}=e^{2}n_{\mathrm{s}}v\,l_{ss}(\xi\Delta k)^{p}\bm{y}E_{y},$ (49)

where $n_{s}=N_{s}/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

$\bm{j}_{c}=ev\hat{\bm{x}}n_{c}(\mu_{+}-\mu_{-})/2=e^{2}v\hat{\bm{x}}n_{c}\chi
l\,E/2,$ (50)

where $n_{c}=N_{c}/V$ is the density of chiral bulk states. Adding all the
current contributions, we obtain

$\displaystyle\frac{\sigma^{(i)}_{xx}}{\sigma_{0}}$
$\displaystyle=1+\frac{3}{2}\frac{N_{c}}{N_{b}}\chi,$ (51)
$\displaystyle\frac{\sigma^{(i)}_{yy}}{\sigma_{0}}$
$\displaystyle=1+3\frac{N_{s}}{N_{b}}\frac{l_{ss}(\xi\Delta k)^{p}}{l}\sim
1+\frac{1}{W}(\xi\Delta k)^{p},$ (52)

where we again used (31). This rough estimate is in qualitative agreement with
the numerical results. Note that $N_{b}$ increases proportional to $W$, which
explains the $W$ dependence of $\sigma_{xx}$ and $\sigma_{yy}$ in Fig. 5(b).
We now understand that the enhancement of $\sigma_{xx}$ and $\sigma_{yy}$ 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\ll W$ (we consider the opposite
limit $l\gg 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\gg W$, or in the crossover regime
$l\sim W$.

###  VI.4 Conductivity, 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

$\frac{\sigma^{(ii)}_{yy}}{\sigma_{0}}=1+\frac{N_{s}}{N}\frac{\Lambda_{s}}{l},$
(53)

where $N_{s}\approx N$ is the number and $\Lambda_{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 $\Lambda_{s}$ is set by the surface-bulk
scattering probability, $\Lambda_{s}\sim l_{sb}$. In Section IV.2 we found
$l_{sb}/l\sim WN_{s}/N_{c}$—Eq. (31)—which leads to

$\frac{\sigma^{(ii)}_{yy}}{\sigma_{0}}\sim W\frac{N_{s}}{N_{c}}.$ (54)

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 $N_{s}$ and $N_{c}$ are $W$ independent in this case.

For large $\xi$, 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 VI.1, the latter should elongate the
full transport length by an additional $l_{ss}(\xi\Delta k)^{p}$. For the
considered parameters, $l_{sb}$ is much larger than this additional part since
$(\xi\Delta k)^{p}\ll WN_{s}/N_{c}$, which explains the weak $\xi$ dependence
in Fig. 4(f).

The transition from (ii) to (i) spoils the strong enhancement of $\sigma_{yy}$
in two ways: First, the ration $N_{s}/N_{c}$ changes to $N_{s}/(N_{c}+N_{b})$
and thus becomes smaller and second, according to Eq. (31), $l_{sb}/l\sim
WN_{s}/N_{b}\to 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 $v_{x}$ averages to zero at each surface separately. Since
$l_{ss}/l\sim 1$, the conductivity in the $x$ direction is not significantly
enhanced.

##  VII 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 $\xi$). Our work highlights the remarkable property of Weyl
semimetals to realize valleys of opposite chirality that are well separated in
momentum space ($\Delta 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[(\xi\Delta k)^{2}]$ and can thus easily reach
macroscopic length scales at realistic values of cone separation $\Delta k$
and inpurity-potential widths $\xi$.

The valley polarization and Fermi-arc surface states lead to a conductivity
enhancement which increases with an increasing width of the disorder potential
and and a decreasing width of the slab ($\sigma\propto 1/W$). Moreover, if the
Fermi energy is reduced towards charge neutrality such 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 ($\gtrsim
100$) number of layers and consequently 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.

###### Acknowledgements.

This research was supported by the European Union Horizon 2020 research and
innovation 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.

## References

  * Adler [1969] S. L. Adler, Axial-Vector Vertex in Spinor Electrodynamics, Phys. Rev. 177(5), 2426 (1969).
  * Bell and Jackiw [1969] J. S. Bell and R. W. Jackiw, A PCAC Puzzle: $\pi^{0}\to\gamma\gamma$ in the $\sigma$-model, Nuovo Cim. 60, 47 (1969).
  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Weyl and Dirac Semimetals in Three Dimensional Solids, Rev. Mod. Phys. 90, 015001 (2018).
  * Yan and Felser [2017] B. Yan and C. Felser, Topological Materials: Weyl Semimetals, Annu. Rev. Condens. Matter Phys. 8, 337 (2017).
  * Rycerz  _et al._ [2007] A. Rycerz, J. Tworzydło, and C. W. Beenakker, Valley filter and valley valve in graphene, Nat. Phys. 3, 172 (2007).
  * Parameswaran  _et al._ [2014] 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).
  * Bednik  _et al._ [2020] G. Bednik, K. S. Tikhonov, and S. V. Syzranov, Magnetotransport and internodal tunnelling in Weyl semimetals, Phys. Rev. Res. 2, 023124 (2020).
  * Nielsen and Ninomiya [1983] H. B. Nielsen and M. Ninomiya, The Adler-Bell-Jackiw anomaly and Weyl fermions in a crystal, Phys. Lett. B 130(6), 389 (1983).
  * Fukushima  _et al._ [2008] K. Fukushima, D. E. Kharzeev, and H. J. Warringa, Chiral magnetic effect, Phys. Rev. D 78, 074033 (2008).
  * Son and Yamamoto [2013] D. T. Son and N. Yamamoto, Kinetic theory with Berry curvature from quantum field theories, Phys. Rev. D 87, 1 (2013).
  * Burkov [2014a] A. A. Burkov, Chiral anomaly and diffusive magnetotransport in weyl metals, Phys. Rev. Lett. 113, 247203 (2014a).
  * Spivak and Andreev [2016] B. Z. Spivak and A. V. Andreev, Magnetotransport phenomena related to the chiral anomaly in Weyl semimetals, Phys. Rev. B 93, 085107 (2016).
  * Behrends and Bardarson [2017] J. Behrends and J. H. Bardarson, Strongly angle-dependent magnetoresistance in Weyl semimetals with long-range disorder, Phys. Rev. B 96, 060201(R) (2017).
  * dos Reis  _et al._ [2016] 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 negative longitudinal magnetoresistance, New J. Phys. 18, 085006 (2016).
  * Xiong  _et al._ [2015] 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 Na$_3$Bi, Science 350, 413 (2015).
  * Hirschberger  _et al._ [2016] 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).
  * Liang  _et al._ [2018] 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 Na3Bi and GdPtBi, Phys. Rev. X 8, 031002 (2018).
  * Zhang  _et al._ [2017] C. Zhang, E. Zhang, W. Wang, Y. Liu, Z.-G. Chen, S. Lu, S. Liang, J. Cao, X. Yuan, L. Tang, Q. Li, C. Zhou, T. Gu, Y. Wu, J. Zou, and F. Xiu, Room-temperature chiral charge pumping in Dirac semimetals, Nat. Commun. 8, 13741 (2017).
  * Ma  _et al._ [2017] Q. Ma, S.-Y. Xu, C.-K. Chan, C.-L. Zhang, G. Chang, Y. Lin, W. Xie, T. Palacios, H. Lin, S. Jia, P. A. Lee, P. Jarillo-Herrero, and N. Gedik, Direct optical detection of Weyl fermion chirality in a topological semimetal, Nat. Phys. 13, 842 (2017).
  * Chan  _et al._ [2017] C.-K. Chan, N. H. Lindner, G. Refael, and P. A. Lee, Photocurrents in Weyl semimetals, Phys. Rev. B 95, 041104(R) (2017).
  * Sbierski  _et al._ [2014] B. Sbierski, G. Pohl, E. J. Bergholtz, and P. W. Brouwer, Quantum Transport of Disordered Weyl Semimetals at the Nodal Point, Phys. Rev. Lett 113, 026602 (2014).
  * Buchhold  _et al._ [2018] M. Buchhold, S. Diehl, and A. Altland, Vanishing Density of States in Weakly Disordered Weyl Semimetals, Phys. Rev. Lett. 121, 215301 (2018).
  * Syzranov and Radzihovsky [2018] 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).
  * Wilson  _et al._ [2020] J. H. Wilson, D. A. Huse, S. D. Sarma, and J. H. Pixley, Avoided quantum criticality in exact numerical simulations of a single disordered Weyl cone, Phys. Rev. B 102, 100201 (2020).
  * Chen  _et al._ [2015] 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), arXiv:1507.00128 .
  * Gorbar  _et al._ [2016] E. V. Gorbar, V. A. Miransky, I. A. Shovkovy, and P. O. Sukhachov, Origin of dissipative Fermi arc transport in Weyl semimetals, Phys Rev B 93, 235127 (2016).
  * Igarashi and Koshino [2017] A. Igarashi and M. Koshino, Magnetotransport in Weyl semimetal nanowires, Phys Rev B 95, 195306 (2017).
  * Resta  _et al._ [2018] 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).
  * Bovenzi  _et al._ [2018] 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).
  * Gorbar  _et al._ [2018] E. V. Gorbar, V. A. Miransky, I. A. Shovkovy, and P. O. Sukhachov, Hydrodynamic electron flow in a Weyl semimetal slab: The role of Chern–Simons terms, Phys. Rev. B 97, 205119 (2018).
  * Sukhachov  _et al._ [2019] P. O. Sukhachov, M. V. Rakov, O. M. Teslyk, and E. V. Gorbar, Fermi arcs and DC transport in nanowires of Dirac and Weyl semimetals, Ann. der Phys. 532, 1900449 (2019).
  * Kaladzhyan and Bardarson [2019] V. Kaladzhyan and J. H. Bardarson, Quantized Fermi-arc-mediated transport in Weyl semimetal nanowires, Phys. Rev. B 100, 085424 (2019).
  * Behrends  _et al._ [2019] J. Behrends, R. Ilan, and J. H. Bardarson, Anomalous conductance scaling in strained Weyl semimetals, Phys. Rev. Res. 1, 032028(R) (2019).
  * Breitkreiz and Brouwer [2019] M. Breitkreiz and P. W. Brouwer, Large contribution of Fermi arcs to the conductivity of topological metals, Phys. Rev. Lett 123, 066804 (2019).
  * Zhang  _et al._ [2019] C. Zhang, Z. Ni, J. Zhang, X. Yuan, Y. Liu, Y. Zou, Z. Liao, Y. Du, A. Narayan, H. Zhang, T. Gu, X. Zhu, L. Pi, S. Sanvito, X. Han, J. Zou, Y. Shi, X. Wan, S. Y. Savrasov, and F. Xiu, Ultrahigh conductivity in Weyl semimetal NbAs nanobelts, Nat. Mater. 18, 482 (2019).
  * Breitkreiz [2020] M. Breitkreiz, Parabolic Hall effect due to Co-Propagating Surface Modes, Phys. Rev. Res. 2, 012071(R) (2020).
  * Burkov [2014b] A. A. Burkov, Anomalous Hall Effect in Weyl Metals, Phys. Rev. Lett 113, 187202 (2014b).
  * Suzuki  _et al._ [2016] 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).
  * Li and Ouyang [2019] H. Li and G. Ouyang, Nonlocal transport in superconducting heterostructures based on Weyl semimetals, Phys. Rev. B 100, 085410 (2019).
  * Perez-Piskunow  _et al._ [2021] 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” at https://zenodo.org/record/4668624 (2021).
  * Kohn and Luttinger [1957] W. Kohn and J. Luttinger, Quantum Theory of Electrical Transport Phenomena, Phys. Rev. 108, 590 (1957).
  * Yang  _et al._ [2011] 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).
  * Lorensen and Cline [1987] W. E. Lorensen and H. E. Cline, Marching cubes: A high resolution 3D surface construction algorithm, SIGGRAPH Comput. Graph 21, 163 (1987).
  * Xiao  _et al._ [2010] D. Xiao, M.-C. Chang, and Q. Niu, Berry phase effects on electronic properties, Rev. Mod. Phys. 82, 1959 (2010).

Supplementary Material

## Derivation of scattering amplitudes

We consider Gaussian-type static impurity potentials,

$\displaystyle V$
$\displaystyle=\sum_{\alpha}U_{\alpha}\phi(\bm{r}-\bm{r}_{\alpha}),$
$\displaystyle\phi(\bm{r})$ $\displaystyle=e^{-|\bm{r}|^{2}/2\xi^{2}},$ (S1)

where the sum runs over impurities with a characteristic width $\xi$, random
and uncorrelated potential magnitudes $U_{\alpha}\in[-\delta,\delta]$, and
random positions $\bm{r}_{\alpha}$.

For our transport consideration we consider the scattering rate
$Q({\bm{\kappa}},{\bm{\kappa^{\prime}}})=2\pi\delta(E_{\bm{\kappa}}-E_{F})q({\bm{\kappa}},{\bm{\kappa^{\prime}}})$
between energy eigenstates $|\psi_{\bm{\kappa}}\rangle$ and
$|\psi_{\bm{\kappa^{\prime}}}\rangle$, which we calculate using Fermi’s Golden
Rule,

$q({\bm{\kappa}},{\bm{\kappa^{\prime}}})=\left\llangle|\langle\psi_{\bm{\kappa}}|V|\psi_{\bm{\kappa^{\prime}}}\rangle|^{2}\right\rrangle,$
(S2)

where the disorder average is defined as

$\displaystyle\llangle\dots\rrangle$
$\displaystyle=\prod_{\alpha}\int_{-\delta}^{\delta}\frac{dU_{\alpha}}{2\delta}\int\frac{d\bm{r}_{\alpha}}{V}(\dots).$
(S3)

We write the normalized wavefunctions as

$\displaystyle\langle\bm{r}\rvert\psi_{\bm{\kappa}}\rangle$
$\displaystyle=\;\frac{1}{\sqrt{L_{x}L_{y}}}e^{i\bm{k}\cdot\bm{\rho}}\,\psi_{\bm{\kappa}}(z),$
(S4)

where $\bm{\rho}=(x,y)$ is the in-plane position and $L_{x}L_{y}$ the in-plane
volume and $\psi_{\bm{\kappa}}(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

$\displaystyle\langle\psi_{\bm{\kappa^{\prime}}}\lvert
V\rvert\psi_{\bm{\kappa}}\rangle$
$\displaystyle=\frac{2\pi\xi^{2}}{L_{x}L_{y}}e^{-\xi^{2}(\bm{k}-\bm{k}^{\prime})^{2}/2}\sum_{\alpha}U_{\alpha}\,\sum_{z}\phi(z-z_{\alpha})\,\psi_{\bm{\kappa^{\prime}}}^{\dagger}(z)\psi_{\bm{\kappa}}^{\phantom{\dagger}}(z).$
(S5)

Inserting into (S2) and using (S3), we obtain

$\displaystyle q({\bm{\kappa}},{\bm{\kappa^{\prime}}})$
$\displaystyle=\frac{4\pi^{2}\xi^{4}}{L_{x}^{2}L_{y}^{2}}e^{-\xi^{2}(\bm{k}-\bm{k}^{\prime})^{2}}\left\llangle\sum_{\alpha,\beta}U_{\alpha}U_{\beta}\,\sum_{z,z^{\prime}}\int\frac{dz_{\alpha}}{W}\phi(z-z_{\alpha})\phi(z^{\prime}-z_{\beta})\,\psi_{\bm{\kappa^{\prime}}}^{\dagger}(z)\psi_{\bm{\kappa}}^{\phantom{\dagger}}(z)\psi_{\bm{\kappa}}^{\dagger}(z^{\prime})\psi_{\bm{\kappa^{\prime}}}^{\phantom{\dagger}}(z^{\prime})\right\rrangle$
(S6)
$\displaystyle=\frac{4\pi^{2}\xi^{4}\delta^{2}}{3L_{x}^{2}L_{y}^{2}}e^{-\xi^{2}(\bm{k}-\bm{k}^{\prime})^{2}}\sum_{\alpha}\,\sum_{z,z^{\prime}}\,\psi_{\bm{\kappa^{\prime}}}^{\dagger}(z)\psi_{\bm{\kappa}}^{\phantom{\dagger}}(z)\psi_{\bm{\kappa}}^{\dagger}(z^{\prime})\psi_{\bm{\kappa^{\prime}}}^{\phantom{\dagger}}(z^{\prime})\,\int\frac{dz_{\alpha}}{W}\phi(z-z_{\alpha})\phi(z^{\prime}-z_{\alpha})$
$\displaystyle=\frac{4\pi^{2}\xi^{4}\delta^{2}n_{i}}{3L_{x}L_{y}}e^{-\xi^{2}(\bm{k}-\bm{k}^{\prime})^{2}}\sum_{z,z^{\prime}}\,\psi_{\bm{\kappa^{\prime}}}^{\dagger}(z)\psi_{\bm{\kappa}}^{\phantom{\dagger}}(z)\psi_{\bm{\kappa}}^{\dagger}(z^{\prime})\psi_{\bm{\kappa^{\prime}}}^{\phantom{\dagger}}(z^{\prime})\int
dz_{i}\,\phi(z-z_{i})\phi(z^{\prime}-z_{i})$
$\displaystyle=\frac{(2\pi\xi^{2})^{3}\delta^{2}n_{i}}{3L_{x}L_{y}}e^{-\xi^{2}(\bm{k}-\bm{k}^{\prime})^{2}}\sum_{z,z^{\prime}}\,\psi_{\bm{\kappa^{\prime}}}^{\dagger}(z)\psi_{\bm{\kappa}}^{\phantom{\dagger}}(z)\psi_{\bm{\kappa}}^{\dagger}(z^{\prime})\psi_{\bm{\kappa^{\prime}}}^{\phantom{\dagger}}(z^{\prime})M(z,z^{\prime},\xi,W),$

where $n_{i}=\sum_{\alpha}/L_{x}L_{y}W$ is the impurity concentration,
furthermore we used

$\left\llangle
U_{\alpha}U_{\beta}\dots\right\rrangle=\delta_{\alpha\beta}\frac{\delta^{2}}{3}\llangle\dots\rrangle,$
(S7)

and defined the function

$\displaystyle M(z,z^{\prime},\xi,W)$ $\displaystyle=\frac{1}{2\pi\xi^{2}}\int
dz_{i}\,\phi(z-z_{i})\phi(z^{\prime}-z_{i})$ (S8)
$\displaystyle=\frac{e^{-(z-z^{\prime})^{2}/4\xi^{2}}}{4\sqrt{\pi}\xi}\left[\mathrm{erf}\left(\frac{W+z+z^{\prime}}{2\xi}\right)+\mathrm{erf}\left(\frac{W-z-z^{\prime}}{2\xi}\right)\right],$
(S9)

plotted in Fig. 3; the error function is defined as
$\mathrm{erf}\left(x\right)=(2/\sqrt{\pi})\int_{0}^{x}e^{-t^{2}}\mathrm{d}t$.
Note that

$\int dz^{\prime}\,M(z,z^{\prime},\xi,W)$ (S10)

is a weak function of $z$, equal to $1$ in the middle of the slab, and going
down to $\approx 0.4$ at the edges in the range $\xi$.


* Papadakis et al. (2007) Papadakis, I. E., Villata, M., & Raiteri, C. M. 2007, A&A, 470, 857, doi: 10.1051/0004-6361:20077516
  * Peterson et al. (1998) Peterson, B. M., Wanders, I., Horne, K., et al. 1998, PASP, 110, 660, doi: 10.1086/316177
  * Peterson et al. (2004) Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004, ApJ, 613, 682, doi: 10.1086/423269
  * Piner & Edwards (2018) Piner, B. G., & Edwards, P. G. 2018, in Fourteenth Marcel Grossmann Meeting - MG14, ed. M. Bianchi, R. T. Jansen, & R. Ruffini, 3074–3079, doi: 10.1142/9789813226609_0389
  * Pollack et al. (2016) Pollack, M., Pauls, D., & Wiita, P. J. 2016, ApJ, 820, 12, doi: 10.3847/0004-637X/820/1/12
  * Rafle et al. (2012) Rafle, H., Webb, J. R., & Bhatta, G. 2012, Journal of the Southeastern Association for Research in Astronomy, 7, 33
  * Raiteri et al. (2010) Raiteri, C. M., Villata, M., Bruschini, L., et al. 2010, A&A, 524, A43, doi: 10.1051/0004-6361/201015191
  * Raiteri et al. (2013) Raiteri, C. M., Villata, M., D’Ammando, F., et al. 2013, MNRAS, 436, 1530, doi: 10.1093/mnras/stt1672
  * Raiteri et al. (2021) Raiteri, C. M., Villata, M., Carosati, D., et al. 2021, MNRAS, 501, 1100, doi: 10.1093/mnras/staa3561
  * Robertson et al. (2015) Robertson, D. R. S., Gallo, L. C., Zoghbi, A., & Fabian, A. C. 2015, MNRAS, 453, 3455, doi: 10.1093/mnras/stv1575
  * Romero et al. (1999) Romero, G. E., Cellone, S. A., & Combi, J. A. 1999, A&AS, 135, 477, doi: 10.1051/aas:1999184
  * Scargle (1982) Scargle, J. D. 1982, ApJ, 263, 835, doi: 10.1086/160554
  * Schwarz (1978) Schwarz, G. 1978, Annals of Statistics, 6, 461
  * Sergison et al. (2020) Sergison, D. J., Naylor, T., Littlefair, S. P., Bell, C. P. M., & Williams, C. D. H. 2020, MNRAS, 491, 5035, doi: 10.1093/mnras/stz3398
  * Shablovinskaya et al. (2023) Shablovinskaya, E., Malygin, E., & Oparin, D. 2023, MNRAS, 519, 3798, doi: 10.1093/mnras/stac3775
  * Sikora et al. (1994) Sikora, M., Begelman, M. C., & Rees, M. J. 1994, ApJ, 421, 153, doi: 10.1086/173633
  * Sikora et al. (2009) Sikora, M., Stawarz, Ł., Moderski, R., Nalewajko, K., & Madejski, G. M. 2009, ApJ, 704, 38, doi: 10.1088/0004-637X/704/1/38
  * Simonetti et al. (1985) Simonetti, J. H., Cordes, J. M., & Heeschen, D. S. 1985, ApJ, 296, 46, doi: 10.1086/163418
  * Singh & Meintjes (2020) Singh, K. K., & Meintjes, P. J. 2020, Astronomische Nachrichten, 341, 713, doi: 10.1002/asna.202013731
  * Steineke et al. (2020) Steineke, R., Waller, L., Reinhart, D., et al. 2020, The Astronomer’s Telegram, 13958, 1
  * Stetson (1987) Stetson, P. B. 1987, PASP, 99, 191, doi: 10.1086/131977
  * Stetson (1992) Stetson, P. B. 1992, in Astronomical Society of the Pacific Conference Series, Vol. 25, Astronomical Data Analysis Software and Systems I, ed. D. M. Worrall, C. Biemesderfer, & J. Barnes, 297
  * Sun et al. (2018) Sun, M., Grier, C. J., & Peterson, B. M. 2018, PyCCF: Python Cross Correlation Function for reverberation mapping studies.  http://ascl.net/1805.032
  * Tavecchio et al. (1998) Tavecchio, F., Maraschi, L., & Ghisellini, G. 1998, ApJ, 509, 608, doi: 10.1086/306526
  * Timmer & Koenig (1995) Timmer, J., & Koenig, M. 1995, A&A, 300, 707
  * Urry & Padovani (1995) Urry, C. M., & Padovani, P. 1995, PASP, 107, 803, doi: 10.1086/133630
  * Urry et al. (1997) Urry, C. M., Treves, A., Maraschi, L., et al. 1997, ApJ, 486, 799, doi: 10.1086/304536
  * Vaughan (2005) Vaughan, S. 2005, A&A, 431, 391, doi: 10.1051/0004-6361:20041453
  * Vaughan (2010) —. 2010, MNRAS, 402, 307, doi: 10.1111/j.1365-2966.2009.15868.x
  * Vermeulen et al. (1995) Vermeulen, R. C., Ogle, P. M., Tran, H. D., et al. 1995, ApJ, 452, L5, doi: 10.1086/309716
  * Villata et al. (2002) Villata, M., Raiteri, C. M., Kurtanidze, O. M., et al. 2002, A&A, 390, 407, doi: 10.1051/0004-6361:20020662
  * Villata et al. (2003) Villata, M., Raiteri, C. M., Kurtanidze, O. M., et al. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 299, High Energy Blazar Astronomy, ed. L. O. Takalo & E. Valtaoja, 221
  * Villata et al. (2004a) —. 2004a, A&A, 421, 103, doi: 10.1051/0004-6361:20035895
  * Villata et al. (2004b) Villata, M., Raiteri, C. M., Aller, H. D., et al. 2004b, A&A, 424, 497, doi: 10.1051/0004-6361:20040439
  * Villata et al. (2009) Villata, M., Raiteri, C. M., Larionov, V. M., et al. 2009, A&A, 501, 455, doi: 10.1051/0004-6361/200912065
  * Wagner & Witzel (1995) Wagner, S. J., & Witzel, A. 1995, ARA&A, 33, 163, doi: 10.1146/annurev.aa.33.090195.001115
  * Weaver et al. (2020) Weaver, Z. R., Williamson, K. E., Jorstad, S. G., et al. 2020, ApJ, 900, 137, doi: 10.3847/1538-4357/aba693
  * Webb et al. (2021) Webb, J. R., Arroyave, V., Laurence, D., et al. 2021, Galaxies, 9, 114, doi: 10.3390/galaxies9040114
  * White & Peterson (1994) White, R. J., & Peterson, B. M. 1994, PASP, 106, 879, doi: 10.1086/133456
  * Wierzcholska et al. (2015) Wierzcholska, A., Ostrowski, M., Stawarz, Ł., Wagner, S., & Hauser, M. 2015, A&A, 573, A69, doi: 10.1051/0004-6361/201423967
  * Xiong et al. (2020) Xiong, D., Bai, J., Fan, J., et al. 2020, ApJS, 247, 49, doi: 10.3847/1538-4365/ab789b
  * Xu et al. (2019) Xu, J., Hu, S., Webb, J. R., et al. 2019, ApJ, 884, 92, doi: 10.3847/1538-4357/ab3e50
  * Ye & Fan (2021) Ye, X.-H., & Fan, J.-H. 2021, PASJ, 73, 775, doi: 10.1093/pasj/psab039
  * Zhang et al. (2020) Zhang, L., Chen, S., Xiao, H., Cai, J., & Fan, J. 2020, ApJ, 897, 10, doi: 10.3847/1538-4357/ab9180
  * Zibecchi et al. (2017) Zibecchi, L., Andruchow, I., Cellone, S. A., et al. 2017, MNRAS, 467, 340, doi: 10.1093/mnras/stx054
  * Zola et al. (2021) Zola, S., Kouprianov, V., Reichart, D. E., Bhatta, G., & Caton, D. B. 2021, in Revista Mexicana de Astronomia y Astrofisica Conference Series, Vol. 53, Revista Mexicana de Astronomia y Astrofisica Conference Series, 206–214, doi: 10.22201/ia.14052059p.2021.53.40