Inhomogeneous Weyl and Dirac semimetals: Transport in axial magnetic fields and Fermi arc surface states from pseudo Landau levels

Topological Dirac and Weyl semimetals have an energy spectrum that hosts Weyl nodes appearing in pairs of opposite chirality. Topological stability is ensured when the nodes are separated in momentum space and unique spectral and transport properties follow. In this work we study the effect of a space dependent Weyl node separation, which we interpret as an emergent background axial vector potential, on the electromagnetic response and the energy spectrum of Weyl and Dirac semimetals. This situation can arise in the solid state either from inhomogeneous strain or non-uniform magnetization and can also be engineered in cold-atomic systems. Using a semiclassical approach we show that the resulting axial magnetic field $\mathbf{B}_{5}$ is observable through an enhancement of the conductivity as $\sigma\sim \mathbf{B}_{5} ^{2}$ due to an underlying chiral pseudo magnetic effect. We then use two lattice models to analyze the effect of $\mathbf{B}_5$ on the spectral properties of topological semimetals. We describe the emergent pseudo-Landau level structure for different spatial profiles of $\mathbf{B}_5$, revealing that (i) the celebrated surface states of Weyl semimetals, the Fermi arcs, can be reinterpreted as $n=0$ pseudo-Landau levels resulting from a $\mathbf{B}_5$ confined to the surface (ii) as a consequence of position-momentum locking a bulk $\mathbf{B}_5$ creates pseudo-Landau levels interpolating in real space between Fermi arcs at opposite surfaces and (iii) there are equilibrium bound currents proportional to $\mathbf{B}_{5}$ that average to zero over the sample, which are the analogs of bound currents in magnetic materials. We conclude by discussing how our findings can be probed experimentally.


I. INTRODUCTION
Electronic and lattice degrees of freedom are inevitably intertwined in solid state physics [1]. With the advent of graphene [2] a remarkable effect was soon acknowledged: elastic deformations of the lattice originating from strain can couple to the low-energy Dirac quasiparticles of graphene as a pseudo or axial magnetic vector potential [3][4][5]. In order to preserve time-reversal symmetry, the axial vector potential couples to the two valleys with an opposite sign. Spatially inhomogeneous strain generates an effective axial-magnetic field B 5 , and gives rise to a pseudo-Landau level spectrum at low energies [6]. Strain-induced pseudo-Landau levels have been observed with scanning tunneling microscopy both in real [7] and artificial graphene [8,9]. The realization of Landau levels without real external magnetic fields results in effective fields as high as 300T, and is a direct consequence of the Dirac nature of the carriers, exemplifying the unique response of these class of systems to strain.
In this respect, Weyl and Dirac semimetals in three dimensions (3D) [10,11] are expected to behave as 3D cousins of graphene and host similar effects. A Weyl semimetal is a state with pairs of band touching points with linear dispersion, also called Weyl nodes. The Weyl nodes are sources and sinks of Berry curvature, i.e., Berry curvature monopoles, in momentum space. The charge of a monopole defines the chirality of the corresponding node, and the two partners of a pair of nodes must have opposite chirality. Since monopoles can only annihilate in pairs, Weyl nodes are topologically protected as long as they are separated in momentum space by a vector b.
Dirac semimetals can be then regarded as special cases of Weyl semimetals, where nodes of opposite chiralities are located at the same momentum but additional symmetries constrain the system to remain gapless.
The vector b can be alternatively interpreted as an axial gauge field since it couples with an opposite sign to Weyl nodes of opposite chirality [12][13][14][15][16][17]. This interpretation suggests that if b is varied in space it will generate a nonzero axial magnetic field B 5 = ∇ × b that couples to fermions of opposite chirality with an opposite sign. This observation, together with the recent theoretical evidence that effective gauge fields can emerge in strained Weyl semimetals [18][19][20][21] are two important motivations of our study.
The main purpose of this work is therefore to address the effect of an axial magnetic field B 5 arising in any type of Weyl and Dirac semimetal. We will first discuss that the physical mechanism for the emergence of B 5 depends on the presence or absence of time reversal symmetry. When a Weyl semimetal preserves time reversal symmetry, each pair of nodes has a time reversed partner (hence the minimal number of nodes is four). In this case, B 5 can be generated by an inhomogeneous strain profile and has an opposite sign for each time reversed pair of nodes such that time reversal symmetry is preserved. In contrast, time reversal breaking Weyl semimetals are not subjected to this constraint and B 5 can emerge either from an inhomogeneous magnetization or strain [13,18]. Then, a richer set of phenomena can occur, with the example of a finite angular momentum for the electronic states [22].
Here we naturally unify these mechanisms under a sin-gle framework by considering any Weyl semimetal in the presence of a spatially varying Weyl node separation. We discuss the effect of inhomogeneous nodal separation on (i) the changes in the spectrum due to the emergence of pseudo-Landau levels and (ii) transport, and in particular, the effect on the conductivity. In exploring (i), we will show that the Fermi arcs, the characteristic surface states of topological semimetals [23,24], can be reinterpreted as the chiral n = 0 pseudo-Landau level associated with a large B 5 confined to the surface. On the other hand, the presence of a finite B 5 in the bulk also creates pseudo-Landau levels. The n = 0 pseudo-Landau level of the bulk is then smoothly connected to the surface arcs. We find that this follows from a known phenomenon in the two dimensional quantum Hall effect, sometimes referred to as position-momentum locking; the average center position of a Landau level wavefunction in one planar direction, say y , is determined by its momentum in the perpendicular direction, k x . Such a concept, considered recently in some detail for real magnetic fields in Weyl semimetals [25,26] will also be of use here. We will find that these results significantly depart from the naive expectation gained from studies of strained graphene [3][4][5] from which the bulk 0-th pseudo-Landau levels is expected to have opposite chirality with respect to the Fermi arcs [13,27]. Regarding (ii), we will show, based on both field theoretic arguments and a rigorous semiclassical calculation, that a finite bulk B 5 enhances the conductivity of a Weyl semimetal as σ ∼ B 2 5 . This response is related to a charge anomaly much like the topological negative magnetoresistance σ ∼ B 2 is rooted in the chiral anomaly [28].
The advances made on the experimental front in recent years allows us to explore the feasibility and detectability of these effects while considering a wide range of possible platforms. While Dirac and Weyl semimetals material growth thrives in condensed matter , they could also be engineered with cold atoms [69]. This provides access to both inversion breaking semimetals that are currently the standard in the solid state (with notable potential exceptions [67,68]), as well as to time reversal breaking semimetals in cold atoms. We end this paper with an estimate of the magnitude of these effects based on real material parameters and discuss how the phenomena described here might be detected.
The structure of this paper is as follows. In section II we review how inhomogeneous strain (magnetization) in time-reversal invariant (broken) Weyl and Dirac semimetals leads to an effective axial magnetic field, or alternatively, a space dependent Weyl node separation. In section III we predict the enhancement of the conductivity of Weyl and Dirac semimetals due to a uni-form axial magnetic field, and relate it to the underlying associated chiral pseudo-magnetic effect. The latter leads to bound currents flowing within the material and along the boundary, even in equilibrium, as well as an enhanced bulk longitudinal conductivity. In section IV we introduce two lattice models that illustrate the spectral breakdown into pseudo-Landau levels, and explicitly calculate the bound current distribution within a finite sample. We describe position momentum locking and show that any equilibrium bulk current is compensated by surface currents. As an example of the relevance of our findings to different contexts, we discuss the implications of our results on the situation where two Weyl semimetals with different Weyl node separation are brought into proximity. Finally, in section V we provide a discussion and proposals to experimentally observe the discussed phenomena.

