Magnetic breakdown spectrum of a Kramers-Weyl semimetal

We calculate the Landau levels of a Kramers-Weyl semimetal thin slab in a perpendicular magnetic field $B$. The coupling of Fermi arcs on opposite surfaces broadens the Landau levels with a band width that oscillates periodically in $1/B$. We interpret the spectrum in terms of a one-dimensional superlattice induced by magnetic breakdown at Weyl points. The band width oscillations may be observed as $1/B$-periodic magnetoconductance oscillations, at weaker fields and higher temperatures than the Shubnikov-de Haas oscillations due to Landau level quantization. No such spectrum appears in a generic Weyl semimetal, the Kramers degeneracy at time-reversally invariant momenta is essential.


I. INTRODUCTION
Kramers-Weyl fermions are massless low-energy excitations that may appear in the Brillouin zone near time-reversally invariant momenta (TRIM). Their gapless nature is protected by Kramers degeneracy, which enforces a band crossing at the TRIM. Crystals that support Kramers-Weyl fermions have strong spin-orbit coupling and belong to one of the chiral point groups, without reflection or mirror symmetry, to allow for a linear rather than quadratic band splitting away from the TRIM. The materials are called topological chiral crystals or Kramers-Weyl semimetals -to be distinguished from generic Weyl semimetals where Kramers degeneracy plays no role. Several candidates were predicted theoretically [1,2] and some have been realized in the laboratory [3][4][5][6][7].
These recent developments have motivated the search for observables that would distinguish Kramers-Weyl fermions from generic Weyl fermions [8][9][10]. Here we report on the fundamentally different Landau level spectrum when the semimetal is confined to a thin slab in a perpendicular magnetic field.
Generically, Landau levels are dispersionless: The energy does not depend on the momentum in the plane perpendicular to the magnetic field B. In contrast, we have found that the Landau levels of a Kramers-Weyl semimetal are broadened into a Landau band. The band width oscillates periodically in 1/B, producing an oscillatory contribution to the magnetoconductance.
The phenomenology is similar to that encountered in a semiconductor 2D electron gas in a superlattice potential [11][12][13][14][15]. In that system the dispersion is due to the drift velocity of cyclotron orbits in perpendicular electric and magnetic fields. Here the surface Fermi arcs provide for open orbits, connected to closed orbits by magnetic breakdown at Weyl points (see Fig. 1).
No open orbits appear in a generic Weyl semimetal [16,17], because the Weyl points are closely separated inside the first Brillouin zone, so the Fermi arcs are short and do not cross the Brillouin zone boundaries (a prerequisite for open orbits). The Landau band dispersion therefore directly ties into a defining property [1] of a Kramers-Weyl semimetal: surface Fermi arcs that span the entire Brillouin zone because they connect TRIM at zone boundaries.
In the next two sections II and III we first compute the spectrum of a Kramers-Weyl semimetal slab in zero magnetic field, to obtain the equi-energy contours that govern the orbits when we apply a perpendicular field. The resonant tunneling between open and closed orbits via magnetic breakdown is studied in Sec. IV. With these preparations we are ready to calculate the dispersive Landau bands and the magnetoconductance oscillations in Secs. V and VI. The analytical calculations are then compared with the numerical solution of a tight-binding model in Secs. VII and VIII. We conclude in Sec. IX.

