Interacting Majorana fermions in strained nodal superconductors
IInteracting Majorana fermions in strained nodal superconductors
Emilian M. Nica and Onur Erten Department of Physics Box 871504 Arizona State University Tempe, Arizona 85287-1504 ∗ (Dated: February 20, 2019)Landau levels (LL) have been predicted to emerge in systems with Dirac nodal points under applied non-uniform strain. We consider 2D, d xy singlet (2D-S) and 3D p ± ip equal-spin triplet (3D-T) superconductors(SCs). We demonstrate the spinful Majorana nature of the bulk gapless zeroth-LLs. Strain along certain direc-tions can induce two topologically distinct phases in the bulk, with zeroth LLs localized at the the interface.These modes are unstable toward ferromagnetism for 2D-S cases. Emergent real-space Majorana fermions in3D-T allow for more exotic possibilities. Strain engineering has emerged in recent years as a promis-ing way to access topologically non-trivial phases. The initialproposal of Landau levels (LLs) in strained graphene [1] wasfollowed by a number of extensions, most notably to nodalDirac superconductors (SCs) [2, 3] and Weyl semi-metals andSCs [4, 5]. These intriguing proposals raise questions regard-ing the nature and the stability of the strain-induced gaplessmodes. Such issues have not been addressed in detail ex-cept in the case of strained graphene[6]. Here, we study thebulk gapless LLs which emerge under applied non-uniformstrain in prototypical models which describe a large classof well-studied 2D d-wave SCs, as exemplified by high- T c cuprates [7]. Furthermore, we consider possible instabilitiesdriven by small residual interactions in the low-energy sector.We likewise analyze 3D time-reversal symmetric equal-spintriplet SCs.We first elucidate the general spinful Majorana nature of thebulk zeroth LLs which satisfy γ s ( k ) = (cid:88) s (cid:48) M ss (cid:48) γ † s (cid:48) ( − k ) , (1)where s = ↑ / ↓ is a spin index with M = − σ y for 2D sin-glet (2D-S) SCs with d xy pairing. Similarly, M = − iσ z for3D triplet (3D-T) with p ± ip pairing. Subsequently, we showthat the topological properties of the strained bulk depend onthe direction of the uniaxial strain. When non-uniform strainis applied along an axis of the Brillouin Zone (BZ), a singlezeroth LL per spin s naturally emerges as a consequence of aphase transition in the bulk which separates two topologicallydistinct gapped regions. Whenever two independent zerothLLs per spin are present, as in the case of 2D-S SCs with strainalong a diagonal of the BZ, the bulk is topologically trivial.Furthermore, we show that the instabilities of these strainedsystems depend on the Majorana properties of the zerothLLs. In 2D-S systems, the halving of the degrees-of-freedom(DOF) associated with the zeroth LLs effectively excludes anyvalley-polarization instabilities, unlike the case in graphene[6]. Instead, the zeroth LLs in flat bands are generically unsta-ble to ferromagnetism induced by small residual short-rangeinteractions. By contrast, we show that the strained 3D-T SCsexhibit emergent real-space Majorana gapless modes in the ∗ Corresponding author: [email protected] 𝑘 𝑥 𝑘 𝑦 𝑘 𝑧 𝑥 𝑘 𝑥 𝑥 𝑘 𝑦 (1) (1’) (1) (1’) (2) (2’) 𝑥′ Trivial SC
Unstrained
Trivial SC
Strain
Topological SC th LL 𝑥 (a) (b) (c) FIG. 1. Nodal Dirac spectrum and Fermi surfaces for (a) 2D d xy pairing and (b) 3D ( p ± ip ) . Note the labeling of the valleys. (c) Thegreen and red allows indicate two directions of applied uniaxial strainalong the axes and diagonal directions, respectively. (c) Emergenceof the zeroth LLs at the boundary between strain-induced topologi-cally distinct phases in the bulk. bulk. Generic density-density interactions are suppressed inthese cases, but non-local spin-exchange four-fermion termsof the form U ( R , R (cid:48) ) c P,σ ( R ) c P, ¯ σ ( R ) c P,σ ( R (cid:48) ) c P, ¯ σ ( R (cid:48) ) areallowed and can lead to more exotic instabilities. The real-space Majorana fermions c P,σ ( R ) arise upon projection ontothe zeroth LL sector. Our conclusions are supported by de-tailed analytical calculations, which are available in the Sup-plementary Materials (SM), as well as by numerical resultspresented further down. Spinful Majorana fermions -
The pairing Hamiltonian underapplied uniaxial strain along the x direction can be written as H = (cid:88) k Ψ † ( k ) (cid:104) ˆ h ( k ) + ˆ∆( k ) (cid:105) Ψ ( k ) , (2)where Ψ T ( k ) = (cid:0) c σ ( x i , k ) , ( iσ y ) σσ (cid:48) c † σ ( x i , − k ) (cid:1) . The index i runs over all N x lattice sites along the x direction. The con-served momentum k is perpendicular to the direction of theapplied strain. The normal part given by the spin-independent ˆ h ∼ σ τ z contains the effects of strain. The pairing part isgiven by ˆ∆ = ˆ∆ S σ τ x , ˆ∆ T,x σ x τ x + ˆ∆ T,y σ y τ x , for 2D-S and a r X i v : . [ c ond - m a t . s up r- c on ] F e b σ and τ represent Paulimatrices in spin and Nambu spaces, respectively. The Hamil-tonian represents a set of effective 1D chains defined for eachconserved momentum k . For more details, we refer the readerto the SM.In the absence of any strain, the low-energy solutions of theHamiltonian determine four nodal points on a generic Fermisurface (FS) for the 2D-S and two nodal points for 3D-T cases.The resulting spectra are illustrated in Figs. 1 (a) and (b), re-spectively.By analogy with graphene [1], weak and non-uniform strainplays the role of an effective vector potential in the low-energytheory and the spectrum exhibits LLs. As pointed out inRef. 2, the strained system is time-reversal symmetric. Conse-quently, the pseudo-magnetic fields are not subject to Meiss-ner screening, in contrast to the case of genuine magneticfields where LLs are suppressed [8].In either 2D-S or 3D-T cases, the Bogoliubov-de Gennes(BdG) Hamiltonian can be decoupled into two independentspin sectors H s ( k ) . To illustrate, in the 2D-S case the two sec-tors involve c ↑ ( x i , k ) , c † ↓ ( x i , − k ) and their spin-flipped coun-terparts, respectively. In the 3D-T case, the two sectors in-volve fermions of equal spin. The BdG equations for eachspin sector are H s ( k )Ψ ns ( k ) = E ns Ψ ns ( k ) , where Ψ Tns =( u ns ( x i , k ) , v ns ( x i , k )) is a BdG spinor, s is the spin of theBdG quasiparticles, and n is a LL index. We define particle-hole (p-h) transformations [9] for singlet and triplet pairing as C S = ( iτ y ) ⊗ ( ) N × N K and C T = τ x ⊗ ( ) N × N K , re-spectively, where K represents complex-conjugation. Since CH s ( k ) C − = − H s ( − k ) , it follows that H s ( k ) has eigen-states Ψ sn ( k ) and C Ψ sn ( − k ) at energies E and − E , respec-tively.Whenever H s ( k ) exhibits a single zeroth LL, only one ofthe two states Ψ s ( k ) and C Ψ s ( − k ) is a non-trivial solution.In practice, only half of all momenta for each spin sector areassociated with independent DOF. We take this as a defini-tion of spinful Majorana modes. Equivalently, the analyticalsolutions for all cases (Sec. II C, III B, and IV B of the SM)indicate that the BdG quasiparticles in the zeroth LLs are pair-wise related via Eq. 1. In order to reconcile this redundancywith the BdG transformation, we constrain the zeroth LL sec-tor to half of the allowed momenta. This procedure preservesthe correct normalization of the c σ ( x i , k ) operators, as shownin Sec. II D of the SM. Topological origin of zeroth LLs -
Whenever strain inducesa single zeroth LL of spin s and momentum k y in the bulk,it is located at the interface of two fully-gapped and topologi-cally distinct phases. To illustrate, we consider an effective 1Dchain which is gapped in the pristine 2D-S case. Non-uniformstrain along one of the axes locally closes the gap in the bulkof the 1D chain, and induces a transition between two topo-logically distinct phases which gives rise to the lowest LL.This mechanism is illustrated in Figs. 1 (c) and by our numer-ical results presented further down. A similar picture emergesin the 3D T case. Whenever two independent zeroth LLs perspin s are induced in the bulk by applied strain, the gappedregions on either side are topologically trivial. This occurs inthe 2D-S case with strain along a diagonal (Fig.1(b)). -0.05-0.025 0 0.025 0.05 1.75 1.775 1.8 1.825 1.85(a) E k y S p e c t r a l W e i g h t Position
FIG. 2. (a) Spectrum for the 2D-S case about valley (1) as a functionof conserved momentum and strain along the x direction for the spinup sector. See text for the parameters of the calculation. (b) Evolutionof the spectral weight with momentum for the two degenerate statesat zero-energy. Note the formation of a broad LL state in the bulkaccompanied by an edge state. Thus the zeroth LLs are closely related to topologicallynontrivial edge states. Indeed, under applied strain, our ef-fective 1D and 2D models belong to the DIII class of theAltland-Zirnbauer classification [10]. The latter predicts Ma-jorana Kramers doublets [10] for point defects in 1D, whileline defects in 2D lead to helical Majorana fermions, both as-sociated with a non-trivial bulk Z invariant. This is preciselywhat we obtain in the 2D-S case with strain along an axis ofthe BZ and in the 3D-T case, respectively. Similarly, the pres-ence of two zeroth LLs in the 2D-S case with strain along adiagonal of the BZ signals topologically trivial bulk phases.We stress that the topological transition in the bulk via ap-plied strain does not imply destruction of the pairing ground-state. We argue by analogy to vortex states in standard Type-IISCs, where the condensate likewise survives [11]. Results -
We first consider the case of a 2D-S cases. Refer-ring to Fig. 1 (a), the unstrained system exhibits Dirac pointson the axes of the BZ, located at k = (0 , ± K F ) and labeledby valley indices (1) and (1 (cid:48) ) . The pair of valleys located at k = ( ± K F , are labeled by (2) and (2 (cid:48) ) , respectively. Weintroduce uniaxial strain along the x direction (green arrow inFig. 1 (a)) via a slowly-varying correction to the hopping term δt ( x i ) = t(cid:15)x i , where t is the pristine hopping coefficient. Asdiscussed in Sec. V, the strain parameter (cid:15) is proportional tothe gradient of the applied strain. The detailed form of the lat-tice Hamiltonian is given in Sec. II A of the SM. Next we solvethis system analytically in the continuum limit about each ofthe four valleys. Details of the calculation are presented inSec. II B and C of the SM. We find that strain induces LLs atvalleys (1) and (1 (cid:48) ) . Furthermore, the zeroth LL states obeyEq. 1. Therefore, the gapless modes of opposite valley andspin indices are not independent. Likewise, the numerical so-lutions for momenta approaching the Dirac points are con-sistent with the analytical results. To illustrate, in Fig. 2 (a) -18-15-12-9-6-3 00.0x10 -5 -5 -4 -4 -4 -4 -4 (a) l o g ( E / E ) ε -2
0 100 200 300 400 500(b) S p e c t r a l W e i g h t Position ε =0 ε =2.0x10 -5 ε =1.0x10 -4
0 100 200 300 400 500(c) Position ε =1.2x10 -4 ε =2.0x10 -4 ε =2.8x10 -4 FIG. 3. 2D-S case with strain along axis: (a) suppression of the bulkgap E for k y = 1 . as a function of increasing strain parameter (cid:15) .Beyond (cid:15) ≈ . × − the gap effectively vanishes. (b) Evolution ofthe spatially-resolved spectral weight with increasing strain indicatedin (a). Although the vanishing strain cannot close the gap, spectralweight accumulates at the edge. (c) Beyond a threshold strain (cid:15) =1 . × − , the gap suppression indicated in (a) is accompanied bythe emergence of the bulk zeroth LL and edge states. In between thetwo localized states, the chain is in a topologically non-trivial phase,as illustrated in Fig. 1 (c). The parameters of the calculation are thesame as in Fig. 2. we show the spectrum of H s = ↑ about valley (1) . The resultswere determined using a chain of 500 sites with lattice spac-ing a = 1 , pairing amplitude ∆ = 0 . , and strain parameter (cid:15) = 9 × − . We can distinguish the two zero-energy statesvia their respective spectral weights. In Fig. 2 (b) we showthe spectral weight | u ( x i , k y ) | + | v ( x i , k y ) | as a function ofposition along the effective 1D chain, for the several decreas-ing values of k y indicated in Fig. 2 (a). As we approach thecenter of the flat zeroth LL band, the spectral weight showswell-defined peaks in the bulk and at the edge, correspondingto a zeroth LL and an edge state, respectively. This confirmsthe topological origin of the zeroth LL modes as pointed outin Fig. 1 (c). Indeed, the bulk of the strained sample now con-tains two topologically-distinct phases, separated from eachother by the zeroth LL. Likewise, an edge state emerges at theboundary between the vacuum and the non-trivial sector. Aswe move away from the node, the zeroth LL merges with anedge state (orange points), separating the topologically trivialvacuum from the non-trivial bulk. These results also illustratethat there is only one bulk, zeroth LL per spin sector, confirm-ing it’s Majorana nature. A similar result is obtained for theopposite spin sector and about valley (1 (cid:48) ) .Our hypothesis can be further supported by considering theevolution of the gapped system at k y = 1 . , as a function ofapplied strain, as shown in Fig. 3 (a). Here E ≈ O (10 − ) is the value of the gap in the absence of strain. As the strainparameter increases, the spectral weight is increasingly local-ized at the edge, as shown in Fig. 3 (b). Upon reaching a strainparameter (cid:15) ≈ . × − , the gap first decreases dramati- cally then stabilizes close to zero. The corresponding spectralweights in Fig. 3 (c) exhibit two sharp peaks corresponding tozeroth LL and edge states, respectively.Next, we discuss the many-body instabilities of the spinfulMajorana fermions. We consider residual short range inter-actions of the Hubbard type in the paired state [12] (see SMSec. II E). Provided that the cyclotron frequency ω c ∼ √ B associated with the pseudo-magnetic field B is much smallerthan the interaction strength U’, we project the lattice oper-ators onto the zeroth LL sector. Since only half of the zeromodes are independent (Eq.1), we can formally eliminate thevalley (1 (cid:48) ) DOF (See SM Sec. II D). Consequently, we ob-tain an effective Hubbard interaction involving quasiparticlesof both spins at valley (1) . At mean-field level, this modelpredicts an instability toward ferromagnetism with a flat gapof magnitude n U (cid:48) / . Here, n ∼ l − B is the density of zerothLLs and l B is a pseudo-magnetic length. The remaining pos-sible instabilities are either energetically unfavored, as in thecase of singlet pairing, or are excluded due to the halving ofthe DOF in the case of a spin-density wave.Note that LLs do not emerge at valleys (2) and (2 (cid:48) ) . Asshown in Sec. II C of the SM, in the continuum limit, theeffects of the strain can be incorporated via a global phase.The resulting states remain gapless and exhibit an approxi-mate linear-in- q dispersion. Using standard RPA arguments,we see that coupling the projected zeroth LLs to the gaplessmodes at valleys (2) and (2 (cid:48) ) does not qualitatively changeour conclusion.We now consider 2D-S SCs with strain applied along thediagonal of the BZ (red arrow in Fig. 1 (a)). These differ fromthe analogous cases of strain along an axis of the BZ in severalways. First, there are only two valleys at k ≈ ± K F , wherethe conserved momentum is perpendicular to the diagonal ofthe BZ. Secondly, the 1D effective models are defined on twosublattices. Moreover, the Hamiltonians are invariant under areflection H s ( k ) = H s ( − k ) . Together with the p-h symmetrydiscussed previously, it ensures that zero-energy modes of thesame spin and at the same valley emerge in pairs with BdGcoefficients Ψ Ts ( k ) = ( u s ( x i , k ) , v s ( x i , k )) and C Ψ s ( k ) =( − v s ( x i , k ) , u s ( x i , k )) . Note the absence of the minus signfor the momentum of the second spinor. For more details,please consult Sec. III A of the SM.A consequence of the combined mirror and p-h symmetry isthat there are two independent zeroth LLs per valley and spin,as illustrated by analytical solutions of the BdG equations inthe continuum limit (Sec. III B of the SM). The numericalsolutions of the effective 1D lattice model likewise supportthis claim. In Fig. 4, we present our numerical results for achain of 500 sites with pairing amplitude ∆ = 10 − and strainparameter (cid:15) = 5 × − . The BdG coefficients of the twodegenerate zero-modes in panel (a) are precisely of the formimposed by the enhanced symmetry, as indicated in panels (b)and (c). This figure also illustrates the absence of any zero-energy edge states, thereby confirming our previous statementthat the bulk is in a topologically trivial phase.The pairs of independent zeroth LLs at opposite valleys arerelated via an analog of Eq. 1. In turn, each pair of indepen-dent zeroth LLs determine the projections of operators defined -0.1-0.05 0 0.05 0.1 0.7 0.8 0.9 1 1.1 1.2(a) E k y -0.15-0.075 0 0.075 0.15 0 100 200 300 400 500(b) x i u II ( x i , k y ) v I ( x i , k y ) 0 100 200 300 400 500(c) x i v II ( x i , k y ) u I ( x i , k y ) FIG. 4. 2D-S case with strain along the diagonal of the BZ. (a) Spec-trum as a function of conserved momentum. The degenerate states atzero energy are both in the bulk. See text for the parameters of thecalculation. (b) Real BdG coefficients as functions of position alongthe 1D chain for the two states indicated in (a). The two sets of co-efficients are related via u II = − v I , v II = u I due to the mirrorsymmetry as discussed in the main text. for the two sublattices. Upon inclusion of a repulsive Hubbardinteraction, as in the case of strain applied along an axis of theBZ, we find that the leading instability is also toward ferro-magnetism. A more detailed discussion is given in Sec. III Cand D of the SM. The ferromagnetic state shows a full gapin the bulk. This state would exhibit exponential activation intemperature dependence of the specific heat or the penetrationdepth to name a few.We discuss the case of 3D-T case next. In the absence ofstrain, we choose a Γ − representation of the tetragonal D h group [13] for the pairing. This leads to two nodal points at k z = ± K F for a spherical FS, as illustrated in Fig. 1 (b). Thistype of pairing can be thought of time reversal symmetric ana-log of the pairing in He-3A [14, 15]. As discussed in Sec. IV Aand B of the SM, under uniaxial non-uniform strain along the x -direction, and for fixed k z ≈ ± K z , the low-energy spec-trum of H s = ↑ along k y consists of single branches of right-moving bulk modes at each valley (1) and (1 (cid:48) ) . For fixed k y ≈ , a flat dispersion emerges along k z . Likewise, theopposite-spin sector H s = ↓ exhibits two left-moving modes ateach of the two valleys. Moreover, we find that for each spinsector, the modes at opposite valleys are related via Eq. 1.These solutions are consistent with the general discussion ofRef. 10, which predicts helical Majorana fermions associatedwith a line defect in a gapped 2D system belonging to the DIIIclass.The analytical results are confirmed by our lattice calcula- -0.2-0.1 0 0.1 0.2-0.15 -0.1 -0.05 0 0.05 0.1 0.15Spin up(a) E k y -0.2-0.1 0 0.1 0.2-0.15 -0.1 -0.05 0 0.05 0.1 0.15(c) Spin down E k y S p e c t r a l W e i g h t Position 0 0.03 0.06 100 200 300 400 500(b) S p e c t r a l W e i g h t Position
FIG. 5. 3D-T case (a) spectrum of spin up sector under uniaxialstrain along x for k z = K F corresponding to valley (1) in Fig. 1 (b),as a function of k y . The parameters of the calculation are discussedin the main text. (b) Spectral weights of the two counter-propagatingmodes indicated in (a) as functions of position along x . Only theright-moving modes are bulk zeroth LL states. (c) and (d) are thesame as (a) and (b) for the spin down sector. Here, the left-movingmode is in the bulk. tions. In Fig. 5 (a), we show the spin up spectrum as a functionof k y at valley (1) corresponding to k z ≈ K F . These resultswere obtained for a chain of 500 sites with a pairing ampli-tude of ∆ = 10 − and strain parameter (cid:15) = 9 × − . Asillustrated by Fig. 5 (b), only one of the two branches corre-sponds to bulk states, while the other is localized at the edge,confirming the scenario set out in Fig. 1 (c). Furthermore, wesee that the two branches are interchanged for the spin-downsector at the same valley, as indicated by Figs. 5 (c) and (d).The Fermi fields in the unpaired basis, projected onto thezeroth LL for each spin sector obey a real-space Majoranacondition (see SM Sec. IV C) c P,σ ( r ) = e iK Fz z c (1) P,σ ( r ) + e − iK Fz z c (1 (cid:48) ) P,σ ( r )=( iσ z ) σσ (cid:48) c † P,σ ( r ) , (3)since the projected operators at each valley are related via c (1 (cid:48) ) P,σ = (cid:80) σ (cid:48) ( iσ z ) σσ (cid:48) c (1) , † P,σ (cid:48) . Note that a similar relation doesnot hold in the singlet case, where zero-modes of oppositespin are related to each other.In contrast to the 2D-S case, the real-space Majorana na-ture of the projected fields in the 3D-T case imposes strongconstraints on possible residual interactions. It is well-known [16] that the product of identical real-space Ma-jorana operators reduces to a constant. This reflects thefact that ground-states of Majorana systems do not con-serve particle number [16]. It follows that any local den-sity c † σ ( R i ) c σ ( R i ) → c P,σ ( R i ) c P,σ ( R i ) = Constant.As such, density-density interactions of any range be-tween projected Majorana operators are dynamically triv-ial. However non-local, four-fermion exchange terms ofthe form U ( R , R (cid:48) ) c P,σ ( R ) c P, ¯ σ ( R ) c P,σ ( R (cid:48) ) c P, ¯ σ ( R (cid:48) ) are al-lowed. These are unlike spinless Majorana fermions where aminimal interaction is defined on a plaquette [17–20]. To ourknowledge, interactions between spinful Majorana fermionshave not been proposed in a condensed matter setting. We an-ticipate that they can lead to more exotic phases which requirefurther investigation. Discussion -
We showed that non-uniform uniaxial strainleads to the formation of spinful Majorana fermions in SCswith a Dirac nodal spectrum. We also demonstrated that theapplied strain can induce topological transitions in the bulk ofthese gapless SCs. The resulting zeroth LLs are localized atthe boundary between two topologically distinct phases withinthe strained sample. In many ways, these bulk states are anal-ogous to the better-known cases of Majorana fermions local-ized in vortex cores of topological SCs [21]. We note thata similar correspondance between edge states and the zerothLLs in a strained Weyl semi-metal was pointed out in Ref. 4.To our knowledge, the more general mechanism involving a local closing of the gap has not been clearly identified as such.We estimate a pseudo-magnetic length l B ≈ a , where a is the inter-atomic distance. We obtain a separation in energyof the LLs E c ≈ . meV, corresponding to a characteristictemperature scale of 1.2 K. These estimates were determinedusing parameters typical of the cuprate class of SCs and arediscussed in detail in Sec. V of the SM.As we have shown, the direct observation of these gap-less states is precluded by various instabilities. In strained2D d-wave SCs, the flat zeroth-LL band is unstable towardferromagnetism. One immediate consequence is the coexis-tence of superconductivity and ferromagnetism in the bulk ofa strained d-wave SC. The finite magnetic moment in the bulkis screened via the Meissner effect. However, signatures ofthe resulting screening currents could be detected by SQUID[22] or NMR experiments. For 3D triplet pairing, the low-energy gapless states are emergent real-space spinful Majo-rana fermions, which are likely to host more exotic phases. Acknowledgements
We thank Joel Moore for fruitful dis-cussions. This work is supported by ASU startup grant. [1] F. Guinea, M. I. Katsnelson, and A. K. Geim, Nat. Phys. , 30(2009).[2] E. M. Nica and M. Franz, Phys, Rev. B , 024520 (2018).[3] G. Massarelli, G. Wachtel, J. Y. T. Wei, and A. Paramekanti,Phys. Rev. B , 224516 (2017).[4] A. Grushin, J. W. F. Venderbos, A. Vishwanath, and R. Ilan,Phys. Rev. X , 041046 (2016).[5] T. Liu, D. I. Pikulin, and M. Franz, Phys. Rev. B , 041201(R)(2017).[6] P. Ghaemi, J. Cayssol, D. N. Sheng, and A. Vishwanath, Phys.Rev. Lett. , 266801 (2012).[7] M. Hashimoto, I. M. Vishik, R.-H. He, T. P. Devereaux, andZ.-X. Shen, Nat. Phys. , 483 (2014).[8] M. Franz and Z. Te˘sanovi´c, Phys. Rev. Lett. , 554 (2000).[9] M. Sato and Y. Ando, Rep. Prog. Phys. , 076501 (2017).[10] J. C. Y. Teo and C. L. Kane, Phys. Rev. B , 115120 (2010).[11] P. G. de Gennes, Superconductivity of Metals and Alloys (West-view, Boulder, 1999).[12] A. C. Potter and P. A. Lee, Phys. Rev. Lett. , 117002 (2014).[13] M. Sigrist and K. Ueda, Rev. Mod. Phys. , 239 (1991).[14] A. J. Leggett, Rev. Mod. Phys. , 331 (1975).[15] G. E. Volovik, The Universe in a Helium Droplet (Oxford Uni-versity Press, NY, 2003).[16] S. R. Elliott and M. Franz, Rev. Mod. Phys. , 137 (2015).[17] A. Rahmani, X. Zhu, M. Franz, and I. Affleck, Phys. Rev. B ,235123 (2015).[18] I. Affleck, A. Rahmani, and D. Pikulin, Phys. Rev. B , 125121(2017).[19] K. Wamer and I. Affleck, Phys. Rev. B , 245120 (2018).[20] A. Rahmani, D. Pikulin, and I. Affleck, Phys. Rev. B , 085110(2019).[21] L. Fu and C. L. Kane, Phys. Rev. Lett. , 096407 (2008).[22] S. Frolov, M. J. A. Stoutimore, T. A. Crane, D. J. Van Harlin-gen, V. A. Oboznov, V. V. Ryazanov, A. Ruosi, C. Granata, andM. Russo, Nat. Phys. , 32 (2008).[23] A. Auerbach, Interacting Electrons and Quantum Magnetism (Springer-Verlag, N.Y., 1994).[24] A. H. Castro, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. , 109 (2009).[25] A. Messiah, Quantum Mechanics (Dover, N. Y., 1999) p. 492.[26] R. Jackiw and C. Rebbi, Phys. Rev. D , 3398 (1976).[27] A. V. Balatskii, G. E. Volovik, and V. A. Konyshev, Zh. Eksp.Teor. Fiz. , 2038 (1986).[28] M. Vozmediano, M. Katsnelson, and F. Guinea, Phys. Rep. ,109 (2010).[29] M. M. Korshunov, V. A. Gavrichkov, S. G. Ovchinnikov,D. Manske, and I. Eremin, Phys. C , 365 (2004).[30] M. Mouallem-Bahout, J. Gaud´e, G. Calvarin, J.-R. Gavarri, andC. Carel, Mater. Lett. , 181 (1994). Supplementary Materials for “Interacting Majorana fermions in strained nodal superconductors”
The Supplementary Materials contain detailed derivations of our results and discussions and in support of the argumentspresented in the main text. Sec. I reviews important aspects of the pairing Hamiltonian and the Bogoliubov-de Gennes (BdG)equations. It provides general arguments in support of the Majorana nature of the zeroth Landau Levels (LLs) and of thetruncation of their Hilbert space, as discussed in the main text. Sec. II considers the case of 2D d-wave superconductors (SCs)under varying uniaxial strain along one axis of the Brillouin Zone (BZ). We introduce the effective lattice models and presenttheir analytical solutions in the continuum limit. This subsection also describes the gapless solutions which do not correspond toLandau quantization, and briefly elaborates on the connection with the Jackiw-Rebbi equations. We also illustrate the projectionof the lattice operators to the zeroth LL sector, and discuss the ensuing ferromagnetic instability. Sec. III considers the case of2D d-wave SCs with strain along the direction of a diagonal of the BZ. The major points of Sec. II are also covered here. Sec. IVfocuses on the case of uniaxial strain in 3D equal-spin triplet superconductors. In particular, it presents the effective latticemodels and their analytical solutions in the continuum limit. It also discusses the projection of the lattice operators onto thezeroth LL sector and illustrates their emergent real-space Majorana nature. Sec.V presents our estimates of the pseudo-magneticlength and LL spacing in energy for realistic values of the non-uniform strain in high- T c compounds. Unless stated otherwise,we work in units where (cid:126) = 1 . CONTENTS
References 5I. General Aspects of the Bogoliubov-de Gennes Hamiltonian 6II. Two dimensional d-wave pairing with strain along one nodal axis 8A. 1D Lattice Hamiltonian in the presence of uniaxial strain 8B. Low-energy effective Hamiltonian under applied uniaxial strain 9C. Landau levels from the low-energy Hamiltonian 11D. Projection of lattice operators onto the zeroth LL modes 13E. Ferromagnetism at mean-field level 15III. Two dimensional d-wave pairing with strain off the nodal axes 16A. 1D Lattice Hamiltonian 16B. Landau levels in the continuum limit 18C. Comparison to numerical solution of the lattice model 20D. Projection onto zeroth LLs 21E. Ferromagnetism at mean-field level 22IV. Helical Majoranas for 3D equal-spin triplet SCs 22A. Lattice Hamiltonian 23B. Landau levels in the continuum limit 23C. Projection of fields onto the zeroth LLs 26V. Estimates of strain 27A. General Formulae 27B. Numerical Estimates 28
I. GENERAL ASPECTS OF THE BOGOLIUBOV-DE GENNES HAMILTONIAN
In this section, we discuss the general properties of spinful Majorana fermions which occur for time-reversal invariant sin-glet and equal spin-triplet pairing Hamiltonians under applied strain. The general discussion in this section is exemplified insubsequent sections which treat specific cases in detail.We define a lattice Hamiltonian in D = 2 , dimensions under applier strain as H Pairing = (cid:88) k c † ( k ) H ( k ) c ( k ) , (S4)where spinors are given by c ( k ) = c ↑ ( x i , k ) c ↓ ( x i , k ) c † ↓ ( x i , − k ) − c † ↑ ( x i , − k ) . (S5)The position x i , i ∈ { , . . . , N } indicates the position along the direction of the applied uniaxial strain, while k represents the D − component conserved momentum along the directions which are perpendicular to the applied strain. The Hamiltoniancan be written as H ( k ) = (cid:18) h ( x i ; x j ; k ) σ ∆ ( x i ; x j ; k )∆ † ( x i ; x j ; k ) − h ∗ ( x i ; x j ; k ) σ (cid:19) . (S6)Due to the presence of the strain, each block is a matrix which depends on both position indices, as well as the momentum.In either singlet or equal-spin triplet pairing cases, we can separate Hamiltonian into two independent spin sectors as H ( k ) = (cid:18) H s = ↑ ( k ) 00 H s = ↓ ( k ) (cid:19) , (S7)corresponding to spin ( ↑↓ ) , ( ↓↑ ) sectors in the singlet case, and spin ( ↑↑ ) and ( ↓↓ ) sectors in the equal-spin triplet cases.The solutions of each H s are then obtained via a standard Bogoliubov-de Gennes (BdG) ansatz [16] γ s,n ( k ) = (cid:88) x i ,σ u ∗ sσ,n ( x i , k ) c σ ( x i , k ) + v ∗ sσ,n ( x i , k ) c † σ ( x i , − k ) , (S8)for energies E ≥ . We likewise define the inverse transformation as c σ ( x i , k ) = (cid:88) n,s u sσ,n ( x i , k ) γ s,n ( k ) + v ∗ sσ,n ( x i , − k ) γ † s,n ( − k ) . (S9)The index s also labels the conserved spin of the BdG quasiparticles. For the case of singlet pairing u sσ = δ sσ u s and v sσ =( iσ y ) sσ v s . In the case of equal-spin triplet pairing u has an identical form while v sσ = − ( σ y ) sσ v s . The index n represents allother labels. In the vicinity of the nodal points, it labels the bulk Landau Levels (LLs).As discussed below and in the following sections, we find that only half of the degrees-of-freedom for the zeroth LL sector areindependent since γ s, ( k ) = M ss (cid:48) γ † s (cid:48) , ( − k ) , where M = σ y for 2D-S, and M = − iσ z for singlet and triplet cases, respectively.Note that the inverse BdG transformation in Eq. S9 implicitly assumes that the two components are independent of each other.Therefore, in order to preserve a well-defined BdG transformation, we restrict the u, v indices to half of the allowed momenta.This procedure is illustrated in Sec. II D.The Majorana nature of the zeroth LLs follows from the general p-h symmetry of the Hamiltonian. Indeed, each H s sectorobeys the BdG eqs. H s ( k ) Φ s,n ( k ) = E n ( k )Φ s,n ( k ) , (S10)where the BdG spinor is Φ s,n ( k ) = (cid:18) u sσ,n ( x i , k ) v sσ (cid:48) ,n ( x i , k ) (cid:19) . (S11)We define a particle-hole (p-h) transformation C Singlet = ( − iτ y ) N × N K, (S12)for singlet and C Triplet = τ x N × N K, (S13)for triplet pairing, respectively. The N × N identity matrix is defined for the position indices and K represents complexconjugation. The two different forms are due to the distinct parities of the two types of pairing. The p-h transformation acts onthe Hamiltonians for each of the two spin sectors as [9] CH s ( k ) C − = − H s ( − k ) . (S14)As is well-know, this implies that the eigenstates of H s ( k ) come in pairs of opposite energy: H s ( k ) [ C Φ s,n ( − k )] = − E n ( − k ) [ C Φ s,n ( − k )] . (S15)When strain induces a single, zeroth LL mode, p-h symmetry imposes an additional constraint. Since only one of either Φ s, ( k ) or C Φ s, ( − k ) is a non-trivial solution of H σσ (cid:48) ( k ) , we require C Φ s, ( − k ) = (cid:18) ± v ∗ s, ( x i , − k ) u ∗ s, ( x i , − k ) (cid:19) (S16) =0 . (S17)This necessarily implies that the Hilbert space of the zeroth LLs must be truncated to half of all allowed momenta k .The above argument excludes the case where the spinor is invariant under the p-h transformation: C Ψ s ( − k ) ∝ Ψ s ( k ) . Thisis confirmed by explicit solutions presented in Sec. II C.Finally, we discuss the case where zero-energy LLs occur in pairs, as discussed in detail in Sec. III B and III C. Consider thecase where H s are invariant under a mirror plane involving the conserved momentum k : H s ( k ) = H s ( − k ) . (S18)Together with the p-h transformation − H s ( − k ) [ C Φ s,n ( k )] = E n ( k ) [ C Φ s,n ( k )] , (S19)this implies H s ( k ) [ C Φ s,n ( k )] = − E n ( k ) [ C Φ s,n ( k )] . (S20)Therefore, zero-energy solutions occur in pairs and are determined by Ψ s, and C Ψ s, ( k ) . II. TWO DIMENSIONAL D-WAVE PAIRING WITH STRAIN ALONG ONE NODAL AXIS
In this section, we discuss the case of a 2D d xy SC with strain along one of the nodal axes of the pairing. In Sec. II A,we discuss the effective 1D lattice Hamiltonians which are appropriate under applied uniaxial strain. An effective low-energyHamiltonian is derived in Sec. II B and it’s detailed LL solutions are presented in Sec. II C. Sec. II D discuss the resultingprojection of the lattice operators onto the zeroth LL sector. The last subsection is devoted to the discussion of a short-rangeHubbard interaction projected onto the zeroth LLs and the ensuing ferromagnetic instability.
A. 1D Lattice Hamiltonian in the presence of uniaxial strain
Using the convention of Eq. S5, we consider a simple lattice model for a 2D d-wave consisting of nearest-neighbor (NN)hopping together with next-nearest-neighbor (NNN) d xy pairing on a square lattice: H = H T B + H ∆ (S21) H T B = (cid:88) x i ,y i ,j,σ t ( x i ) (cid:2) c † σ ( x i , y i ) c σ ( x i + δ x,j , y i + δ y,j ) + c † σ ( x i + δ x,j , y i + δ y,j ) c σ ( x i , y i ) (cid:3) − µc † σ ( x i , y i ) c σ ( x i , y i ) (S22) H Pair = (cid:88) x i ,y i ,j,σ,σ (cid:48) ( iσ y ) σσ (cid:48) ∆( δ (cid:48) j ) (cid:2) c σ ( x i , y i ) c σ (cid:48) ( x i + δ (cid:48) x,j , y i + δ (cid:48) y,j ) + c σ ( x i , y i ) c σ (cid:48) ( x i − δ (cid:48) x,j , y i − δ (cid:48) y,j ) (cid:3) + H.c. . (S23)Note that the pairing amplitude is independent of spin indices σ, σ (cid:48) . We define the NN and NNN vectors δ = ( a, , δ = (0 , a ) and δ (cid:48) = ( a, a ) , δ (cid:48) = ( − a, a ) , where a is the NN distance. The sign change of the d xy pairing is induced via ∆( δ (cid:48) / ) = ± ∆ .Without any strain, we have t ( x i ) = t and the resulting dispersion shows gapless excitations around the four momenta ( ± K F , and (0 , ± K F ) , which denote the positions of four valleys denoted by (2) , (2 (cid:48) ) and (1) , (1 (cid:48) ) , respectively.We allow for the effects of strain by introducing a position dependent hopping along the x direction t ( x i ) = t + δt ( x i ) (S24) δt ( x i ) = t(cid:15)x i , (S25)where (cid:15) is approximately proportional to the gradient of the strain along x (Sec. V A). The modification to the unperturbedhopping t is due to the non-uniform stretching of the NN bonds. For a microscopic derivation of this term, we refer the readerto the supplementary material of Ref. 2. We can in principle allow an analogous variation in the pairing amplitudes. However,as shown in the supplementary material of Ref. 2, in the low-energy limit, such a contribution can be eliminated via a gaugetransformation and does not change our results qualitatively.Our central assumption is that the contribution to the hopping due to strain is small compared to any of the parameters of theunperturbed system over the finite length of the sample i.e. max ( δt ( x i )) (cid:28) ∆ , t . As a consequence, we assume that the appliedstrain does not destroy the d-wave superconductor.We assume periodic boundary conditions along the y direction and apply a Fourier transform: c σ ( x i , y i ) = 1 N y (cid:88) k y e ik y y i c σ ( x i , k y ) , (S26)where the 1D Brillouin Zone (BZ) is defined by k y ∈ [ − π/a, π/a ] . The Hamiltonian becomes H T B = (cid:88) i,k y ,σ t ( x i ) (cid:2) c † σ ( x i , k y ) c σ ( x i + a, k y ) + c † σ ( x i + a, k y ) c σ ( x i , k y ) (cid:3) − µc † σ ( x i , k y ) c σ ( x i , k y )+2 t ( x i ) cos ( k y a ) c † σ ( x i , k y ) c σ ( x i , k y ) (S27) H Pair = (cid:88) x i ,k y ,σ,σ (cid:48) ( − i ∆)( iσ y ) σσ (cid:48) sin ( k y a ) [ c σ ( x i , − k y ) c σ (cid:48) ( x i − a, k y ) − c σ ( x i , − k y ) c σ (cid:48) ( x i + a, k y )] + H.c. (S28)Using the ansatz for singlet pairing in Eq. S8, we obtain the BdG Eqs for each H s sector as t ( x i ) [ u s,n ( x i +1 , k ) + u s,n ( x i − , k )] − µu s,n ( x i , k ) + 2 t ( x i ) cos ( ka ) u s,n ( x i , k )+ i ∆ sin ( ka ) [ v s,n ( x i − , k ) − v s,n ( x i +1 , k )] = E n ( k ) u s,n ( x i , k ) − i ∆ sin ( ka ) [ u s,n ( x i +1 , k ) − u s,n ( x i − , k )] − (cid:26) t ( x i ) [ v s,n ( x i +1 , k ) + v s,n ( x i − , k )] − µv s,n ( x i , k ) + 2 t ( x i ) cos ( k y a ) v s,n ( x i , k ) (cid:27) = E n ( k ) v sσ,n ( x i , k ) , (S29)where we dropped the y index. B. Low-energy effective Hamiltonian under applied uniaxial strain
We consider the low-energy limit of the lattice BdG eqs. As we are interested in solutions in the bulk of a sample, we consideran infinite system and ignore the effects of the boundary along the direction of the strain. Our numerical results justify thisapproximation.Note that the lattice and continuum fields are related via the Wannier states [23] c σ ( x i , y i ) = (cid:90) dxdyφ ∗ ( x, y ; x i , y i )Ψ σ ( x, y ) (S30) Ψ σ ( x, y ) = (cid:88) x i ,y j φ ( x, y ; x i , y j ) c σ ( x i , y i ) . (S31)A partial Fourier transform along the y direction gives c σ ( x i , k ) = (cid:90) dx dk (cid:48) π φ ∗ ( x, k (cid:48) ; x i , k )Ψ σ ( x, k (cid:48) ) . (S32)For sufficiently large systems, we allow the Fourier transform of the Wannier function to be sharply peaked φ ∗ ( x, k (cid:48) ; x i , k ) ≈ πδ ( k − k (cid:48) ) ˜ φ ( x ; x i ) . (S33)Using these expression in Eq. S8, we obtain γ s,n ( k ) = (cid:88) σ (cid:90) dx (cid:104) u ∗ sσ,n ( x, k )Ψ σ ( x, k ) − sgn (¯ σ ) v ∗ sσ,n ( x, k )Ψ † ¯ σ ( x, − k ) (cid:105) . (S34)0where we defined the continuum version of the BdG coefficients as u ( α ) , ∗ σ,n ( x, k ) = (cid:88) x i ˜ φ ∗ ( x ; x i ) u ( α ) , ∗ σ,n ( x i , k ) (S35) v ( α ) , ∗ σ,n ( x, k ) = (cid:88) x i ˜ φ ( x ; x i ) v ( α ) , ∗ σ,n ( x i , k ) (S36)We assume for simplicity that the Wannier states can be chosen to be real. We can formally switch to the continuum bymultiplying the each one of eqs. S29 by the Wannier state ˜ φ and summing over all i . We subsequently drop the index i .In the low-energy limit and for momenta in the vicinity of each valley s.t. k ≈ K ( α ) F y + q y , we employ the approximation u s,n ( x, K ( α ) F y + q ) ≈ e iK ( α ) Fx x u ( α ) s,n ( x, q ) (S37) v s,n ( x, K ( α ) F y + q ) ≈ e iK ( α ) Fx x v ( α ) s,n ( x, q ) , (S38)where α ∈ { , (cid:48) , , (cid:48) } is a valley index and ( u, v ) ( α ) are envelope functions which are assumed to vary slowly on the scale ofthe lattice. Such an approximation is well-known in the context of graphene [24]. We introduce a cutoff scale | q | ≤ Λ .Our ansatz allows for a separation of the solutions in terms of valley index. We simplify common phase factors e iK x x andexpand the envelope functions s.t. e iK ( α ) Fx a u ( α ) s,n ( x + a, q ) + e − iK ( α ) Fx a u ( α ) s,n ( x − a, q ) ≈ (cid:16) K ( α ) F x a (cid:17) u ( α ) s,n ( x, q ) + 2 ai sin (cid:16) K ( α ) F x a (cid:17) ∂ x u ( α ) s,n ( x, q ) (S39) e − iK Fx a v ( α ) s,n ( x − a, q ) − e iK Fx a v ( α ) s,n ( x + a, q ) ≈ − i sin (cid:16) K ( α ) F x a (cid:17) v ( α ) s,n ( x, q ) − a cos (cid:16) K ( α ) F x a (cid:17) ∂ x v ( α ) s,n ( x, q ) . (S40)We can further expand cos ( ka ) ≈ cos (cid:16) K ( α ) F y a (cid:17) − sin (cid:16) K ( α ) F y a (cid:17) qa (S41) sin ( ka ) ≈ sin (cid:16) K ( α ) F y a (cid:17) + cos (cid:16) K ( α ) F y a (cid:17) qa. (S42)We keep only terms in Eqs. S29 which are first order in either of the three small quantities δt ( x ) , ∂ x , q . We note that the zerothorder terms vanish either due to presence of a Fermi surface or to the vanishing of the pairing at each valley. To first order, thefirst Eq. of S29 becomes (cid:104) δt ( x ) cos (cid:16) K ( α ) F x a (cid:17) + 2 ta sin (cid:16) K ( α ) F x a (cid:17) i∂ x − ta sin (cid:16) K ( α ) F y a (cid:17) q + 2 δt ( x ) cos (cid:16) K ( α ) F y a (cid:17)(cid:105) u ( α ) s,n ( x, q ) (cid:104) − a sin (cid:16) K ( α ) F y a (cid:17) cos (cid:16) K ( α ) F x a (cid:17) i∂ x + 2∆ a cos (cid:16) K ( α ) F y a (cid:17) sin (cid:16) K ( α ) F x a (cid:17) q (cid:105) v ( α ) s,n ( x, q ) = E n ( k ) u ( α ) s,n ( x i , k ) . (S43)We define v F =2 ta sin ( K F a ) (S44) v ∆ =2∆ a sin ( K F a ) (S45) A ( x ) =2 δt ( x ) [cos ( K F x a ) + cos ( K F y a )] , (S46)where K F = max ( | K ( α ) F x | , | K ( α ) F y | ) is independent of α and positive. The BdG Eqs can be recast as v F (cid:20) sgn (cid:16) K ( α ) F x (cid:17) i∂ x − sgn (cid:16) K ( α ) F y (cid:17) q + Av F (cid:21) u ( α ) s,n ( x i , q ) + v ∆ (cid:104) − sgn (cid:16) K ( α ) F y (cid:17) i∂ x + sgn (cid:16) K ( α ) F x (cid:17) q (cid:105) v ( α ) s,n ( x i , q )= E n ( q ) (cid:16) K ( α ) F y + q (cid:17) u ( α ) s,n ( x, q ) v ∆ (cid:104) − sgn (cid:16) K ( α ) F y (cid:17) i∂ x + sgn (cid:16) K ( α ) F x (cid:17) q (cid:105) u ( α ) s,n ( x, q ) − v F (cid:20) sgn (cid:16) K ( α ) F x (cid:17) i∂ x − sgn (cid:16) K ( α ) F y (cid:17) q + Av F (cid:21) v ( α ) s,n ( x, q )= E n ( q ) (cid:16) K ( α ) F y + q (cid:17) v ( α ) s,n ( x, q ) . (S47)1 C. Landau levels from the low-energy Hamiltonian
We first focus on the (1), (1’) valleys which correspond to K F = (0 , ± K F ) . The strain-induced vector potential induces LLsabout the nodal momenta. By contrast, LLs do not emerge in the vicinity of valleys (2), (2’), as shown further down below.The BdG Eqs. in the vicinity of valley (1) are given by v F (cid:18) − q + Av F (cid:19) u (1) s,n ( x, q ) + v ∆ ( − i∂ x ) v (1) s,n ( x, q ) = E n ( K (1) F y + q ) u (1) s,n ( x, q ) v ∆ ( − i∂ x ) u (1) s,n ( x, q ) − v F (cid:18) − q + Av F (cid:19) v (1) s,n ( x, q ) = E n ( K (1) F y + q ) v (1) s,n ( x, q ) . (S48)We apply a unitary transformation U = 1 √ σ + iσ x ) , (S49)where σ x acts on the space of the BdG spinor. The new Bdg Eqs. for valley (1) read v ∆ ( − i∂ x ) − iv F (cid:16) − q + e ˜ Av F (cid:17) v ∆ ( − i∂ x ) + iv F (cid:16) − q + e ˜ Av F (cid:17) (cid:32) ˜ u (1) s,n ( x, q )˜ v (1) s,n ( x, q ) (cid:33) = E n ( q ) (cid:32) ˜ u (1) s,n ( x, q )˜ v (1) s,n ( x, q ) (cid:33) , (S50)where we introduced the effective vector potential ˜ A = (cid:18) , A ( x ) e (cid:19) . (S51)Recalling that ˜ A ( x ) = Bx by construction , we define a pseudo-magnetic field as B = ∇ × A . (S52)We now clearly see that the BdG equations correspond to Dirac fermions in the presence of a uniform magnetic field in theLandau gauge [24]. The well-known solutions for n (cid:54) = 0 are given by the Landau levels (cid:32) ˜ u (1) s,n ( x i , q )˜ v (1) s,n ( x i , q ) (cid:33) = 1 √ (cid:18) ψ n − ( x, q ) ± ψ n ( x, q ) (cid:19) , (S53)with energies E n = ± ω c √ n. (S54)The cyclotron frequencies and magnetic lengths are defined by ω c = √ v ∆ l B (S55) l B = (cid:114) v ∆ eB . (S56) Ψ n ( x, q ) represent harmonic oscillator wavefunctions [25] ψ n ( ξ ) = (cid:114) l B u n ( ξ ) , (S57) u n = 1 π / √ n n ! e − ξ / H n ( ξ ) , (S58) (cid:88) n u n ( ξ ) u n ( ξ (cid:48) ) = δ ( ξ − ξ (cid:48) ) , (S59) (cid:90) ∞−∞ dξu n ( ξ ) u m ( ξ ) = δ nm . (S60)2where H n are Hermite polynomials and the dimensionless variable is ξ = ( x/l B ) − ( v F /v ∆ ) l B q . The ψ n functions obey thenormalization condition (cid:90) dxψ n ( x, q ) ψ m ( x, q ) = (cid:90) dxλl B u n ( ξ ) u m ( ξ )= (cid:90) dξu n ( ξ ) u m ( ξ )= δ nm . (S61)Referring to the BdG Eqs. (S47), for valley (1’) we have v F (cid:18) q + Av F (cid:19) u (1 (cid:48) ) s,n ( x, q ) + v ∆ ( i∂ x ) v (1 (cid:48) ) s,n ( x, q ) = E n ( q ) u (1 (cid:48) ) s,n ( x, q ) v ∆ ( i∂ x ) u (1 (cid:48) ) s,n ( x, q ) − v F (cid:18) q + Av F (cid:19) v (1 (cid:48) ) s,n ( x, q ) = E n ( q ) v (1 (cid:48) ) s,n ( x, q ) . (S62)These can be obtained from Eqs. S48 for valley (1) via complex conjugation and sending q → − q . This procedure amounts to atime-reversal operation and illustrates the that this symmetry is preserved by the applied strain.In the un-rotated basis, the solutions for the non-zeroth LLs are (cid:32) u ( α ) s,n ( x, q ) v ( α ) s,n ( x, q ) (cid:33) = 12 (cid:18) ψ n − ( x, sgn ( α ) q ) − i sgn ( α ) ψ n ( x, sgn ( α ) q ) − i sgn ( α ) ψ n − ( x, sgn ( α ) q ) + ψ n ( x, sgn ( α ) q ) (cid:19) , (S63)where sgn ( α ) = ± for α ∈ { , (cid:48) } .The zeroth LL in the original basis is determined by (cid:32) u ( α ) s, ( x i , q ) v ( α ) s, ( x i , q ) (cid:33) = ψ ( x, sgn ( α ) q ) √ (cid:18) − i sgn ( α )1 (cid:19) . (S64)The continuum version of the BdG ansatz in Eq. S8 is given by γ ( α ) s,n ( q ) = (cid:90) dx (cid:104) u ( α ) , ∗ sσ,n ( x, q )Ψ ( α ) σ ( x, q ) + v ( α ) , ∗ sσ,n ( x, q )Ψ (¯ α ) , † σ ( x, − q ) (cid:105) . (S65)In the spin-singlet case u sσ = δ sσ and similarly for v sσ = ( iσ y ) sσ v s . The explicit expressions are γ (1) ↑ , ( q ) = (cid:90) dx ψ ( x, q ) √ (cid:104) i Ψ (1) ↑ ( x, q ) + Ψ (1 (cid:48) ) , † ↓ ( x, − q ) (cid:105) (S66) γ (1) ↓ , ( q ) = (cid:90) dx ψ ( x, q ) √ (cid:104) i Ψ (1) ↓ ( x, q ) − Ψ (1 (cid:48) ) , † ↑ ( x i , − q ) (cid:105) (S67) γ (1 (cid:48) ) ↑ , ( q ) = (cid:90) dx ψ ( x, − q ) √ (cid:104) − i Ψ (1 (cid:48) ) ↑ ( x, q ) + Ψ (1) , † ↓ ( x, − q ) (cid:105) (S68) γ (1 (cid:48) ) ↓ , ( q ) = (cid:90) dx ψ ( x, − q ) √ (cid:104) − i Ψ (1 (cid:48) ) ↓ ( x, q ) − Ψ (1) , † ↑ ( x, − q ) (cid:105) . (S69)It follows that γ (1 (cid:48) ) ↑ , ( q ) = iγ (1) , † ↓ , ( − q ) (S70) γ (1 (cid:48) ) ↓ , ( q ) = − iγ (1) , † ↑ , ( − q ) . (S71)We now comment on the solutions at valleys (2), (2’). Based on Eqs. S47, the first of these is subject to the following eqs. v F (cid:18) i∂ x + Av F (cid:19) u (2) s,n ( x, q ) + v ∆ qv (2) σ,n ( x, q ) = E n ( q ) u (2) s,n ( x, q ) v ∆ qu ( α ) s,n ( x, q ) − v F (cid:18) i∂ x + Av F (cid:19) v (2) s,n ( x, q ) = E n ( q ) v (2) s,n ( x, q ) . (S72)3These do not result in the emergence of LLs. The vector potential can be taken into account by incorporating a suitable globalphase: u (2) s,n ( x, q ) = e i (cid:82) xx dx (cid:48) A ( x (cid:48) ) vF ˜ u (2) s,n ( x, q ) (S73) u (2) s,n ( x, q ) = e i (cid:82) xx dx (cid:48) A ( x (cid:48) ) vF ˜ u (2) s,n ( x, q ) . (S74)The solutions correspond to linearly-dispersing, counter-propagating modes centered around K F y = 0 .We finally comment on the topological nature of the zeroth
LLs. To illustrate, we consider the zero-energy solutions ofEqs. S50. By canceling factors of i these eqs. can be reduced to [ v ∆ ( − iσ y ) ∂ x + m ( x ) σ x ] Φ (1)0 ( x, q ) = 0 , (S75)where m ( x ) = v F q − A ( x ) . These are identical to the well-known Jackiw-Rebbi eqs. [26]. The latter host a topologically-protected domain-wall solution centered about the point where m ( x ) changes sign. Going beyond the low-energy approximation,the zeroth LL modes are expected to persist as topologically-protected modes centered around line-defects, in accordance withthe findings of Ref. 10. D. Projection of lattice operators onto the zeroth LL modes
We consider solutions deep in the bulk of a strained sample and thus ignore contributions from any edge states. As before, wecan pass onto the continuum via an expansion using the Wannier states. The continuum version of the inverse BdG transformationis Ψ ( α ) σ ( x, q ) = (cid:88) n,s u ( α ) sσ,n ( x, q ) γ ( α ) s,n ( k ) + v (¯ α ) , ∗ sσ,n ( x, − q ) γ (¯ α ) , † s,n ( − q ) . (S76)Note that for the zeroth LL contribution, only two out of the four states in Eqs. S66-S69 are independent. This necessarilyimplies that the inverse BdG transformation must be restricted to half of the the zeroth LL states. Without loss of generality, wewrite the expansion of the fields as Ψ ( α ) σ ( x, q ) = Ψ ( α ) P,σ + 12 (cid:88) | n |≥ [ ψ n − ( x, sgn ( α ) q ) − i sgn ( α ) ψ n ( x, sgn ( α ) q )] γ ( α ) σ,n ( q ) − sgn ( σ ) [ i sgn (¯ α ) ψ n − ( x, − sgn (¯ α ) q ) + ψ n ( x, − sgn (¯ α ) q )] γ (¯ α ) , † ¯ σ,n ( − q ) , (S77)where we introduced the projection onto the zeroth LL Ψ (1) P, ↑ ( x, q ) = − iψ ( x, q ) √ γ (1) ↑ , ( q ) (S78) Ψ (1) P, ↓ ( x, q ) = − iψ ( x, q ) √ γ (1) ↓ , ( q ) (S79) Ψ (1 (cid:48) ) P, ↑ ( x, q ) = − ψ ( x, − q ) √ γ (1) , † ↓ , ( − q ) (S80) Ψ (1 (cid:48) ) P, ↓ ( x, q ) = ψ ( x, − q ) √ γ (1) , † ↑ , ( − q ) . (S81)We chose to keep both positive and negative values of q , but retained only those contributions from valley (1) . These expressionsare consistent with Eqs. S66-S69. Using the completeness of the Hermite polynomials (Eq. S59) it is straightforward to verifythat the fields in Eq. S77 satisfy canonical anti-commutation relations.Via Eqs. S70-S71, and S76, one can verify that, without any restriction on the Hilbert space of the zeroth LLs, the inverseBdG transformation is indeed singular. The projections onto the zeroth LL are themselves not independent since Ψ (1 (cid:48) ) P,σ ( x, q ) = i sgn ( σ )Ψ (1) , † P, ¯ σ ( x, − q ) . (S82)We now illustrate that Eq. S82 is independent of the choice of Hilbert space for the zeroth LL. Consider retaining both valleyindices, but restricting the momenta q to positive values. The contribution of the zeroth LLs is given by4 Ψ (1) P, ↑ ( x, q ) = ψ ( x, q ) √ (cid:104) − iθ ( q ) γ (1) ↑ , ( q ) − θ ( − q ) γ (1 (cid:48) ) , † ↓ , ( − q ) (cid:105) (S83) Ψ (1) P, ↓ ( x, q ) = ψ ( x, q ) √ (cid:104) − iθ ( q ) γ (1) ↓ , ( q ) + θ ( − q ) γ (1 (cid:48) ) , † ↑ , ( − q ) (cid:105) (S84) Ψ (1 (cid:48) ) P, ↑ ( x, q ) = ψ ( x, − q ) √ (cid:104) iθ ( − q ) γ (1 (cid:48) ) ↑ , ( q ) − θ ( q ) γ (1) , † ↓ , ( − q ) (cid:105) (S85) Ψ (1 (cid:48) ) P, ↓ ( x, q ) = ψ ( x, − q ) √ (cid:104) iθ ( − q ) γ (1 (cid:48) ) ↓ , ( q ) + θ ( q ) γ (1) , † ↑ , ( − q ) (cid:105) . (S86)Using Eqs. S70-S71 one can verify that these satisfy the relations in Eqs. S82.The Fourier transforms of the projected operators, defined by Ψ ( α ) P,σ ( x, y ) = (cid:90) Λ − Λ dq π e iqy Ψ ( α ) P,σ ( x, q ) , (S87)likewise obey Ψ (¯ α ) P,σ ( x, y ) = i sgn ( σ )Ψ ( α ) , † P, ¯ σ ( x, y ) . (S88)In addition, they are subject to the anti-commutation relations { Ψ ( α ) σ , Ψ ( α (cid:48) ) , † σ (cid:48) } = − i sgn ( σ (cid:48) ) { Ψ ( α ) σ , Ψ ( ¯ α (cid:48) )¯ σ (cid:48) } = δ σ,σ (cid:48) (cid:90) Λ − Λ (cid:18) dq π (cid:19) e iq ( y − y (cid:48) ) (cid:26) ψ ( x, q ) ψ ( x (cid:48) , q ) (cid:27) = δ σ,σ (cid:48) l B √ π (cid:90) Λ − Λ (cid:18) dq π (cid:19) e iq ( y − y (cid:48) ) e − (cid:16) xlB − l B λq (cid:17) / e − (cid:16) x (cid:48) lB − l B λq (cid:17) / = δ σ,σ (cid:48) l B √ π (cid:90) Λ − Λ (cid:18) dq π (cid:19) e − l B λ q + λq ( x + x (cid:48) )+ iq ( y − y (cid:48) ) e − x l B e − x (cid:48) l B ≈ δ σ,σ (cid:48) λl B √ π (cid:114) πl B π e [ ( x + x (cid:48) )+ i ( y − y (cid:48) ) ] l Bλ e − x l B e − x (cid:48) l B ≈ δ σ,σ (cid:48) πλl B e − ( x − x (cid:48) )24 λ l B e − ( y − y (cid:48) )24 λ l B e i ( x + x (cid:48) )( y − y (cid:48) )2 λl B (S89) { Ψ ( α ) σ , Ψ ( α ) σ (cid:48) } = 0 . (S90)where λ = v F /v ∆ . In the fourth line we approximated the integral by a Gaussian in the limit Λ l B → ∞ . This approximation isjustified for weak pseudo-magnetic fields considered here. We identify the coefficient as half of the density of LLs n = 12 πλl B , (S91)by using the well-known formula for the degeneracy of LLs N LL = L x L y πl B , (S92)where L x , L y are the lengths of the system in either direction and we accounted for the anisotropy between the Fermi velocitiesvia the factor λ .The field operator projected onto the zeroth LLs gives Ψ P,σ ( x, y ) = e iK F y Ψ (1) P,σ ( x, y ) + e − iK F y Ψ (1 (cid:48) ) P,σ ( x, y ) . (S93)5 E. Ferromagnetism at mean-field level
In order to investigate many body instabilities of the zeroth LL, we consider a repulsive, short-range Hubbard interaction inthe continuum limit H U = U (cid:90) d r Ψ † ↑ ( r )Ψ ↑ ( r )Ψ † ↓ ( r )Ψ ↓ ( r ) . (S94)Using Eq. S93, we project this term onto the zeroth LLs: H U,P = U (cid:90) d r (cid:104) ρ (1) P, ↑ ρ (1) P, ↓ + ρ (1) P, ↑ ρ (1 (cid:48) ) P, ↓ + ρ (1 (cid:48) ) P, ↑ ρ (1) P, ↓ + ρ (1 (cid:48) ) P, ↑ ρ (1 (cid:48) ) P, ↓ + Ψ (1) , † P, ↑ Ψ (1 (cid:48) ) P, ↑ Ψ (1 (cid:48) ) , † P, ↓ Ψ (1) P, ↓ + Ψ (1 (cid:48) ) , † P, ↑ Ψ (1) P, ↑ Ψ (1) , † P, ↓ Ψ (1 (cid:48) ) P, ↓ (cid:105) . (S95)Note that we have only kept terms which lack rapidly oscillating phase factors such as e ± iK F x . We anticipate that in thelong-wavelength limit considered here, the terms we kept are the most relevant. We have also ignored any contributions fromboundary terms.Using Eq. S88, we can formally eliminate the degrees of freedom at valley (1 (cid:48) ) in favor of their analogoues at valley (1) : ρ (1) P ↑ ρ (1 (cid:48) ) P ↓ → Ψ (1) , † P ↑ Ψ (1) P ↑ Ψ (1) P ↑ Ψ (1) , † P ↑ (S96) ρ (1 (cid:48) ) P ↑ ρ (1) P ↓ → Ψ (1) P ↓ Ψ (1) , † P ↓ Ψ (1) , † P ↓ Ψ (1) P ↓ (S97) ρ (1 (cid:48) ) P ↑ ρ (1 (cid:48) ) P ↓ → Ψ (1) P ↓ Ψ (1) , † P ↓ Ψ (1) P ↑ Ψ (1) , † P ↑ (S98) Ψ (1) , † ↑ Ψ (1 (cid:48) ) ↑ Ψ (1 (cid:48) ) , † ↓ Ψ (1) ↓ → − Ψ (1) , † P ↑ Ψ (1) , † P ↓ Ψ (1) P ↑ Ψ (1) P ↓ (S99) Ψ (1 (cid:48) ) , † P ↑ Ψ (1) P ↑ Ψ (1) , † P ↓ Ψ (1 (cid:48) ) P ↓ → − Ψ (1) P ↓ Ψ (1) P ↑ Ψ (1) , † P ↓ Ψ (1) , † P ↑ . (S100)We take advantage of the fact that the anti-commutators of the projected operators is finite and proportional to half of thedensity of zeroth LLs per spin (Eq. S89). This way the Hamiltonian reduces to H U,P = U (cid:90) d r (cid:20) ρ (1 P ↑ ρ (1) P ↓ − n ( ρ (1) P ↑ + ρ (1) P ↓ ) + 12 n (cid:21) . (S101)(S102)We write ρ (1) ↑ = 12 (cid:104) S (1) z + ρ (1) (cid:105) (S103) ρ (1) ↓ = − (cid:104) S (1) z − ρ (1) (cid:105) , (S104)where S (1) z = 12 (cid:16) ρ (1) ↑ − ρ (1) ↓ (cid:17) (S105) ρ (1) = ρ (1) ↑ + ρ (1) ↓ . (S106)The Hamiltonian can be re-written as H U,P = U (cid:90) d r (cid:20) − (cid:16) S (1) P,z (cid:17) + (cid:16) ρ (1) P − n (cid:17) − n (cid:21) (S107)The Hamiltonian can be decoupled at mean-field level in the ferromagnetic channel. For repulsive interactions, the pairinginstability is energetically unfavorable. A possible competing spin-density order is not allowed, since we have eliminated thedegrees of freedom at valley (1 (cid:48) ) .Taking the local spin-density and total density as variational parameters, we take the expectation value of H U,P . From theexpression above, it is clear that a half-filled ground state corresponding to ρ (1) P = n / is energetically favored.Using S87, the spin-part is given by6 (cid:104) H U,P,
Spin (cid:105) = − (cid:32) √ πUλl B (cid:33) (cid:90) (cid:18) dk π (cid:19) (cid:18) dq π (cid:19) (cid:18) dQ π (cid:19) e − l B (cid:26) [ Q +( k − q )]2+ Q (cid:27) (cid:104) γ (1) , † k + Q ↑ γ (1) k ↑ − γ (1) , † k + Q ↓ γ (1) k ↓ (cid:105) (cid:104) γ (1) , † q − Q ↑ γ (1) q ↑ − γ (1) , † q − Q ↓ γ (1) q ↓ (cid:105) . (S108)The Gaussian term ensures that the dominant contributions are due to momenta k, q, Q which are within all within a scale set by /l B . At half-filling, it is clear that this term is maximally negative at ferromagnetic alignment.In conclusion, a ferromagnetic ground-state is favored at mean-field level. The resulting Stoner spectrum consists of flat bandsat E q ≈ (cid:32) √ πUλl B (cid:33) π ) (cid:90) ∞−∞ dke − l B ( k − q )22 = ± n U . (S109) III. TWO DIMENSIONAL D-WAVE PAIRING WITH STRAIN OFF THE NODAL AXES
In this section we discuss a 2D d-wave SC with strain applied off the nodal axes of the pairing. In Sec. III A, we introduce theeffective 1D lattice Hamiltonians under uniaxial strain. We also discuss their enhanced symmetry and sublattice structure, whichare important differences w.r.t. the models of Sec. II A. In Sec. III B, we introduce the low-energy, continuum limit of the BdGeqs. and obtain their solutions. We find two degenerate zeroth LLs in the vicinity of each of the two valleys with BdG coefficientswhich are subject to the constraint imposed by the combined mirror and p-h symmetry discussed in the main text and below.Sec. III C presents a comparison of the analytical solutions with the numerical results. Sec. III D is devoted to a discussion of theprojection of the lattice operators onto the zeroth LL section, while Sec. III E illustrates the ensuing ferromagnetic instability.
A. 1D Lattice Hamiltonian
We consider the tight-binding Hamiltonian in Eq. S23 without strain and rotate by π/ . In order to preserve translationalsymmetry along the rotated y direction at the edges we choose a two-site unit cell. This lattice is illustrated in Fig. S6.We subsequently introduce a spatially-varying tight-binding parameter t ( x i ) , where x is along the diagonal of the un-rotatedsystem. Applying a Fourier transform along the y direction we obtain H = H T B + H ∆ (S110) H T B = N x (cid:88) i =1 ,σ t ( x i ) cos (cid:18) k y a √ (cid:19) c † A,σ ( x i , k y ) c B,σ ( x i , k y ) + 2 t (cid:18) x i + a √ (cid:19) cos (cid:18) k y a √ (cid:19) c † B,σ ( x i , k y ) c A,σ ( x i +1 , k y ) + H.c. − µc † A,σ ( x i , k y ) c A,σ ( x i , k y ) − µc † B,σ ( x i , k y ) c B,σ ( x i , k y ) H ∆ = N x (cid:88) i =1 (cid:88) σσ (cid:48) ( iσ y ) σσ (cid:48) ∆ (cid:2) c † σ ( A, x i , k y ) c † A,σ (cid:48) ( x i +1 , − k y ) + c † A,σ ( x i +1 , k y ) c † A,σ (cid:48) ( x i , − k y ) − (cid:16) k y a √ (cid:17) c † A,σ ( x i , k y ) c A,σ (cid:48) ( x i , − k y ) (cid:3) + ( A ↔ B ) + H.c. (S111)The position index runs over all even values while the momentum is along a folded 1D BZ. Without strain, the effective chainmodel is periodic by translation via a √ or i → i + 1 . Subsequently, we distinguish between two inequivalent sublattices whichwe label by A, B .We consider the ansatz in Eq. S8 modified to take into account the sublattice structure: γ s,n ( k ) = (cid:88) x i ,σ (cid:26) u ∗ sσ,A,n ( x i , k ) c A,σ ( x i , k ) + v ∗ sσ,A,n ( x i , k ) c † A,σ ( x i , − k )+ u ∗ sσ,B,n ( x i , k ) c B,σ ( x i , k ) + v ∗ sσ,B,n ( x i , k ) c † B,σ ( x i , − k ) (cid:27) . (S112)Using the convention in Eq. S4 where u sσ = δ sσ u s and v sσ = ( iσ y ) sσ v s , the singlet pairing Hamiltonian becomes independentof the spin indices σ , which we subsequently drop for simplicity.7 A B x y FIG. S6. Lattice model for 2D-S with strain along the diagonal of the BZ. The x, y directions are rotated wrt the axes of the BZ. In order topreserve translational symmetry in the y direction along the edges, we use a two-site unit cell, illustrated by the dashed line. The strain is alongthe x direction. The BdG eqs. for sublattice B are given by (cid:18) k y a √ (cid:19) (cid:104) t ( x i ) u s,A,n ( x i , k y ) + t ( x i + a/ √ u s,A,n ( x i +1 , k y ) (cid:105) − µu s,B,n ( x i , k y )+∆ [ v s,B,n ( x i +1 , k y ) + v s,B,n ( x i − , k y )] −
2∆ cos (cid:16) k y a √ (cid:17) v s,B,n ( x i , k y ) = E n ( k y ) u s,B,n ( x i , k y ) (S113) ∆ [ u s,B,n ( x i +1 , k y ) + u s,B,n ( x i − , k y )] −
2∆ cos (cid:16) k y a √ (cid:17) u s,B,n ( x i , k y ) − (cid:26) (cid:18) k y a √ (cid:19) (cid:104) t ( x i ) v s,A,n ( x i , k y ) + t ( x i + a/ √ v s,A,n ( x i +1 , k y ) (cid:105) − µv s,B,n ( x i , k y ) (cid:27) = E n ( k y ) v s,B,n ( x i , k y ) . (S114)The corresponding eqs. for sublattice A can be obtained by simply exchanging A, B indices.Finally, we note that the BdG equations for each sublattice are invariant under k → − k , as discussed in the main text.8 B. Landau levels in the continuum limit
In the limit of zero strain, the bonding solution of Eqs. S113-S114 exhibits nodes at four momenta ( ± K F , ± K F ) , which welabel by valley indices (1) , (2) , (1 (cid:48) ) , (2 (cid:48) ) in a clockwise direction. The anti-bonding solution is generically gapped and will beignored in the following.We extend the ansatz for the pristine case to the strained system by writing the wavefunctions as u σ,A,n ( x i )( K F y + q y ) = e iK ( α ) Fx x i u ( α ) σ,A,n ( x i , q y ) + e iK ¯( α ) Fx x i u (¯ α ) σ,A,n ( x i , q y ) (S115) u σ,B,n ( x i )( K F y + q y ) = e iK ( α ) Fx ( x i + a/ √ u ( α ) σ,A,n ( x i + a/ √ , q y ) + e iK ¯( α ) Fx ( x i + a/ √ u (¯ α ) σ,A,n ( x i + a/ √ , q y ) , (S116)in the vicinity of k y ≈ K F y + q y with q y small. A similar expansion holds for v A/B . Note that we are considering bondingsolutions where the envelope functions u ( α ) A = u ( α ) B , while keeping track of the phase difference between the two sublattices.Explicitly, for positive K F y , α = 1 and ¯ α = 2 (cid:48) , and K ( α ) = K F x while K ( ¯ α ) = − K F x . A similar expansion holds for valleys(2), and (1’) with opposite signs for both K F x and K F y . This form for the ansatz is clearly justified for zero strain, where thepairs of valleys (1) and (2 (cid:48) ) and (1 (cid:48) ) and (2) are decoupled. Any amount of small strain breaks translational symmetry along the x direction and can in principle mix the valleys in pairs. The ansatz continues to hold, as indicated by numerical results whichdo not assume this decomposition (see main text).Following the procedure in Sec. II B, we formally pass to the continuum by expanding in terms of Wannier states. We alsodrop the i , y and the sublattice index. Our ansatz becomes γ ( α ) s,n ( q ) = (cid:88) σ (cid:90) dxe − iK ( α ) Fx x (cid:26) (cid:104) u ( α ) , ∗ sσ,n ( x, q )Ψ ( I ) Aσ ( x, q ) + v ( α ) , ∗ sσ,n ( x, q )Ψ ( II ) , † Aσ ( x, − q ) (cid:105) + e − iK ( α ) Fx a/ √ (cid:104) u ( α ) , ∗ sσ,n ( x + a/ √ , q )Ψ ( I ) B,σ ( x, q ) + v ( α ) , ∗ sσ,n ( x + a/ √ , q )Ψ ( II ) , † Bσ ( x, − q ) (cid:105) (cid:27) , α ∈ { , (cid:48) } . (S117)We introduced new field operators Ψ ( I/II ) , † σ ( x, q ) in the vicinity of k y ≈ ± K F y + q . For valleys (cid:48) , we must interchange I, II indices. This decomposition is equivalent to γ s,n (cid:16) K ( α ) F y + q (cid:17) = γ ( α ) s,n ( q ) + γ (¯ α ) s,n ( q ) , (S118)With this ansatz, the BdG equations can be separated into valley sectors.To illustrate the general solution, we consider valley (1) in particular. We expand the envelope functions in real space and thetrigonometric functions in momentum space about the nodal points. We furthermore approximate t ( x i + a/ √ by t ( x i ) sincethe effect of the strain is assumed to vary slowly on the scale of the lattice. We keep only terms which are first order in the smallterms δt ( x ) , ∂ x , q . To zeroth order we obtain: (cid:20) t cos (cid:18) K F x a √ (cid:19) cos (cid:18) K F y a √ (cid:19)(cid:21) u (1) s,n ( x, q ) − µu (1) s,n ( x, q ) = 0 . (S119)To first order, we obtain (cid:20) δt ( x ) cos (cid:18) K F x a √ (cid:19) cos (cid:18) K F y a √ (cid:19) − ta √ (cid:18) K F x a √ (cid:19) sin (cid:18) K F y a √ (cid:19) q + 2 ta √ (cid:18) K F y a √ (cid:19) sin (cid:18) K F x a √ (cid:19) i∂ x (cid:21) u (1) s,n ( x, q )+ (cid:104) √ a sin (cid:16) K F x a √ (cid:17) i∂ x + 2∆ a √ (cid:16) K F y a √ (cid:17) q (cid:105) v (1) s,n ( x, q ) = E n ( q ) u (1) s,n ( x, q ) (S120)We define v F =2 ta √ (cid:18) K F a √ (cid:19) sin (cid:18) K F a √ (cid:19) (S121) = ta √ (cid:16) K F a √ (cid:17) (S122) v ∆ =2∆ a √ (cid:16) K F a √ (cid:17) , (S123) A =4 t cos (cid:18) K F a √ (cid:19) , (S124) A ( x ) =4 δt ( x ) cos (cid:18) K F a √ (cid:19) (S125)9where K F = | K F x | = | K F y | ≥ .Thd BdG eqs. simplify to v F (cid:20) i∂ x − q + Av F (cid:21) u (1) s,n ( x, q ) + v ∆ [ i∂ x + q ] v (1) s,n ( x, q ) = E n ( q ) u (1) s,n ( x, q ) v ∆ [ i∂ x + q ] u (1) s,n ( x, q ) − v F (cid:20) i∂ x − q + Av F (cid:21) v (1) s,n ( x, q ) = E n ( q ) v (1) s,n ( x, q ) . (S126)We now focus on the zero-energy solutions. We transform τ z → τ y and obtain (cid:26) v ∆ [ i∂ x + q ] − iv F (cid:20) i∂ x − q + A v F (cid:21)(cid:27) ˜ v (1) s,a, ( x, q ) = 0 (S127) (cid:26) v ∆ [ i∂ x + q ] + iv F (cid:20) i∂ x − q + A v F (cid:21)(cid:27) ˜ u (1) s,a, ( x, q ) =0 . (S128)These equations are very similar to those for the zeroth LL obtained in Sec. II. By introducing the quantity v F ± iv ∆ = ve ± iθ , (S129)we can write (cid:20) ∂ x + ie − iθ q − iv F e − iθ A ( x ) v (cid:21) ˜ v s, ( x, q ) =0 (S130) (cid:20) ∂ x + ie iθ q − iv F e iθ A ( x ) v (cid:21) ˜ u s, ( x, q ) =0 . (S131)The solutions are given by ˜ v s, ( x, q ) = C v exp (cid:104) sin(2 θ ) q − sin( θ ) v F A ( x ) v (cid:105) v F v ∆ sin( θ ) ∂ x A ( x ) exp i (cid:104) cos(2 θ ) q − cos( θ ) v F A ( x ) v (cid:105) v F v ∆ cos( θ ) ∂ x A ( x ) (S132) ˜ u s, ( x, q ) = C u exp − (cid:104) sin(2 θ ) q − sin( θ ) v F A ( x ) v (cid:105) v F v ∆ sin( θ ) ∂ x A ( x ) exp i (cid:104) cos(2 θ ) q − cos( θ ) v F A ( x ) v (cid:105) v F v ∆ cos( θ ) ∂ x A ( x ) (S133)We recall that A ( x ) ∼ x by construction. Therefore, only one solution is normalizable since it contains a decaying Gaussianform. This is typical of the zeroth LL wavefunctions for Dirac systems in a Landau gauge.Without loss of generality, we consider the case where ˜ u is the allowed non-trivial solution: (cid:32) u (1) s, v (1) s, (cid:33) = 1 √ u ( x, q ) (cid:18) − i (cid:19) . (S134)The equations for valley (2’) can be obtained from Eqs. S126 by changing the sign of ∂ x . The analogues of Eqs. S127-S128involve complex conjugation and interchanging ˜ u and ˜ v : (cid:32) u (2 (cid:48) ) s, ( x, q ) v (2 (cid:48) ) s, ( x, q ) (cid:33) = e iφ (2 (cid:48) ) s (cid:32) − v (1) , ∗ s, ( x, q ) u (1) , ∗ s, ( x, q ) (cid:33) . (S135)We introduced a global phase term, which must be chosen to ensure consistency with other symmetries. We can readily checkthat these are consistent with the mirror symmetry given by Eq. S20. Explicitly, we have (cid:32) u (2 (cid:48) ) s, ( x, q ) v (2 (cid:48) ) s, ( x, q ) (cid:33) = sgn ( s ) √ u ∗ ( x, q ) (cid:18) − i (cid:19) , (S136)where we fixed the global phase.0The BdG eqs. at valleys (1’) can be obtained from Eqs. S126 by changing the signs of both ∂ x and q . The solutions are givenby complex-conjugating, interchanging ˜ u and ˜ v and changing the sign of q . These steps are nothing but the p-h transformationin Eq. S12. The explicit solutions read (cid:32) u (1 (cid:48) ) s, ( x, q ) v (1 (cid:48) ) s, ( x, q ) (cid:33) = i √ u ∗ ( x, − q ) (cid:18) − i (cid:19) , (S137)and (cid:32) u (2) s, ( x, q ) v (2) s, ( x, q ) (cid:33) = − i sgn ( s ) √ u ( x, − q ) (cid:18) − i (cid:19) . (S138)Note that we allowed for additional global phases i to ensure time-reversal symmetric solutions at opposite valleys. This doesnot violate the p-h correspondence since solutions are determined modulo a global phase.To sum up the discussion thus far, we obtain two independent solutions for k y ≈ K F y + q y , where K F y is positive: γ (1) s, ( q ) = (cid:88) σ (cid:90) dxe − iK Fx x (cid:26) ˜ u ∗ ( x, q ) (cid:104) δ sσ Ψ ( I ) Aσ ( x, q ) + i ( iσ y ) sσ Ψ ( II ) , † Aσ ( x, − q ) (cid:105) + e − iK Fx a/ √ ˜ u ∗ ( x + a/ √ , q ) (cid:104) δ sσ Ψ ( I ) B,σ ( x, q ) + i ( iσ y ) sσ Ψ ( II ) , † Bσ ( x, − q ) (cid:105) (cid:27) (S139) γ (2 (cid:48) ) s,n ( q ) = sgn ( s ) (cid:88) σ (cid:90) dxe iK Fx x (cid:26) ˜ u ( x, q ) (cid:104) iδ sσ Ψ ( I ) Aσ ( x, q ) + ( iσ y ) sσ Ψ ( II ) , † Aσ ( x, − q ) (cid:105) + e iK Fx a/ √ ˜ u ( x + a/ √ , q ) (cid:104) iδ sσ Ψ ( I ) B,σ ( x, q ) + ( iσ y ) sσ Ψ ( II ) , † Bσ ( x, − q ) (cid:105) (cid:27) (S140) γ (1 (cid:48) ) s, ( q ) = (cid:88) σ (cid:90) dxe iK Fx x (cid:26) ˜ u ( x, − q ) (cid:104) δ sσ Ψ ( II ) Aσ ( x, q ) − i ( iσ y ) sσ Ψ ( I ) , † Aσ ( x, − q ) (cid:105) + e iK Fx a/ √ ˜ u ( x + a/ √ , − q ) (cid:104) δ sσ Ψ ( II ) B,σ ( x, q ) − i ( iσ y ) sσ Ψ ( I ) , † Bσ ( x, − q ) (cid:105) (cid:27) (S141) γ (2) s, ( q ) = sgn ( s ) (cid:88) σ (cid:90) dxe − iK Fx x (cid:26) ˜ u ∗ ( x, − q ) (cid:104) − iδ sσ Ψ ( II ) Aσ ( x, q ) + ( iσ y ) sσ Ψ ( I ) , † Aσ ( x, − q ) (cid:105) + e − iK Fx a/ √ ˜ u ∗ ( x + a/ √ , − q ) (cid:104) − iδ sσ Ψ ( II ) B,σ ( x, q ) + ( iσ y ) sσ Ψ ( I ) , † Bσ ( x, − q ) (cid:105) (cid:27) (S142)These satisfy γ (1 (cid:48) ) s (cid:48) , ( q ) = (cid:88) s ( σ y ) s (cid:48) s γ (1) , † s, ( − q ) (S143) γ (2) s (cid:48) , ( q ) = (cid:88) s − ( σ y ) s (cid:48) s γ (2 (cid:48) ) , † s, ( − q ) . (S144)These imply that at each valley there are two independent solutions. However, solutions at opposite valleys are not independent,in analogy with the results of Sec. II C. C. Comparison to numerical solution of the lattice model
In order to compare with the results of the numerical calculation presented in Fig. 4 of the main text, we consider the linearcombinations of the analytical solutions found in the previous section1 γ ( I ) s, , ( q ) = γ (1) s, ( q ) − i sgn ( s ) γ (2 (cid:48) ) s, ( q ) (S145) γ ( I ) , † s, , ( q ) = − iγ (1) s, ( q ) + sgn ( s ) γ (2 (cid:48) ) s, ( q ) , (S146)(S147)which imply γ ( I )1 ,s, ( q ) = (cid:90) dx (cid:104) u ( I )1 ,A,sσ Ψ ( I ) Aσ ( x, q ) + v ( I )1 ,A,sσ Ψ ( II ) † Aσ ( x, − q ) (cid:105) + (cid:104) u ( I )1 ,B,sσ Ψ ( I ) Bσ ( x, q ) + v ( I )1 ,B,sσ Ψ ( II ) † Bσ ( x, − q ) (cid:105) (S148) γ ( I )2 ,s, ( q ) = (cid:90) dx (cid:104) u ( I )2 ,A,sσ Ψ ( I ) Aσ ( x, q ) + v ( I )2 ,A,sσ Ψ ( II ) † Aσ ( x, − q ) (cid:105) + (cid:104) u ( I )2 ,B,sσ Ψ ( I ) Bσ ( x, q ) + v ( I )2 ,B,sσ Ψ ( II ) † Bσ ( x, − q ) (cid:105) , (S149)where u ( I )1 ,A,sσ =2 Re (cid:2) e iKx ˜ u ( x, q ) (cid:3) δ sσ (S150) v ( I )1 ,A,sσ =2 Im (cid:2) e iKx ˜ u ( x, q ) (cid:3) ( iσ y ) sσ (S151) u ( I )1 ,B,sσ =2 Re (cid:104) e iK ( x + a/ √ ˜ u ( x, q ) (cid:105) δ sσ (S152) v ( I )1 ,B,sσ =2 Im (cid:104) e iK ( x + a/ √ ˜ u ( x, q ) (cid:105) ( iσ y ) sσ (S153) u ( I )2 ,A,sσ = − Im (cid:2) e iKx ˜ u ( x, q ) (cid:3) δ sσ (S154) v ( I )2 ,A,sσ =2 Re (cid:2) e iKx ˜ u ( x, q ) (cid:3) ( iσ y ) sσ (S155) u ( I )2 ,B,sσ = − Im (cid:104) e iK ( x + a/ √ ˜ u ( x, q ) (cid:105) δ sσ (S156) v ( I )2 ,B,sσ =2 Re (cid:104) e iK ( x + a/ √ ˜ u ( x, q ) (cid:105) ( iσ y ) sσ . (S157)Similar linear combinations can be taken wrt valley (II) operators. These can be obtained via a time-reversal operation (cid:88) s (cid:48) ( − iσ s ) ss (cid:48) γ ( II )1 ,s (cid:48) , ( − q ) = θγ ( I )1 ,s, ( q ) , (S158)which gives (cid:32) u ( II )1 / ,A/B,s ( q ) v ( II )1 / ,A/B,s ( q ) (cid:33) = ( − iσ y ) ss (cid:48) (cid:32) u ( I )1 / ,A/B,s (cid:48) ( − q ) v ( I )1 / ,A/B,s (cid:48) ( − q ) , (cid:33) (S159)where u ( α ) sσ = δ sσ u ( α ) s and v ( α ) sσ = ( iσ y ) sσ v ( α ) s . Hence, the BdG coefficients can be made real.One can see that these obey the constraints imposed by the mirror symmetry discussed in the main text. Furthermore, theseanalytical solutions are consistent with the numerical results of Fig. 4 of the main text. D. Projection onto zeroth LLs
The inverse BdG eqs. can be obtained via generalizing Eqs. S76 to the case of two suballtices: Ψ ( I ) A,σ ( x, q ) = (cid:88) n,s e iK ( α ) Fx x u ( α ) A,sσ,n ( x, q ) γ ( α ) s,n ( q ) + e − iK (¯ α ) Fx x v (¯ α ) , ∗ A,sσ,n ( x, − q ) γ (¯ α ) , † s,n ( − q ) , α ∈ { , (cid:48) } (S160) Ψ ( I ) B,σ ( x, q ) = (cid:88) n,s e iK ( α ) Fx x u ( α ) B,sσ,n ( x, q ) γ ( α ) s,n ( q ) + e − iK ( α ) Fx x v (¯ α ) , ∗ B,sσ,n ( x, − q ) γ (¯ α ) , † s,n ( − q ) , α ∈ { , (cid:48) } (S161)(S162) Ψ ( II ) A,σ ( x, q ) = (cid:88) n,s e iK ( α ) Fx x u ( α ) A,sσ,n ( x, q ) γ ( α ) s,n ( q ) + e − iK (¯ α ) xFx v (¯ α ) , ∗ A,sσ,n ( x, − q ) γ (¯ α ) , † s,n ( − q ) , α ∈ { (cid:48) , } (S163) Ψ ( II ) B,σ ( x, q ) = (cid:88) n,s e iK ( α ) Fx x u ( α ) B,sσ,n ( x, q ) γ ( α ) s,n ( q ) + e − iK ( α ) Fx x v (¯ α ) , ∗ B,sσ,n ( x, − q ) γ (¯ α ) , † s,n ( − q ) , α ∈ { (cid:48) , } . (S164)2As discussed in Sec. II D, we truncate the Hilbert space by considering only (1) and (2 (cid:48) ) solutions. This ensures that the BdGtransformation on the zeroth LL sector is well-defined are well-defined. The projected operators are given by Ψ ( I ) P,A ↑ ( x, q ) = e iK F x u (1) A, ↑ , ( x, q ) γ (1) ↑ , ( q ) + e − iK F x u (2 (cid:48) ) A, ↑ , ( x, q ) γ (2 (cid:48) ) ↑ , ( q ) (S165) Ψ ( I ) P,A ↓ ( x, q ) = e iK F x u (1) A, ↓ , ( x, q ) γ (1) ↓ , ( q ) + e − iK F x u (2 (cid:48) ) A, ↓ , ( x, q ) γ (2 (cid:48) ) ↓ , ( q ) (S166) Ψ ( II ) P,A ↑ ( x, q ) = − e − iK F x v (1) , ∗ A, ↑ , ( x, − q ) γ (1) , † ↓ , ( − q ) − e iK F x v (2 (cid:48) ) , ∗ A, ↑ , ( x, − q ) γ (2 (cid:48) ) , † ↓ , ( − q ) (S167) Ψ ( II ) P,A ↓ ( x, q ) = e − iK F x v (1) , ∗ A, ↓ , ( x, − q ) γ (1) , † ↑ , ( − q ) + e iK F x v (2 (cid:48) ) , ∗ A, ↓ , ( x, − q ) γ (2 (cid:48) ) ↑ , ( − q ) , (S168)with similar expressions for B . Using Eqs. S134, S136, it is straightforward to verify that Ψ ( II ) P,Aσ ( x, q ) =( σ y ) σσ (cid:48) Ψ ( I ) , † P,Aσ ( x, − q ) . (S169)Similar relations hold for sublattice B . These expressions are the analogues of Eq. S82.Following the procedure set out in Sec. II D, we consider the real space fields on each sublattice: Ψ ( I ) A/B,σ ( x, y ) = (cid:90) (cid:18) dq π (cid:19) e iqy Ψ ( I ) A/Bσ ( x, q ) (S170) Ψ ( II ) A/B,σ ( x, y ) = (cid:90) (cid:18) dq π (cid:19) e iqy Ψ ( II ) A/Bσ ( x, q ) . (S171)Finally, the real-space fields are determined by Ψ A/Bσ ( x, y ) = e iK Fy a Ψ ( I ) A/B,σ ( x, y ) + e − iK Fy a Ψ ( II ) A/B,σ ( x, y ) . (S172)Note that the overall phases for momenta along the y direction. E. Ferromagnetism at mean-field level
As for the case of strain along a nodal axis, we allow small but finite residual Hubbard interactions. In contrast to the previouscase, there are now two independent states per spin at each valley, related to each other by Eqs. S143, S144. As discussed inSec. III A the real-space system is described by a two-site unit cell. The two independent zero-energy LLs per spin at each valleycan be used to determine the projections of the two fields Ψ P,Aσ and Ψ P,B,σ defined on each sublattice, as shown in the previoussection.Consider a continuum version of the Hubbard interaction on each sublattice: H U,P = U (cid:48) (cid:90) d r (cid:104) Ψ † A ↑ ( r )Ψ A ↑ ( r )Ψ † A ↓ ( r )Ψ A ↓ ( r ) + Ψ † B ↑ ( r )Ψ B ↑ ( r )Ψ † B ↓ ( r )Ψ B ↓ ( r ) (cid:105) . (S173)Using Eqs. S169, each sublattice sector maps onto the single-valley case studied in Sec. II E. Hence, each sublattice is expectedto order ferromagnetically at mean-field level.The large degeneracy due to possible relative orientations of the resulting sublattice moments can be lifted by arbitrarily smallinter-sublattice interactions. In view of expected short range, repulsive nature of the residual interactions in a SC, and in theabsence of negligible dispersion, we conclude that the most likely instability is still ferromagnetic. IV. HELICAL MAJORANAS FOR 3D EQUAL-SPIN TRIPLET SCS
In this Section we consider the emergence of helical Majorana modes in 3D equal-spin triplet superconductors under uniaxialstrain. Sec. IV A introduces a family of effective 1D lattice models. In Sec. IV B, we obtain helical Majorana states in thecontinuum limit. In Sec. IV C we show that the Fermi fields projected onto the zeroth LLs are emergent real-space Majoranafermions.3
A. Lattice Hamiltonian
We consider a 3D lattice model with NN hopping along either x, y, z directions and strain along the x directions. The pairingis chosen to belong to the Γ − representation of the tetragonal D h point group [13]. Without strain, this can be written as H = (cid:88) k c † k ↑ c † k ↓ c − k ↓ − c − k ↑ T h k − ∆ k h k ∆ (cid:48) k
00 ∆ (cid:48) , ∗ k − h k − ∆ ∗ k − h k c k ↑ c k ↓ c † − k ↓ − c † − k ↑ , (S174)where h k =2 t (cos( k x a ) + cos( k y a ) + cos( k z a )) − µ (S175) ∆ k =∆ ( − sin( k x a ) + i sin( k y a )) (S176) ∆ (cid:48) k =∆ (sin( k x a ) + i sin( k y a )) . (S177)The spectrum is given by E k = ± (cid:113) h k + ∆ (sin( k x a ) + sin( k y a ) ) , (S178)with Dirac point nodes at K / (cid:48) = (0 , , ± K F ) .The lattice Hamiltonian can be chosen as H = H T B + H Pair (S179)where H T B = (cid:88) R i (cid:88) δ j (cid:88) σ t ( x i ) c † σ ( R i ) c σ ( R i + δ j ) + H.c. − µc † σ ( R i ) c σ ( R i ) . (S180)where δ j = x,y,z equals ( a, , , (0 , a, , (0 , , a ) , respectively. The pairing part is given by H Pair = ∆2 (cid:88) x i ,y j ,z l (cid:88) σ (cid:26) i sgn ( σ ) (cid:2) c † σ ( x i , y j , z l ) c † σ ( x i +1 , y j , z l ) − c † σ ( x i , y j , z l ) c † σ ( x i − , y j , z l ) (cid:3) + (cid:2) c † σ ( x i , y j , z l ) c † σ ( x i , y j +1 , z l ) − c † σ ( x i , y j , z l ) c † σ ( x i , y j − , z l ) (cid:3) (cid:27) + H.c. (S181)Note the difference in sign due to the convention in Eq. S174.We apply Fourier transforms along the y and z directions. As in the pristine case, the Hamiltonian, and consequently the BdGequations separate into two independent sectors: t ( x i ) [ u s,n ( x i +1 , k ) + u s,n ( x i − , k )] + { t ( x i ) [cos( k y a ) + cos( k z a )] − µ } u s,n ( x i , k ) − i ∆2 [( v s,n ( x i +1 , k ) − v s,n ( x i − , k )) + 2 sgn ( σ ) sin( k y a ) v s,n ( x i , k )] = E n ( k ) u s,n ( x i , k ) (S182) i ∆2 [( u s,n ( x i − , k ) − u s,n ( x i +1 , k )) + 2 sgn ( s ) sin( k y a ) u s,n ( x i , k )] − t ( x i ) [ v s,n ( x i +1 , k ) + v s,n ( x i − , k )] − { t ( x i ) [cos( k y a ) + cos( k z a )] − µ } v s,n ( x i , k ) = E n ( k ) v s,n ( x i , k ) . (S183)In contrast to the singlet cases of Sec. II and III, the BdG eqs. for each H s retain a dependence on spin index s . B. Landau levels in the continuum limit
We proceed as in Sec. II B and consider ansatze u s,n ( x i , k y , k z ) ≈ e iK ( α ) F z l u ( α ) s,n ( x i , k y , k z ) + e iK (¯ α ) F z l u (¯ α ) s,n ( x i , k y , k z ) , (S184)4for the wavefunctions in the vicinity of the two valleys located at K = (0 , ± K F ) and labeled by (1) and (1 (cid:48) ) . A similarexpansion is expected to hold for the v ’s. With this ansatz, the BdG eqs. for each H s further separate into two valley sectors.We proceed by taking the continuum limit. Since the envelope functions in Eq. S184 are assumed to vary slowly on the scale ofthe lattice, we can expand u ( α ) s,n (cid:16) x + a, q y , K ( α ) F + q z (cid:17) + u ( α ) s,n (cid:16) x − a, q y , K ( α ) F + q z (cid:17) ≈≈ u ( α ) s,n ( x, q y , q z ) + a∂ x u ( α ) s,n ( x, q y , q z ) + u ( α ) s,n ( x, q y , q z ) − a∂ x u ( α ) s,n ( x, q y , q z )= 2 u ( α ) s,n ( x, q y , q z ) . (S185)Similarly, we have v ( α ) s,n ( x + a, q y , q z ) − v ( α ) s,n ( x − a, q y , q z ) ≈ a∂ x v ( α ) s,n ( x, q y , q z ) . (S186)We can likewise expand to first order in q y , q z cos (( K F y + q y ) a ) =1 (S187) cos (( K F z + q z ) a ) = cos ( K F a ) − a sgn ( K F z ) sin ( K F ) q z (S188) sin ( K F y + q y ) = aq y , (S189)since K F y = 0 .Allowing for t ( x ) = t + δt ( x ) , we write Eq. S182 by retaining terms to first order in either q y , ∂ x , δt ( x ) : { δt ( x ) [4 + 2 cos ( K F a )] − ta sgn ( K F z ) sin ( K F ) q z } u ( α ) s,n ( x, q y , q z ) + [ − ia ∆ ∂ x − ia ∆ sgn ( s ) q y ] v ( α ) s,n ( x, q y , q z )= E n ( q y , q z ) u ( α ) s,n ( x, q y , q z ) , (S190)with a similar expression for Eq. S183. We define A ( x ) = δt ( x ) [4 + 2 cos ( K F a )] e (S191) v F =2 ta sin ( K F ) (S192) v ∆ = a ∆ . (S193)The BdG eqs. can be written as v F (cid:20) − sgn ( K F z ) q z + eAv F (cid:21) u ( α ) s,n ( x, q y , q z ) − iv ∆ [ ∂ x + sgn ( s ) q y ] v ( α ) s,n ( x, q y , q z ) = E n ( q y , q z ) u ( α ) s,n ( x, q y , q z ) (S194) − iv ∆ [ ∂ x − sgn ( s ) q y ] u ( α ) s,n ( x, q y , q z ) − v F (cid:20) − sgn ( K F z ) q z + eAv F (cid:21) v ( α ) s,n ( x, q y , q z ) = E n ( q y , q z ) v ( α ) s,n ( x, q y , q z ) , (S195)or H ( α ) s Ψ ( α ) σ,n = E n ( q y , q z )Ψ ( α ) s,n , (S196) H ( α ) s = v F (cid:20) − sgn (cid:16) K ( α ) F z (cid:17) q z + eAv F (cid:21) τ z + ( − i∂ x ) τ x + v ∆ sgn ( s ) q y τ y (S197)Let us consider the solution at valley (1) . We assume that δt ( x ) is chosen to generate a positive pseudo-magnetic field. As inSec. II C, we apply the transformation in Eq. S49 s.t. σ z → σ y . The BdG Hamiltonian becomes ˜ H (1) s = v F (cid:20) − q z + eAv F (cid:21) τ y + v ∆ ( − i∂ x ) τ x − v ∆ sgn ( s ) q y τ z = − v ∆ sgn ( s ) q y τ z + v ∆ (cid:26) ( − i∂ x ) τ x + (cid:20) − λq z + eAv ∆ (cid:21) τ y (cid:27) . (S198)Following Ref. 27, we can choose the positive-energy solutions as5 E n ( q y , q z ) = (cid:113)(cid:0) ω c √ n (cid:1) + ( v ∆ q z ) , n = 0 , , . . . (S199) ˜Ψ (1) s,n = sgn ( s ) (cid:18) − iα n ψ n − [ x, λq z ] β n ψ n [ x, λq z ] (cid:19) , (S200)where ψ ’s are the normalized harmonic oscillator wavefunctions and the coefficients are given by α n = (cid:115) E n − v ∆ sgn ( s ) q y E n (S201) β n = (cid:115) E n + v ∆ sgn ( s ) q y E n . (S202)The overall phase sgn ( s ) is added to ensure that the solutions obey time-reversal symmetry. In the un-rotated basis the solutionsare Ψ (1) s,n = sgn ( s ) 1 √ (cid:18) − iα n ψ n − [ x, λq z ] − iβ n ψ n [ x, λq z ] − α n ψ n − [ x, λq z ] + β n ψ n [ x, λq z ] (cid:19) . (S203)For the zeroth LL, we have E ( q y , q z ) = v ∆ sgn ( q y ) q y , (S204)and we require β = 1 . The solutions are given by (cid:32) u (1) s, ( x, q y , q z ) v (1) s, ( x, q y , q z ) (cid:33) = sgn ( s ) ψ ( x, q ) √ (cid:18) − i (cid:19) , sgn(s) q y ≥ . (S205)Therefore, the chiral spin-up zeroth LL mode is right-moving while it’s spin-down counterpart is left-moving.We apply a similar procedure for valley (1 (cid:48) ) . Namely, we apply the inverse transformation which maps σ z → − σ y to get ˜ H (1 (cid:48) ) s = − v F (cid:20) q z + eAv F (cid:21) τ y + v ∆ ( − i∂ x ) τ x + v ∆ sgn ( s ) q y τ z = v ∆ sgn ( s ) q y τ z + v ∆ (cid:26) ( − i∂ x ) τ x + (cid:20) − λq z − eAv ∆ (cid:21) τ y (cid:27) . (S206)Note the opposite sign for the vector potential wrt valley (1) , in accordance with time-reversal symmetry. From Ref. 27, wewrite the solutions as ˜Ψ (1 (cid:48) ) s,n = sgn ( s ) (cid:18) − iα n ψ n [ x, λq z ] β n ψ n − [ x, λq z ] (cid:19) , (S207)with α n = (cid:115) E n + v ∆ sgn ( s ) q y E n (S208) β n = (cid:115) E n − v ∆ sgn ( s ) q y E n . (S209)In the un-rotated basis the solutions are Ψ (1 (cid:48) ) s,n = sgn ( s ) √ (cid:18) − iα n ψ n [ x, λq z ] + iβ n ψ n − [ x, λq z ] α n ψ n [ x, λq z ] + β n ψ n [ x, λq z ] (cid:19) . (S210)Formally, the solutions can be converted to the initial sign of the vector potential by sending q z → − q z and taking into accountthe parity of the Hermite polynomials.For the zeroth LL, β = 0 and α = 1 and the coefficients are given by6 Ψ (1 (cid:48) ) s,n = sgn ( s ) ψ [ x, − λq z ] √ (cid:18) − i (cid:19) , (S211)with energies E ( q y , q z ) = v ∆ sgn ( q y ) q y . (S212)The zeroth LLs at valley (1 (cid:48) ) have the same chiralities as their valley (1) correspondents.The explicit zeroth LL solutions are γ (1) ↑ , ( q y , q z ) = (cid:90) dx ψ ( x, q z ) √ (cid:104) i Ψ (1) ↑ ( x, q y , q z ) − Ψ (1 (cid:48) ) , † ↑ ( x, − q y , − q z ) (cid:105) , (S213) γ (1) ↓ , ( q y , q z ) = (cid:90) dx ψ ( x, q z ) √ (cid:104) − i Ψ (1) ↓ ( x, q y , q z ) − Ψ (1 (cid:48) ) , † ↓ ( x, − q y , − q z ) (cid:105) (S214) γ (1 (cid:48) ) ↑ , ( q y , q z ) = (cid:90) dx ψ ( x, − q z ) √ (cid:104) i Ψ (1 (cid:48) ) ↑ ( x, q y , q z ) − Ψ (1) , † ↑ ( x, − q y , − qz ) (cid:105) (S215) γ (1 (cid:48) ) ↓ , ( q y , q z ) = (cid:90) dx ψ ( x, − q z ) √ (cid:104) − i Ψ (1 (cid:48) ) ↓ ( x, q y , q z ) − Ψ (1) , † ↓ ( x, − q y , − q z ) (cid:105) . (S216)As for the case of singlet pairing, only two out of the above four states are independent since γ (1 (cid:48) ) ↑ , ( q y , q z ) = − iγ (1) , † ↑ , ( − q y , − q z ) (S217) γ (1 (cid:48) ) ↓ , ( q y , q z ) = iγ (1) , † ↓ , ( − q y , − q z ) . (S218) C. Projection of fields onto the zeroth LLs
We follow the arguments of Sec. II D to obtain the Fermi fields via the inverse BdG transformation: Ψ ( α ) σ ( x, q ) = (cid:88) n,s u ( α ) sσ,n ( x, q ) γ ( α ) s,n ( k ) + v (¯ α ) , ∗ sσ,n ( x, − q ) γ (¯ α ) , † s,n ( − q ) . (S219)In the spin-triplet case u ( α ) sσ = δ sσ u s while v ( α ) sσ = ( − σ z ) sσ v ¯ s . Only two out of the four degrees-of-freedom in Eqs. S213-S216are independent. Therefore we must truncate the Hilbert space of the zeroth LL to ensure a well-defined BdG transformation: Ψ (1) P, ↑ ( x, q ) = − iψ ( x, q z ) √ γ (1) ↑ , ( q ) (S220) Ψ (1) P, ↓ ( x, q ) = iψ ( x, q z ) √ γ (1) ↓ , ( q ) (S221) Ψ (1 (cid:48) ) P, ↑ ( x, q ) = − ψ ( x, − q z ) √ γ (1) , † ↓ , ( − q ) (S222) Ψ (1 (cid:48) ) P, ↓ ( x, q ) = − ψ ( x, − q z ) √ γ (1) , † ↑ , ( − q ) . (S223)In turn, the projections onto the zeroth LL are not independent: Ψ (1 (cid:48) ) P,σ ( x, q ) = i sgn ( σ )Ψ (1) , † P,σ ( x, − q ) . (S224)As in the case of singlet pairing, these relations are matched by the real-space operators Ψ (1 (cid:48) ) P,σ ( x, y, z ) = (cid:90) (cid:90) (cid:18) dq y π (cid:19) (cid:18) dq z π (cid:19) e iq y y e iq z z Ψ (1 (cid:48) ) P,σ ( x, q y , q z )= (cid:90) (cid:90) (cid:18) dq y π (cid:19) (cid:18) dq z π (cid:19) e − iq y y e − iq z z Ψ (1 (cid:48) ) P,σ ( x, − q y , − q z )= i sgn ( σ ) (cid:90) (cid:90) (cid:18) dq y π (cid:19) (cid:18) dq z π (cid:19) e − iq y y e − iq z z Ψ (1) , † P,σ ( x, q y , q z )= i sgn ( σ )Ψ (1) , † P,σ ( x, y, z ) . (S225)7We can map the valley (1 (cid:48) ) operators to those at valley (1) . From Sec. II D, these satisfy { Ψ (1) σ ( r ) , Ψ (1) , † σ (cid:48) ( r (cid:48) ) } = ≈ δ σ,σ (cid:48) λl B e − ( x − x (cid:48) )24 λ l B e − ( z − z (cid:48) )24 λ l B e i ( x + x (cid:48) )( z − z (cid:48) )2 λl B δ ( y − y (cid:48) ) (S226) { Ψ (1) σ , Ψ (1) σ (cid:48) } =0 (S227) { Ψ (1) , † σ , Ψ (1) , † σ (cid:48) } =0 . (S228)The projection of the field Ψ σ ( r ) with contributions from both valleys, onto the zeroth LL is given by Ψ P,σ ( r ) = e iK Fz z Ψ (1) P,σ ( r ) + e − iK Fz z Ψ (1 (cid:48) ) P,σ ( r )= e iK Fz z Ψ (1) P,σ ( r ) + i sgn ( σ ) e − iK Fz z Ψ (1) , † P,σ ( r ) . (S229)We deduce that Ψ † P,σ ( r ) = e − iK Fz z Ψ (1) , † P,σ ( r ) − i sgn ( σ ) e iK Fz z Ψ (1) P,σ ( r )= − i sgn ( σ )Ψ P,σ ( r ) . (S230)This relation clearly signals emergent real-space Majorana fermions. Indeed, these projected fields obey { Ψ P,σ ( r ) , Ψ P,σ (cid:48) ( r (cid:48) ) } = cos (cid:18)(cid:20) K F + ( x + x (cid:48) )2 λl B (cid:21) ( z − z (cid:48) ) (cid:19) (cid:34) √ λl B e − ( x − x (cid:48) )22 λ l B (cid:35) (cid:34) √ λl B e − ( z − z (cid:48) )24 λ l B (cid:35) δ σ,σ (cid:48) δ ( y − y (cid:48) ) . (S231)The first term represents fast oscillations. We recognize the remaining terms as the contributions from the zeroth LL harmonicoscillator wavefunctions. The last term is a Dirac delta function which is typical of free Majorana fields [16].The real-space Majorana nature of the projected fields leading to Eq. S231 was illustrated in the continuum. However, thisproperty is not intrinsic to this limit. Indeed, it follows from the general p-h symmetry and topology of the bulk Hamiltonianunder strain as illustrated by the general arguments of Sec. I. We can therefore extend the Majorana nature to operators definedon the lattice such that c † P,σ ( R ) = − i sgn ( σ ) c P, ¯ σ ( R ) , as mentioned in the main text. These are expected to obey a lattice versionof Eq. S230. V. ESTIMATES OF STRAINA. General Formulae
In the continuum limit, the effect of uniaxial strain of the hopping coefficients of our models is given by [2, 24] δt ( x ) = tβu xx , (S232)where t is the pristine value of the hopping. β is given by β = d ln tlna , (S233)where a is the nearest-neighbor distance. It is typically associated with the Gr¨uneisen parameter [28]. u xx is a diagonalcomponent of the strain tensor.In the low-energy theory, δt ( x ) determines an effective vector potential via ˜ A y ( x ) = δt ( x ) e = tβu xx ( x ) e , (S234)where e is the bare electronic charge. For illustration purposes, we chose the vector potential along the y direction, but general-ization to other cases of uniaxial strain are obvious.In order to produce an approximately uniform pseudo-magnetic field, we require8 ∂ x ˜ A y ( x ) = (cid:18) tβe (cid:19) ∂ x u xx = B z . (S235)From the previous expression, it is clear that we require a strain which, to leading order, varies linearly with distance in thebulk of the sample: u xx ≈ ρx, (S236)where ρ is the gradient of the strain. It is convenient to define a dimensionless parameter ˜ ρ = ρa. (S237)Note that the strain parameter (cid:15) discussed in the main text is given by (cid:15) = βρ. (S238)All of the relevant scales of the problems can be written in terms of the parameters t, ∆ , β, ˜ ρ , representing hopping coefficient,pairing amplitude, Gr¨uneisen parameter, and rate of strain increase per unit cell, respectively. The pseudo-magnetic length isdetermined by l B = (cid:114) v ∆ eB ≈ a (cid:115)(cid:18) ∆ t (cid:19) (cid:18) β ˜ ρ (cid:19) . (S239)The LL separation in energy is given by E c = (cid:126) ω c = √ al B =∆ (cid:115) β ˜ ρ ∆ /t , (S240)where we re-introduced factors of (cid:126) in Eq. S45 with v ∆ ≈ ∆ a/ (cid:126) . B. Numerical Estimates
Using the formulae of the previous section, we can estimate typical pseudo-magnetic lengths and LL separation in energyusing parameters for high- T c materials: ∆ ≈
30 meV [7], a ≈ A , t ≈ β ≈ . [30].Considering a maximum strain of 0.05 achieved over a typical length scale approximately L = ˜ ρ as ˜ ρ =0 . (cid:16) aL (cid:17) =3 . × − . (S241)This gives l B ≈ a (S242) E c ≈ . meV (1 ..
30 meV [7], a ≈ A , t ≈ β ≈ . [30].Considering a maximum strain of 0.05 achieved over a typical length scale approximately L = ˜ ρ as ˜ ρ =0 . (cid:16) aL (cid:17) =3 . × − . (S241)This gives l B ≈ a (S242) E c ≈ . meV (1 .. K ) ..