II. EMERGENCE OF SPACE DEPENDENT NODE SEPARATION
In this section we review how a non-uniform magnetization or strain can lead to a space dependent node separation. This discussion of the physical origin of axial magnetic fields in solid states systems will set the stage for our general study presented in the next sections, by introducing the main concepts and notation.

A. Time-reversal symmetric Weyl semimetal
The existence of Weyl nodes relies on breaking either time-reversal (T ) or inversion symmetry (I). If both are present, all energy bands are manifestly twofold degenerate, and, in particular, point nodes must have a degenerate partner with opposite chirality. If T is broken, the two nodes of a pair are allowed to be separated in momentum space, with the separation determined by b. If I is broken but T is preserved, each pair must have a time-reversed partner, such that the minimum number of Weyl nodes in I-broken Weyl semimetals is four. In an I-broken Weyl semimetal, the Hamiltonian of the low-energy Weyl fermions for each pair is given by where v i are the Fermi velocities, σ i are Pauli matrices spanning the orbital space of the Weyl fermions, and s = ± denotes the chirality of the nodes. It is clear from Eq. (1) that the vector b = (b x , b y , b z ) quantifies the separation of the two nodes of a pair in momentum space.
In general, in a crystalline material a uniform strain of strength g alters the overlaps between wavefunctions in different orbitals [18,21]. In Ref. 21, for instance, which considered HgTe with spin-orbit induced band inversion in the presence of strain, found that at low energies the spectrum consists of four pairs of Weyl nodes. Due to the presence of strain, the Hamiltonian of each pair then takes the form From (2), we observe the effect of strain is twofold. Firstly it modifies the Fermi velocities through v i (g) ∼ √ g. Secondly, it controls the Weyl node separation through the vector b with magnitude |b(g)| ∼ √ g.
Both of these effects are analogous to those predicted in graphene [3,5,70].
Ref. [21] analyzed the effect of strain as an homogeneous perturbation characterized by g. However in practice, strain can have a non uniform profile g(r). This can occur, for instance, when there is a lattice mismatch between a substrate and the semimetal sample or alternatively via chemically induced strain. In the former scenario, for thin enough samples, it can be expected that the strain profile relaxes smoothly as a function of the distance from the substrate (quantitative estimates will be provided in section V). Let us then choose the direction of the strain inhomogeneity to be y such that g(r) = g(y), and assume that g(y) changes slowly with y, on a length scale that is much smaller than any other microscopic length scale. In such a perturbative regime, g(y) g 0 + δg(y) with δg(y)/g 0 1, it is reasonable to expand Eq. (2) to lowest order in (3) We first note that the effect of the space dependent Fermi velocity due to strain [70] is of order δg(y) g0 ∂ i , so it can be neglected in the low-energy linear approximation. An inhomogeneous velocity has been predicted recently to result in an anomalous Hall signal transverse to the direction in which the strain is applied [71,72]. However, these works neglected the effect of emergent inhomogeneous Weyl node separation that, according to the above analysis, it is of lower order in δg(x) and therefore dominant.
Motivated by the above effective model, we consider the following real space Hamiltonian where b(r) ≡ b (0) − δb(r) encodes the strain profile. From such a model with a single pair of nodes, the Dirac semimetal and the inversion breaking Weyl semimetal can be constructed by adding time-reversal partners of Eq. (4) in the appropriate location in momentum space in order to restore time reversal symmetry. We note that in the literature, the vector b has been referred to, depending on the context, as chiral, axial or pseudo vector potential. In this work, we will refer to b as the axial vector potential as custom in high-energy physics where this field was first discussed.

B. Time-reversal breaking Weyl semimetal
In T -breaking Weyl semimetals, the axial vector potential b that separates the Weyl nodes in momentum space can have its physical origin in a finite magnetization. A hallmark example of a Weyl semimetal of this kind is the proposal by Burkov and Balents [73], which is based on a topological insulator-trivial magnetic insulator heterostructures. Another example is the Weyl semimetal state induced by magnetic order at a quadratic band crossing point on pyrochlore iridates [74]. In addition, a time-reversal breaking Weyl semimetal has been argued to be consistent with ARPES measurements on YbMnBi 2 [67]. Other promising recent proposals [68] suggest, relying on ab-initio calculations, that magnetic Heusler compounds can host Weyl quasiparticles. Remarkably, alloys of the latter kind could realize the minimal model of a T -breaking Weyl semimetal with only two nodes, and a node separation that spans a significant fraction of the Brillouin zone.
Similar to other magnetic materials, the magnetization that induces the Weyl semimetal state by separating two nodes of opposite chirality in momentum space, may itself be non-uniform and become a function of position [75]. As a result, the Weyl node separation becomes a function of position, generating an axial magnetic field strength. To lowest order in the spatial variation of the Weyl node separation (i.e., magnetization), the effective model (4) will describe this state. We note as well that interfaces between magnetic domains, either abrupt or continuous could be modeled by the spatial profiles of b that we consider in this work. In addition, it is plausible to expect that coating one of the surfaces with a ferromagnetic layer or enough surface magnetic dopants will also result in a spatially varying magnetization. The successful realization of magnetically doped topological insulators [76] suggests that similar techniques could be applied in Weyl samples as well.
In addition to an inhomogeneous magnetization, a nonuniform strain can generate an axial magnetic field in Tbreaking Weyl semimetals. Once T is broken and a finite b exists, strain can allow to effectively change b from its bare value, thus changing the separation of Weyl nodes in momentum space. Note, however, that strain on its own cannot break time-reversal symmetry and thus cannot induce a finite b from b = 0.

C. Effective realization in cold atomic systems
While the above sections discuss realizations of inhomogeneous Dirac and Weyl semimetals in the solid state, our discussion is also relevant for cold atomic sys-tems [69]. Two appealing aspects of the latter is that different topological properties can be extracted easily from experimental data [77][78][79] and models that respect or break time reversal could be engineered.
There are two further advantages of considering cold atomic systems that are relevant to probe the effects we discuss. Firstly, engineering a space dependent Weyl node separation is experimentally feasible. A natural way to achieve this is to smoothen the confining trapping potential that contains the lattice [80]. Although plausible, it does not offer the degree of control that will expose all of the salient features we are after, but only a subset of them. A more controlled and realistic approach is to create domain walls between topological states [81] as already described in Ref. [82] for Chern insulators. We elaborate on more particulars of this proposal in the concluding section.
The second advantage of cold-atomic systems is the controlled experimental access to transport properties that was recently discussed in the context of Chern insulators [83] [84]. The key observable is the center-of-mass velocity v c.m. since it is simple to measure experimentally and is connected to the current density j and particle density n through [83] v c.m. = j n .
Although in general n also depends on external fields, in our particular case it acts as an unimportant constant factor as we describe below. Thus, via Eq. (5) it is possible to directly probe the longitudinal conductivity.

