Bulk-boundary quantum oscillations in inhomogeneous Weyl semimetals

Weyl fermions in an external magnetic field exhibit the chiral anomaly, a non-conservation of chiral fermions. In a Weyl semimetal, a spatially inhomogeneous Weyl node separation causes similar effect by creating an intrinsic pseudo-magnetic field with an opposite sign for nodes of opposite chirality. In the present work we study the interplay of external and intrinsic fields. In particular, we focus on quantum oscillations due to bulk-boundary trajectories. When caused by an external field, such oscillations are a proven experimental technique to detect Weyl semimetals. We show that the intrinsic field leaves hallmarks on such oscillations by decreasing the period of the oscillations in an analytically traceable manner. The oscillations can thus be used to test the effect of an intrinsic field and to extract its strength.

Weyl fermions in an external magnetic field exhibit the chiral anomaly, a non-conservation of chiral fermions. In a Weyl semimetal, a spatially inhomogeneous Weyl node separation causes similar effect by creating an intrinsic pseudo-magnetic field with an opposite sign for nodes of opposite chirality.
In the present work we study the interplay of external and intrinsic fields. In particular, we focus on quantum oscillations due to bulk-boundary trajectories. When caused by an external field, such oscillations are a proven experimental technique to detect Weyl semimetals. We show that the intrinsic field leaves hallmarks on such oscillations by decreasing the period of the oscillations in an analytically traceable manner. The oscillations can thus be used to test the effect of an intrinsic field and to extract its strength.
Band structure of a three-dimensional crystalline material can exhibit non-degenerate band crossings (nodes) in momentum space around which the Berry curvature behaves as a magnetic field emerging from a monopole. These sources or sinks of the "Berry flux" are known as the Weyl nodes, and must come in pairs of opposite sign in the Brillouin zone. Fermions in the states around such nodes can be assigned a quantum number, chirality, determined by the charge of the "Berry monopole". The total chirality in the Brillouin zone must be zero. The fermions around the nodes are relativistic, and their low energy dynamics is described by the Weyl equation. The materials are thus termed Weyl semimetals [1][2][3][4][5][6]. As such, they can realize an old concept from quantum field theory -chiral anomaly [7].
The chiral anomaly is the non-conservation of chiral charge in relativistic quantum field theory. The consequences of the existence of the chiral charge in the field theory are well known and celebrated. In particular, the chiral magnetic effect (CME) is a current response generated along the direction an externally applied magnetic field [7,8]. The same phenomenon is expected to occur in WSM. The observation is complicated by the Neilsen-Ninomiya theorem, which constraints the chiral charges to come in pairs of opposite chirality on the lattice. Therefore the CME cannot exist in a crystal in equilibrium due to a cancellation of contributions from nodes of opposite chiralities [9][10][11][12]. Nevertheless, one can obtain a CME in non-equilibrium conditions, either through dynamics, or through an imbalance of chiral chemical potentials. The latter is possible to generate by applying parallel electric and magnetic fields. The resulting conductivity enhancement has been observed in experiments measuring magnetoresistance [13][14][15][16][17][18][19].
Recent works have discussed novel ways for creating a CME in WSMs that can exist in equilibrium. The trigger for such a CME is an intrinsic pseudo-magnetic field and can be sustained locally while vanishing when integrating over the entire sample [20][21][22]. This phenomenon is linked to a known fact from classical electromagnetism -bound currents as a result of inhomogeneous magne-tization can flow within a medium so long as the total dissipation-less current vanishes when averaged over the volume of the sample. If the Weyl node separation is stemming from the magnetization, then the magnetization drives the CME.
In this work we address an outstanding challenge concerning the physics of WSM by proposing a direct way to detect and quantify the emergence of pseudo-magnetic fields and the resulting pseudo-CME. We achieve this by exploring the effects of inhomogeneities on a known and proven experimental scheme: quantum oscillations due to semiclassical trajectories traversing the bulk and surface of WSM [23][24][25]. This is a striking transport measurement carried out on WSM and exemplifying the existence of bulk-surface trajectories that result in a coherent periodic motion driven solely by the external magnetic field. In the present work we show that pseudo-magnetic fields can bend the bulk quasi-particle trajectories and hence have immediate and quantifiable effects on the period of such oscillations.
The emergence of pseudo-fields in Dirac materials has striking manifestations in graphene and WSMs. In WSMs it has been recently been claimed to play a key role both in the understanding of the physics of Fermi arcs, as well as in unraveling an equilibrium CME. In graphene, lattice deformations couple to the electronic degrees of freedom as gauge potentials that do not break time-reversal symmetry, but nevertheless result in the formation of Landau levels [26,27]. Preserving timereversal symmetry comes about through the coupling of the pseudo-fields with an opposite sign to the two valleys of graphene. The situation in WSM is similar, albeit richer. Pseudo-fields can emerge via lattice deformation as well as inhomogeneous magnetization, and the emergent pseudo-gauge fields couple to fermions of opposite chirality with an opposite sign [20,21,[28][29][30][31][32][33][34]. The resulting Landau level structure from spatial inhomogeneities is similar in many ways to the one resulting from external fields: it generates a lowest Landau level that disperses along the direction of the field near each node. But it is crucially different in one respect: the chirality of the low-arXiv:1802.00512v1 [cond-mat.mes-hall] 1 Feb 2018 est Landau level is indifferent to the sign of the node. The origin of the pseudo-CME in Weyl semimetals is rooted in the fact that while the usual CME is zero in equilibrium due to the cancellation of currents from nodes of opposite chirality, there is no such local cancellation for pseudo-fields. Therefore, for a time-reversal broken WSM with two nodes there will be a net local CME. Inversion-broken systems with a minimum of four nodes have counter-propagating currents stemming from time reversed pair of nodes.
The principle behind the effect of pseudo-fields on quantum oscillations is simple. Quantum oscillations stem from trajectories that traverse the bulk via the dispersion of the lowest Landau level, combined with a semiclassical sliding motion along the arcs at the surface perpendicular to the direction of the field [23,24]. When fixing the direction of the external magnetic field such trajectories are deformed due to intrinsic magnetic fields. Particles in the bulk are forced to move in the direction of the total effective magnetic field felt by the Weyl node, which is a superposition of the two components (external and intrinsic).
Below we derive the relevant formula for the period of the semiclassical oscillations. As we show, deformed trajectories have a strong quantifiable effect on the density of states as well as the frequency of oscillations in experimentally available responses (e.g. conductivity). The striking difference as compared to oscillations emerging purely from the external fields is that while oscillations in the absence of pseudo-fields are periodic and depend only on the total momentum space enclosed by the Fermi arcs, with bulk pseudo-fields the interval between oscillations becomes field-, pseudo-field-, as well as thicknessdependent. We support our predictions with numerical simulations performed using a tight-binding model compatible with the physics of Cd 3 As 2 , a Dirac semimetal on which the original quantum oscillations experiment was performed [25].
To make our discussion concrete, consider a film of a WSM. For simplicity, we consider a WSM with a single pair of Weyl nodes, but the analysis straightforwardly generalizes. We take the Weyl node separation p 0 to be along p y as depicted in Fig 1. Then the low energy Hamiltonian is Where v is the velocity, which we take here to be isotropic (the case of anisotropic velocity is discussed in the supplementary material), and µ 0 is the chemical potential offset with respect to the Dirac point. We note that written in this form, p 0 couples to the Hamiltonian as an axial vector potential. When p 0 is a constant, its importance is in the separation of nodes but beyond that it brings about no other interesting structure since ∇ × p 0 = 0. This is in accordance with the vector potential interpretation. Strain renders the Weyl node sep- aration space-dependent [20,21,31,35], and may make ∇ × p 0 non-zero. For simplicity we consider a strain profile that makes p 0 depend linearly on the z coordinate.
As we show below, such profile corresponds to a physical strain configuration. In such case we can write the Weyl nodes separation as p 0 (z) = (b 0 − B 5 z)ŷ Now, taking the curl of p 0 (z) we can define which is a pseudo-magnetic field that couples to Weyl nodes of opposite chirality with an opposite sign. Therefore, this position-dependent Weyl node separation and, as a result, bulk strain leads to intrinsic pseudo-magnetic fields. The magnetic field breaks the linear Dirac spectrum around each node into Landau levels, with a dispersion relation Here and later B(B5) = 1 √ eB(B5) , = 1. As the pseudomagnetic field acts oppositely in the two Weyl nodes, for the given node configuration the chirality of the lowest Landau level is the product of the sign of B 5 and corresponding Weyl node chirality. We now turn to the influence of B 5 on the bulk trajectories. Assuming an externally applied field in the z direction, B = Bẑ the total field experienced by a parti-cle with chirality s (with s = ±1) is The intrinsic field thus tilts the bulk trajectories from the direction of the external field by an angle θ = tan −1 B 5 /B. Hence, the bulk path traversed is of length In order to determine the period of quantum oscillations, we follow the analysis in Ref. [24], and derive the phase space quantization condition for the closed quasiparticle trajectories c p · dr = 2π(n + γ) where γ is a constant offset. According to the discussion above, the integral for the mixed bulk-surface trajectories is broken into two pieces, due to the presence of the intrinsic and external vector potentials, namely the integral is taken over four segments of the trajectories, including the two arcs, and two bulk branches linking the top and the bottom surfaces. See Fig. 1 for depiction of the trajectories in the mixed real-momentum and purely real spaces. For the arcs, the integral yields Where Φ is the total flux enclosed by the real space orbit in the surface plane. If the surface encircled is S, then with S k the momentum space area enclosed by the arcs. At small chemical potential this area is approximately given by S k = k 0 (µ + µ 0 )/v, where k 0 is the total length of the arcs, µ is the chemical potential measured from the Weyl nodes, µ 0 is the chemical potential offset as discussed in Ref. [24] and v the Fermi velocity at the surface which we take to be equal to that of the bulk. Note that in principle, k 0 may depend on B 5 : the presence of B 5 in the bulk necessarily means the length of the two arcs on opposite surfaces is inequivalent. Here, we will analyze the simplest case in which strain enhances the arc length on one surface by the same amount it shortens the arc on the opposite surface. Then, the total length of the surface trajectory is not modified by B 5 , although B 5 is finite in the bulk. This corresponds to the physical strain, corresponding to bending the Weyl semimetal field, discussed in [36]. In the supplementary material we discuss other cases where changes in the two arcs do not compensate one another.
In the bulk, the trajectory of the particles is parallel to the total magnetic field, so that [37] p · dr = L 1 + (B 5 /B) 2 (2µ/v).
FIG. 2. Density of states of the single spin block of particlehole symmetric Cd3As2 under both external magnetic field and stress-produced pseudo-magnetic field. Parameters of the simulations are: es = ep = 0.0574eV, m s⊥ = m p⊥ = 9.014eV nm 2 , m s = m p = 6.407eV nm 2 , A = 1.21287eV nm. The simulation is performed on a cubic lattice with lattice constant a = 8nm, the thickness of the material is 240nm, and the width of the stripe is 480nm. External magnetic field is B = 5T. Horizontal scale shows the effective pseudomagnetic field, and vertical is the energy as measured from the Weyl nodes. Striped lines are expressions from (10) without a free parameter Defining L ef f = 2L 1 + (B 5 /B) 2 and summing the two contributions together we have From equation (9) we can obtain our first testable prediction. The positions of the bulk-boundary energy levels, given by: 10) and are strongly affected by B 5 . Increasing B 5 makes the levels more dense. Furthermore, we can consider quantum oscillations as a function of B or B 5 . At B 5 = 0, the oscillations' period is ∆(1/B) ≈ 2πe/S k . As a small B 5 is introduced, a correction is added to the denominator, S k → S k + 2µeLB 2 5 /vB, making the oscillations nonperiodic, as the separation between peaks becomes magnetic field dependent. Moreover, as opposed to the case of a purely external magnetic field, the separation between peaks is now thickness dependent. In the opposite limit B 5 B, we obtain ∆(1/B) = 2πe/(S k +2µLB 5 /v), from which it is clear that that while oscillations are periodic in 1/B, B 5 decreases the period of oscillations, and makes it depend on the sample thickness.
Thus we obtain our main experimental predictions: closed bulk-boundary trajectories produce peaks in DOS at energies corresponding to the solutions of the Eq. (10). These can be observed in conductance (Shubnikovde Haas, SdH) and magnetization (de Haas-van Alphen, dHvA). Numerical tests. To confirm the validity of the results above and their applicability to realistic materials and conditions we performed numerical simulations of a Hamiltonian applicable to Cd 3 As 2 and Na 3 Bi Dirac semimetals. In these semimetals we can neglect the spinorbit coupling, thus we use the basis of a single spin, |s ↑, p ↑ : where: p ± = p y ± ip z , p = (p y , p z ), and the parameters used are summarized in Fig. 2. Note that we use the particlehole symmetric version of the model for simplicity (E p = −E s ). For this model we identify: distance between the Weyl nodes p 0 = ( es m s⊥ , 0, 0), and velocities around the Weyl points, v ⊥ = 2 √ e s m s⊥ , and v = A. For the purpose of our simulations we set v ⊥ = v by changing A. This makes comparison to (10) straightforward. We use the same procedure as in [36] to introduce the B 5 field according to the displacement vector: where α controls the strength of the strain. From this we compute the elements of the symmetric strain tensor u ij = (∂ i u j + ∂ j u i )/2. Then u 13 = 2αz, and u 11 = 2αx and correspondingly the pseudomagnetic field generated by the strain. In this model the u 31 has much smaller contribution to the pseudomagnetic field than u 11 due to a small prefactor 1/(ap 0 ) 2 , where a is the lattice constant of the material. In Cd 3 As 2 this prefactor is ≈ 1/57 [36]. We thus only use u 11 , which gives uniform pseudo-magnetic field in y direction of strength B 5 = 2α c ea cot ap 0 . Such strain corresponds to the hopping modification according to: Such modification makes the distance between the Weyl nodes, set by hopping in z direction, position-dependent, in accordance with the definition of B 5 we used above.
To introduce a real magnetic field we use the standard Peierls substitution which produces a real magnetic field in the z direction.
With both real and pseudo-magnetic field present only x direction remains infinite in the simulations. Thus, even though the obtained agreement with the theory seen in Fig. 2 and Fig. 3 is very good, we could not get rid of the finite-size effects completely.
In our numerical results we show DOS of a slab of the Cd 3 As 2 for a fixed B(B 5 ), while varying B 5 (B) correspondingly. This allows us to model the two experimental scenarios. We imagine putting a sample into fixed external field and continuously bending it to create the pseudo-magnetic field (see Fig. 3 for the change in DOS, corresponding to this scenario). Alternatively, one can fix the bend of the sample and change the external field (see Fig. 2 for similar results in this case). We show the result of the equation (10) without fitting parameters together with the numerically computed DOS. There is visible disagreement for small B 5 regime seen in Fig. 2, as the traverse of the Fermi arc is the relatively large part of the trajectory. The linear dependence on the chemical potential is a simplistic approximation for the motion along Fermi arc, thus causing discrepancy. The good agreement otherwise shows reliability of our model for predicting the influence of the external and pseudofields. Thus our prediction enable extraction of the values of B 5 as a function of strain applied to material by applying external magnetic field and measuring SdH or dHvA quantum oscillations.
We stress that results presented here apply both to time-reversal-and inversion-broken Weyl semimetals, since one can think of the latter as two time-reversed copies of the former. While locally the pseudo-CME might add up to a zero net contribution in time reversal symmetric systems due to the cancellation between time reversed pairs of nodes, the trajectories are still modified by them, and the effect on quantum oscillations should still be present.
The case of Dirac semimetals is more subtle: it is known that the strain can develop spin-orbit coupling gapping out the Dirac semimetals like Cd 3 As 2 , as the symmetry protecting the cones is broken [38,39]. Nevertheless, we predict that small B 5 is still accessible in the experiment in the limit of high magnetic field or high chemical potential with respect to Dirac point compared to spin-orbit gap. In the first case the two Weyl cones corresponding to the same Dirac cone have opposite spins, and are shifted in energy and momentum due to Zeeman term. Thus, we predict two sets of quantum oscillations corresponding to the two spin sectors to be present. In the second case the gap near the Dirac points does not influence the physics at high chemical potential and our predictions remain intact.

