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

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.

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 −2 B W , where l B = ħ h/eB ≈ 26 nm 1/B[T] is the magnetic length, and W the width of the slab. Taking an experimentally realistic value of k F = 0.01 Å −1 , 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][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 ∼ 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 nonequilibrium 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.

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.

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 eE · r e st with the position operator r and an adiabatic time dependence e st with s → 0 + , 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 where p is the equilibrium density matrix and g the non-equilibrium correction. Inserting into the von Neumann equation for the density matrix, and expanding to first order in E we obtain 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 |κ〉, where κ combines the quantum numbers. Off-diagonal and diagonal elements read, respectively, where the notation is A κκ = 〈κ|A|κ 〉, A κ = A κκ , E κ = 〈κ|H 0 |κ〉. The field-dependent term (8) expands in powers of V starting with the zeroth order, 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 Inserting into (7), taking the adiabatic limit s → 0 + , and applying disorder averaging ⟪. . .⟫ we obtain 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 where v κ = ∂ κ E κ is the velocity. Making the ansatz (11) simplifies to 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 Particle conservation however requires κ g κ = 0, which fixes the constant to c = 0. The current-density expectation value reads where V is the system volume. The conductivity tensor σ, defined as j = σE, becomes, using (13), 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.

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 κκ , 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 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, While it also vanishes upon averaging on the left-hand side (lhs) of (18), the new term on the rhs becomes In this sum there are terms that are proportional to |V κκ | 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, 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.

Weyl-semimetal slab model
We consider a minimal lattice model of a Weyl semimetal [42], where m κ = t(2 + cos β − cos k x − cos k y − cos k z ) 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 k z → (δ i, j+1 + δ i, j−1 )/2, where i = 0, 1, . . . , (W − 1) is the site number, corresponding to the discrete position in z, 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 diagonalization 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 Fermilevel 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 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.

Disorder potential
We model the disorder by static Gaussian potentials, 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(κ, 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 where ψ κ (z) = 〈z|ψ κ 〉 and n i is the impurity concentration. A detailed derivation can be found in the Appendix. The function M (z, z , ξ, W ) is given by where the error function is defined as erf (x) = (2/ π) x 0 e −t 2 dt. As illustrated in Fig. 3, M (z, z , ξ, W ) is mainly the (z − z ) 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 Despite the rather complex form, the impurity scattering is fully determined by two parameters -the real-space impurity width ξ and the overall impurity strength (set by δ 2 n i ), 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 | ξ −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

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 where the Fermi-surface average is given by 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 where the sum runs over the type j of states and the averaging is analogous to (26) but only over type i of states ( κ → i κ ).
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 N s and (ii) the number of states being dominated by surface states, N ∼ 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 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 respectively. Consequently, for W λ (which we always consider), we can estimate the scaling of the scattering rate between the different types of states as From this, (28), and (27) the width dependence of the scattering probabilities summarizes to 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.

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.

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

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, E =x E. The average occupation of the valley at k x = ±β is given by ± κ g κ /(N b /2), where the sum runs over all bulk states at the valley ±, and N b /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, where in the second line we used Eq. (13) and defined the average over valley bulk states 〈. . .
. 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 ξ 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 width above l e (ξ∆k) 2 , which becomes exponentially large for ξ∆k 1.

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.

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 where M i j = δ i j k q ik − q i j is given by the scattering rates between states, q i j . We assume scattering only between the nearest neighbors, with the rate q nn , and direct scattering between countermovers with the rate q d , in which case the matrix M becomes where we used that the direct-scattering part of q i j cancels when multiplied with Λ. For the nearest-neighbor part M nn of M there is a left pseudo-inverse, so that P M nn Λ = Λ (note i Λ i = 0 due to particle conservation). With its help, the full solution of (34) becomes We first ignore direct scattering. From the form of P, it is clear that the solution can depend 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 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 where the exponent p depends on the velocity arrangement. In the opposite limit we instead obtain 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/(q d + q 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Λ This explains the general power-law conductivity enhancement with ξ in Fig. 4, except for σ y y in Fig. 4(f), which we discuss separately.

Valley polarization
The averaged occupation of the two valleys ± can be expressed as The occupation difference is related to the valley polarization χ defined in Eq. (5.2), 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, where we used N b N c . 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 Together with (44) we obtain 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 N b /N s so that 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 ∝ W smaller than the total number of bulk states N b . 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 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 χ would enter in the form l ss (ξ∆) p /l ∝ 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 ξ, 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.

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 ∼ 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 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), 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 where n c = N c /V is the density of chiral bulk states. Adding all the current contributions, we obtain 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 σ 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 .

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) where N s ≈ 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 ∼ l sb . In Section 4.2 we found l s b /l ∼ W N s /N c -Eq. (31)-which leads to 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 ξ, 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 l ss (ξ∆k) p . For the considered parameters, l s b is much larger than this additional part since (ξ∆k) p W N s /N c , 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 N s /N c changes to N s /(N c + N b ) and thus becomes smaller and second, according to Eq. (31), l s b /l ∼ W N s /N b → 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 ∼ 1, the conductivity in the x direction is not significantly enhanced.

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

Derivation of scattering amplitudes
We consider Gaussian-type static impurity potentials, 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(κ, κ ) = 2πδ(E κ − E F ) × q(κ, κ ) between energy eigenstates |ψ κ 〉 and |ψ κ 〉, which we calculate using Fermi's Golden Rule, q(κ, κ ) = ⟪|〈ψ κ |V |ψ κ 〉| 2 ⟫ , where the disorder average is defined as We write the normalized wavefunctions as 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 〈ψ κ |V |ψ κ 〉 = 2πξ 2 L x L y e −ξ 2 (k−k ) 2 /2 Inserting into (S2) and using (S3), we obtain where n i = α /L x L y W is the impurity concentration, furthermore we used and defined the function plotted in Fig. 3; the error function is defined as erf (x) = (2/ π) x 0 e −t 2 dt. Note that dz M (z, z , ξ, 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 ξ.