III. ENHANCED TOPOLOGICAL LONGITUDINAL CONDUCTANCE
Many peculiar transport properties have been attributed to topological semimetals [10,11], while not many of them have proven easy to measure. One particularly striking feature that has been confirmed in several experiments 85] is the appearance of a non-saturating negative magnetoresistance [86,87]. Conductance is one of the most accessible experimental observables to characterize the properties of a material. Therefore, we aim here to explore the effect of an inhomogeneous Weyl node separation on this quantity. We show that a space dependent axial vector potential b enhances the longitudinal conductance of a Weyl semimetal through its corresponding axial magnetic field B 5 = ∇ × b. We begin by introducing a coupling to an external electromagnetic field and discussing the role of gauge invariance. We follow this by extending a quantum field theory argument previously used in Ref. [40] to include pseudo gauge fields. We then support our claims with a more rigorous semiclassical calculation by solving the semiclassical kinetic equation.
A. Coupling to an external electromagnetic field and the role of gauge invariance In order to calculate the conductivity we couple the Hamiltonian (4) to an external electromagnetic field A µ = (A 0 , A) using the minimal substitution The first important property of this Hamiltonian is that left handed (s = 1) and right handed (s = −1) Weyl fermions are decoupled, experiencing a chirality dependent gauge field, a s µ = A µ + sb µ , where b µ = (0, b(r)). For what follows it is useful to define the left and right handed field strengths f µν s = ∂ µ a s ν − ∂ ν a s µ . The continuity equation for the current for each species j µ s = (n s , j s ) is determined by [13,88] ∂ µ j µ s = s which we now discuss in some depth. The right hand side of Eq.
µ , and f µν s , respectively. The second term on the right hand side is a scattering term that acts to equilibrate the number density n s of the left and right chiralities [89,90] with a typical scattering time τ [91].
The conservation laws for the current j µ = s j µ s and axial current j µ 5 = s sj µ s can be obtained from the addition or subtraction of the two equations composing Eq. (7). A first inspection reveals a seemingly striking feature; neither of the two currents is conserved. The non-conservation of j µ 5 is not forbidden, and it is referred to as the chiral anomaly [11,88,92]. If both E 5 = 0 and B 5 = 0 then one recovers the celebrated chiral anomaly that schematically reads ∂ µ j µ 5 ∼ E · B [28,88]. In this case the vector current satisfies ∂ µ j µ = 0 and thus charge is conserved. However, it seems that charge is not conserved when either E 5 = 0 or B 5 = 0. Indeed, from Eq. (7) we find that the vector current satisfies [13,17].
In fact, there is no contradiction with current conservation, which can be understood in two different ways. Firstly, the axial vector potential b is an observable and thus must be single valued, with a zero vacuum expectation value [11]. Thus the total flux of the resulting B 5 must be zero over a surface enclosing the entire sample; regions with B 5 and −B 5 compensate each other by generating an equal number of left and right handed fermions separated in real space. Then, upon applying an electric field charge flows from one region to the other, respecting global charge conservation.
Second, at the quantum field theory level, the nonconservation of charge is fixed by defining a consistent current that preserves gauge invariance. The procedure has been described both in the context of high energy physics [93] and Weyl semimetals [33]: the current calculated above (the covariant current) is complemented by Chern Simons currents (the Bardeen-Zumino polynomials) of the form − 1 4π 2 µνρσ b ν F µν , which exactly cancel the anomaly and define the consistent current. This procedure effectively imposes a boundary condition for the spectral flow at the energy of the cut-off that bounds the field theory. In a lattice system, such a cut-off is a natural quantity, and the spectral flow is bounded by construction. Note also that this procedure applies wherever there is a finite axial magnetic field B 5 and therefore includes the boundaries where b jumps from zero to a finite value or viceversa.