II. BOUNDARY CONDITION FOR KRAMERS-WEYL FERMIONS
The first step in our analysis is to characterize the surface Fermi arcs in a Kramers-Weyl semimetal, which requires a determination of the boundary condition for Kramers-Weyl fermions. This is more strongly constrained by time-reversal symmetry than the familiar boundary condition on the Dirac equation [18]. In that case the confinement by a Dirac mass V µ = µ(n · σ) generates a boundary condition Ψ = (n ⊥ ×n ) · σΨ. (2.1) The unit vectorsn andn ⊥ are parallel and perpendicular to the boundary, respectively. Although σ → −σ upon time reversal, the Dirac mass may still preserve time-reversal symmetry if the Weyl fermions are not at a time-reversally invariant momentum (TRIM). For example, in graphene a Dirac mass +µ at the K-point in the Brillouin zone and a Dirac mass −µ at the K -point preserves time-reversal symmetry. In contrast, for Kramers-Weyl fermions at a TRIM the V µ term in the Hamiltonian is incompatible with timereversal symmetry. To preserve time-reversal symmetry the boundary condition must couple two Weyl cones, it cannot be of the single-cone form (2.1).
FIG. 1: Electron orbits in a thin slab geometry perpendicular to a magnetic field (along the x-axis), for a generic Weyl semimetal [16,17] (at the left) and for a Kramers-Weyl semimetal (at the right). In each case we show separately a front view (in the x-y plane, to show how the orbits switch between top and bottom surfaces of the slab) and a top view (in the y-z plane, to indicate the magnetic flux enclosed by the orbits). The Kramers-Weyl semimetal combines open orbits (red arrows) with closed orbits enclosing either a large flux Φ or a small flux δΦ. Open and closed orbits are coupled by a periodic chain of magnetic breakdown events, spaced by l 2 m /a0 (with a0 the lattice constant and lm = /eB the magnetic length). The open orbits broaden the Landau levels into a band, the band width varies from minimal to maximal when δΦ is incremented by h/e. Because δΦ ∝ Bl 4 m ∝ 1/B, the band width oscillations are periodic in 1/B.
In App. A we demonstrate that, indeed, pairs of Weyl cones at a TRIM are coupled at the boundary of a Kramers-Weyl semimetal. Relying on that result, we derive in this section the time-reversal invariant boundary condition for Kramers-Weyl fermions.
We consider a Kramers-Weyl semimetal in a slab geometry, confined to the y-z plane by boundaries at x = 0 and x = W . In a minimal description we account for the coupling of two Weyl cones at the boundary. To first order in momentum k, measured from a Weyl point, the Hamiltonian of the uncoupled Weyl cones is (2. 2) The ± sign indicates whether the two Weyl cones have the same chirality (+) or the opposite chirality (−). The two Weyl points need not be at the same energy, we allow for an offset ε. We also allow for anisotropy in the velocity components v α . The σ α 's are Pauli matrices acting on the spin degree of freedom. We will also use τ α Pauli matrices that act on the Weyl cone index, with σ 0 and τ 0 the corresponding 2 × 2 unit matrix. We can then write The current operator in the x-direction is j + = v x σ x τ 0 for H + and j − = v x σ x τ z for H − . The time-reversal operation T does not couple Weyl cones at a TRIM, it only inverts the spin and momentum: An energy-independent boundary condition on the wave function Ψ has the general form [18] Ψ = M ± · Ψ, M ± = M † ± , M 2 ± = 1, (2.5) in terms of a Hermitian and unitary matrix M ± . The matrix M ± anticommutes with the current operator j ± perpendicular to the surface, to ensure current conservation. Time-reversal symmetry further requires that These restrictions reduce M ± to the single-parameter form The angle φ has a simple physical interpretation in the case H + , M + case of two coupled Weyl cones of the same chirality: It determines the direction of propagation of the helical surface states (the Fermi arcs). We will take φ = 0 at x = 0 and φ = π at x = W . This produces a surface state that is an eigenstate of τ y σ y with eigenvalue +1 on one surface and eigenvalue −1 on the opposite surface, so a circulating surface state in the ±ydirection. (Alternatively, if we would take φ = ±π/2 the state would circulate in the ±z-direction.) Notice that these are helical rather than chiral surface states: The eigenstates Ψ of τ y σ y with eigenvalue +1 contain both right-movers (σ y Ψ = +Ψ) and left-movers (σ y Ψ = −Ψ). This is the key distinction with surface states in a magnetic Weyl semimetal, which circulate unidirectionally around the slab [19][20][21][22].
In the case H − , M − that the coupled Weyl cones have the opposite chirality there are no helical surface states and the physical interpretation of the angle φ in Eq. (2.7) is less obvious. Since our interest here is in the Fermi arcs, we will not consider that case further in what follows.

