Deconfinement of Majorana vortex modes produces a superconducting Landau level
M. J. Pacholski, G. Lemut, O. Ovdat, İ. Adagideli, C. W. J. Beenakker
DDeconfinement of Majorana vortex modes produces a superconducting Landau level
M. J. Pacholski, G. Lemut, O. Ovdat, ˙I. Adagideli,
2, 3 and C. W. J. Beenakker Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands Faculty of Engineering and Natural Sciences, Sabanci University, Orhanli-Tuzla, Istanbul, Turkey MESA+ Institute for Nanotechnology, University of Twente, 7500 AE Enschede, The Netherlands (Dated: January 2021)A spatially oscillating pair potential ∆( r ) = ∆ e i K · x with momentum K > ∆ / (cid:126) v drives a decon-finement transition of the Majorana bound states in the vortex cores of a Fu-Kane heterostructure(a 3D topological insulator with Fermi velocity v , on a superconducting substrate with gap ∆ , ina perpendicular magnetic field). In the deconfined phase at zero chemical potential the Majoranafermions form a dispersionless Landau level, protected by chiral symmetry against broadening dueto vortex scattering. The coherent superposition of electrons and holes in the Majorana Landaulevel is detectable as a local density of states oscillation with wave vector (cid:112) K − (∆ / (cid:126) v ) . Thestriped pattern also provides a means to measure the chirality of the Majorana fermions. Introduction —
Deconfinement transitions in physicsrefer to transitions into a phase where particles can existas delocalized states, rather than only as bound states.Unlike thermodynamic phase transitions, the deconfine-ment transition is not associated with a spontaneouslybroken symmetry but with a change in the momentumspace topology of the ground state [1]. A prominent ex-ample in superconductors is the appearance of a Fermisurface for Bogoliubov quasiparticles when a supercon-ductor becomes gapless [2–5]. Such a Bogoliubov Fermisurface has been observed very recently [6].Motivated by these developments we consider here thedeconfinement transition for Majorana zero-modes in thevortex core of a topological superconductor. We willdemonstrate, analytically and by numerical simulations,that the delocalized phase at zero chemical potential re-mains a highly degenerate zero-energy level — a super-conducting counterpart of the Majorana Landau level ina Kitaev spin liquid [7, 8]. Unlike a conventional elec-tronic Landau level, the Majorana Landau level has anon-uniform density profile: quantum interference of theelectron and hole components creates spatial oscillationswith a wave vector set by the Cooper pair momentumthat drives the deconfinement transition.The system of Ref. 6 is shown in Fig. 1. It is a thin layerof topological insulator deposited on a bulk superconduc-tor, such that the proximity effect induces a pairing gap∆ in the surface states. A superflow with Cooper pairmomentum K lowers the excitation energy for quasipar-ticles with velocity v by the Doppler shift v · K , closingthe gap when vK exceeds ∆ . Following Fu and Kane[9], we add a perpendicular magnetic field B to confine aMajorana zero-mode to the core of each h/ e vortex thatpenetrates the superconductor. We seek to characterizethe deconfined phase that emerges when vK > ∆ . Confined phase —
To set the stage we first investigatethe confined phase for vK < ∆ . Electrons on the two-dimensional (2D) surface of a 3D topological insulatorhave the Dirac Hamiltonian v k · σ − µ , with µ the chem-ical potential, v the energy-independent Fermi velocity, k = ( k x , k y ) the momentum operator in the x – y surfaceplane, and σ = ( σ x , σ y ) two Pauli spin matrices. (The FIG. 1: Schematic of the Fu-Kane heterostructure [9], a topo-logical insulator with induced superconductivity (gap ∆ ) ina perpendicular magnetic field B . Vortices (red) bind midgapstates known as Majorana zero-modes. Here we study thedeconfinement transition in response to an in-plane super-current (blue arrows, momentum K ). When vK > ∆ thezero-modes delocalize into a Majorana Landau level. × σ is implicit when the Hamiltonian con-tains a scalar term.) Application of a perpendicular mag-netic field B (in the z -direction), adds an in-plane vectorpotential A = ( A x , A y ) to the momentum, k (cid:55)→ k − e A .The electron charge is + e and for ease of notation we willset v and (cid:126) both equal to unity in most equations.The superconducting substrate induces a pair potential∆ = ∆ e iφ . The phase field φ ( r ) winds by ± π aroundeach vortex, at position R n , as expressed by ∇ × ∇ φ ( r ) = ± π ˆ z (cid:80) n δ ( r − R n ) , ∇ φ = 0 . (1)The pair potential couples electrons and holes in the 4 × H = (cid:18) Kσ x + ( k − e A ) · σ ∆ e iφ ∆ e − iφ Kσ x − ( k + e A ) · σ (cid:19) , (2)at zero chemical potential, including a superflow momen-tum field K ≥ x -direction [10]. The superflowcan be a screening current in response to a magnetic fieldin the y -direction [6], or it can result from an externallyimposed flux bias or current bias. The Zeeman energy a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n from an in-plane magnetic field has an equivalent effect[3] (although it was estimated to be negligible relative tothe orbital effect of the field in the experiment [6]).The top and bottom surfaces of the topological insu-lator layer decouple if the layer thickness W is largerthan the penetration depth of the surface electrons (ofthe order of the atomic separation a ). Both surfacesmay still have induced superconductivity if W is less thanthe superconducting coherence length ξ = (cid:126) v/ ∆ . For vK < ∆ a pair of Majorana zero-modes will appear ineach vortex core, one at the top surface and one at thebottom surface. We can consider these separately.Setting ∆( r ) = ∆ ( r ) e ± iθ , in polar coordinates ( r, θ )for a ± π phase vortex at the origin, we need to solvethe zero-mode equation H ± Ψ ± = 0 with H ± = (cid:18) Kσ x − ( i ∇ + e A ) · σ ∆ ( r ) e ± iθ ∆ ( r ) e ∓ iθ Kσ x + ( i ∇ − e A ) · σ (cid:19) . (3)The pair potential amplitude ∆ ( r ) increases from 0 at r = 0 to a value ∆ > r (cid:38) ξ .When K = 0 this is a familiar calculation [11], whichis readily generalized to K >
0. The Majorana zero-mode has a definite chirality C , meaning that its four-component wave function Ψ ± is an eigenstate of thechirality operator Λ = diag (1 , − , − ,
1) with eigen-value C = ±
1. One has Ψ + = ( iψ + , , , ψ + ), Ψ − =(0 , iψ − , ψ − ,
0) with [12] ψ ± ( r ) = e ∓ Ky e ∓ χ ( r ) exp (cid:18) − (cid:90) r ∆ ( r (cid:48) ) dr (cid:48) (cid:19) , (4a) χ ( r ) = e π (cid:90) d r (cid:48) B ( r (cid:48) ) ln | r − r (cid:48) | . (4b)The factor e ∓ χ ( r ) is a power law for large r , so the zero-mode is confined exponentially to the vortex core as longas K < ∆ . When K > ∆ the solution (4) is no longernormalizable, it diverges exponentially along the y -axis.This signals a transition into a deconfined phase, whichwe consider next. Deconfined phase —
In Fig. 2 we show results from anumerical simulation of the deconfinement transition forthe model Hamiltonian described below. The left panelshows zero-modes confined to a pair of vortex cores for
K < ∆ , the right panel shows the deconfined state for K > ∆ . The decay | Ψ | ∝ e − Ky e − ∆ r in the confinedphase is anisotropic, with a decay rate ∆ along the x -axis and two different decay rates ∆ ± K in the ± y -direction. The direction into which the zero-mode decaysmore slowly is set by the chirality [14]: Fig. 2 shows C = +1 with a slow decay in the − y direction, for C = − y direction.In the deconfined phase the zero-mode density profilehas a pronounced periodic modulation in the x -direction,parallel to the superflow, with bifuration points at thevortex cores. This striped pattern is unexpected for aLandau level. We present an analytical description. Chiral symmetry protected Majorana Landau level —
The chiral symmetry of the Hamiltonian (2) plays a key
FIG. 2: Intensity profile | Ψ( x, y ) | of a Majorana zero-mode inthe vortex lattice [13]. The left panel shows the confined phase( K < ∆ ), the right panel the deconfined phase ( K > ∆ ).The dotted square indicates the unit cell containing a pairof h/ e vortices. These plots are for Majorana fermions ofpositive chirality, for negative chirality the density profile isinverted y (cid:55)→ − y . role in our analysis of the Majorana Landau level, sim-ilar to the role it plays for Landau level quantization ingraphene [15, 16] and in a Weyl superconductor [17]. Chi-ral symmetry means that H at µ = 0 anticommutes withΛ. The Hamiltonian then becomes block-off-diagonal inthe basis of eigenstates of Λ, U † HU = (cid:18) † (cid:19) , U = , (5a)Ξ = (cid:18) k − − eA − + K ∆ e iφ ∆ e − iφ − k + − eA + + K (cid:19) , (5b)where we have abbreviated k ± = k x ± ik y , A ± = A x ± iA y .A zero-mode is either a wave function ( u,
0) of positivechirality with Ξ † u = 0, or a wave function (0 , u ) of neg-ative chirality with Ξ u = 0. The difference between thenumber of normalizable eigenstates of either chirality iscalled the index of the Hamiltonian. It is topologicallyprotected, meaning insensitive to perturbations [18].Vortices are strong scatterers [19], completely obscur-ing the Landau level quantization in a nontopological su-perconductor [20]. Here chiral symmetry ensures thatthe vortices cannot broaden the zeroth Landau level. Helmholtz equation for the Majorana Landau level —
Let us focus on the Landau level of positive chirality,described by the equation Ξ † u = 0. This 2 × u ( r ) = e − Ky − q ( r ) e iφ ( r ) σ z ˜ u ( r ) , (6)with ∂ x q = − ∂ y φ + eA y , ∂ y q = ∂ x φ − eA x , (7) ⇒ (cid:18) − i∂ x + ∂ y ∆ ∆ i∂ x + ∂ y (cid:19) ˜ u = 0 . (8)The fields A , φ , and K no longer appear explicitly in thedifferential equation (8) for ˜ u , but they still determinethe solution by the requirements of normalizability andsingle-valuedness of the zero-mode u .Outside of the vortex core the spatial dependence ofthe pair potential amplitude ∆ may be neglected andone further simplification is possible: Substitution of ˜ u =( f, g ) gives g = ∆ − ( i∂ x − ∂ y ) f and a scalar second-orderdifferential equation for f , ∇ f = ∆ f. (9)In the context of classical wave equations this is theHelmholtz equation with imaginary wave vector.Eq. (6) requires that ˜ u and hence f have an expo-nential envelope e Ky in the y -direction. The Helmholtzequation (9) then ties that to a plane wave ∝ e ± iQx inthe x -direction, with wave vector Q = (cid:112) K − ∆ . Thisalready explains the striped pattern in the numerical sim-ulations of Fig. 2. For a more detailed comparison weproceed to a full solution of the Helmholtz equation. Analytical solution of the Majorana Landau level wavefunction —
The solutions of Eq. (9) for f are con-strained by the requirements of normalizability andsingle-valuedness of u . To determine the normalizabil-ity constraint we use that the field q ( r ) defined in Eq.(7) has the integral representation [21] q ( r ) = 12Φ (cid:90) d r (cid:48) B ( r (cid:48) ) ln | r − r (cid:48) | − (cid:88) n ln | r − R n | . (10)We consider N vortices (each of +2 π vorticity) in a re-gion S enclosing a flux Φ = N Φ , with Φ = h/ e thesuperconducting flux quantum [22]. If we set B → S , the field q ( r ) → (Φ / Φ − N ) ln r = 0 for r → ∞ . In view of Eq. (6), normalizability requires that e − Ky f is square integrable for r → ∞ . Near a vortexcore e − q f ∝ | r − R n | / f must be square integrable [23].Concerning the single-valuedness, the factor e iφ/ inEq. (6) introduces a branch cut at each vortex position R n , across which the function f should change sign —to ensure a single-valued u . This is a local constraint:branch cuts can be connected pairwise, hence there is nosign change in f on a contour encircling a vortex pair.We have obtained an exact analytical solution [24] ofthe Helmholtz equation in the limit that the separationof a vortex pair goes to zero. We place the two vorticesat the origin of a disc of radius R , enclosing a flux h/e ,with zero magnetic field outside of the disc. The enve-lope function then equals e − q ( r ) = r min e − r / R , with r min = min( r, R ).The two independent solutions are given by ˜ u =( f , f ) and ˜ u (cid:48) = σ x ˜ u ∗ , with f n = 2 i n e − inθ K n (∆ r ) − (cid:90) Q − Q dp C n ( p ) e ixp + y √ ∆ + p ,C n ( p ) = ∆ − n (∆ + p ) − / (cid:0) p − (cid:113) ∆ + p (cid:1) n . (11)The vortex pair is at the origin, with x + iy = re iθ , andK n is a Bessel function. FIG. 3: Dispersion relation of the topological superconductor,calculated from the model Hamiltonian (14) for zero magneticfield (black dashed lines, chemical potential µ = 0) and in thepresence of the magnetic vortex lattice (colored flat bandsat charge ± q eff e , for two values of µ ). For both data sets K = 2∆ = 20 (cid:126) v/d . The corresponding zero-modes follow from Eq. (6), u = e − q ( r ) e − Ky ( e iθ f , e − iθ f ) , u (cid:48) = σ x u ∗ . (12)For small r the zero-modes tend to a constant (the fac-tor 1 /r from K is canceled by the factor r from e − q ).The large- r asymptotics follows upon an expansion of theintegrand around the extremal points ± Q , giving f n → ( − n e Ky ∆ n (cid:18) ( K + Q ) n e − iQx iKx − Qy − ( K − Q ) n e iQx iKx + Qy (cid:19) . (13)The zero-modes decay as e − Ky f n ∝ /r for r (cid:29) R ,which needs to be regularized for a square-integrablewave function [25–27]. In a chain of vortices (spacing b ), the superposition of the solution (13) decays expo-nentially in the direction perpendicular to the chain [24].The decay length is λ = bK/Q or λ = bQ/K for a chainoriented along the x -axis or y -axis, respectively. Numerical simulation —
For a numerical study of thedeconfinement transition we represent the topological in-sulator layer by the low-energy Hamiltonian [28, 29] H ( k ) = ( v/a ) (cid:80) j = x,y σ j sin k j a + σ z M ( k ) − µ,M ( k ) = M − ( M /a ) (cid:80) j = x,y (1 − cos k j a ) , (14)in the basis Ψ = 2 − / ( ψ ↑ upper + ψ ↑ lower , ψ ↓ upper − ψ ↓ lower )of spin-up and spin-down states on the upper and lowersurfaces [30]. The atomic lattice constant is a , the Fermivelocity is v , and µ is the chemical potential. Hybridiza-tion of the states on the two surfaces introduces the massterm M ( k ). We set M = 0, to avoid the opening of agap at k = 0, but retain a nonzero M = 0 . a v in orderto eliminate the fermion doubling at a k = ( π, π ).In the corresponding BdG Hamiltonian the electronblock H ( k − e A + K ) is coupled to the hole block − H ( k + e A + K ) by the s -wave pair potential ∆ e iφ — which we take the same for both layers. We assumea strong type-II superconductor, for which we can take auniform magnetic field B and uniform pair potential am-plitude ∆ . The +2 π vortices are positioned on a square FIG. 4: Left panel: Numerically calculated intensity profile | Ψ( x, y ) | of the zeroth Landau level in a vortex lattice witha pair of h/ e vortices at the center of the unit cell ( K =2∆ = 40 (cid:126) v/d , µ = 0). Right panel: Analytical result fromthe solution of the Helmholtz equation (9) for a single h/e vortex [35]. lattice (lattice constant d = 302 a ) with two vorticesper unit cell.The spectrum is calculated using the Kwant tight-binding code [31, 32]. In Fig. 3 we show the dispersion-less Landau levels, both for chemical potential µ = 0and for nonzero µ . The zeroth Landau level has en-ergy E = ± q eff µ , with q eff e the charge expectationvalue. For the model Hamiltonian (2) we have [33] q eff = Q/K = (cid:112) − ∆ /K . The numerics at K = 2∆ gives a value 0.85, within 2% of (cid:112) / . E = E L ± q eff µ with E L = √ πq eff (cid:126) v/d , again in very good agreementwith the numerics. Notice that the flatness of the disper-sion persists at nonzero µ — even though the topologicalprotection due to chiral symmetry [34] is only rigorouslyeffective at µ = 0.In Fig. 4 we compare numerical and analytical resultsfor the case that the two h/ e vortices are both placedat the center of the unit cell. The agreement is quite sat-isfactory, given the different geometries (a vortex latticein the numerics, a single h/e vortex in the analytics). Striped local density of states —
The striped patternof the Majorana Landau level is observable by tunnelingspectroscopy, which measures the local density of states ρ ( r ) = (cid:80) k (cid:2) | ψ e ( r ) | f (cid:48) ( E − eV ) + | ψ h ( r ) | f (cid:48) ( E + eV ) (cid:3) , (15)averaged over the 2D magnetic Brillouin zone, (cid:80) k =(2 π ) − (cid:82) dk x dk y , weighted by the derivative of the Fermifunction. If E is much larger than temperature, the signof the bias voltage V determines whether the electroncomponent ψ e or the hole component ψ h contributes, sothese can be measured separately.As shown in Fig. 5, the oscillations are most pro-nounced for the hole component when µ > µ < V = ± E is an addi-tional experimental signature for the identification of the FIG. 5: Electron and hole contributions to the local density ofstates in the zeroth Landau level, along a line parallel to the x -axis which passes close through a vortex core at x = y =3 d /
4. The curves are plots of (cid:80) k | ψ e,h ( x, y ) | normalizedto unit peak height at the vortex core. The parameters are K = 2∆ = 40 (cid:126) v/d , µ = 0 . (cid:126) v/d . The expected oscillationperiod of π (cid:126) /Q = 0 . d is indicated. effect. Conclusion —
We have shown that the pair densitywave in a topological insulator with induced supercon-ductivity, observed in a recent experiment [6], can inducea deconfinement transition of Majorana vortex modes.For a superflow momentum K in the x -direction the con-finement is suppressed in either the + y or − y direction,depending on the chirality of the Majorana fermions.When (cid:126) vK exceeds the induced gap ∆ the zero-modesfrom a vortex lattice hybridize into a deconfined flatband. This Majorana Landau level is characterized by astriped interference pattern in the local density of states,with wave number Q = (cid:112) K − (∆ / (cid:126) v ) . These signa-tures of chirality and electron-hole interference charac-teristic of a Majorana fermion should be accessible in ascanning probe experiment.The Majorana Landau level provides a realization ofa flat band with extended wave functions, in which in-teraction effects are expected to be enhanced due tothe quenching of kinetic energy. Interacting Majoranafermions in a Fu-Kane superconductor have been studiedby placing vortices in close proximity inside a quantumdot [36]. The deconfinement transition provides a meansto open up the system and obtain a fully 2D flat bandwith widely separated vortices. An intriguing topic forfurther research is to investigate how the exchange ofvortices operates on this highly degenerate manifold. Acknowledgements —
This project has received fund-ing from the Netherlands Organization for Scientific Re-search (NWO/OCW) and from the European ResearchCouncil (ERC) under the European Union’s Horizon2020 research and innovation programme. [1] G. E. Volovik,
Quantum phase transitions from topologyin momentum space , Lect. Notes Phys. , 31 (2007).[2] D. F. Agterberg, P. M. R. Brydon, and C. Timm,
Bogoli-ubov Fermi surfaces in superconductors with broken time-reversal symmetry , Phys. Rev. Lett. , 127001 (2017).[3] Noah F. Q. Yuan and Liang Fu,
Zeeman-induced gaplesssuperconductivity with a partial Fermi surface , Phys. Rev.B , 115139 (2018).[4] S. Autti, J. T. M¨akinen, J. Rysti, G. E. Volovik, V. V.Zavjalov, and V. B. Eltsov, Exceeding the Landau speedlimit with topological Bogoliubov Fermi surfaces , Phys.Rev. Res. , 033013 (2020).[5] J. M. Link and I. F. Herbut, Bogoliubov-Fermi sur-faces in non-centrosymmetric multi-component supercon-ductors , Phys. Rev. Lett. , 237004 (2020).[6] Zhen Zhu, Micha(cid:32)l Papaj, Xiao-Ang Nie, Hao-Ke Xu,Yi-Sheng Gu, Xu Yang, Dandan Guan, Shiyong Wang,Yaoyi Li, Canhua Liu, Jianlin Luo, Zhu-An Xu, HaoZheng, Liang Fu, and Jin-Feng Jia,
Discovery of seg-mented Fermi surface induced by Cooper pair momentum ,arXiv:2010.02216. For a commentary on this experiment,see DOI: 10.36471JCCM October 2020 01[7] S. Rachel, L. Fritz, and M. Vojta,
Landau levels of Ma-jorana fermions in a spin liquid , Phys. Rev. Lett. ,167201 (2016).[8] B. Perreault, S. Rachel, F. J. Burnell, and J. Knolle,
Majo-rana Landau-level Raman spectroscopy , Phys. Rev. B ,184429 (2017).[9] Liang Fu and C. L. Kane, Superconducting proximity ef-fect and Majorana fermions at the surface of a topologicalinsulator , Phys. Rev. Lett. , 096407 (2008).[10] The term Kσ x in the BdG Hamiltonian (2) is equivalent,upon a gauge transformation, to a gradient Kx in φ .[11] R. Jackiw and P. Rossi, Zero modes of the vortex-fermionsystem , Nucl. Phys. B , 681 (1981).[12] To understand how the solution (4) relates to the K = 0 solution in Ref. [11], note the (non-unitary)transformation e Ky Λ H ± e Ky Λ = H ± + Kσ x , with Λ =diag (1 , − , − , ± is an eigenstate of Λwith eigenvalue ±
1, so if H ± Ψ ± = 0 for K = 0, then H ± e ± Ky Ψ ± = 0 for K (cid:54) = 0.[13] The data in Fig. 2 is obtained from the tight-bindingHamiltonian (14) of the topological insulator layer. Theparameters are ∆ = 20 (cid:126) v/d , d = 302 a , B = h/ed , µ = 0, M = 0, M = 0 . /a . The vortex pair in a unit cellis at the positions ( x, y ) = ( d / ,
1) and ( d / , K equals 0 . /v in the leftpanel and 2 ∆ /v in the right panel.[14] The anisotropic decay of the Majorana zero-mode in theleft panel of Fig. 2 can be understood as the effect of theMagnus force which the superflow momentum K = K ˆ x exerts on the axial spin S = C ˆ z of the Majorana fermions(as determined by their chirality C = ± K × S .[15] M. I. Katsnelson and M. F. Prokhorova, Zero-energystates in corrugated bilayer graphene , Phys. Rev. B ,205424 (2008).[16] J. Kailasvuori, Pedestrian index theorem `a la Aharonov-Casher for bulk threshold modes in corrugated multilayergraphene , EPL , 47008 (2009). [17] M. J. Pacholski, C. W. J. Beenakker, and I. Adagideli, Topologically protected Landau level in the vortex latticeof a Weyl superconductor , Phys. Rev. Lett. , 037701(2018).[18] Y. Aharonov and A. Casher,
Ground state of a spin-1/2charged particle in a two-dimensional magnetic field , Phys.Rev. A , 2461 (1979).[19] A. S. Mel’nikov, Quantization of the quasiparticle spec-trum in the mixed state of d-wave superconductors , J.Phys. Condens. Matter , 4219 (1999).[20] M. Franz and Z. Teˇsanovi´c, Quasiparticles in the vortexlattice of unconventional superconductors: Bloch waves orLandau levels? , Phys. Rev. Lett. , 554 (2000).[21] The integral equation (10) for q ( r ) follows from the def-inition (7), which implies that ∇ q ( r ) = ˆ z · ∇ × ( e A − ∇ φ ) = eB − π (cid:80) n δ ( r − R n ). The Green function of this2D Poisson equation is (2 π ) − ln | r − r (cid:48) | . Also note thatΦ ≡ π/e in units where (cid:126) ≡ S . Ifthe number of vortices is odd, a zero-energy edge statealong the perimeter of S will ensure that the total numberof Majorana zero-modes remains even.[23] This normalization requirement at the vortex core tiesthe chirality of the Majorana zero-modes to the sign ofthe vorticity. If we would have chosen − π vortices thefield q ( r ) would tend to + ln | r − R n | near a vortex core,and the product e − q f ∝ | r − R n | − / f would not havebeen square integrable.[24] Details of the solution of the Helmholtz equation aregiven in Apps. A and B of the Supplemental Material.[25] The 1 /r decay of the deconfined Majorana zero-mode im-plies a density of states peak which decays slowly ∝ / ln L as a function of the system size L . There is a formal simi-larity here with the zero-modes originating from vacanciesin a 2D bipartite lattice [26, 27].[26] B. Sutherland, Localization of electronic wave functionsdue to local topology , Phys. Rev. B , 5208 (1986).[27] V. M. Pereira, F. Guinea, J. M. B. Lopes-dos Santos,N. M. R. Peres, and A. H. Castro-Neto, Disorder inducedlocalized states in graphene , Phys. Rev. Lett. , 036801(2006).[28] Wen-Yu Shan, Hai-Zhou Lu, and Shun-Qing Shen, Effec-tive continuous model for surface states and thin films ofthree-dimensional topological insulators , New J. Phys. ,043048 (2010).[29] Song-Bo Zhang, Hai-Zhou Lu, and Shun-Qing Shen, Edgestates and integer quantum Hall effect in topological insu-lator thin films , Scientif. Rep. , 13277 (2015).[30] In the basis Ψ = ( ψ ↑ upper , ψ ↓ upper , ψ ↑ lower , ψ ↓ lower ) the4 × H = t (cid:80) j = x,y τ z σ j sin k j a + τ x σ M ( k ) − µ , with Pauli matrix τ z acting on the layer index. A unitary transformationblock-diagonalizes the Hamiltonian. One of the 2 × M replaced by − M .[31] C. W. Groth, M. Wimmer, A. R. Akhmerov, and X.Waintal, Kwant: A software package for quantum trans-port , New J. Phys. , 063065 (2014).[32] Details of the method of numerical simulation, with sup-porting data, are given in App. C of the Supplemental Material.[33] The renormalized charge q eff in the Majorana Landaulevel is calculated in App. D of the Supplemental Mate-rial. That calculation also gives the renormalized Fermivelocity v eff = √ v x vy = √ q eff v that appears in the Lan-dau level energy E L .[34] The chiral symmetry at µ = 0 is broken by the massterm M ( k ) in the Hamiltonian (14). This residual chiralsymmetry breaking is visible in Fig. 3 as a very smallsplitting of the µ = 0 Landau levels (green flat bands).[35] The comparison between numerics and analytics in Fig.4 involves no adjustable parameters. To compare the samestate in the degenerate zeroth Landau level we choose thestate with left-right reflection symmetry. There are two ofthese, the other is compared in App. E of the SupplementalMaterial.[36] D. I. Pikulin and M. Franz, Black hole on a chip: proposalfor a physical realization of the SYK model in a solid-statesystem , Phys. Rev. X , 031006 (2017). Appendix A: Solution of the Helmholtz equation forthe Majorana Landau level
The general solution of the 2D Helmholtz equation ∇ f = ∆ f that governs the Majorana Landau level isa superposition of waves e ipx ± y √ p +∆ . Which super-position we need is determined by the requirement that e − Ky − q ( r ) f ( x, y ) is square integrable in the x – y plane,with K > ∆ >
0. We denote Q = (cid:112) K − ∆ . For easeof notation we will set ∆ ≡ q ( r ) = (cid:15)r − N ln min( r, , N = 1 , , . . . , (A1)corresponding to 2 N vortices, each of vorticity +2 π , atthe origin. The positive infinitesimal (cid:15) > r → ∞ . The restriction toan even number of overlapping vortices means that thebranch cut which connects vortices pairwise can be ig-nored. (We have not succeeded in finding an analyticalsolution that incorporates the branch cut, but of coursein the numerics this is not a limitation.)The superposition of elementary solutions e ipx ± y √ p +1 that cancels the exponential growth factor e − Ky has thegeneral form f = (cid:82) | p | >Q dp C ( p ) e ipx + y √ p +1 if y < , − (cid:82) | p | . (A2)(We can use the symbol C twice without loss of generalitybecause the integration ranges do not overlap.)The solution should be continuously differentiable at r (cid:54) = 0, which is satisfied if f ( x, y ) and ∂ y f ( x, y ) are con-tinuous functions of y at y = 0, x (cid:54) = 0. The continuityrequirement is that the Fourier transform (cid:82) · · · e ipx dp of C ( p ) equals the Fourier transform of D ( p ) for x (cid:54) = 0,which means that C ( p ) and D ( p ) differ by a polynomial L ( p ) of p . [Recall that the Fourier transform of a poly-nomial is given by derivatives of δ ( x ).] Similarly, the re-quirement of a continuous derivative is that (cid:112) p + 1 C ( p )and − (cid:112) p + 1 D ( p ) differ by a polynomial T ( p ). Theunique solution of these two requirements is C ( p ) = T ( p ) (cid:112) p + 1 − L ( p ) ,D ( p ) = T ( p ) (cid:112) p + 1 + L ( p ) . (A3)We are free to choose a convenient basis for the poly-nomials T ( p ) and L ( p ), we will choose one for which theintegral over D ( p ) has a closed-form expression. The ba-sis polynomials T n ( p ) and L n ( p ), n = 0 , , , . . . are T n ( p ) = (cid:16) p + (cid:112) p + 1 (cid:17) n + (cid:16) p − (cid:112) p + 1 (cid:17) n ,L n ( p ) = (cid:16) p + (cid:112) p + 1 (cid:17) n (cid:112) p + 1 − (cid:16) p − (cid:112) p + 1 (cid:17) n (cid:112) p + 1 . (A4)This choice of basis is related to a basis of Chebyshevpolynomials T n , via the identities T n ( p ) = 2( − i ) n T n ( ip ) ,L n ( p ) = 2( − i ) n − n − (cid:88) m =0 T m − n +1 ( ip ) . (A5)Note that T − n ( p ) = ( − n T n ( p ) , L − n ( p ) = − ( − n L − n ( p ) . (A6)A complete basis for the pairs of polynomials T ( p ) , L ( p ) istherefore given by the two sets { T n , L n }∪{ T n , − L n } with n = 0 , , , . . . , or equivalently by the single set { T n , L n } with n = 0 , ± , ± , . . . . The corresponding basis of thefunctions C ( p ) and D ( p ) in Eq. (A3) is C n ( p ) = T n ( p ) (cid:112) p + 1 − L n ( p ) = (cid:16) p − (cid:112) p + 1 (cid:17) n (cid:112) p + 1 ,D n ( p ) = T n ( p ) (cid:112) p + 1 + L n ( p ) = (cid:16) p + (cid:112) p + 1 (cid:17) n (cid:112) p + 1 , (A7)with n = 0 , ± , ± , . . . .We next use the Bessel function identities [S1]K n ( r ) = (cid:40) i n e inθ (cid:82) ∞−∞ dp D n ( p ) e ipx − y √ p +1 if y ≥ , i n e inθ (cid:82) ∞−∞ dp C n ( p ) e ipx + y √ p +1 if y ≤ , (A8)where r = (cid:112) x + y and e iθ = ( x + iy ) /r , to write thesolution (A2) in the form f n ( x, y ) = − (cid:90) Q − Q dp (cid:16) p − (cid:112) p + 1 (cid:17) n (cid:112) p + 1 e ixp + y √ p +1 + 2 i n e − inθ K n ( r ) , (A9)which is Eq. (11) in the main text (upon restoring theunits of ∆ ).The function f n is the first component of the spinor˜ u = ( f, g ), the second component is g n = ( i∂ x − ∂ y ) f n = f n − . (A10)We now obtained an infinite countable set of solutions˜ u n = ( f n , f n − ), n = 0 , ± , ± , . . . of the Helmholtzequation, such that e − Ky e − (cid:15)r ˜ u n is square integrable atinfinity. The condition that r N ˜ u is square integrable atthe origin (containing 2 N overlapping vortices) selects afinite subset. For r → f n (cid:39) r −| n | if n (cid:54) = 0 and f (cid:39) ln r . Normalizability requires that both | n | ≤ N and | n − | ≤ N , hence there are 2 N allowed values of n ∈ {−N + 1 , −N + 2 , . . . N − , N } .All of this was for zero-modes Ψ = ( f, g, ,
0) of posi-tive chirality, in a lattice of +2 π vortices. Alternatively,we can consider zero-modes Ψ = (0 , , f, g ) of negativechirality in a lattice of − π vortices. The differentialequations for f and g remain the same, but now the ex-ponential factor that needs to be canceled is e Ky ratherthan e − Ky . The sign change gives the negative chiralitysolution f n ( x, y ) = − (cid:90) Q − Q dp (cid:16) p − (cid:112) p + 1 (cid:17) n (cid:112) p + 1 e ixp − y √ p +1 + 2 i n e inθ K n ( r ) , (A11a) g n = ( i∂ x − ∂ y ) f n = − f n +1 . (A11b) The 2 N zero-modes are now labeled by the index n ∈{−N , −N + 1 , . . . N − , N − } . Appendix B: Chain of vortices
The regularization at infinity by the (cid:15) term in Eq. (A1)is not needed if we have a periodic lattice of vortices. Wedemonstrate this by considering a linear chain of vorticesat positions R (cid:96) , spaced by b at an angle ϑ ∈ [0 , π/ x -axis. We take a linear superposition of thesolutions e − Ky f n ( r − R (cid:96) ) from Eq. (A9), with complexweights, F n ( r ) = ∞ (cid:88) (cid:96) = −∞ e i(cid:96)κ e (cid:96)Kb sin ϑ e − Ky f n ( r − R (cid:96) ) . (B1)We do not include the envelope e − q , because it tends tounity for large r if we set (cid:15) ≡
0. The Bloch phase κ isarbitrary.We substitute the large- r expansion (13), F n → ( − n ∞ (cid:88) (cid:96) = −∞ e i(cid:96)κ (cid:18) ( K + Q ) n e − iQ ( x − (cid:96)b cos ϑ ) iK ( x − (cid:96)b cos ϑ ) − Q ( y − (cid:96)b sin ϑ ) − ( K − Q ) n e iQ ( x − (cid:96)b cos ϑ ) iK ( x − (cid:96)b cos ϑ ) + Q ( y − (cid:96)b sin ϑ ) (cid:19) . (B2)We seek the decay of F n in the direction perpendicular to the chain, so for large | ρ | when ( x, y ) = ( − ρ sin ϑ, ρ cos ϑ ).We thus need to evaluate an infinite sum of the form[S2] S ( α, z ) = ∞ (cid:88) (cid:96) = −∞ e i(cid:96)α z + (cid:96) , α ∈ (0 , π ) , z ∈ C \ Z , (B3a) S ( α, z ) = 2 πie iαz − e i ( α − π ) z . (B3b)In the limit | Im z | → ∞ this tends to S ( α, z ) → (cid:40) − πie − (2 π − α )Im z if Im z → ∞ , πie α Im z if Im z → −∞ . (B4)Substitution of Eq. (B3) into Eq. (B2) gives, for x = − ρ sin θ , y = ρ cos θ , F n → ( − n ( K + Q ) n e iQρ sin ϑ Qb sin ϑ − iKb cos ϑ S ( α + , z − )+ ( − n ( K − Q ) n e − iQρ sin ϑ Qb sin ϑ + iKb cos ϑ S ( α − , z + ) , (B5)where we abbreviated α ± = κ ± Qb cos ϑ mod 2 π,z ± = ρb sin 2 ϑ ± iKQK − sin ϑ . (B6)Provided that α ± (cid:54) = 0 mod 2 π , the decay is exponen-tial: | F n | (cid:39) e − c | ρ | /λ , with (reinserting the units of ∆ ) λ = b K − ∆ sin ϑK (cid:112) K − ∆ (B7)and c a coefficient of order unity that depends on the signof ρ , c = (cid:40) min( α + , π − α − ) if ρ > , min( α − , π − α + ) if ρ < . (B8)For a chain oriented along the x -axis or y -axis we have λ equal to bK/Q or bQ/K , respectively. Appendix C: Details of the numerical simulation
The model Hamiltonian we consider is H = (cid:18) H ( k − e A + K ) ∆ e iφ ∆ e − iφ − H ( k + e A − K ) (cid:19) , (C1a) H ( k ) = ( v/a ) σ x sin a k x + ( v/a ) σ y sin a k y + σ z M ( k ) − µ, (C1b) M ( k ) = M − ( M /a )(2 − cos a k x − cos a k y ) . (C1c)The Hamiltonian acts on a spinor with the four compo-nentsΨ( k ) = 1 √ [ ψ ↑ upper + ψ ↑ lower ]( k )[ ψ ↓ upper − ψ ↓ lower ]( k ) − i [ ψ ↓ upper + ψ ↓ lower ] ∗ ( − k ) i [ ψ ↑ upper − ψ ↑ lower ] ∗ ( − k ) , (C2)for spin-up and spin-down electrons on the upper andlower surface of the topological insulator layer. The massterm M ( k ) is chosen with M = 0, M = 0 . a v , suchthat H has a single gapless Weyl point at k = 0. Nearthis Weyl point the upper and lower surface are uncou-pled, so the eigenstate can equivalently be written in thesingle-surface basis ( ψ ↑ , ψ ↓ , − iψ ∗↓ , iψ ∗↑ ).The Hamiltonian is discretized on a square lattice (lat-tice constant a ) with nearest neighbor hopping (hoppingenergy v/a ). The magnetic field B is uniform in the z -direction, vector potential A = − By ˆ x . The superflowmomentum is K = K ˆ x . The amplitude ∆ of the pairpotential is taken as a constant, the phase φ ( x, y ) windsby 2 π around each vortex.We take a square vortex lattice, with lattice constant d = N a . The flux through each magnetic unit cell is h/e , so it contains a pair of h/ e vortices. The integer N determines the magnetic field via B = ( N a ) − h/e .The vortices are placed on the diagonal of the magneticunit cell, at the positions ( x, y ) = ( N a / ,
1) and(
N a / , N twice an odd integer,we ensure that the singularity in the phase field at thevortex core does not coincide with a lattice point. Thephase field is discretized along the lines set out in App.B of Ref. S3. The eigenvalues and eigenfunctions of H are calculated using the Kwant tight-binding code [S4].Here we collect some additional results to those shownin the main text. In the confined phase vK < ∆ we showin Fig. 6 the anisotropic decay rates of the Majorana zero-modes bound to a vortex core, as in the left panel of Fig.2. FIG. 6: Decay rate of the Majorana mode confined to a vortexcore. The data from the numerical simulation (colored points,∆ = 20 v/d ) closely follows the analytical prediction | Ψ | ∝ e − Ky e − (∆ /v ) r (dashed lines). –X XMY Γ FIG. 7: Dispersion relation in zero magnetic field (blackdashed lines) and in the presence of the magnetic vortex lat-tice (green solid lines, the right panel shows the magneticBrillouin zone). Both band structures are for µ = 0, and thesame parameters as in Fig. 3. The red dots indicate the Weylpoints at k = ( ± Q,
0) in zero magnetic field. The Landaulevels are at ±√ n E L , n = 0 , ,
2, with E L = √ πq eff (cid:126) v/d . In the deconfined phase vK > ∆ we show in Fig. 7the Landau levels in the vortex lattice (complementingFig. 3). Fig. 8 shows the local density of states in thezeroth Landau level. This shows the variation over theentire unit cell of the vortex lattice, to complement theline cut through a vortex core shown in Fig. 5 of the maintext. Appendix D: Renormalized charge in the MajoranaLandau level
The charge expectation value of the deconfined zero-mode can be calculated by means of the block diago-nalization approach of Ref. S3. Starting from the BdGHamiltonian (2) we first make the gauge transformation
FIG. 8: Local density of states in the unit cell of the vor-tex lattice, at the energy E > µ > (cid:80) k | ψ e,h ( x, y ) | , summed over themagnetic Brillouin zone, normalized to unit maximum value.The white dotted line indicates the cut shown in Fig. 5 of themain text, at the same parameters. The electron contributionto the local density of states (right panel) and the hole con-tribution (left panel) can be measured separately by tunnelspectroscopy at voltages V = E and V = − E , respectively. H (cid:55)→ U † HU with U = (cid:18) e iφ
00 1 (cid:19) , resulting in H = (cid:18) ( k + a + q ) · σ − µ ∆ ∆ − ( k + a − q ) · σ + µ (cid:19) , a = ∇ φ, q = ∇ φ − e A + K ˆ x. (D1)We have included the chemical potential µ .For K > ∆ in zero magnetic field there are gaplessWeyl points at k = ( k x , k y ) = ( ˜ K,
0) with˜ K = ± κK, κ = (cid:113) − ∆ /K . (D2)To focus on the effect of a magnetic field on states near˜ K we set k x = ˜ K + δk x and consider δk x small.A unitary transformation H (cid:55)→ V † HV with V = (cid:18) σ cos( α/ σ x sin( α/ − σ x sin( α/ σ cos( α/ (cid:19) , (D3)tan α = − ∆ / ˜ K, cos α = − (1 + ∆ / ˜ K ) − / = − κ, approximately block-diagonalizes the Hamiltonian; the2 × δk x , a , q , and µ . The 2 × k =( κK,
0) is given by H + = κµ − ( κδk x + κa x − q x ) σ x +( k y + a y − κq y ) σ y , (D4)while the electron-like states near k = ( − κK,
0) are de-scribed by H − = − κµ + ( κδk x + κa x + q x ) σ x − ( k y + a y + κq y ) σ y . (D5) FIG. 9: Comparison between numerical and analytical inten-sity profiles | Ψ( x, y ) | , normalized to unit maximal value, forone of the two reflection-symmetric states in the zeroth Lan-dau level. The parameter values are the same as in Fig. 4,which compared the other state. The block diagonalization removes any interference be-tween the electron and hole blocks, so this approxima-tion cannot describe the striped density of states of Fig.2 — for that we need the Helmholtz equation consid-ered in the main text. Because the charge operatorˆ Q = − e∂H ± /∂µ = ∓ κe commutes with H ± , the ex-pectation value is given simply by (cid:104) ˆ Q (cid:105) = ∓ κe ⇒ q eff = κ. (D6)The Fermi velocity in the x -direction is renormalizedby the same factor, v x = κv , while v y is unaffected. Thisaffects the Landau level energy E L = √ π (cid:126) v eff /d of theanisotropic Dirac cone, via v eff = √ v x v y = √ κv . Appendix E: Comparison of numerics and analytics
In order to compare the analytic solution (11) of theHelmholtz equation with the numerical results from thetight-binding Hamiltonian (14) we proceed as follows.For the analytic solution we take a single pair of vor-tices located at r = 0, in a uniform magnetic field withtotal flux h/e in a large disc centered at the origin. Thereare then two independent zero-modes u , u (cid:48) given by Eq.(12) with q ( r ) = − ln r .For the numerical calculation we consider an infinitelattice of vortices, with pairs of vortices positioned atpoints R n = d n , n ∈ Z , in a uniform magnetic field B = ( h/e ) d − , vector potential A = − B ( y, T n = (cid:18) e ihn y x/d e − ihn y x/d (cid:19) T n ,T n r T † n = r + d n . (E1)(The 2 × e i k · n of the eigenstates defines0the magnetic momentum k ∈ [0 , π ) . At each value of k there are two independent zero-modes.To make sure we are comparing the same state in thedegenerate manifold we consider the operator product P x = (cid:18) e iφ ( r ) e − iφ ( r ) (cid:19) σ x P x (cid:18) e − iφ ( r ) e iφ ( r ) (cid:19) , (E2)with eigenvalues ±
1, which is a symmetry respected bothby the analytic and by the numerical calculation. Theoperator P x is the mirror symmetry operator in the x -direction, P x xP † x = − x , P x yP † x = y . (E3)The magnetic momentum transforms under P x as k x (cid:55)→− k x , k y (cid:55)→ k y . For the comparison we set k = 0, which is invari-ant under the action of P x . Then we can take the twozero-modes obtained numerically to be eigenstates of P x ,and compare them with the corresponding eigenstatesobtained analytically. Those are u ± ( r ) = u ( r ) ± u (cid:48) ( r ) , (E4)which, in view of the fact that f n ( − x, y ) = f ∗ n ( x, y ) (E5)are eigenfunctions of P x with eigenvalues ±
1. Figs. 4and 9 compare the modulus squared of the +1 and − P x respectively, with quite satisfactory cor-respondence. [S1] The identities (A8) follow from the integral representa-tion K n ( r ) = ( r/ n (cid:82) ∞ t − n − exp( − t − r /t ) dt , uponthe substitution p = ( t − /t ).[S2] For a derivation of Eq. (B3b), and its relation to theLerch zeta function, see https://mathoverflow.net/q/379157/11260 .[S3] M. J. Pacholski, C. W. J. Beenakker, and I. Adagideli, Topologically protected Landau level in the vortex latticeof a Weyl superconductor , Phys. Rev. Lett. , 037701(2018).[S4] C. W. Groth, M. Wimmer, A. R. Akhmerov, and X.Waintal,
Kwant: A software package for quantum trans-port , New J. Phys.16