Effect of Zeeman coupling on the Majorana vortex modes in iron-based topological superconductors
Areg Ghazaryan, Pedro L. S. Lopes, Pavan Hosur, Matthew J. Gilbert, Pouyan Ghaemi
EEffect of Zeeman coupling on the Majorana vortex modes in iron-based topologicalsuperconductors
Areg Ghazaryan, P. L. S. Lopes, Pavan Hosur, Matthew J. Gilbert,
4, 5, 6 and Pouyan Ghaemi
7, 8 IST Austria (Institute of Science and Technology Austria), Am Campus 1, 3400 Klosterneuburg, Austria Stewart Blusson Quantum Matter Institute, University of British Columbia, Vancouver, British Columbia, Canada V6T 1Z4 Department of Physics, University of Houston, Houston, TX 77204, USA Micro and Nanotechnology Laboratory, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA Department of Electrical and Computer Engineering,University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA Department of Electrical Engineering, Stanford University, Stanford, California 94305, USA Physics Department, City College of the City University of New York, New York, NY 10031, USA Graduate Center of the City University of New York, NY 10031, USA (Dated: July 26, 2019)In the superconducting regime of FeTe (1 − x ) Se x , there exist two types of vortices which are dis-tinct by the presence or absence of zero energy states in their core. To understand their origin,we examine the interplay of Zeeman coupling and superconducting pairings in three-dimensionalmetals with band inversion. Weak Zeeman fields are found to suppress the intra-orbital spin-singletpairing, known to localize the states at the ends of the vortices on the surface. On the other hand,an orbital-triplet pairing is shown to be stable against Zeeman interactions, but leads to delocal-ized zero-energy Majorana modes which extend through the vortex. In contrast, the finite-energyvortex modes remain localized at the vortex ends even when the pairing is of orbital-triplet form.Phenomenologically, this manifests as an observed disappearance of zero-bias peaks within the coresof topological vortices upon increase of the applied magnetic field. The presence of magnetic im-purities in FeTe (1 − x ) Se x , which are attracted to the vortices, would lead to such Zeeman-induceddelocalization of Majorana modes in a fraction of vortices that capture a large enough number ofmagnetic impurities. Our results provide an explanation to the dichotomy between topological andnon-topological vortices recently observed in FeTe (1 − x ) Se x . Introduction:
To date, one of the major impedimentsin the search for Majorana fermions (MFs) is that a req-uisite topological superconductivity, either intrinsic [1, 2]or induced in a host material via a proximity coupling toa standard s -wave superconductor [3–6]. Of the avail-able materials that possess topology, superconductivityand magnetism, iron-based superconductors are of recentinterest [7–15]. In particular, the iron-based supercon-ductor FeTe . Se . (FTS) has recently been shown tohave strong spin-orbit interactions and band inversionthat result in a helical, topologically-protected, Diraccone on the surface [16–19]. The phenomenology of vor-tices, proliferated in the presence of magnetic fields, isalso noteworthy in FTS [20–23]. The low charge densityin the superconducting phase of this system is experimen-tally advantageous as it results in large Caroli-de Gennes-Matricon (CDM) vortex mode gaps [24] which facilitatesthe spectral detection of zero-energy vortex modes viascanning tunneling microscopy (STM). Intriguingly, vor-tices in FTS show two distinct types of behavior: topo-logical, carrying zero energy states consistent with thepresence of MF, and trivial vortices that lack the zeroenergy state but carry finite energy CDMs.While the evidence of MFs in FTS has been observed, acomprehensive understanding of the salient physics of thevortex composition is lacking. More precisely, the energyspectra of the vortices follow different hierarchies relativeto the CDM vortex gap, δ = ∆ µ , where ∆ is the bulk su- perconducting gap and µ is the chemical potential. Thetrivial vortex energy spectrum scales as ( n + 1 / δ , (with n ∈ Z ) which does not include the zero mode, whereasthe topological vortices follow nδ . Additionally, the per-centage of vortices with zero-energy modes decreases asthe perpendicular magnetic field increases, despite thefact that the inter-vortex distances are generically largerthan the superconducting coherence length [20, 21, 25].It was also shown that the distribution of vortices withand without zero mode has no correlation with the chargedisorder on the surface of the material [21]. These fea-tures, on aggregate, suggest that the properties that dis-tinguish the two classes of vortices stem from the bulkproperties of individual vortex rather than those of thesurface states. In particular the effects of magnetic field,beyond generation of vortices, might crucially affect theproperties of superconducting state and the vortices.Motivated by experimental observations, we examinethe effects of Zeeman coupling on the vortex modes inFTS. As the topological properties of FTS are driven byband inversion, we eschew more complex band modelsand utilize a simple model of a doped 3D time-reversalsymmetric (TRS) topological insulator (TI) which hasbeen used as an appropriate toy model to investigate theproperties of vortices in FTS [26, 27]. An essential prop-erty of the electronic band structure in TRS TI is thepresence of degenerate Fermi surfaces. Due to the strongspin-orbit coupling, Zeeman field splits the the degen- a r X i v : . [ c ond - m a t . s up r- c on ] J u l erate Fermi surface into two helical Fermi surfaces withopposite helicity (see supplementary material). In thisletter, we show that the split Fermi surfaces prefer anorbital-triplet superconducting pairing which delocalizesthe MF modes at the ends of the vortices on the surface.To make direct connection with the dichotomy of vor-tices in FTS we note that the Zeeman field may resultfrom the magnetic impurities along the vortex core [28].Interestingly, such magnetic Fe impurities are known toexists in FTS [29, 30] and they are attracted to the vor-tices [28]. In addition, increase of magnetic field naturallyleads to enhancement of Zeeman coupling which furtherdestabilize the topological vortices, as is experimentallyobserved. We should also note that neutron scatteringmeasurements have shown the evidence of ferromagneticclusters of Fe atoms in FTS [29]. Observation of theclusters of vortices, with and without zero energy vortexmodes, would further support our theory. Model Hamiltonian:
The two-orbital model of a 3DTRS TI is represented by the tight-binding Hamiltonian H k = (cid:80) k ψ † k [ τ x d k . σ + m k τ z − µ ] ψ k where Pauli matri-ces σ i and τ i act on spin and orbital space respectively, d k i = 2 t sin( k i ), m k = M + m (cid:80) i cos( k i ) and M , m , t are parameters of the model and µ is the chemical poten-tial. By varying the parameters, this model Hamiltoniancould represent both strong and week TRS TIs [31]. Thetight-binding model is used for the numerical calculationswhile our analytical results are based on the low-energyeffective model H = (cid:126) v F τ x σσσ · k + m k τ z , (1)where m k = m + (cid:15)k is the effective mass term and k is themomentum relative to the center of the Brillouin zone.The trivial (topological) insulator corresponds to m(cid:15) > < (cid:126) v F = 1.With the Hamiltonian defined, we begin our analysisin the metallic phase, when µ > | m | . The model dis-plays two degenerate Fermi surfaces that split by a Zee-man field ∆ Z . Since the bulk band structure gap is largecompared with superconducting gap, we use the effectiveHamiltonian resulting from projecting the Hamiltonianin Eq. (1) into the states at the two Fermi surfaces: H = (cid:90) d k f † k [( E k − µ ) ν + d k · ν ] f k . (2)Here ν i are the identity or Pauli matrices acting onthe the space of two Fermi surfaces. The vector d k = ∆ Z (cid:16) − m k E k k x + k y | k | , , k z k (cid:17) presents the Zeeman field and f k is the Fermion fields ψ k projected onto the Fermi sur-faces (see supplementary materials). The two fermi spin-split fermi surfaces are identified by diagonalizing theprojected Hamitonian (2) in the ν space.Previous analysis [32] identified two types ofsuperconductivity in doped TIs which are en-ergetically favorable: intra-orbital spin singlet, (cid:82) d r ψ † ∆ iτ σ y ψ † T + H.c. , and inter-orbital orbital-triplet spin-singlet, (cid:82) d k ψ † k iτ x σ y ψ † T − k + H.c. . Hence-forth, we will refer to these superconducting pairingsas intra-orbital singlet and inter-orbital triplet pairings,respectively.Upon projection onto the Fermi surfaces the supercon-ducting pairing potentials assume the following form,∆ α (cid:90) d k f † k e − iφ k ν α f † T − k , (3)where for the intra-orbital singlet pairing α = 1 and forinter-orbital triplet pairing α = 0. Examination of thesuperconducting pairing term in Eq. (3) in comparisonwith the kinetic Hamiltonian in Eq. (2) shows that theintra-orbital singlet pairing, ∆ , corresponds to pairingelectrons between different Zeeman split fermi surfaces.In contrast, the inter-orbital triplet pairing, ∆ , pairselectrons solely within each of the Zeeman split Fermisurfaces. Thus, the effect of Zeeman coupling is to breakthe TRS within the model and imbalance the two pairingpotentials in the favor of ∆ which couples the electronssolely within each Zeeman split Fermi surface [33].To examine the outlined effect of Zeeman coupling onthe dominant form of superconducting pairing, we utilizea linear gap equation to determine the critical tempera-tures of the two superconducting pairings [32, 34–36](seesupplementary materials). The corresponding U − V model, in which U and V are the intra and inter orbitalinteraction, leads to the equation :det (cid:12)(cid:12)(cid:12)(cid:12) U ¯ χ − U χ V χ V χ − (cid:12)(cid:12)(cid:12)(cid:12) = 1 , V χ = 1 . (4)Here ¯ χ , χ and χ are the superconducting suscepti-bilities characterizing intra-orbital spin-singlet pairingand χ describes the inter-orbital spin triplet pairing.¯ χ = − (cid:82) w D − w D D ( ξ ) tanh ( ξ/ T ) / ξdξ is the standard s -wave susceptibility, D ( ξ ) is the density of states and w D is the Debye frequency.By numerically solving the U − V Eq. (S74), we ob-tain the critical temperature T c for each pairing channel.Fig. 1 shows the resulting phase boundaries that delin-eate the regions where either ∆ or ∆ correspond tohigher critical temperature and so is the dominant formof pairing. In Fig. 1, we plot the phase boundary as wevary ∆ Z and the ratio U/V . It is evident in Fig. 1 thatthe inclusion of the Zeeman effect results in the enhance-ment of the triplet (∆ ) pairing and the suppression ofthe singlet one (∆ ) at a given chemical potential. Vortex Modes:
Having established the phase diagramof the superconducting pairing, we proceed to study theinternal structure of the vortices as we insert a π -flux-tube into the pairing potential: ∆ α ( r ) = | ∆ α ( r ) | e iθ [37,38].We fix k z = 0 as we are interested in points wherethe topological Z index changes and this only occurs at (cid:2) (cid:1) (cid:6) (cid:2) (cid:2) (cid:1) (cid:6) (cid:4) (cid:2) (cid:1) (cid:6) (cid:5) (cid:2) (cid:1) (cid:6) (cid:7) (cid:2) (cid:1) (cid:6) (cid:8)(cid:2) (cid:1) (cid:2)(cid:2) (cid:1) (cid:4)(cid:2) (cid:1) (cid:5)(cid:2) (cid:1) (cid:7)(cid:2) (cid:1) (cid:8)(cid:3) (cid:1) (cid:2) (cid:1) (cid:1) (cid:6) (cid:8) (cid:13) (cid:4) (cid:11) (cid:1) (cid:9) (cid:11) (cid:3) (cid:6) (cid:13) (cid:2) (cid:7)(cid:13) (cid:11) (cid:6) (cid:10) (cid:7) (cid:4) (cid:13) D (cid:1) (cid:3) (cid:1) (cid:2) (cid:1) (cid:4) (cid:6) (cid:8) (cid:13) (cid:11) (cid:2) (cid:1) (cid:9) (cid:11) (cid:3) (cid:6) (cid:13) (cid:2) (cid:7)(cid:12) (cid:6) (cid:8) (cid:5) (cid:7) (cid:4) (cid:13) FIG. 1. Phase boundaries between regions with supercon-ducting order parameters ∆ or ∆ as function of the ratio U/V of interaction strengths of each channel and the fieldmagnitude, ∆ Z (Red solid curve). The blue dashed line isguide for the eye of the case when ∆ Z = 0. The parametersof the Hamiltonian are m = − . (cid:15) = 0 .
5. IncreasedZeeman coupling results in larger regions where inter-orbitaltriplet pairing is the ground state. k z = 0 or π . Since the vortex modes stem from statesclose to the Fermi energy, their wave functions can beexpressed in terms of a superposition of the TI conduc-tion band eigenstates in cylindrical coordinates χ νl,k (seesupplementary material).Here l s are angular momentum quantum numbers and ν = ± correspond to the two energy bands of Hamilto-nian (1) which are split by Zeeman coupling.Given the rotational symmetry of the vortex profile,the vortex modes with different l s are not hybridized andthe vortex Hamiltonian decouples into sectors associatedwith each l . Also for inter-orbital triplet pairing, thevortex does not mixes the ν ’s. On the other hand, thetranslation symmetry in the plane perpendicular to thevortex is broken and different radial momenta k , as wellas Nambu particle-hole states are mixed by the vortex.The effective vortex Hamiltonian acting in the radial mo-mentum and Nambu particle-hole spaces takes the formof a 1D Jackiw-Rebbi model [39] (see supplementary ma-terial) H νl = Π z ∆ λ ν ( k ) ξ + Π y ∆ ξ ( i∂ k ) + Π x E νl, ( k ) . (5)The Π matrices act on a Nambu space, ξ is thesuperconducting coherence length and ∆ λ ν ( k ) /ξ = k + hm k E k ν − µ . The Jackiw-Rebbi lowest-energy solu-tions are localized at the Fermi surface where the co-efficient λ ν ( k νF ) changes sign. These states have theform e − (cid:82) kνF + kkνF − k dk (cid:48) λ ν ( k (cid:48) ) (1 , T in Nambu particle-hole ba-sis and has energy E νl, ( k F ). Notice that vortex-modewave functions are exponentially localized around Fermi Energy
C h e m i c a l P o t e n t i a l ( a ) ( b )
C h e m i c a l P o t e n t i a l
FIG. 2. Dependence of the vortex modes energies on thechemical potential µ obtained from the 3D lattice model withperiodic boundary conditions along the z -direction for k z = 0.(a) Intra-orbital singlet pairing ∆ = 0 . Z = 0. (b)Inter-orbital triplet pairing ∆ = 0 . Z = 0 .
2. Theparameters for the model are: M = 4 . m = − . t = 1 . ×
48 lattice inthe x − y plane. We observe a zero-energy mode state existsfor all chemical potentials inside the conduction band whenthe pairing is inter-orbital triplet type. wavevector. Since the BdG Hamiltonian in equation (5)is in the bases of TI conduction band states, the full wavefunction of the two vortex modes (each associated withone Zeeman split Fermi surface) takes the formΨ νV ( l, r ) ≈ χ νl,k F ( r ) (1 , T Π (6)where (1 , T Π is the spinor in Nambu particle-hole space.Previously, a similar result was obtained for the intra-orbital spin-singlet case [37, 38]. In that case, solutionswhere again centered at the metallic phase Fermi surface,with corresponding energies E νl, = ∆ k F ξ (2 πl + π + νφ B ). φ B is a Berry-phase-like term which permits zero-modeswhenever φ B = π . In contrast, for the inter-orbitaltriplet case, the energies of the vortex modes are E νl, = ∆ k νF ξ (2 πl ). They noticeably lack the Berry phase termpresent in the intra-orbital singlet case and the l = 0states always has zero energy. Therefore, in the inter-orbital triplet case, a zero-energy channel always exists along the vortex. Its presence delocalizes the zero-modesat the end of the vortex, on the sample surface. De-localization of the zero mode at the ends of the vortexresults in a suppression of the zero-bias signal in STMmeasurements.In Fig. 2 we verify the analytic results using a 3D lat-tice model with periodic boundary conditions along thez-direction. Fig. 2 (a) shows the spectrum of the vortexmodes for the intra-orbital singlet pairing where the vor-tex gap closes solely when φ B = π . In contrast, Fig. 2 (b)shows that for the inter-orbital triplet pairing, once thechemical potential is in the conduction band the systemdevelops vortex zero modes which remain gapless for allchemical potentials. Thus, increasing the Zeeman cou-pling results in a shift of the Fermi surfaces that destabi-lize the intra-orbital singlet pairing in favor of the inter- FIG. 3. Spatial profile of the lowest energy vortex mode in3D slab geometry for different values of parameter α , whichcontrols the amplitude of the inter-orbital triplet pairing ∆ = α ∆ , intra-orbital s-wave ∆ = (1 − α )∆ and Zeeman field∆ Z = α ∆ Z around the vortex. The calculation is performedon a 24 × ×
24. The parameters used are: µ = 1 .
1, ∆ = 0 . = 0 .
4, ∆ Z = 0 . M = 4 . m = − . t = 1 . orbital triplet pairing leading to the omnipresence of azero mode in the vortex core.In Fig. 3, we examine the manifestation of the changein superconductive pairing by numerically inserting avortex in the tight-binding lattice model. We set µ =1 . t , where t is hoping amplitude in the lattice model.For small Zeeman couping, the intra-orbital singlet chan-nel is dominant and the vortex states, both topologicaland trivial, are localized at the ends of the vortex on thesurface. As the Zeeman splitting is increased, we desta-bilize the intra-orbital singlet pairing channel in favor ofthe inter-orbital triplet pairing allowing MFs, localized atthe ends of the vortex string, to penetrate into the bulk.At sufficiently high Zeeman coupling, the intra-orital sin-glet pairing is fully suppressed and MFs on the surfacedelocalize through the vortex modes and extend into thebulk of the superconductor. We may also understandthis mechanism, starting from the case where the pairingis solely inter-orbital triplet type where the vortex hoststwo zero modes extended through the bulk. Intra-orbitalsinglet pairing hybridizes these two modes and generatea gap for the vortex modes extending through the bulkand localizes the MFs at the ends of the vortex on thesurface. Interestingly, as will be shown below, this mech-anism would only delocalize the Majorana zero modesand could not do the same for finite energy vortex modesat the ends of the vortex. It is then well consistent withthe experimental results which shows that finite energyvortex modes are intact in all vortices, whereas the MFsare absent in some of the vortices. Effective 1D Vortex Model:
Using the wavefunctionsfrom Eq. (6), we now derive an effective Hamiltonian forthe modes along the vortex line connecting the surfaces (cid:1)(cid:3) (cid:2)(cid:6) (cid:1)(cid:3) (cid:2)(cid:5) (cid:1)(cid:3) (cid:2)(cid:4) (cid:3) (cid:2)(cid:3) (cid:3) (cid:2)(cid:4) (cid:3) (cid:2)(cid:5) (cid:3) (cid:2)(cid:6)(cid:1)(cid:3) (cid:2)(cid:6)(cid:1)(cid:3) (cid:2)(cid:5)(cid:1)(cid:3) (cid:2)(cid:4)(cid:3) (cid:2)(cid:3)(cid:3) (cid:2)(cid:4)(cid:3) (cid:2)(cid:5)(cid:3) (cid:2)(cid:6) (cid:1)(cid:3) (cid:2)(cid:6) (cid:1)(cid:3) (cid:2)(cid:5) (cid:1)(cid:3) (cid:2)(cid:4) (cid:3) (cid:2)(cid:3) (cid:3) (cid:2)(cid:4) (cid:3) (cid:2)(cid:5) (cid:3) (cid:2)(cid:6)(cid:1)(cid:3) (cid:2)(cid:5)(cid:1)(cid:3) (cid:2)(cid:4)(cid:3) (cid:2)(cid:3)(cid:3) (cid:2)(cid:4)(cid:3) (cid:2)(cid:5) (cid:1) (cid:4) (cid:1) (cid:1) (cid:2) (cid:1) (cid:3) (cid:1) (cid:1) (cid:1) (cid:4) (cid:2)(cid:5) (cid:3)(cid:6) (cid:1) (cid:3) (cid:2) (cid:1) (cid:1) (cid:4) (cid:1) (cid:1) (cid:2) (cid:1) (cid:3) (cid:1) (cid:1) (cid:1) (cid:4) (cid:2)
FIG. 4. k z momentum dependence of the vortex modes for(a) intra-orbital s -wave pairing ∆ = 0 . = 0 . Z = 0 . µ = 1 .
1. Additional parameters of themodel are the same as in Fig. 2. [40]. To this end, we calculate the matrix elements of the k z -dependent terms in Hamiltonian, Eq. (1), within thespace of the two vortex zero modes presented in Eq. (6).Without lose of generality we consider m > (cid:15) < η i as Pauli matrices acting on the space of thetwo zero modes, the effective vortex Hamiltonian havethe general form H lV = E + l + E − l η (7)+ (cid:18) E + l − E − l − ˜ (cid:15)∂ z (cid:19) η z + ˜ v z ( ∂ z ) η x . Here, E νl are the energies of l -th vortex modes from Fermisurface corresponding to ν = ± . In Eq. (7), both ˜ (cid:15) and ˜ v z are parameters of the k z -dependent terms. The effectiveone dimensional Hamiltonian in Eq. (7) corresponds tothe effective vortex Hamiltonian and supports localizedstates at the two ends if (cid:0) E + l − E − l (cid:1) ˜ (cid:15) < λ ν ( k νF ), we notice k + F < k − F andconsequently, E + l − E − l >
0. Therefore, the condition forthe presence of localized states at the end of vortex withfinite energy (cid:0) E + l − E − l (cid:1) ˜ (cid:15) < l (cid:54) = 0 and finite energies of E + l + E − l are localized atthe ends of the vortex lines. In contrast, for the two zeroenergy modes ( l = 0) with (cid:0) E +0 − E − (cid:1) ˜ (cid:15) = 0. As a result,the MFs can not stay localized at the ends of the vortexwhen the bulk pairing is of inter-orbital triplet form. Asshown in Fig. 4(b), the effective vortex Hamiltonian forthe zero energy states is linearly dispersing with momen-tum along the vortex and, due to their Nambu particle-hole character, protected against back-scattering whichfurther supports the delocalization of zero modes uponformation of inter-orbital triplet pairing. Conclusion:
We have examined the effect of Zeemancoupling on the structure of vortex modes in FTS. Wefind that the inter-orbital triplet paring and intra-orbitalsinglet pairing compete as a function of Zeeman cou-pling, resulting in dramatically different vortex struc-tures. Intra-orbital singlet pairing leads to the presenceof Majorana vortex modes localized on the surface atthe end of vortex strings. On the other hand, inter-orbital triplet paring supports localized finite-energy triv-ial vortex modes but destabilize the zero-energy Majo-rana modes at the end of vortex strings. This delocal-ization is through the formation of zero modes which ex-tend along the vortex through the bulk. Such extendedstates have small coupling with STM probe and do notshow strong signals, in comparison with the case wherethe zero modes are localized at the end of vortex, on thesample surface. Our results shed light in the existenceof two types of vortices experimentally observed in FTSand support the topological nature of vortex modes inthis system.The main property of the band structure of FTS tofacilitate our models is the presence of degenerate Fermisurfaces that are split into two helical Fermi surfaces bythe Zeeman field. The size of the superconducting gapin FTS is of the order of 2 . . µ B [28, 29]. Given the average distance of the Fe im-purity atoms, their associated Zeeman coupling is of theorder of ∆ Z ≈ .
84 meV. It is then evident that the Zee-man coupling can affect the form of the superconductingparing. Our results then demonstrate that the nature ofthe vortices in FTS is inextricably linked to the effect ofZeeman coupling, which determines the form of super-conducting pairing. It also indicates that suppression ofZeeman coupling, by reduction of magnetic impurities,stabilizes the vortex MFs in FTS.
Acknowledgments.
This research was supported un-der National Science Foundation Grants EFRI-1542863and PSC-CUNY Award, jointly funded by The Profes-sional Staff Congress and The City University of NewYork (PG). PLSL acknowledges support by the CanadaFirst Research Excellence Fund. AG acknowledges sup-port from the European Unions Horizon 2020 researchand innovation program under the Marie Skodowska-Curie grant agreement No 754411. PH was supportedby the Department of Physics and the College of Natu-ral Sciences and Mathematics at the University of Hous-ton. M.J.G. acknowledges financial support from the Na-tional Science Foundation (NSF) under Grant No. DMR-1720633, CAREER Award ECCS-1351871 and the Officeof Naval Research (ONR) under grant N00014-17-1-3012. [1] Y. Li and Z.-A. Xu, Advanced Quantum Technologies ,1800112.[2] M. Sato and Y. Ando, Reports on Progress in Physics ,076501 (2017).[3] L. Fu and C. L. Kane, Phys. Rev. Lett. , 096407(2008). [4] R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev.Lett. , 077001 (2010).[5] V. Mourik, K. Zuo, S. Frolov, S. Plissard, E. Bakkers, andL. Kouwenhoven, Science , 1003 (2012).[6] S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon,J. Seo, A. H. MacDonald, B. A. Bernevig, and A. Yazdani,Science , 602 (2014).[7] N. Hao and J. Hu, National Science Review , 213 (2018).[8] P. Dai, Rev. Mod. Phys. , 855 (2015).[9] R. P. Day, G. Levy, M. Michiardi, B. Zwartsenberg,M. Zonno, F. Ji, E. Razzoli, F. Boschini, S. Chi, R. Liang,P. K. Das, I. Vobornik, J. Fujii, W. N. Hardy, D. A. Bonn,I. S. Elfimov, and A. Damascelli, Phys. Rev. Lett. ,076401 (2018).[10] M. H. Christensen, J. Kang, and R. M. Fernandes, Phys.Rev. B , 014512 (2019).[11] V. Cvetkovic and O. Vafek, Phys. Rev. B , 134510(2013).[12] A. Lau and C. Timm, Phys. Rev. B , 024517 (2014).[13] C. Youmans, A. Ghazaryan, M. Kargarian, andP. Ghaemi, Phys. Rev. B , 144517 (2018).[14] E. J. K¨onig and P. Coleman, Phys. Rev. Lett. ,207001 (2019).[15] R.-X. Zhang, W. S. Cole, X. Wu, and S. D.Sarma, “Higher order topology and nodal topological su-perconductivity in fe(se,te) heterostructures,” (2019),arXiv:1905.10647.[16] T. Hanaguri, S. Niitaka, K. Kuroki, and H. Takag, Sci-ence , 474 (2010).[17] Z. Wang, P. Zhang, G. Xu, L. K. Zeng, H. Miao, X. Xu,T. Qian, H. Weng, P. Richard, A. V. Fedorov, H. Ding,X. Dai, and Z. Fang, Phys. Rev. B , 115119 (2015).[18] P. Zhang, K. Yaji, T. Hashimoto, Y. Ota, T. Kondo,K. Okazaki, Z. Wang, J. Wen, G. D. Gu, H. Ding, andS. Shin, Science , 182 (2018).[19] G. Xu, B. Lian, P. Tang, X.-L. Qi, and S.-C. Zhang,Phys. Rev. Lett. , 047001 (2016).[20] D. Wang, L. Kong, P. Fan, H. Chen, S. Zhu, W. Liu,L. Cao, Y. Sun, S. Du, J. Schneeloch, R. Zhong, G. Gu,L. Fu, H. Ding, and H.-J. Gao, Science , 333 (2018).[21] T. Machida, S. P. Y. Sun, Y. K. S. Takeda, T. Hanaguri,T. Sasagawa, and T. Tamegai, “Zero-energy vortex boundstate in the superconducting topological surface state offe(se,te),” (2018), arXiv:1812.08995.[22] L. Kong, S. Zhu, M. Papaj, L. Cao, H. Isobe, W. Liu,D. Wang, P. Fan, H. Chen, Y. Sun, S. Du, J. Schneeloch,R. Zhong, G. Gu, L. Fu, H.-J. Gao, and H. Ding, “Obser-vation of half-integer level shift of vortex bound states inan iron-based superconductor,” (2019), arXiv:1901.02293.[23] S. Zhu, L. Kong, L. Cao, H. Chen, S. Du, Y. Xing,D. Wang, C. Shen, F. Yang, J. Schneeloch, R. Zhong,G. Gu, L. Fu, Y.-Y. Zhang, H. Ding, and H.-J. Gao,ArXiv:1904.06124.[24] C. Caroli, P. D. Gennes, and J. Matricon, Physics Letters , 307 (1964).[25] C.-K. Chiu, Y. H. T. Machida, T. Hanaguri, and F.-C. Zhang, “Scalable majorana vortex modes in iron-basedsuperconductors,” (2019), arXiv:1904.13374.[26] S. Qin, L. Hu, X. Wu, X. Dai, C. Fang, F. chun Zhang,and J. Hu, “Topological vortex phase transitions in iron-based superconductors,” (2019), arXiv:1901.03120.[27] P. Zhang, Z. Wang, X. Wu, K. Yaji, Y. Ishida,Y. Kohama, G. Dai, Y. Sun, C. Bareille, K. Kuroda,T. Kondo, K. Okazaki, K. Kindo, X. Wang, C. Jin, J. Hu, R. Thomale, K. Sumida, S. Wu, K. Miyamoto, T. Okuda,H. Ding, G. D. Gu, T. Tamegai, T. Kawakami, M. Sato,and S. Shin, Nature Physics , 41 (2019).[28] K. Jiang, X. Dai, and Z. Wang, Phys. Rev. X , 011033(2019).[29] V. Thampy, J. Kang, J. A. Rodriguez-Rivera, W. Bao,A. T. Savici, J. Hu, T. J. Liu, B. Qian, D. Fobes, Z. Q.Mao, C. B. Fu, W. C. Chen, Q. Ye, R. W. Erwin, T. R.Gentile, Z. Tesanovic, and C. Broholm, Phys. Rev. Lett. , 107002 (2012).[30] J.-X. Yin, Z. Wu, J.-H. Wang, Z.-Y. Ye, J. Gong, X.-Y.Hou, L. Shan, A. Li, X.-J. Liang, X.-X. Wu, J. Li, C.-S.Ting, Z.-Q. Wang, J.-P. Hu, P.-H. Hor, H. Ding, and S. H.Pan, Nature Physics , 543 EP (2015).[31] P. Hosur, S. Ryu, and A. Vishwanath, Phys. Rev. B ,045120 (2010).[32] L. Fu and E. Berg, Phys. Rev. Lett. , 097001 (2010).[33] Y. Kim, T. M. Philip, M. J. Park, and M. J. Gilbert,Phys. Rev. B , 235434 (2016).[34] S. Nakosai, Y. Tanaka, and N. Nagaosa, Phys. Rev. Lett. , 147003 (2012).[35] T. Hashimoto, S. Kobayashi, Y. Tanaka, and M. Sato,Phys. Rev. B , 014510 (2016).[36] P. Hosur, X. Dai, Z. Fang, and X.-L. Qi, Phys. Rev. B , 045130 (2014).[37] P. Hosur, P. Ghaemi, R. S. K. Mong, and A. Vishwanath,Phys. Rev. Lett. , 097001 (2011).[38] C.-K. Chiu, M. J. Gilbert, and T. L. Hughes, Phys. Rev.B , 144507 (2011).[39] R. Jackiw and C. Rebbi, Phys. Rev. D , 3398 (1976).[40] P. L. e. S. Lopes and P. Ghaemi, Phys. Rev. B , 064518(2015).[41] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. , 1057(2011). Effect of Zeeman coupling on the Majorana vortex modes in iron-basedtopological superconductors - Supplementary Material
In this supplementary material, we present (a) the derivation of the vortex spectrum in doped topological insulatorwith intra-orbital singlet and inter-orbital triplet superconductivity, (b) a derivation of effective vortex Hamiltonianwhich captures the the properties of localized stats at the ends of the vortex (c) the self-consistent calculations todetermine the effect of Zeeman coupling on the type of superconducting state.
VORTEX MODES IN DOPED TOPOLOGICAL INSULATORSModel Hamiltonian
The effective Hamiltonian from Eq. (1) of the main text is H = (cid:82) d rψ † ( r ) ( H − µ ) ψ ( r ) where ψ ( r ) =( ψ ↑ A , ψ ↓ A , ψ ↑ B , ψ ↓ B ) T . Here ψ σJ is the destruction operator for electron in orbital J with spin σ . In momentumspace, the Hamiltonian reads H = (cid:90) d k / (2 π ) ψ † k ( H k − µ ) ψ k , (S1)where H k = α · k + m k β (S2)with m k = m − (cid:15)k , α = τ x σ and β = τ z σ . τ and σ matrices act on the A, B orbitals and ↑↓ spins, respectively.The dispersions of the four bands follow E k = E k = − E k = − E k ≡ E k = (cid:112) k + m k . We expand the Fermionoperators in an eigenbasis as ψ k = (cid:88) a =1 ϕ a k f k ,a , (S3)where the eigenstates ϕ a k ’s are given by ϕ a k = e iφ k / e − iσ z φ k / e − iσ y θ k / e − iτ y σ z α k / ˆ e a , (S4)with tan φ k = k y k x cos θ k = k z k cos α k = m k E k , (S5)where ˆ e a is a basis of R corresponding to the four bands. The phases are chosen so that the wavefunctions are singlevalued upon a 2 π evolution of the azimuthal momentum-angle variable. The angle-variables transform according to k → − k as φ − k = φ k − π and θ − k = π − θ k , while α − k = α k . Notice α k depends only on the modulus of k . Explicitly,the two upper-band wavefunctions are ϕ k = e iφ k / e − iσ z φ k / e − iσ y θ k / e − iτ y σ z α k / = cos θ k cos α k e iφ k sin θ k cos α k cos θ k sin α k e iφ k sin θ k sin α k (S6)and ϕ k = e iφ k / e − iσ z φ k / e − iσ y θ k / e − iτ y σ z α k / = − sin θ k cos α k e iφ k cos θ k cos α k sin θ k sin α k − e iφ k cos θ k sin α k . (S7)We are going to assume a positive chemical potential, tuned inside the upper two bands. This allows neglecting thelower bands with negative energy as ψ k = (cid:80) a =1 ϕ a k f k ,a ≈ (cid:80) a =1 , ϕ a k f k ,a . This leads to H ≈ (cid:88) a =1 , (cid:90) k ( E k − µ ) f † k ,a f k ,a . (S8)where f k = ( f k , , f k , ) T contains the two positive energy modes only. Notice that the conduction band consists oftwo degenerate bands. Zeeman coupling
The Zeeman magnetization is incorporated through the operator H Z = ∆ Z (cid:82) d rψ † ( r ) σ z ψ ( r ). Projecting in thetwo upper bands, we have to consider only the matrix elements for ˆ e and ˆ e . Introducing Pauli matrices ν torepresent the two Fermi surfaces in the two degenerate conduction bands, the projected Zeeman term takes the form H Z ≈ (cid:90) d k f † k [ b k · ν ] f k , (S9)where b k = ∆ Z (cid:0) − cos α k sin θ k , , cos θ k (cid:1) .The signature of a vortex phase transition is the appearance of zero-energy vortex modes with zero momen-tum along the vortices. In what follows, we will thus consider k z = 0 ⇒ θ k = π/ b k = ∆ Z (cid:0) − cos α k , , (cid:1) ≡ − b x k ˆ x . The Hamiltonian, including the projected Zeeman coupling then reads[ H + H Z ] k z =0 ≈ (cid:90) d k f † k [( E k − µ ) ν − b x k ν x ] f k . (S10)It is evident that the Zeeman field mixes the modes of the two upper bands and can be diagonalized by the unitarytransformation d k = (cid:0) d + k , d − k (cid:1) T = e iν y π/ f k . Hamiltonian (S10) then becomes[ H + H Z ] k z =0 ≈ (cid:90) d k d † k [( E k − µ ) ν − b x k ν z ] d k . (S11) Effect of Zeeman coupling on intra-orbital singlet pairing
The uniform intra-orbital s-wave pairing has the form: H SC, = (cid:90) d r (cid:16) ψ †↑ A ∆ ψ †↓ A + ψ †↑ B ∆ ψ †↓ B (cid:17) + H.c. = 12 (cid:90) d r ψ † ∆ iσ y ψ † T + H.c.. (S12)In momentum space and projecting in the upper Fermi surfaces, H SC, = ∆ (cid:90) d k ψ † k iσ y ψ † T − k + H.c. ≈ ∆ (cid:88) a,b =1 , (cid:90) d k f † k ,a (cid:16) ϕ a † k iσ y ϕ b ∗− k (cid:17) f †− k ,b + H.c.. (S13)For a, b = 1 , k → − k results in H SC, = − ∆ (cid:90) d k (cid:16) f † k e − iφ k ν z f † T − k (cid:17) + H.c.. (S14)Rotating to the basis that also diagonalizes Zeeman coupling term we get[ H + H Z ] k z =0 ≈ (cid:90) d k d † k [( E k − µ ) ν − b x k ν z ] d k [ H SC, ] k z =0 ≈ − ∆ (cid:90) d k d † k (cid:0) e − iφ k ν x (cid:1) d † T − k + H.c.. (S15)The chemical potential will be crossing both conduction bands and leads to two split Fermi surfaces which haveopposite relative helicity. The intra-orbital singlet pairing can only pair states on different Fermi surfaces. Thissuggests, as we confirm below, that the intra-orbital singlet pairing is suppressed by Zeeman couplings.
Uniform orbital-triplet pairing
The different types of pairings allowed by symmetries in regular centro-symmetric topological insulators werepreviously studied in Ref. [1]. In particular, a triplet-type pairing deserves special attention here, which consists ofinter-orbital orbital-triplet cooper pairs. In our basis choice, it reads as H SC, = ∆ (cid:90) d k ψ † k iτ x σ y ψ † T − k + H.c. ≈ ∆ (cid:88) a,b =1 , (cid:90) d k f † k ,a (cid:16) ϕ a † k iτ x σ y ϕ b ∗− k (cid:17) f †− k ,b + H.c.. (S16)Using the unitary transformation ϕ a † k iτ x σ y ϕ b ∗− k = − e − iφ k sin α k ˆ e Ta ˆ e b , where sin α k = | k | √ k + m k , we get H SC, = − ∆ (cid:90) d k e − iφ k sin α k (cid:16) f † k ν f † T − k (cid:17) + H.c.. (S17)This pairing term commutes with the Zeeman term and so can be simultaneously diagonalized, returning[ H + H Z ] k z =0 ≈ (cid:90) d k d † k [( E k − µ ) ν − b x k ν z ] d k [ H SC, ] k z =0 ≈ − ∆ (cid:90) d k e − iφ k sin α k (cid:16) d † k ν d † T − k (cid:17) + H.c.. (S18)It is clear that, contrary to the intra-orbital singlet pairing, the inter-orbital triplet superconductivity pairs solelythe states on each Zeeman-split helical Fermi surface, exclusively. As a result, the intra-orbital triplet pairing is notsuppressed by the Zeeman coupling.0
Vortex Hamiltonian in inter-orbital triplet paired state
Having in hands a type of pairing that is not suppressed by the Zeeman-induced splitting of the Fermi surfaces,we can derive the the spectrum of vortex bound states. We keep the axial symmetry of the problem and introduce avortex along the z direction, still focusing on k z = 0. We have (cid:2) H vrtxSC, (cid:3) k z =0 = 12 (cid:90) d r ψ † ∆ ( r ) iτ x σ y ψ † T + H.c., (S19)where ∆ ( r ) = ∆ ξ ( x + iy ) . (S20)In momentum space and in polar coordinates (cid:18) ∂ k x ∂ k y (cid:19) = (cid:18) cos φ k − sin φ k sin φ k cos φ k (cid:19) (cid:18) ∂ k k ∂ φ k (cid:19) . (S21)Projecting the vortex Hamiltonian into the conduction bands we get (cid:2) H vrtxSC, (cid:3) k z =0 = i ξ (cid:90) d k ψ † k (cid:0) ∂ k x + i∂ k y (cid:1) iτ x σ y ψ † T − k + H.c. ≈ i ξ (cid:88) a,b =1 , (cid:90) d k f † k ,a (cid:20) e iφ k ϕ a † k (cid:18) ∂ k + i k ∂ φ k (cid:19) iτ x σ y ϕ b ∗− k (cid:21) f †− k ,b + H.c.. (S22)The matrix elements associated with the Bloch wave-functions on each band read ϕ a † k (cid:18) ∂ k + i k ∂ φ k (cid:19) iτ x σ y ϕ b ∗− k (S23)= (cid:2) i ( A k ) ab + i ( A φ k ) ab (cid:3) + ϕ a † k iτ x σ y ϕ b ∗− k (cid:18) ∂ k + ik ∂ φ k (cid:19) where i [ A k ] ab ≡ ϕ a † k iτ x σ y (cid:0) ∂ k ϕ b ∗− k (cid:1) i [ A φ k ] ab ≡ ϕ a † k iτ x σ y (cid:18) ik ∂ φ k ϕ b ∗− k (cid:19) . (S24)Since we solely consider states close to the Fermi surfaces in conduction bands a, b = 1 ,
2, the Berry connection-liketerms read i ( A k ) ab = − e − iφ k ∂ k sin α k (cid:0) ˆ e Ta ˆ e b (cid:1) i ( A φ k ) ab = − e − iφ k sin α k k (cid:0) ˆ e Ta ˆ e b (cid:1) . (S25)The vortex Hamiltonian then takes the form (cid:2) H vrtxSC, (cid:3) k z =0 = − ∆ ξ (cid:90) d k f † k iν sin α k (cid:20) ∂ k + A k + i∂ φ k + A φ k k (cid:21) f † T − k + H.c. (S26)with A k = ∂ k (cid:16) log (cid:112) sin α k (cid:17) ν A φ k = ν . (S27)1Notice that the radial component of the Berry connection is a pure gauge and can be neglected. Furthermore, theangular part vanishes by Fermi statistics. Going to the Zeeman diagonal basis is again trivial and the Hamiltonianreduces to (cid:2) H vrtxSC, (cid:3) k z =0 = − ∆ ξ (cid:90) d k d † k iν sin α k (cid:20) ∂ k + i∂ φ k k (cid:21) d † T − k + H.c.. (S28)
Spectrum of vortex modes
We are now ready to fully determine the spectrum of the Caroli-de Gennes modes, as well as their correspondingwavefunctions. Introducing Nambu spinors as Ψ k = (cid:18) d k d † T − k (cid:19) , (S29)the Hamiltonian at k z = 0 reads H k z =0 T ISC = 12 (cid:90) k Ψ † k ( E k − µ ) ν − h k ν z − i ∆ ξ sin α k ν (cid:104) ∂ k + i∂ φ k k (cid:105) − i ∆ ξ sin α k ν (cid:104) ∂ k − i∂ φ k k (cid:105) − ( E k − µ ) ν + h k ν z Ψ k . (S30)We can preform a mode expansion as Ψ k = (cid:88) n U n,k ˜Ψ n,k (S31)where U n,k = e − inφ k √ π i i . (S32)The Hamiltonian then takes the form H = 12 (cid:88) n,n (cid:48) (cid:90) dk π ˜Ψ † n,k (cid:20)(cid:90) dφ k π kU † n,k h BdG U n (cid:48) ,k (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) ˜ h BdG δ nn (cid:48) ˜Ψ n (cid:48) ,k , (S33)where ˜ h BdG = (cid:32) ( E k − µ ) ν − b x k ν z ∆ ξ sin α k ν (cid:0) ∂ k + nk (cid:1) ∆ ξ sin α k ν (cid:0) − ∂ k + nk (cid:1) − ( E k − µ ) ν + b x k ν z (cid:33) . (S34)We have two decoupled ν sectors: ˜Ψ n,k = (cid:88) ν = ± Θ νn,k a νn,k , (S35)where Θ + n,k = u + n ( k )0 v + n ( k )0 , Θ − n,k = u − n ( k )0 v − n ( k ) (S36)2and ˜ h νBdG Θ νn,k = E νn Θ νn,k . (S37)The Nambu constraint enforces u νn ( k ) a νn,k = − i ( − n v ν ∗− n ( k ) a ν †− n,k . (S38)Introducing Π Pauli matrices in the Nambu space, the BdG Hamiltonian for each set of modes reads˜ h νBdG = Π z ( k − b x k ν − µ ) (cid:124) (cid:123)(cid:122) (cid:125) ≡ ∆ λ ν ( k ) /ξ +Π y ∆ ξ ( i∂ k ) + Π x ∆ ξ (cid:16) nk (cid:17) (S39)and we can solve for the lowest energy eigenstates by a Jackiw-Rebbi argument for each ν . The first two terms abovehave a zero-energy mode when k νF − h k νF ν = µ so that we expect to have the lowest energy states localized close toeach Fermi surface. For these modes E νn = ∆ ξ k νc n, n ≥ . (S40)A zero-energy solution always exists for each ν . The negative n modes are not independent from the positive n (whichcompensates the factor of 1 / H ).Explicitly, the wavefunctions are fixed by[Π z λ ν ( k ) + Π y ( i∂ k )] Θ νn,k = 0 ⇒ [ ∂ k + Π x λ ν ( k )] Θ νn,k = 0 , (S41)with solutions Θ νn,k ≈ e − (cid:82) kνF + kkνF − k dk (cid:48) λ ν ( k (cid:48) ) √ N (cid:18) (cid:19) , (S42)where N is a normalization factor from the radial momentum integration. In other words, u νn ( k ) = v νn ( k ) = e − (cid:82) kνF + kkνF − k dk (cid:48) λ ν ( k (cid:48) ) √ N (S43)for the low energy modes. Also note the exponential localization of the wavefunctions at k νF . Generically we willapproximate the solutions by evaluating them at the Fermi momenta. Notice that, since the wavefunctions for thelow energy modes are independent of n and real, the Nambu constraint reduces to a ν − n,k = − i ( − n a ν † n,k ⇒ a ν ,k = − ia ν † ,k . (S44)Absorbing a π/ n = 0 modes corresponds to Majorana fermions.We arrive at the expansion d ν k = (cid:88) n e − inφ k √ πk u νn ( k ) a νn,k ≈ f ν ( k ) (cid:88) n (cid:0) e − inφ k a νn,k (cid:1) , (S45)with f ν ( k ) ≡ e − (cid:82) kνF + kkνF − k dk (cid:48) λ ν ( k (cid:48) ) √ π N k . ψ k ≈ (cid:88) a =1 ϕ a k (cid:16) e − iν y π/ d k (cid:17) a ≈ √ (cid:2) ϕ k (cid:0) d + k − d − k (cid:1) + ϕ k (cid:0) d + k + d − k (cid:1)(cid:3) . (S46)Using ϕ a k at k z = 0 and the expansion of the d ν k operators, we obtain ψ k ≈ (cid:88) n e − i ( n − φ k cos α k e − inφ k sin α k f + ( k ) a + n,k − e − inφ k cos α k e − i ( n − φ k sin α k f − ( k ) a − n,k . (S47)We finish by going to real space by Fourier transforming, ψ ( r ) = (cid:82) k e i r · k ψ k . The angular integrals introduce Besselfunctions; using the definitions of the α k angles and the exponential localization of the wavefunctions at the Fermisurfaces, we finally arrive at ψ ( r ) ≈ (cid:88) ν =1 , (cid:88) l c νl χ νl ( r ) a νl,k F , where χ + l ( r ) = 1 (cid:113) N + k + F e − i ( l − θ J l − (cid:0) k + F r (cid:1) (cid:16) m k + F + E k + F (cid:17) e − ilθ J l (cid:0) k + F r (cid:1) k + F χ − l ( r ) = 1 (cid:113) N − k − F e − ilθ J l (cid:0) k − F r (cid:1) k − F e − i ( l − θ J l − (cid:0) k − F r (cid:1) (cid:16) m k − F − E k − F (cid:17) , (S48)and c νl = ( − l k νF f ν ( k νF ).
1D VORTEX HAMILTONIAN
In order to explicitly examine the structure of vortex modes localized at the end of the vortices on the surface, weneed to derive the effective vortex Hamiltonian which can be applied for k z (cid:54) = 0. With the wavefunctions derived inlast section, finite k z may be considered perturbatively. Considering the linear k z term in H k z ≈ − iα z ∂ z , we obtain (cid:90) d rχ − l (cid:48) ( r ) τ x σ z χ + l ( r )= δ l,l (cid:48) (cid:90) rdr (cid:104) J l (cid:0) k + F r (cid:1) J l (cid:0) k − F r (cid:1) k + F k − F − J l − (cid:0) k + F r (cid:1) J l − (cid:0) k − F r (cid:1) (cid:16) m k + F + E k + F (cid:17) (cid:16) m k − F − E k − F (cid:17)(cid:105) ≡ ˜∆ l (cid:54) = 0 . (S49)So, H k z ≈ − i ˜∆ η x ∂ z . (S50)where µ i are the Pauli matrices acting in the space of the two states given in equation (S48). As outlined in the maintext, this linear gapless Hamiltonian does not support localized states at the ends of the vortex at zero energy.To derive the vortex Hamiltonian for higher energy vortex modes, we have to consider the other matrix elementassociated with the term (cid:15)k z β . The term is diagonal in the bases of states in equation (S48) and its matrix elementsare given by4 (cid:15) + = (cid:90) d rχ + l (cid:48) ( r ) τ z σ χ + l ( r )= δ l,l (cid:48) (cid:90) rdr (cid:104) J l − ( k + F r )( m k + F + E k + F ) − J l ( k + F r ) k + F (cid:105) (S51) (cid:15) − = (cid:90) d rχ − l (cid:48) ( r ) τ z σ χ − l ( r )= δ l,l (cid:48) (cid:90) rdr (cid:104) − J l − ( k − F r )( m k − F − E k − F ) + J l ( k − F r ) k − F (cid:105) (S52)Note that (cid:15) + > (cid:15) − . In fact, for chemical potential close to regions where m k F ≈ (cid:15) + ≈ − (cid:15) − >
0. As a result, the second derivativeterm, leads to a a term of the approximate form − ¯ (cid:15)η z ∂ z (S53)where ¯ (cid:15) has the same sign as (cid:15) .The the energy of vortex modes with finite energy add naturally to the 1D model as E + n + E − n η + E + n − E − n η z . Thestability of finite energy states at the end of vortex would then dependent on the sign of E + n − E − n . This can be readilychecked from equation (S40), which gives an inverse proportionality of finite energies with corresponding Fermi wavevector. The Fermi wave vector is given by k νF = µ + h k νF ν where h k F ν = − ∆ z m kνF µ . It is then readily clear that for m k > h µk F < k + F < k − F , leading to E + n > E − n . As is outlined in the main letter, this relationship leads tostability of finite energy states. EFFECT OF ZEEMAN COUPLING ON THE TYPE OF SUPERCONDUCTING STATE
For a careful study of the energetics of the problem and which pairing symmetry is favored in the presence of aZeeman coupling, we perform a mean-field self-consistent analysis. For this, we start with the wavefunctions at finite k z and in the presence of Zeeman coupling. This is slightly more involved.Now we consider the full TI plus Zeeman Hamiltonian with all momenta H = H T I + H Z . (S54)Using the usual squaring trick on H ≡ H k + ∆ Z σ z , we may find the spectrum to be ± E ± , k = ± (cid:115) E k + (cid:18) ∆ Z (cid:19) ± | ∆ Z | (cid:113) k z + m k = ± (cid:115) k ⊥ + (cid:18) ξ k ± ∆ Z (cid:19) , (S55)where we identified the TI energy spectrum E k ≡ (cid:112) k + m k but also defined ξ k = (cid:112) k z + m k and rewrote theexpression in a way that may suggest us some hints on how to proceed. First let us rotate away the in-planemomenta, H (cid:48) = e iσ z φ k / He − iσ z φ k / = τ x σ x k ⊥ + τ x σ z k z + m k τ z σ + ∆ Z τ σ z . (S56)Then we rotate k z and mass terms, τ x σ z k z + m k τ z σ = ξ k τ z σ e iτ y σ z β k tan β k = k z m k (S57)So H (cid:48)(cid:48) = e iτ y σ z β k / H (cid:48) e − iτ y σ z β k = τ x σ x k ⊥ + ξ k τ z σ + ∆ Z τ σ z . (S58)5Now since we rotated away k z , H (cid:48)(cid:48) now commutes with the mirror operation M = τ z σ z . (S59)So we can project with P ν = 1 + ν M , (S60)where ν = ± . We obtain H (cid:48)(cid:48) ν ≡ P ν H (cid:48)(cid:48) P ν = (cid:18) ν M (cid:19) (cid:18) τ x σ x k ⊥ + ξ k τ z σ + ∆ Z τ σ z (cid:19) (cid:18) ν M (cid:19) = (cid:18) τ x σ x − ντ y σ y (cid:19) k ⊥ + (cid:18) ξ k + ν ∆ Z (cid:19) (cid:18) τ z σ + ντ σ z (cid:19) ≡ Λ νx k ⊥ + Λ νz (cid:18) ξ k + ν h (cid:19) where we introduced the Λ ν × H (cid:48)(cid:48) ν = Λ νz (cid:20) Λ νz Λ νx k ⊥ + (cid:18) ξ k + ν ∆ Z (cid:19)(cid:21) = E ν, k Λ νz e i Λ νy γ ν k , (S61)where tan γ ν k ≡ k ⊥ ξ k + ν h (S62)and, by definition, Λ νz Λ νx = i Λ νy = i τ y σ x + ντ x σ y ) . (S63)The wavefunctions now read ϕ a ν ,ν k = e − iσ z φ k / e − iτ y σ z β k / P ν e − i (cid:16) τyσx + ντxσy (cid:17) γ ν k e a ν . (S64)Notice that the a variable here is enslaved by ν . For ν = +, a + = 1 , ν = − , a − = 2 , − a ν +1 E ν, k .Considering only the the positive energy bands we define the corresponding projection operators P ν k = 14 ( τ σ + ν cos β k τ z σ z + ν sin β k τ x σ + cos β k cos γ ν k τ z σ + sin β k cos γ ν k τ x σ z + ν cos γ ν k τ σ z + cos φ k sin γ ν k τ x σ x + sin φ k sin γ ν k τ x σ y − ν cos φ k cos β k sin γ ν k τ y σ y + ν sin φ k cos β k sin γ ν k τ y σ x + ν cos φ k sin β k sin γ ν k τ σ x + ν sin φ k sin β k sin γ ν k τ σ y ) . (S65)The single particle Green’s function in the normal state for two positive bands can be written as G ( iω n , k ) = P + k iω n − (cid:15) + , k + P − k iω n − (cid:15) − , k , (S66)where ω n = (2 n + 1) πT are the Matsubara frequencies and (cid:15) ν, k = E ν, k − µ ( µ is the chemical potential). In order towrite the linearized gap equation for the superconducting phase we define irreducible susceptibility as χ ij = − TN (cid:88) ω n , k Tr (cid:34) ˆ∆ i ∆ G ( iω n , k ) ˆ∆ j ∆ G ( − iω n , k ) (cid:35) , (S67)6where N is the number of unit cells, ˆ∆ i / ∆ defines the orbital and spin structure of different pairing potentials.Following the work of Ref. [1], we consider only intra-orbital s-wave pairing and inter-orbital time-reversal invariantpairing channels. We treat the Zeeman field as a small perturbation compared to critical temperatures, so the phasediagram will still involve only these two phases, as was the case for ∆ Z = 0. For these, the orbital and spin structurehave the following form ˆ∆ a = ∆ τ σ , and ˆ∆ b = ∆ τ z σ (S68)ˆ∆ = ∆ τ x σ . (S69)Here, ˆ∆ a and ˆ∆ b are two different intra-orbital components for the s-wave channel and ˆ∆ is the structure ofinter-orbital pairing.In order to evaluate the susceptibilities we make the following approximations: 1) that despite the presence of theZeeman term the Fermi surface is approximately isotropic, 2) that the Zeeman field strength is small and expand theform of the energies of two Fermi surfaces (S55) E ± , k ≈ E k ± ∆ Z ξ k E k , (S70)where E k = (cid:112) k ⊥ + ξ k is the energy of the bands when Zeeman term is zero. Now the integration with k of (S67)can be carried for the case with Zeeman term. The only difference is that the Fermi energy where the integral ispeaked will correspond to different momenta and different density of states for Zeeman split bands. Assuming thatthe density of states is the same for both bands and using the isotropy condition of Fermi surfaces, the Fermi momentacan be determined from the energy equation (S55). The only subtlety here arises for the cross term case in (S67),which involves components from both bands. In that case we use the energy expansion (S70) in the denominator anddisregard the Zeeman term in the numerator. After that the sum for that case is similar to the case without Zeemanfield, so the momenta where it is peaked is determined from E k = µ . It should be noted that both the expansion ofthe energy (S70) and the neglecting of Zeeman terms in the numerator of the cross term are valid for ∆ Z (cid:28) µ and∆ Z (cid:28) T . After this, the susceptibilities related to ∆ and ∆ can be written as χ ≡ χ a b = ¯ χ (cid:20) m + k µ (cid:18) Z ξ + k (cid:19) + m − k µ (cid:18) − ∆ Z ξ − k (cid:19)(cid:21) (S71) χ ≡ χ b b = ¯ χ ( m + k ) µ (cid:18) Z ξ + k (cid:19) + ( m − k ) µ (cid:18) − ∆ Z ξ − k (cid:19) + (cid:18) k ⊥ ξ k (cid:19) − µ − ∆ Z (cid:16) µ + hξ k µ (cid:17) (cid:16) µ − ∆ Z ξ k µ (cid:17) (S72) χ ≡ χ = ¯ χ (cid:18) k + z ξ + k (cid:19) + (cid:18) k − z ξ − k (cid:19) + (cid:18) m k ξ k (cid:19) − (cid:0) ξ k + ∆ Z (cid:1) (cid:0) ξ k − ∆ Z (cid:1) − ( k ⊥ ) (cid:16) µ + ∆ Z ξ k µ (cid:17) (cid:16) µ − ∆ Z ξ k µ (cid:17) . (S73)The superscripts + , − and 0 denote that the corresponding quantities are evaluated at the energy E + , k = µ , E − , k = µ and E k = µ , respectively. ¯ χ = − (cid:82) w D − w D D ( ξ ) tanh ( ξ/ T ) / ξdξ is the standard s -wave susceptibility, D ( ξ ) is the densityof states and w D is the Debye frequency. Plugging these forms of susceptibilities into linearized gap equations fors-wave and inter-orbital triplet pairings [1–3]det (cid:12)(cid:12)(cid:12)(cid:12) U ¯ χ − U χ V χ V χ − (cid:12)(cid:12)(cid:12)(cid:12) = 1 , V χ = 1 , (S74)we obtain the resulting phase boundary in this case (see Fig. 1 of the main text). In (S74) U and V describe intra-and inter-orbital interactions of electrons, respectively [1–3]. [1] L. Fu and E. Berg, Phys. Rev. Lett. , 097001 (2010).[2] S. Nakosai, Y. Tanaka, and N. Nagaosa, Phys. Rev. Lett. , 147003 (2012).[3] T. Hashimoto, S. Kobayashi, Y. Tanaka, and M. Sato, Phys. Rev. B94