A. Dispersion relation
We calculate the energy spectrum of H + with boundary condition M + from Eq. (2.7) along the lines of Ref. 23. Integration in the x-direction of the wave equation H ± Ψ = EΨ with k x = −i ∂/∂x relates the wave amplitudes at the top and bottom surface via Ψ(W ) = e iΞ Ψ(0), with As discussed in Sec. II we impose the boundary condition Ψ = M + (0)Ψ on the x = 0 surface and Ψ = M + (π)Ψ on the x = W surface. The round-trip evolution then gives the determinantal equation with the definitions In the zero-offset limit ε = 0 Eq. (3.4) reduces to the more compact expression (3.6) which is a squared Weiss equation [23,24].
The dispersion relation E(k y , k z ) which follows from Eq. (3.4) is plotted in Fig. 2. The surface states (indicated in red) are nearly flat as function of k z , so they propagate mainly in the ±y direction. In the limit ε → 0 the bands cross at k z = 0, this crossing is removed by the energy offset.

B. Fermi surface topology
The equi-energy contours E(k y , k z ) = E F are plotted in Fig. 3 for several values of W . The topology of the Fermi surface changes at a critical width At W = W c the surface bands from upper and lower surface touch at the Weyl point k y = k z = 0, and for larger widths the upper and lower surface bands decouple from a bulk band, in the interior of the slab. For ε = 0 the surface and bulk bands intersect at k z = 0 when W > W c . The gap δk y which opens up for nonzero ε is For later use we also record the area S 0 enclosed by the bulk band,

IV. RESONANT TUNNELING BETWEEN OPEN AND CLOSED ORBITS IN A MAGNETIC FIELD
Upon application of a magnetic field B in the xdirection, perpendicular to the slab, the Lorentz force causes a wave packet to drift along an equi-energy contour. Becausek = eṙ × B the orbit in real space is obtained from the orbit in momentum space by rotation over π/2 and rescaling by a factor /eB = l 2 m (magnetic length squared).
Inspection of Fig. 3 shows that for W > W c closed orbits in the interior of the slab coexist with open orbits on the surface. The open and closed orbits are coupled via tunneling through a momentum gap δk y (magnetic breakdown [25,26]), with tunnel probability T MB = 1 − R MB given by the Landau-Zener formula (4.1) In the expression for the breakdown field B c a numerical prefactor of order unity is omitted [26,27].
The real-space orbits are illustrated in Fig. 4: An electron in a Fermi arc on the top surface switches to the bottom surface when the Fermi arc terminates at a Weyl point [16]. The direction of propagation (helicity) of the surface electron may change as a consequence of the magnetic breakdown, which couples a right-moving electron on the top surface to a left-moving electron on the bottom surface. This backscattering process occurs with reflec- tion probability The phase shift φ accumulated in one round trip along the closed orbit is determined by the enclosed area S 0 in momentum space, with ν ∈ [0, 1) a magnetic-field independent offset. Resonant tunneling through the closed orbit, resulting in R = 1, happens when φ is an integer multiple of 2π. We thus see that the resonances are periodic in 1/B, with period (We have substituted the small-ε expression (3.9) for S 0 .) The Shubnikov-de Haas (SdH) oscillations due to Landau level quantisation are also periodic in 1/B. Their period is determined by the area Comparison with Eq. (4.4) shows that the period of the SdH oscillations is smaller than that of the magnetic breakdown oscillations by a factor k F a 0 (W/W c − 1) 3/2 , which is typically 1.