Changes in the total arc length
In this section we consider the generic case in which the total length of the Fermi arcs is modified by the existence of a bulk B 5 . This can occur, for example, in strained samples where strain is applied on one surface and gradually relaxes to zero away from that surface such that the opposite side of the sample maintains the original unstrained value.
Getting back to equation (9), we note that both L ef f and S k depend on B, B 5 , and L. Writing out S k explicitly 2π(n + γ) =µ(2L /v + k 0 (B 5 )/evB) where k 0 (B 5 ) = (k 0 (0, B 5 ) + k 0 (L, B 5 ))/2, k 0 (z, B 5 ) is the Weyl node separation as a function of position in z direction and B 5 . This expression transforms to eq. (9) when k 0 (0) = k 0 (L) = k 0 . Let us now estimate the change in the overall arc length as follows. We assume that B 5 is uniform in the bulk of the sample, i.e the Weyl node separation changes linearly from one surface to another in the z direction. Hence the total arc length is given by where z 0 denotes the position in which the nodal separation is unperturbed. We set the sample position between z = 0 and z = L, so that k 0 (0, B 5 ) = k 0 − B 5 z 0 and k 0 (L, B 5 ) = k 0 + B 5 (L − z 0 ) and the change is the total length of the arcs can be estimated as B 5 (L − 2z 0 ) and is linear in B 5 and L. Consequently, while we expect oscillations, they will not be at constant intervals as a function of B or B 5 . We can consider the two limits of B 5 /B 1 and B 5 /B 1. For the latter, L ≈ L and we get that ∆(1/B) = 2πev/k 0 (B 5 )(µ 0 + µ) (19) the expression is of a similar form to the one appearing in [24], but now k 0 changes linearly in B 5 . For B 5 B we have that L ≈ LB 5 /B so that 2π (n+γ) = µ(2LB 5 /vB+k 0 (B 5 )/evB)+k 0 (B 5 )µ 0 /evB (20) and therefore ∆(1/B) = 2πev [2eµLB 5 + k 0 (B 5 )(µ + µ 0 )] −1 (21) And the denominator again changes linearly in both B 5 , L.

Unisotropic Fermi velocity
In equation (1) the Fermi velocity was taken to have a constant and isotropic value, v. In practice, the Fermi velocity might have a different value depending on the direction, hence we now extend the calculation to the case where v = (v ⊥ , v ⊥ , v z ). Since we choose the direction of the external field to be in the z direction and perpendicular to the surface, we take the surface velocity to be v ⊥ as well. The modification introduced by the anisotropy