B. Quantum field theory approach
Our aim is to use the chiral anomaly Eq. (7) to find the longitudinal conductance in the presence of a finite axial magnetic field B 5 . To promote the derivation in Ref. [40] we first use that for three-dimensional Weyl fermions n s = µ 3 s /(6 π 2 3 v 3 ), which we employ to rewrite the steady state form of Eq. (7) as We now recall that the component of the current parallel to B and B 5 is given by The first term is the chiral magnetic effect [29,30] and has been thoroughly studied in the context of Weyl semimetals [14][15][16][31][32][33]. This term must be zero in equilibrium [36,37,[94][95][96] since it is proportional to a chemical potential imbalance µ 5 = 1 2 s sµ s . The second term on the left hand side is the key to our results: it represents the analog of the chiral anomaly in the presence of an axial magnetic field B 5 and can be finite in equilibrium. In the context of dense relativistic matter it was shown that it is possible to generate an axial current j 5 from a vector field B with a conductivity proportional to the chemical potential µ = 1 2 s µ s [97]. It follows that there is a contribution to the vector current j generated by an axial field with the same coefficient [98]. We refer to this term as the chiral pseudo-magnetic effect.
In order to obtain the longitudinal conductivity we combine Eqs. (9) and (8), Assuming that and for the case where B = Bẑ and B 5 = B 5ẑ , B s = (B + sB 5 )ẑ and E 5 = 0, we obtain to linear order in E. The first term is a transport (or free) current. For the case when B 5 = 0 it reproduces the chiral anomaly enhanced magneto-conductivity [86]. The central result here is that for B 5 = 0 there is a contribution to the bulk longitudinal conductivity coming entirely from the spatial dependence of b(r): inhomogeneous strain or magnetization enhances conductance in a Weyl or a Dirac semimetal. We now provide a physical interpretation of the second term as a magnetization (or bound) current. Recall that b(r) is analogous to a finite magnetization M(r). Thus local bound currents given by j b ∝ ∇ × M(r) ∝ B 5 are allowed, consistent with the form of the second term in Eq. (13). Note that in the related context of the quantum Hall effect, a finite bound current proportional to µ and the curl of the magnetization is expected on general grounds [99]. In that case, the bound current has its origin on the edge states, while here they can be associated with the bulk.
We remark that the second term in Eq. (13) highlights an important difference between semimetals that respect time reversal symmetry and those that do not. In the former the sum over time reversed pairs of nodes will cancel out the µB 5 term. In Section III C we provide further arguments that support the interpretation of this term as a bound current, while in Section IV we corroborate these findings by numerically studying the bound current profile within specific lattice models.

C. Boltzmann equation approach
In this section we outline a different and more rigorous derivation of Eq. (13) that relies on solving the Boltzmann equation. Within this approach, the semiclassical equations of motion are extended to include an anomalous velocity term that arises due to the existence of a non-zero Berry curvature [100]. Typically, for a Weyl semimetal, the equations of motion are written for a single flavor of chiral fermions accounting for the physics in momentum space centered around a particular Weyl node [86]. Both chiralities feel the same electric and magnetic fields. The key difference in the present analysis is that the effective external fields are now chirality dependent due to the axial vector potential; both chiralities are still decoupled but feel different effective fields.
More precisely, the effect of the effective magnetic fields B s = (B + sB 5 ) on the left (s = +1) and right (s = −1) chiralities can be incorporated by promoting the semiclassical equations of motion [100] tȯ Here E s p = ε s p − m s p · B s and Ω s p are, respectively, the dispersion relation, which includes a correction due to the magnetic orbital moment m s p , and the corresponding Berry curvature for each chirality s. For each chirality the unperturbed dispersion relation of the upper band is ε s p = sv|p|. The Berry curvature and magnetic orbital moment in this case take the simple form Ω s p = s 1 2|p| 2p and m s p = −ev|p|Ω s p . We emphasize that, consistent with the discussion in Sec. II, we neglect effects of the higher order corrections due to the inhomogeneous Fermi velocity which can enter through a space dependent Berry curvature [71,72,101].
The distribution function for each chirality f s p satisfies a semiclassical kinetic equation where I coll {f } is the collision integral. Using Eqs. (14) and (15) we can writė where D p,s = (1 + e c B s · Ω s p ) and v s p = ∂E s p /∂p is the perturbed velocity. We find that the contribution of a single chirality to the current density is given by Within the relaxation time approximation [102] and f 0 p,s (E k,s ) is the equilibrium distribution function to be evaluated at the modified dispersion relation E k,s . We are interested in a stationary and homogeneous solution to the kinetic equation, which takes the forṁ Expanding the left hand side of the previous equation to lowest order in the fields and rearranging terms, we obtain (20) Inserting the first term in Eq. (20) into Eq. (18) and ignoring the contribution transverse to electric field leads to To leading order in the magnetic field, we can make the replacement E s p → ε s p in both the equilibrium distribution function and the definition of the velocity. The integral results in with µ the chemical potential, consistent with earlier findings [19,32]. Summing Eq. (22) over chiralities results in the second term in Eq. (13).
Remarkably, we have also performed a rate of entropy production calculation in the spirit of Ref. [86] and find this term to be absent. The reason for this absence is that the rate of entropy production is determined exclusively by the free current density rather than the total current density. This result implies, as anticipated above, that Eq. (22) should be physically interpreted as a bound current proportional to the curl of the magnetization in the system and as such does not contribute to transport.
Next, we insert the second term in Eq. (20) into Eq. (18) to calculate the correction to the current corresponding to the deviation δf s p = f s p − f 0 p,s up to first order in E and second order in B s , which reads Using that for the lower band Ω s p = −s 1 2|p| 2p and m s p = −ev|p|Ω s p and after tedious but straightforward manipulations we obtain where ν(ε) = ε 2 /2π 2 3 v 3 is the unperturbed density of states of a Weyl node at energy ε. The first term encodes the conductivity due to a finite Fermi surface [103,104]. The second term is novel to this work and upon summing chiralities leads to a longitudinal contribution to the conductivity that reads for B = 0 and represents the main result of this section.
Up to the numerical coefficient, this contribution has exactly the same form as the first term of Eq. (13) confirming that an inhomogeneous strain or magnetization contributes to increase the conductivity. The numerical differences between Eqs. (13) and (25) can be traced back to the expansion of the modified density of states D p,s appearing in the denominator of Eq. (23) and the inclusion of the magnetic moment in the semiclassical calculation. We note as well that, as with conventional magnetorresistance, we expect that the term Eq. (24) is also supplemented by a Fermi surface contribution due to the Lorentz force when the Fermi surface is anisotropic [105].
We conclude this section by computing the center-ofmass velocity in this approach, an observable quantity in cold atomic experiments. From Eq. (5), to obtain the center-of-mass velocity and relate it to the conductivity we must compute the electron density n that to leading order reads The relevant center-of-mass velocity component is obtained by inserting Eqs. (24) and (26) into Eq. (5) to obtain which is valid in the absence of magnetic field (B = 0). Therefore, through Eq. (27) a cold atomic experiment can probe the anomalous longitudinal conductivity σ ∼ B 2 5 that we predict by monitoring the center-of-mass motion of an atomic cloud.

IV. SPECTRAL PROPERTIES OF INHOMOGENEOUS SEMIMETALS
Based on our field-theoretic and semiclassical discussion regarding the axial-gauge field coupling in Weyl semimetals, we now study the effects of magnetization or strain inhomogeneities within a microscopic lattice model realization. This allows us to numerically corroborate the arguments presented in the previous section. We will consider two canonical lattice models whose low-energy long-wavelength limit is given by Eq. (4). For simplicity, we choose lattice realizations of Weyl semimetals with the minimum number of low-energy Weyl fermions, i.e., one Weyl fermion of each chirality, which implies that timereversal symmetry is broken in both models. Our results, however, generalize to time-reversal symmetric (but inversion symmetry-breaking) models, and occasionally we will explicitly comment on such generalizations.
The first minimal lattice model we use to describe Weyl fermions coupled to axial-gauge fields, is motivated by the solid state realizations of Weyl and Dirac semimetals. It can be obtained starting from a lattice-regularized version of a topological insulator k · p Hamiltonian [106], such as Bi 2 Se 3 or related materials. A simple generalization of this model has been used to describe the Dirac semimetals Cd 2 As 3 and Na 3 Bi [34,81]. The model has four bands, originating from an orbital degree of freedom A, B and spin ↑, ↓, and the corresponding electron operators are defined as c r = (c rA↑ , c rA↓ , c rB↑ , c rB↓ ) T . The Hamiltonian H 4b of this four band model is the sum of two terms and is given by The first term is the topological insulator Hamiltonian with Γ-matrices Γ j = (σ z s y , σ z s x , σ y s 0 , σ x s 0 ), Γ b = σ y s z , and with the components of the D vector given by (29) Here we have set the kinetic energy scale t = 1. The matrices σ x,y,z and s x,y,z are Pauli matrices acting on orbital and spin degrees of freedom, respectively (σ 0 and s 0 are identity matrices). The Dirac mass M , which describes a hybridization of the orbitals into bonding and anti-bonding states, controls whether the material is on a trivial-or topological-insulator phase (which may be strong or weak). Here, we set M ≡ 3 + m, such that m = 0 corresponds to a Dirac semimetal state with a 3D Dirac point at k = 0, and m > 1 (m < 0) corresponds to a trivial (topological) insulator.
Whereas the first term in Eq. (28) respects both timereversal symmetry (T ) and inversion symmetry (I), a nonzero b = (b x , b y , b z ) breaks T . To see how it is responsible for generating the Weyl semimetal phase, one may expand (28) to linear order in δk around k = 0 and obtain This result shows that b is responsible for the separation of the two nodes in momentum space, and thus couples as an axial gauge field. Indeed, Eq. (30) should be compared to the Weyl Hamiltonian of Eq. (4), where the node separation was identified with the axial gauge field. In Eq. (30), the matrix Γ b takes the role of giving the axial gauge field b a different sign at the two nodes. It should be noted, however, that in Eq. (30) the node separation depends on both b and m. In particular, the node separation is proportional to |b| 2 − m 2 and vanishes (i.e., the system is insulating) for m > |b| [14,81]. In the following, we will study Hamiltonian (28) numerically in the presence of a spatially non-uniform b(r) which gives rise to axial gauge fields B 5 ∇ × b(r) = 0. It is important to note that ∇ × b(r) can only have the meaning of an axial magnetic field coupled to Weyl fermions in regions where m < |b(r)|. We always take m ≥ 0, since we are not interested in the regime where D describes a topological insulator. A full phase diagram of the model defined by Eq. (28) was described in Ref. 81.
The second lattice model we use to verify our results consists of effectively spinless electrons, has only two  bands, and can be regarded as one of the two time reversal partners that compose the model in Ref. [69]. It is therefore motivated by two practical considerations: (i) it is plausible to implement it in the cold atomic context, and (ii) it falls into the class of models tailored for the methods presented in Ref. [82] to generate domain walls in topological systems. In this case, the electron operators are given by c r = (c rA , c rB ), where A, B represent some generalized orbital degree of freedom, and the Hamiltonian takes the form Here, σ z = ±1 again represents the orbital degree of freedom, and the components of d are given by This model breaks time-reversal symmetry; it has a pair of linearly dispersing Weyl cones at k = {0, 0, ± cos −1 (M/t − 2)} for 1 < |M/t| < 3, two pairs of Weyl cones for |M/t| < 1 and is a gapped insulator when |M/t| > 3. The parameter M/t controls the distance between Weyl nodes; interpolating between M/J = 2 and M/t > 3 simulates the boundary between a Weyl semimetal with two nodes and an insulator. Therefore, we promote M → M (y) which sets the Weyl node separation as b(y) by b(y) = (0, 0, 2 arccos(−2 + M (y)/t)).
We note that, although this model separates the Weyl nodes in k z , any other separation direction can be chosen by redefining d appropriately.
From our calculations we find that both models qualitatively exhibit the same behaviour, so we will focus mainly on results obtained for the four-band model of Eq. (28).

A. Lattice pseudo-Landau Level structure of B5
We begin by specifying the spatial profiles of b(r) that generate the axial magnetic configurations we will study. In what follows, we will always take b(r) = b x (y)x, such that it only depends on y and corresponds to an axial magnetic field along z. Since nonzero b x implies a separation of Weyl nodes along the x axis in momentum space, b x (y) describes Weyl nodes whose separation ∆k x depends on the y coordinate. Note that this effectively corresponds to a Landau gauge, and consequently,  (k x , k z ) remain good quantum numbers. Furthermore, it follows that the axial magnetic field ∇ × b is orthogonal to the direction of node separation. The extension of the system in the y direction is L y (similarly for x and z), which we measure in units of the lattice constant a, and we assume b x (1) = b x (L y ).
We first consider the series of profiles b x (y) shown schematically in the bottom row of Fig. 1. The form of b x (y) is such that it increases linearly from b x 1 to b x 2 between y = 1 and y = 1 , then stays flat at b x 2 , and subsequently decreases linearly again from b x 2 back to b x 1 between y = 2 and y = 2 . As a result, between 1 and 1 we have a constant and negative ∇×b B 5 , whereas between 2 and 2 we have a positive constant ∇ × b B 5 , with a magnitude that depends on the ratio ∆b x /∆ 1,2 , with ∆b x = b x 2 − b x 1 and ∆ 1,2 = 1,2 − 1,2 . For fixed ∆b x , the strength of the effective axial magnetic field increases for decreasing ∆ 1,2 . By taking ∆ 1,2 → 0, as shown in the right most configuration of Fig. 1, we create two sharp interfaces between regions of constant b = b x 1x and b = b x 2x . We choose m [see Eqs. (28) and (30)] such that the latter corresponds to a boundary between a Weyl semimetal and an insulator, i.e., m > b x 1 [81]. Thus, the configurations of Fig. 1 interpolate between axial magnetic fields in the 3D bulk of the system, and a sharp 2D boundary between a Weyl semimetal and an insulator. The axial magnetic field is gradually confined to a 2D surface while increasing its strength. The energy spectra corresponding to these axial magnetic field profiles are shown in the top panels of Fig. 1, where we have set b x 1 = 0.0, b x 2 = 1.0, and m = 0.5. Before turning to an analysis of the energy spectra obtained from the lattice model, it is useful to recall the Landau level spectrum associated with axial magnetic fields in the continuum, cf. Eq. (4). We will refer to these Landau levels as pseudo-Landau levels. For Weyl fermions coupled to a uniform axial magnetic field B 5 along the z direction, B z 5 , the n ≥ 1 pseudo-Landau levels have energies E n± (k z ) = ± v k 2 z + 2|B z 5 |n and thus disperse in k z . In addition, there is an n = 0 or zeroth pseudo-Landau level with energy E 0 (k z ) = sgn(B z 5 ) vk z . Importantly, all energy levels are doubly degenerate: one for each of the two Weyl nodes. In particular, the two Weyl nodes, which have opposite chiralities, have the same E 0 (k z ) dispersion of the zeroth Landau level. In contrast, in the presence of a (vector) magnetic field B z , the chirality of the n = 0 branch, i.e., the upward or downward slope as a function of k z , depends on the chirality of the Weyl nodes [28].
The pseudo-Landau level structure of the continuum is reflected in the low-energy part of the lattice energy spectra shown in Fig. 1. In particular, in the left most panel we observe both a flat branch of zero energy states and flat branches of states at higher energies. These can be identified with pseudo-Landau levels since these do not disperse in k x . Moreover, labeling the flat branches of states at nonzero energies by an index n, we find that the energies indeed scale as ∼ √ n. In the right most panel, where the axial magnetic field is confined to the sharp boundary between the Weyl semimetal and the insulator, we still find a branch of zero energy states. The flat branches at higher energies are absent. The zero energy states are simply the Fermi arc surface states, which must exist at such a surface boundary due to the topology of the Weyl nodes and they connect the projections of the bulk Weyl nodes. Figure 1 thus suggests that as 1,2 − 1,2 is gradually taken to zero, confining the axial magnetic field to a narrower region along y while increasing its magnitude, the n = 0 pseudo-Landau level becomes the Fermi arc. The energy of the higher n ≥ 0 pseudo-Landau levels scales as ∼ √ B 5 , and are therefore "pushed" out of the spectrum as B 5 increases, as can be observed from left to right in Fig. 1 .
The key implication of Fig. 1 is that the Fermi arc surface states of a Weyl semimetal can be thought of as an n = 0 pseudo-Landau level corresponding to an axial magnetic field spatially confined to the surface boundary, i.e., ∼ B z 5 δ(y − 1,2 ). The chirality of the Fermi arc states, i.e., their dispersion in k z (discussed in more detail below), corresponds to the chirality of the n = 0 pseudo-Landau level modes and depends on sgn(B z 5 ).

B. Fermi arcs as n = 0 pseudo-Landau Level
To study this correspondence in more detail, we now consider a different set of b x (y) profiles, which are shown in the middle row of Fig. 2. As is schematically demonstrated, in this set of profiles b x (y) increases and decreases linearly from −b x to b x , and we take b x = 0.5. Furthermore, m is set to zero, such that the profiles interpolate between bulk axial magnetic fields on the left, and an interface between two Weyl semimetals with inverted Weyl node separation b x on the right. The corresponding energy spectra are shown in the top panels of Fig. 2.
Focusing first on the right panel, i.e., the sharp interface between two Weyl semimetals with inverted b x , we find that the branch of zero energy states, which connect the projections of the bulk nodes, is fourfold degenerate. These are the Fermi arcs localized at the boundaries between the Weyl semimetals, two for each boundary. This is consistent with the number of arcs mandated by topol-ogy. Similar to Fig. 1, the spectrum in the left panel exhibits the pseudo-Landau level structure at low energies. The flat branch of zero energy states corresponding to the n = 0 pseudo-Landau level is fourfold degenerate, in agreement with the degeneracy of the Fermi arcs (2 × 2 = 4). The pseudo-Landau level degeneracy can be understood by recalling that each Weyl node contributes one Landau level, which is doubled due to the two spatially separated regions of ±B 5 . (The counting in Fig.  1 is more subtle, which we explain below.) From Fig. 2 we again observe that the Fermi arcs are adiabatically connected to n = 0 pseudo-Landau levels as the profile of b x (y) is varied. The n ≥ 1 pseudo-Landau levels are pushed to higher energies (and are eventually absent from the spectrum) due to the increasing axial magnetic field strength.
More insight can be gained by studying the wavefunctions of the Fermi arc states and comparing them to Landau orbital wave functions. Recall that wavefunctions of the lowest n = 0 Landau level orbitals are given by where l b is the magnetic length. The Landau orbitals are centered around y = k x l 2 b , which locks the y coordinate to momentum k x . The distance between the Landau orbitals along y, as well as the spread of the wavefunction, are determined by the magnetic length l b . Since l b ∼ 1/ √ B the magnetic length decreases as the field strength B increases.
In Fig. 2 we plot the support of the wavefunctions of the zero energy states corresponding to different values of k x . The states for which wave functions are shown are indicated by solid black dots in the upper panels. From the bottom left panel of Fig. 2, we see that at each momentum there are indeed four states, with a wave function support similar to that of Landau orbitals. At k x = 0, these are localized at the center of each of the two regions of positive and negative axial magnetic fields, B z 5 = ±B 5 . As k x increases, the support of the wave function is shifted along y, and in opposite directions for Landau orbitals associated with Weyl nodes of opposite chiralities. This result follows from the positionmomentum locking expressed in Eq. (34).
As the profile of b x (y) is changed from left to right, the axial magnetic field strength is increased and thus the effective magnetic length is decreased. This is clearly reflected in the wavefunction support of the zero energy states: the spread becomes narrower and they move closer together. Eventually, when the axial magnetic field is confined to the interfacial boundary between the Weyl semimetals at 1 and 2 , the Landau orbitals are localized at the boundaries, and should be viewed as the wave functions of the Fermi arc states.
An analogous picture arises when we plot the wavefunction support of zero energy states of the spectra in Fig. 1, which is presented in Fig. 3. Figure 3 . 3. Plot of the support of the wavefunction, as a function of y, of selected zero energy states corresponding to the spectra shown in Fig. 1. Panels (a), (b), and (c) correspond to left, middle, and right most panels of Fig. 1, respectively (graphically indicated by the square, triangle and hexagon). Here, in each panel, the seven different curves correspond to different values of kx, which are indicated by solid black dots in Fig. 1 and red arrows in (a). Two pseudo-Landau orbitals correspond to each value of kx, as is most clearly seen in (a). This matches the degeneracy of n = 0 pseudo-Landau level and is consistent with the number of Fermi arcs: one Fermi arc per boundary. Note that as compared to the bottom row of Fig. 2, the Landau orbitals centered in the region where b x (y) < m are absent, which explains the difference in degeneracy (twofold vs. fourfold as discussed in the main text).
position-momentum locking y ∝ k x . Note that there are only half the number of Landau orbitals as compared to Fig. 2, reflecting the different degeneracy of zero energy states: twofold versus fourfold. As is shown in Fig.  3(a), the Landau orbitals in regions where b x (y) < m are absent. The region where b x (y) > m divides into a part where B z 5 is positive and a part where B z 5 is negative. In each region, both Weyl nodes contribute n = 0 pseudo-Landau level orbitals, but of opposite momentum. Figure 3(a) shows the contribution of only one Weyl node (chirality). Similar to Fig. 2, we observe that in Fig. 3, as the magnetic length decreases towards zero from (a) to (c), the pseudo-Landau orbitals become the wave functions of the Fermi arc states localized at the surface boundary between the insulator and Weyl semimetal. From the perspective of Fermi arc surface states, the degeneracy of states is can be inferred from topology: the number of Fermi arcs at each boundary is equal to the Chern number change across the boundary. The Chern number can be defined when Weyl nodes of opposite (or different) Berry monopole charge are separated in momentum space. This suggests that, per the adiabatic connection between Fermi arcs and n = 0 pseudo-Landau levels, the degeneracy of pseudo-Landau levels in the solid state systems has a topological origin.
The pseudo-Landau levels of the continuum are flat as a function of k x , but disperse in k z . Therefore, it is useful to compare the spectra of the lattice model as a function of k z . The energy spectra, obtained for the two b x (y) profiles of the left most and right most panel of Fig. 1, are shown in Fig. 4. The low-energy part of the spectrum in the left panel of Fig. 4 clearly exhibits the pseudo-Landau level structure, with n ≥ 1 pseudo-Landau level energies dispersing in k z as ∼ k 2 z + nB 5 . In addition, there is an n = 0 chiral (E ∝ +k z ) and anti-chiral (E ∝ −k z ) mode, which are localized in different spatial regions corresponding to opposite B z 5 . This is confirmed by the wave function support of states with different k z in the bottom panel. Note that the momentum k z is not locked to the y coordinate. The right panel of Fig. 4 shows the dispersion of the Fermi arc states in k z . The wave functions (bottom panel) are localized at the surface boundaries. Once more, one may think of the Fermi arcs as the n = 0 chiral (and anti-chiral) modes of the n = 0 pseudo-Landau level. The difference in the spread of the wave functions between the bottom left and bottom right panels originates from the difference in effective magnetic length.
At this stage, two remarks regarding the generality of our results are in order. The first concerns the profiles of b x (y) shown in Figs. 1 and 2. In all these cases the increase or decrease of b x (y) is chosen to be linear, giving rise to a constant B z 5 . In solid state materials, axial vector potentials are due to strain or magnetization inhomogeneities, and one may expect the axial magnetic field strength to have a more general dependence on position, i.e., B 5 = B 5 (r). To verify the application of our results to the more general case, we calculate the energy spectrum of Hamiltonian (28) with b x (y) as shown in Fig. 5(d). The corresponding spectrum is shown in Fig.  5(c), which demonstrates two important features. First, the n = 0 pseudo-Landau level remain non-dispersive in k x and at zero energy (for k z = 0). Second, the higher pseudo-Landau levels appear to have acquired a dispersion and shifted in energy. This is in agreement with the energies of the pseudo-Landau levels, E 0 ∝ sgn(B z 5 )k z and E n± ∝ k 2 z + 2|B z 5 |n, respectively. Replacing B z 5 by B z 5 (y) and noting that y ∝ sgn(B z 5 )k x due to Eq. (34), it follows that the n ≥ 1 pseudo-Landau levels should disperse in k x , whereas the energy of the n = 0 pseudo-Landau level should not change. This shows that our results hold for a general axial magnetic field configuration.
The second remark concerns the shape of the Fermi arcs. In Figs. 1, 2, and 4, Fermi arcs states at zero energy  4. Energy spectra, as a function of kz (with kx = 0), and wavefunction support obtained from Eq. (28) in the presence of a bulk axial magnetic field B z 5 (left) and for an interface boundary between an insulator and a Weyl semimetal (right). The spatial profile of b x (y) corresponding to the left and right panels here are given in the most left and most right bottom panels of Fig. 1, respectively. On the left, the low-energy branch of the spectrum has the structure of pseudo-Landau levels. The linearly dispersing chiral (black) and anti-chiral (red) modes are the n = 0 pseudo-Landau levels. The wavefunction support of n = 0 modes is shown for states with kz values indicated by solid black dots. On the right, the linearly dispersing modes correspond to the Fermi arc surface states and are sharply localized at the boundary interfaces 1 and connect the bulk Weyl nodes in a straight line located on the k x axis. This is due to a non-fundamental symmetry of the model (28). In general, the shape of the Fermi arcs is not restricted and the condition E(k x , k z ) = ε node , where ε node is the energy of the bulk nodes (assuming they are both at the same energy), may define a curve connecting the bulk nodes with arbitrary shape. To study the more general case, we add −δt i cos k i aσ 0 s 0 to Hamiltonian (28) and calculate the energy spectra for two different b x (y) profiles. As shown in Figs. 5(a) and (b), both the pseudo-Landau levels (a) and the Fermi arcs (b) disperse in k x , implying that the Fermi arcs trace out a curve of the form shown in the inset of (b). We conclude from Figs. 5(a) and (b) that the correspondence between Fermi arcs and n = 0 pseudo-Landau levels holds for more general Fermi arc shapes.

C. Calculation of the bound current density
As pointed out in section III, the axial vector potential b in time-reversal broken Weyl semimetals physically corresponds to a magnetization vector, and is therefore expected to give rise to bound currents j b . Generally, the total current in an electronic system with a finite magne- tization can be expressed as a sum of bound j b and free current j f densities [107]. The former has the property that it must average to zero over the system's volume, but it is allowed to be non-zero locally. It therefore can be expressed as the curl of a local vector M(r), the magnetization, such that j b = ∇ × M. In the case where M is constant in the bulk, the bound currents only exist as surface currents, j surf b , that are perpendicular to the surface normaln such that j surf b =n × M. We now provide numerical evidence that justifies the interpretation of the second term in Eq. (13), i.e., µB 5 , as a bound current density, by calculating the current density from our microscopic lattice models. To this end, we consider five different linear profiles of the axial vector potential of the form b = b x a (y)x, parametrized by a = 1, 2, 3, 4, 5. These trace out a finite Weyl semimetal with a Weyl node separation that increases linearly in the y direction [see Fig. 6 (b)]. The current density j z (y) is computed through the expression j z (y) = Ĵ z (y) = 1 L y L z n,kx,kz u n kxkz J z (y) u n kxkz f (ε n kxkz ) where J z (y) = ∂H kxkz (y)/∂k z is the current operator, u n kxkz are the single particle eigenstates and f (ε n kxkz ) is the Fermi-Dirac distribution function evaluated at the n−th eigenvalue ε n kxkz . Its real space distribution is shown in Fig. 6 (a) for different values of b x a , with an offset for clarity. From Fig. 6 (a), we observe two main features. Firstly, for the flat profile b x 1 , the current is localized at the boundaries with equal weight but opposite sign. This is entirely consistent with what is expected of a bound current; for a constant magnetization the bound currents are localized at the interface and are normal to it (j surf b explained above). Secondly, as the slope increases (profiles b x a =1 ), we observe that the weight associated with one boundary is transferred to the bulk, but the total current density remains zero overall.
In order to understand both of these features in more detail, we express each of these linear profiles mathematically as which is schematically shown in Fig. 7 (b) by the light blue curve. The corresponding B 5 = ∇ × b profile is set to a constant in the bulk (which is zero for the flat profile outlined by b x 1 and non-zero for b x a =1 ), and to two Dirac delta functions of opposite sign corresponding to each boundary. Setting 1 results in the spectral profile already discussed in Fig. 1 (lower row, right most panel): a Weyl semimetal with two Fermi arc surface states. As chiral surface states, the two Fermi arcs carry two compensating current densities, consistent with what is observed for the b x 1 case shown in Fig. 6 (a). As the difference b x f − b x i increases the profile traces higher values of a in b x a =1 and an increasing and constant B 5 = ∇ × b emerges in the bulk. The relevant energy spectrum for this situation is shown in Fig 7 (c); pseudo-Landau level emerge, and the 0th chiral pseudo-Landau level seems indistinguishable from a Fermi arc.
In fact, it is possible to distinguish the surface and bulk contribution due to position-momentum locking. As was discussed following Eq. (34), the average guiding center position y is tied to the momentum k x such that y ∼ k x . Thus, by tracking the dependence of the wavefunction of the zeroth Landau Level Ψ n=0 (y) as a function of k x we can extract real-space information.
In Fig 7 (d) we show |Ψ n=0 (y)| 2 for two representative values of k x marked in Fig 7 (c). At k x = 0 the wavefunction only has a finite weight at the edges, outlining a purely surface state. However, when b x f < |k x | < b x i the wave function acquires weight in the bulk associated with the 0th pseudo-Landau level emerging in this momentum space region. Importantly both regions are continuously connected as a function of k x and thus, from position-momentum locking, they are also connected in real space. The bulk B 5 creates bulk pseudo-Landau levels that connect to the surface arcs as a consequence   of position-momentum locking. This is a central result of this work and it allows us to understand the different instances of Fig. 6 (a): as B 5 increases, a 0th bulk pseudo-Landau level forms that connects the unbalanced surface states. The current density is sensitive to this fact; the current density lost by one Fermi arc is compensated and transferred to a bulk 0th pseudo-Landau level. In other words the bound current at one surface is compensated by the sum of bound currents in the bulk and the remaining surface.
Taking one step back, our results can be summarized by the identification of the second term in Eq. (13) with a bound current due to the curl of a magnetization, j b = µB 5 = ∇ × M. We have further corroborated this by explicitly checking that for small fields the magnitude of the bound current in the bulk of Fig. 6 (a) increases linearly with B 5 (cf. Fig. 6 (c)) and it is only finite if µ = 0.
All these arguments combined establish the three main points collected in the abstract: (i) the Fermi arcs can be reinterpreted as n = 0 pseudo-Landau levels resulting from a B 5 confined to the surface (ii) a bulk B 5 creates bulk-pseudo Landau levels that connect to the surface arcs as a consequence of position-momentum locking and (iii) there are bound currents proportional to B 5 and the chemical potential that average to zero over the sample, as occurs in magnetic materials.
It is important to stress that the effects we discuss are sensitive to the pseudo-magnetic field direction, in particular whether it is parallel or perpendicular to the Weyl node separation. In the parallel case, the bulk and surface behave as expected; the two bulk 0-th pseudo-Landau Levels have the same chirality, which is opposite to the chirality of the Fermi arcs at the boundaries. The surface spectral weight is evenly distributed between each arc, and compensates the bulk [27]. In contrast, in the perpendicular case considered in this work, the surface states are allowed to have a chirality that matches that of the bulk 0-th pseudo-Landau levels, a fact that is enforced by position-momentum locking. The Fermi arcs at different surfaces do not have the same spectral weight, rendering the rich spectral structure reported here.

D.
Interface between two Weyl semimetals as a finite bulk B5 Given the above discussion, it is interesting to consider the case of two Weyl semimetals of different Weyl node separation that share a boundary. At this point, we can intuitively predict the outcome: if the system abruptly changes the Weyl node separation, this is equivalent to the appearance of a B 5 that is confined to the interface. This case occurs when the change from b x i to b x f is localized at a single point in space, as shown in Fig. 7 (b) dark curve. It represents two Weyl semimetals with two constant axial vectors potentials b x i and b x f brought into contact. From the previous section, we expect that cur-rents bound to the surfaces of the two semimetals have a different magnitude, and their difference is proportional to ∇ × (b i − b f ). Alternatively, the extensions of their Fermi arcs in momentum space are different. Therefore, we expect that there is a finite amount of current bound to the interface, that compensates for that difference.
From a spectral point of view, it is interesting to observe how the Fermi arcs behave in this case. To this end, we consider the profile of b x (y) shown by the dark curve in Fig. 7 (b). Mathematically it is described by The corresponding spectrum is shown in Fig. 7 (a) and is characterized by the appearance of"mini-arcs", fragments of longer arcs that are surface states belonging to the uniform system with the larger Weyl node separation b x i . This situation is quite generic for an interface between two topological phases with the same topological invariant -the surface states hybridize and gap out along the region of contact. For Weyl semimetals we can consider this from the perspective of assigning a Chern number to two-dimensional slices of momentum space that lie between two Weyl nodes [23]. The corresponding chiral edge states hybridize and gap out in the region of overlap in momentum space, which is equal to the length of the shorter axial vector potential b x f . It is now straightforward to understand how the spectrum evolves when the interface between two Weyl semimetals is smeared across the entire sample, which corresponds to Fig. 7 (b) light blue curve. The spectrum in this case, presented in Fig. 7 (c), shows that the smeared interface transforms the interfacial arc states at l s into zeroth pseudo-Landau Levels. Moreover, this conclusion based on the spectral information of Fig. 7 (c) is corroborated by the wavefunction spread in real space shown in Fig. 7 (d). As discussed in the last section, the latter portrays how, as k x is changed, the real space probability density shifts from edge to bulk. We stress that this is another complementary instance of the discussion around Figs. 1 and 2. The surface Fermi arcs blend with the 0th-Landau levels within a single b x (y) profile, Fig. 7 (b) (light blue curve) rather than a sequence of them as in Figs. 1 and 2.
Moreover, the comparison between the light blue and dark blue profiles in Fig. 7 (b) inspires the following physically appealing picture of how the current density profile is distributed for the linear profile studied in Fig. 6 (b). The region between i and f in Fig. 7 (b) (light blue curve) can be approximated by a collection of infinitesimal discontinuities in the spirit of the trapezoidal method for curve integration. At each discontinuity there are Fermi arcs that meet and annihilate with a part of the arcs of the neighboring layer. These bulk arcs carry current density and thus act as a collection of bulk sheets that "leak" current density into the bulk from the surface.

V. DISCUSSION AND CONCLUSION
In this work we studied inhomogeneous Weyl and Dirac semimetals with a space dependent Weyl node separation. We have discussed how such a scenario can arise either from inhomogeneous strain or magnetization in existing solid state systems as well as cold-atomic setups. Underlying our results is an axial magnetic field B 5 that couples to the electronic degrees of freedom with two main experimental consequences for both Dirac and Weyl semimetals, which we have addressed in depth.
The first experimental consequence is a drastic change of the spectral properties of topological semimetals. Using two different lattice models we have established two novel spectral features attributed to the emergence of B 5 : (i) Fermi arcs are secretly 0th pseudo-Landau levels due to a finite and large B 5 at the boundary and (ii) bulk pseudo-Landau levels form due to B 5 and compensate the difference in density of states of inequivalent Fermi arcs at opposite boundaries via position-momentum locking. Point (i) is an inevitable consequence of the boundary supporting a finite B 5 since the Weyl node separation must vanish in vacuum. This allows us to identify the existence of Fermi arcs with the emergence of a boundary 0th pseudo-Landau level due to B 5 , even in the absence of bulk inhomogeneities. Such a correspondence provides a novel perspective on surface physics of Weyl and Dirac semimetals. Point (ii) is particularly relevant for spectral probes like ARPES or STM that are sensitive to a modification of the electronic spectrum [7][8][9]. Remarkably, for the particular case of time-reversal broken realizations of inhomogeneous Weyl semimetals, the axial magnetic field results in an inhomogeneous distribution of bound currents throughout the sample, which exist in equilibrium and average to zero over the entire volume. The appearance and field dependence of these bound currents in the lattice realization that we study is also consistent with our semiclassical treatment of the response.
We envision two plausible routes to probe the effects resulting from the discussed bound currents. Magnetic sensors such as scanning superconducting quantum interference devices (SQUIDs) are a natural way to probe local distributions of bound currents. Such probes were proven very successful at detecting magnetization modulations and inhomogeneous current distributions at interfaces such as LAO/STO heterostructures, and in two-and three-dimensional topological insulators (cf. Refs. [108][109][110]). It is important to emphasize that non-local currents due to an inhomogeneous magnetization and the bound currents due to B 5 are physically analogous phenomena, and thus it is very likely for the latter to be detectable. The difference between bound currents in an ordinary magnetized material and bound currents within a Weyl semimetal is that those in the latter emerge from a unique coupling between the Weyl fermions and the background magnetization as described throughout this paper.
Another promising alternative are torque experiments,  8. A schematic point contact transport set-up to measure the current J (orange arrow) to determine the strain dependent conductivity σ ∼ B 2 5 predicted in this work. Strain is induced by the substrate with a small lattice mismatch that relaxes along the y direction (red dashed arrow). If the Weyl node separation is b x (blue arrow) then B5 ẑ (green arrow). Contacts V1 (V2) probe a surface with (without) Fermi arcs, resulting in an anisotropic contribution to the conductivity.
which have in fact already been conducted with Weyl semimetals [111]. The magnetic torque τ = M × B is a direct measure of the magnetic anisotropy of a crystal, and thus can reveal noncompensation of surface magnetic domains.
Within our semiclassical treatment we have predicted a second experimentally relevant consequence, which is related to transport: inhomogeneities enhance the longitudinal conductivity in Dirac or Weyl semimetals as σ ∼ B 2 5 . We find that this transport property is routed in a chiral pseudo magnetic effect similar to the enhancement of the magneto-conductance due to the existence of the chiral magnetic effect. We stress that this prediction is applicable to all kinds of topological semimetals. Since this enhancement is insensitive to the sign of B 5 , it will be generated by each pair of Weyl nodes, a number that is of order ten in realistic materials.
A typical transport setup with point contact probes, schematically shown in Fig. 8, suffices. It represents a thin film, or nanowire of a topological semimetal that is strained by a substrate with a small lattice mismatch. This setup can be particularly relevant for strained HgTe/CdTe heterostructures discussed in Ref. [21]. As illustrated in Fig. 8, point contact measurements have the additional freedom to choose the voltage probes to be on a surface with or without Fermi arcs (V 1,2 respectively), depending on wether the momentum space projection of the nodes on that particular surface coincides or not. Although for the TaAs class of materials all accessible surfaces host Fermi arcs, Dirac materials like Cd 3 As 2 do present surfaces free of arcs. Such anisotropy highlights the difference between topological metals and more conventional strained semiconductors. While in the latter, conductance can also be enhanced via strain by modifying the band structure, changing the position of the voltage probes is not expected to result in anisotropic measurements if current jetting is negligible [45].
The above considerations suggests a way to probe the surface bound currents. As we have discussed extensively in section IV, surface currents on opposite surfaces in the presence of bulk strain will not compensate each other completely since their difference will be carried by the bulk pseudo-Landau levels. It could be expected that two distinct surfaces that are perpendicular to B 5 will carry different currents. Put differently, the Fermi arcs of two opposite surfaces will not be of equal length in the presence of bulk strain. Thus, measuring an anisotropy in the surface currents could reveal the size of the bulk strain gradients. However, typical surface current signals are small, and thus this detection mechanism is less feasible in practice. Nonetheless we note that the special nature of the Fermi arcs can conspire to make the surface contribution sizable, as in other related situations [112].
In order to assess the significance of the effects we predict we now estimate the magnitude of a bulk axial-field B 5 . For a thin sample, it is expected that strain relaxes linearly with height over several unit cells if the lattice constant mismatch is small. For a typical sample, such as a Cd 3 As 2 thin film, we can consider a height of the order of ∆L ∼ 10 nm, and a conservative value of strain is translated into a change in lattice constant of ∼ 1% but can be as large as ∼ 10%. If we assume that the Weyl node separation spans typically a tenth of the Brillouin zone |b| ∼ 1/10Å −1 , then the effective magnetic field is |B 5 | ∼ ∆b e∆L ∼ e 10 15 m −2 ∼ 4 T. Remarkably, this conservative estimate results in sizable magnetic fields, certainly above the detectable threshold of magnetic loops and SQUIDS. These fields will also induce detectable changes in conductance that can be probed by growing samples with different ∆L that result in different intrinsic B 5 .
Cold atomic systems offer a controlled alternative. For these systems, we have proposed that the conductivity enhancement can be measured by monitoring the centerof-mass velocity v c.m. . This quantity is the quotient between the current density j and the particle density n, which will depend on the external fields in general. For small fields the density reduces to a trivial constant and the center-of-mass position is determined solely by the current density times a constant factor. This experiment could be performed in a cold atomic realization of the two band model used in this work. Moreover, this model, along with the family of two band models that host Weyl fermions [69], are ideal to apply the method proposed in Ref. [82] to study inhomogeneities. It relies on implementing a space dependent offset between neighboring sites. Applying such an offset profile, which is feasible experimentally, results in a space dependent onsite potential (M (y) in our two band model) leading to a finite axial magnetic field, as described in this work. These two considerations render cold atomic systems as a natural platform to engineer and observe the effects presented.
To conclude, we have shown that inhomogeneous strain and magnetization have profound and observable implications on the electronic spectrum and transport properties of Dirac and Weyl semimetals. Moreover, our general analysis provides an alternative angle to explore distinct topological to trivial and topological to topological interfaces, both smooth and abrupt [81]. The fundamental correspondence between Fermi arcs and pseudo-Landau levels and their connection to bound currents departs from naive expectations and can inspire future theoretical and experimental studies of surface effects in Weyl and Dirac semimetals.
Note added: While preparing this manuscript, we became aware of the recent related work by Pikulin et al. [27] discussing the effect of strain generated by torsion in inversion breaking Weyl semimetal nanowires. Their findings are consistent with and complementary to our results. In particular torsion in wires generates an axial magnetic field B 5 that is parallel to the Weyl node separation rather than the perpendicular component described here.