V. DISPERSIVE LANDAU BANDS
Let us now discuss how magnetic breakdown converts the flat dispersionless Landau levels into dispersive bands. The mechanism crucially relies on the fact that the surface Fermi arcs in a Kramers-Weyl semimetal connect Weyl points at time-reversally invariant momenta. Consider two TRIM K and K in the (k y , k z ) plane of the surface Brillouin zone. We choose K = (0, 0) at the zone center and K = (0, π/a 0 ) at the zone boundary, with G = (0, 2π/a 0 ) a reciprocal lattice vector.
In the periodic zone scheme, the Weyl points can be repeated along the k z -axis with period 2π/a 0 , to form an infinite one-dimensional chain (see Fig. 5). The perpendicular magnetic field B induces a flow along this chain in momentum space, which in real space is oriented along the y-axis with period In the weak-field regime l m a 0 the period L of the magnetic-field induced superlattice is large compared to the period a 0 of the atomic lattice. We seek the band structure of the superlattice.
We distinguish the Weyl points at K and K by their different magnetic breakdown probabilities, denoted respectively by T MB = 1 − R MB and T MB = 1 − R MB . We focus on the case that T MB and T MB are close to unity and the areas S 0 and S 0 of the closed orbits are the same -this is the small-ε regime in Eqs. (3.9) and (4.1). (The more general case is treated in App. C.) The phase shift ψ accumulated upon propagation from one Weyl point to the next is gauge dependent, we choose the Landau gauge A = (0, −Bz, 0). For simplicity we ignore the curvature of the open orbits, approximating them by straight contours along the line k y = E/ v y . The phase shift is then given by the same for each segment of an open orbit connecting two Weyl points. The quantization condition for a Landau level at energy E n is 2ψ + φ = 2πn, n = 1, 2, . . ., which amounts to the quantization in units of h/e of the magnetic flux through the real-space area S Σ l 4 m . Since S Σ S 0 the Landau level spacing is governed by the energy depen- dence of ψ, The Landau level spacing increases ∝ B and not ∝ √ B, as one might have expected for massless electrons. The origin of the difference is explained in Fig. 6. The where the phase φ is to be evaluated at E = (n − ν) ω c . Each Landau level is split into two subbands having the same band width The band width oscillates periodically in 1/B with period (4.4).

VI. MAGNETOCONDUCTANCE OSCILLATIONS
The dispersive Landau bands leave observable signatures in electrical conduction, in the form of magnetoconductance oscillations due to the resonant coupling of closed and open orbits. These have been previously studied when the open orbits are caused by an electrostatic superlattice [11][12][13][14][15]. We apply that theory to our setting.
From the dispersion relation (5.4) we calculate the square of the group velocity V = ∂E/ ∂k y , averaged over the Landau band, For weak impurity scattering, scattering rate 1/τ imp ω c , the effective diffusion coefficient [15], and the 2D density of states N 2D = (π v y a 0 ) −1 of the Landau band, determine the oscillatory contribution δσ yy to the longitudinal conductivity via the Drude formula for a 2D electron gas, The magnetoconductance oscillations due to magnetic breakdown (MB) coexist with the Shubnikov-de Haas (SdH) oscillations due to Landau level quantization. Both are periodic in 1/B, but with very different period, see Eqs. (4.4) and (4.5).
The difference in period causes a different temperature dependence of the magnetoconductance oscillations. A conductance measurement at temperature T corresponds to an energy average over a range ∆E ≈ 4k B T (being the full-width-at-half-maximum of the derivative of the Fermi-Dirac distribution). The oscillations become unobservable when the energy average changes the area S 0 or S Σ by more than π/l 2 m . This results in different characteristic energy or temperature scales,

VII. TIGHT-BINDING MODEL ON A CUBIC LATTICE
We have tested the analytical calculations from the previous sections numerically, on a tight-binding model of a Kramers-Weyl semimetal [1]. In this section we describe the model, results are presented in the next section.

A. Hamiltonian
We take a simple cubic lattice (lattice constant a, one atom per unit cell), when the nearest-neighbor hopping terms are the same in each direction α ∈ {x, y, z}. There are two terms to consider, a spin-independent term ∝ t 0 that is even in momentum and a spin-orbit coupling term ∝ t 1 σ α that is odd in momentum, The offset is arbitrarily fixed at −t 0 .
The geometry is a slab, with a normaln in the x-z plane at an angle φ with the x-axis (so the normal is rotated by φ around the y-axis). The boundaries of the slab are constructed by removing all sites at x < 0 and x > W . In the rotated basis aligned with the normal to the slab one has We will work in this rotated basis and for ease of notation omit the prime, writing k x or k ⊥ for the momentum component perpendicular to the slab and (k y , k z ) = k for the parallel momenta.

B. Folded Brillouin zone
The termination of the lattice in the slab geometry breaks the translation invariance in the perpendicular x- Since the Weyl points at Γ and M have the same chirality, for N + M even we are in the situation that the surface of the slab couples Weyl points of the same chirality -which is required for surface Fermi arcs to appear (see Sec. II). For N +M odd, in contrast, the Weyl points at the Γ and X points of opposite chirality are coupled by the surface, since these line up along the k ⊥ axis. Then surface Fermi arcs will not appear. In App. B we present a general analysis, for arbitrary Bravais lattices, that determines which lattice terminations support Fermi arcs and which do not.

VIII. TIGHT-BINDING MODEL RESULTS
We present results for M = N = 1, corresponding to a φ = π/4 rotation of the lattice around the y-axis.The folded and rotated Brillouin zone has a pair of Weyl points of + chirality at K = (0, 0, 0) and a second pair of − chirality at K = (π/a , 0, π/a ) in the rotated coordinates (see Fig. 9, second panel, with a = a √ 2). There is a second pair translated by k y = π/a. Each Weyl point supports a pair of Weyl cones of the same chirality, folded onto each other in the first Brillouin zone. The Weyl cones at K have energy offset ε = |2t 0 |, while those at K have ε = 0. We may adjust the offset by adding a rotational symmetry breaking term δH = δt 0 cos k z a to the tight-binding Hamiltonian (7.1). This changes the offsets into In Fig.10 we show how the Fermi arcs appear in the dispersion relation connecting the Weyl cones at k z = 0 and k z = π/a . This figure extends the local description near a Weyl cone from Fig. 2 to the entire Brillouin zone. The corresponding equi-energy contours are presented in Fig. 11. Increasing the spin-independent hopping term t 0 introduces more bands, but the qualitative picture near the center of the Brillouin zone remains the same as in Fig. 3 for W > W c .
The effect on the dispersion of a magnetic field B, perpendicular to the slab, is shown in Fig. 12 (see also App. D). The field was incorporated in the tight-binding model via the Peierls substitution in the gauge A = (0, −Bz, 0), with coordinate z restricted to |z| < L/2. Translational invariance in the y-direction is maintained, so we have a one-dimensional dispersion E(k y ). The boundaries of the system at z = ±L/2 introduce edge modes, which are visible in panel a as linearly dispersing modes near k y = ± 1 2 L/l 2 m (modulo π/a). Panels b,c,d focus on the , for the same system as in Fig. 10. In panels c and d the spin-independent hopping term t0 is increased by a factor 5 (at the same δt0 = −0.02 t1).
region near k y = 0, where these edge effects can be neglected. The effect on the dispersion of a variation in ε and ε is qualitatively similar to that obtained from the analytical solution of the continuum model, compare the four panels of Fig. 12 with the corresponding panels in Fig. 16. The width δE of the dispersive Landau bands (from maximum to minimum energy) is plotted as a function of 1/B in Fig. 13 and the periodicity ∆(1/B) is compared with the predicted Eq. (4.4) in Fig. 14. To remove the rapid Shubnikov-De Haas (SdH) oscillations we averaged over an energy interval ∆E around E F . This corresponds to a thermal average at effective temperature T eff = ∆E/4k B . From Eq. (6.4), with k F a ≈ 0.2, W/W c ≈ 1.5, we estimate that the characteristic energy scale at which the oscillations average out is five times smaller for the SdH oscillations than for the oscillations due to magnetic breakdown, consistent with what we see in the numerics.

IX. CONCLUSION
In conclusion, we have shown that Kramers-Weyl fermions (massless fermions near time-reversally invariant momenta) confined to a thin slab have a fundamentally different Landau level spectrum than generic massless electrons: The Landau levels are not flat but broadened with a band width that oscillates periodically in 1/B. The origin of the dispersion is magnetic breakdown at Weyl points, which couples open orbits from surface Fermi arcs to closed orbits in the interior of the slab.
The band width oscillations are observable as a slow modulation of the conductance with magnetic field, on which the rapid Shubnikov-de Haas oscillations are superimposed. The periodicities are widely separated because the quantized areas in the Brillouin zone are very different (compare the areas S 0 and S Σ in Fig. 5). This is a robust feature of the band structure of a Kramers-Weyl semimetal, as illustrated in the model calculation of Fig. 11. Since generic Weyl fermions have only the Shubnikov-de Haas oscillations, the observation of two distinct periodicities in the magnetoconductance would provide for a unique signature of Kramers-Weyl fermions. The dispersive Landau band is interpreted as the band structure of a one-dimensional superlattice of magnetic breakdown centra, separated in real space by a distance L = (eBa 0 /h) −1 -which in weak fields is much larger than the atomic lattice constant a 0 . Such a magnetic breakdown lattice has been studied in the past for massive electrons [26], the Kramers-Weyl semimetals would provide an opportunity to investigate their properties for massless electrons.

Acknowledgments
The tight-binding model calculations were performed using the Kwant code [29]. This project has received funding from the Netherlands Organization for Scientific Research (NWO/OCW) and from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme The derivation of the boundary condition for Kramers-Weyl fermions in Sec. II relies on pairwise coupling of Weyl cones at a TRIM by the boundary. Let us demonstrate that this is indeed what happens.
Consider a 3D Bravais lattice and its Brillouin zone. A time-reversally-invariant momentum (TRIM) is by definition a momentum K such that K = −K + G with G a reciprocal lattice vector, or equivalently, K = 1 2 G. Now consider the restriction of the lattice to x > 0, by removing all lattice points at x < 0. Assume that the restricted lattice is still periodic in the y-z plane, with an enlarged unit cell. Fig. 8 shows an example for a cubic lattice.
The enlarged unit cell will correspond to a reduced Brillouin zone, with a new set of reciprocal lattice vectors G. The original set K 1 , K 2 , K 3 , . . . of TRIM is folded onto a new setK 1 ,K 2 ,K 3 , . . . in the reduced Brillouin zone. The folding may introduce degeneracies, such that two different K's are folded onto the sameK. The statement to prove is this: • Each TRIMK in the folded Brillouin zone is either degenerate (because two K's were folded onto the sameK), or there is a second TRIMK along the k x -axis. Fig. 9 illustrates that this statement is true for the cubic lattice. We wish to prove that it holds for any Bravais lattice. Enlargement of the unit cell changes the primitive lattice vectors from a 1 , a 2 , a 3 intoã 1 ,ã 2 ,ã 3 . The two sets are related by integer coefficients n ij , n ij a j , n ij ∈ Z. (A1) The corresponding primitive vectors b,b in reciprocal space satisfy Any momentum k can thus be expanded as A TRIM K α in the first Brillouin zone of the original lattice is given by The index α labels each TRIM, identified by the 8 distinct triples (m α,1 , m α,2 , m α,3 ) ∈ Z 2 ⊗ Z 2 ⊗ Z 2 . Subsitution into the expansion (A3) gives  , for each triple ni1 ni2 ni3 and each triple mα,1 mα,2 mα,3 (both ∈ Z2 ⊗ Z2 ⊗ Z2). If we select any two rows and intersect with any column to obtain an ordered pair of values ν, ν , we can then find a second column with the same ν, ν at the intersection.
We now fold K α →K α into the first Brillouin zone of theb reciprocal vectors, In Table I we list for each TRIM and each choice of (n i1 , n i2 , n i3 ) ∈ Z 2 ⊗ Z 2 ⊗ Z 2 the corresponding value of ν α,i ∈ {0, 1 2 }. We fix the y and z-components ofK α by specifying ν α,2 and ν α,3 ∈ {0, 1 2 } and ask how many choices of α remain, so how many values of α satisfy the two equations Inspection of Table I shows that the number of solutions is even. More specifically, there are • 4 solutions if n 21 , n 22 , n 23 and n 31 , n 32 , n 33 are identical and different from 000 mod 2; • 2 solutions otherwise.
The multiple solutions correspond to pairs K α and K β that are either folded onto the sameK α =K β (if det n = 0 mod 2), or ontoK α andK β that differ only in the xcomponent (if det n = 1 mod 2). These are the TRIM that are coupled by the boundary normal to the x-axis. When the boundary couples only Weyl cones of the same chirality, these persist and give rise to surface Fermi arcs. If, however, opposite chiralities are coupled, then the boundary gaps out the Weyl cones and no Fermi arcs appear. Which of these two possibilities is realized can be determined by using that the parity of m α1 + m α2 + m α3 determines the chirality of the Weyl cone at K α . Table II identifies for each choice of n 21 , n 22 , n 33 and n 31 , n 32 , n 33 how many pairs of Weyl cones of opposite chirality are folded onto the same point of the surface Brillouin zone. We conclude that surface Fermi arcs appear if either • n 2i + n 3i = 1 mod 2 for each i, or    To calculate the effect of the coupling of open and closed orbits on the Landau levels we apply the scattering theory of Refs. 15, 26, 28 to the equi-energy contours shown in Fig. 15. We distinguish the two Weyl points at k z = 0 and k z = π/a 0 by their different magnetic breakdown probability, denoted respectively by T MB = 1 − R MB and T MB = 1 − R MB . The areas of the closed orbits may also differ, we denote these by S 0 and S 0 and the corresponding phase shifts by φ = S 0 l 2 m + 2πν and φ = S 0 l 2 m + 2πν. The coupling of the closed and open orbits at these two Weyl points is described by a pair of scattering matrices, given by for the Weyl point at k z = 0, and similarly for the other Weyl point at k z = π/a 0 (with T MB → T MB , φ → φ ). The coefficients can be rearranged in an energydependent transfer matrix, and similarly for T (with t → t , r → r ). The transfer matrices are energy dependent via the energy dependence of S 0 and hence of φ. Because det T total = 1, the eigenvalues of T total come in inverse pairs λ, 1/λ. The transfer matrix translates the wave function over a period L in real space, so we require that λ = e iqL for some real wave number q, hence λ + 1/λ = e iqL + e −iqL , or equivalently [28] tr T total (E) = 2 cos qL.
(C6) (In the main text we denote q by k y , here we choose a different symbol as a reminder that q is a conserved quan- For T MB and T MB close to unity an analytical solution E n (q) for the dispersive Landau bands can be obtained. We substitute ψ = π(n − ν) − (φ + φ )/4 + πδE/ ω c into Eq. (C5) and expand to second order in δE and to first order in R MB , R MB . Then we equate to 2 cos qL to arrive at E ± n (q) = (n − ν) ω c ± δE(q), (πδE/ ω c ) 2 = ρ + ρ + 2 ρρ cos qL, where φ and φ are evaluated at E = (n − ν) ω c . Corrections are of second order in R MB and R MB and we have assumed that the areas S 0 , S 0 of the closed orbit are small compared to k F /a 0 -so that variations of φ and φ over the Landau band can be neglected relative to the band spacing ω c .
Appendix D: Landau levels from surface Fermi arcs As explained in Fig. 6, the spacing of Landau levels formed out of surface Fermi arcs varies ∝ B -in contrast to the √ B dependence for unconfined massless electrons. In the tight-binding model of Sec. VIII we can test this by setting ε = ε = 0, so that there are only closed orbits and the Landau levels are dispersionless. The expected quantization is E n = (n−ν) ω c , ω c = eBv F a / , n = 0, 1, 2, . . . (D1) with v F the velocity in the surface Fermi arc, connecting Weyl points spaced by π/a . As shown in Fig. 17, this agrees nicely with the numerics.
In an unconfined 2D electron gas, the offset ν equals 1/2 or 0 for massive or massless electrons, respectively. For the surface Fermi arcs we observe that ν depends on the parity of the number of unit cells between top and bottom surface: ν = 0 if W/a is odd, while ν = 1/2 if W/a is even. This parity effect suggests that the coupling of Fermi arc states on opposite surfaces, needed to close the orbit in Fig. 1, introduces a phase shift that depends on the parity of W/a . We are not aware of such a phase shift for generic Weyl semimetals [16,17,30,31], it seems to be a characteristic feature of Kramers-Weyl fermions that deserves further study.