Analytical and semianalytical tools to determine the topological character of Shiba chains
NNew tools to determine the topological character of Shiba chains
Nicholas Sedlmayr, ∗ Vardan Kaladzhyan, and Cristina Bena Institute of Physics, Maria Curie-Sk(cid:32)lodowska University,Plac Marii Sk(cid:32)lodowskiej-Curie 1, PL-20031 Lublin, Poland Department of Physics, University of Basel, Klingelbergstrasse 82, CH-4056 Basel, Switzerland Institut de Physique Th´eorique, Universit´e Paris Saclay,CEA CNRS, Orme des Merisiers, 91190 Gif-sur-Yvette Cedex, France (Dated: February 10, 2021)We introduce three new analytical and semi-analytical tools that allow one to determine thetopological character of impurity Shiba chains. The analytical methods are based on calculatingthe effective Green’s function of an infinite embedded chain using the T-matrix formalism. We thusprovide a solution to the longstanding size-effects problem affecting the only general alternativemethod, the numerical tight-binding analysis. As an example we consider a chain of magnetic im-purities deposited on an s-wave superconducting substrate with Rashba spin-orbit and we calculateits topological phase diagram as a function of the magnetic impurity strength and the chemicalpotential. We find a perfect agreement between all our new techniques and a numerical analysis.
Introduction.
Magnetic impurities in conventional super-conductors give rise to impurity-bound subgap states, theso-called Yu-Shiba-Rusinov states [1–5]. The topologi-cal properties of chains of such magnetic impurities, alsoknown as Shiba chains, have come into focus in recentyears. Since the pioneering theoretical proposals [6–10]these chains became one of the most prominent platformsfor Majorana bound states. Various alternatives havebeen considered theoretically [11–33], while experimen-tal work has mostly focused on magnetic islands [34, 35]and on emulating magnetic Shiba chains on top of su-perconducting substrates [36–42]. In the latter experi-ments, zero-bias peaks observed in the tunneling spectraat the ends of magnetic-atom chains are assumed to cor-respond to Majorana bound states. The true nature ofthese states remains, however, open to question, since thepresence of zero-bias peaks is not an unambiguous proofof the topological character of the underlying system.Previously, the only general versatile tool to deter-mine the topological character of an impurity chain wasa numerical tight-binding analysis. Other theoreticalapproaches proposed were not so generally applicable:among the few alternative proposals one can mentiondescribing the fine-tuned case of dilute impurities cor-responding to Shiba states forming very close to zero en-ergy [7], and approaches based either on Pfaffian invari-ants for effective one-dimensional models [11, 33, 36, 41,43] or on the zeros of the Green’s function in the presenceof impurities [44]. However, the tight-binding method,when applied to a large but finite two-dimensional sys-tem with a modified charge/spin chemical potential alonga finite chain of atoms, exhibits size-effects problems,and thus the computing times and resources necessaryto model the chain accurately are very large. This is dueespecially to the interactions between the states formingat the ends of the chain, as well as due to the leakage ofthe latter into the two-dimensional bulk [19, 45].In what follows we introduce two analytical tools and
Figure 1. Schematics of the systems considered. Panel (a)shows a chain of magnetic impurity atoms (green) on the sur-face of the superconductor (dark blue). For I (see main text)the system is infinite in both directions. For III the surfaceis finite in the x -direction, and infinite in the y -direction. Inpanel (b) the red atom is an additional probe point-impurity,as described in II in the main text. All numerical tight-binding solutions are performed for systems finite in bothdirections. one semi-analytical tool that allow one to bypass the size-effects problem and provide us with at least one fullygeneral and size-independent approach to test if a givenimpurity chain is topological. To exemplify our methodwe focus on a two-dimensional s -wave superconductingsubstrate with Rashba spin-orbit interaction and a chainof magnetic impurities deposited on top.I. We consider an infinite two-dimensional bulk andwe model the Shiba chain as an infinite impurity line de-scribed by a finite magnetic potential along the centralline of atoms (see Fig. 1). We solve this system exactlyusing the T -matrix formalism [46–53]. We thus obtainthe effective Green’s function for the impurity wire. Fi-nally, using the latter we calculate the standard topolog-ical chiral invariant [54] enabling us to predict in whichparameter range the Shiba chain is topologically non-trivial.II. In the second approach we use the effective Green’sfunction obtained within the first approach, but we testthe formation of Majorana bound states by adding an ex-tra probe point-impurity on the wire, see Fig. 1(b). Wemodel the impurity as an infinite-amplitude scalar poten- a r X i v : . [ c ond - m a t . s up r- c on ] F e b tial that effectively cuts the wire in two. We check theformation of Majorana states by calculating the averagezero-energy density of states in the wire in the presence ofthe probe impurity. If such states form, the zero-energydensity of states in the wire is finite, while if the sys-tem is nontopological it is equal to zero. Additionally,one needs to check that no zero-energy states exist in thewire in the absence of the probe impurity.III. We compare these methods with a semi-analyticalapproach based on Ref. [55] for calculating the Z in-variant of quasi-1D systems. We consider a quasi-one-dimensional system, infinite in one direction and finitein the other. The Shiba chain is modeled the same wayas in I, see Fig. 1(a) but now our system is finite in the x direction. We consider systems of different width; forsufficiently wide systems the invariant no longer dependsexplicitly on the width.The first two approaches are fully analytical and arebereft of finite-size problems, rendering the problemmuch more tractable and accurate than tight-bindingsimulations. The effective Green’s function for the wirecan be obtained using at most a numerical integral inmomentum space, and in some cases, such as the onepresented here, a closed-form expression can be derivedwithout resorting to numerical integration.The third approach is semi-analytical in the sense thatwe have a fully analytical formula for the topological in-variant, but its complexity depends on the transverse sizeof the system, and thus for large widths it can be evalu-ated only numerically.The main advantage of our approaches is that theydo not depend on the size of the system, thus provid-ing reliable tools that are not perturbed by finite-size ef-fects. Furthermore, while the first and third approachesrequire some basic symmetry analysis in order to be im-plemented, the second one is fully general and can be ef-fortlessly implemented and straightforwardly applied toany system. Model.
We consider a two-dimensional square latticewith Rashba spin-orbit interaction and s -wave supercon-ductivity, which can be described in real space by thefollowing tight-binding Hamiltonian H = (cid:88) r Ψ † r ( − µ τ z − ∆ τ x ) Ψ r − Ψ † r ( t + i λ σ y ) τ z Ψ r + x + H . c . − Ψ † r ( t − i λ σ x ) τ z Ψ r + y + H . c . , (1)written in the basis Ψ r = ( c r , ↑ , c r , ↓ , c † r , ↓ , − c † r , ↑ ) T , where c ( † ) r ,σ annihilates (creates) an electron with spin σ = ↑ , ↓ atsite r = ( x, y ). We denote by µ the chemical potential, ∆is the superconducting pairing parameter, t the hoppingconstant, and λ the Rashba spin-orbit coupling constant.The Pauli matrices σ and τ act in the spin and particle-hole subspaces, respectively. In momentum space, in the basis Ψ k = ( c k , ↑ , c k , ↓ , c †− k , ↓ , − c †− k , ↑ ) T , the Hamiltonianof this system is given by H k = − [ µ + 2 t (cos k x + cos k y )] τ z (2) − λ [sin k x σ y − sin k y σ x ] τ z − ∆ τ x . The Shiba chain is modeled as a set of onsite potentialscorresponding to local Zeeman fields of the form V = V σ z , H imp = (cid:88) r ∈ C Ψ † r V Ψ r , (3)with C describing the sites of the one-dimensional infi-nite chain. Therefore, the full system in real space isdescribed by H + H imp . Effective Green’s function for the chain.
The effectiveone-dimensional Green’s function G for the chain canbe computed using the technique presented in Ref. [56].Thus, for a Hamiltonian H k we can define the bulk un-perturbed Green’s function G ( ω, k ) = (i ω + H k ) − . Inthe presence an impurity potential V δ ( x ) we can com-pute an effective Green’s function with the help of the T -matrix [46–50, 52]: T = [ I − V G ( ω, k y )] − V , (4)where G ( ω, k y ) ≡ G ( ω, x = 0 , k y ) = (cid:90) π − π dk x π G ( ω, k ) . (5)In what follows we are interested in evaluating the effec-tive Green’s function at x = 0, i.e. along the impuritychain. In the T-matrix formalism this can be written as: G ( ω, k y ) = G ( ω, k y ) + G ( ω, k y ) T G ( ω, k y ) . (6)where G ( ω, k y ) at ω = 0 can be obtained either byperforming the integral in Eq. (5) numerically, or forthe system considered here can also be computed ana-lytically (see Eq. (S49) in the Supplementary Material).Hereinafter we set k y → k where possible. This Green’sfunction defines an effective Hamiltonian for the chain: H − k ≡ G (0 , k ). In what follows we apply the three dif-ferent methods mentioned in the introduction, and cal-culate the topological phase diagram for the magneticShiba chains. I. The chiral invariant from the Green’s function.
Inorder to calculate the appropriate topological invariantfrom the effective Green’s function we first investigatethe symmetries of the problem. The model we con-sider exhibits several symmetries important for the topol-ogy [57, 58]: a particle–hole (PH) symmetry describedby an anti-unitary operator P , such that PH ∗ k P − = −H − k and P = 1; a “time-reversal” (TR) symme-try described by the anti-unitary operator T , where T H ∗ k T − = H − k with T = 1; and finally we have - ν GF / t μ / t ( a ) δρδρ max
0. 1. 2. 3. 4.0.1.2.3.4.V / t μ / t ( b )
0. 1. 3. 4.0.2.3.4. V / t μ / t δ N ( ) δ N max / t μ / t ( c ) ϵ / t / t μ / t ( d ) Figure 2. (a) The chiral invariant, ν , plotted as a function ofthe chemical potential and the magnetic impurity strength,as obtained from the effective wire GF. (b) The zero-energyDOS, ρ , in the presence of a probe point-impurity; shown inthe inset is the total density of states in the absence of theprobe. (c,d) The Majorana polarization C and the energycorresponding to the lowest-energy state calculated using atight-binding model for a chain of 60 sites introduced into atwo-dimensional system of size 21 ×
80. We take λ = 0 . t and∆ = 0 . t the combination of these symmetries, the sublattice or“chiral” symmetry, described by the unitary operator S = PT , with S − H k S = −H k . When all of thesesymmetries are present the Hamiltonian is in the BDIsymmetry class [59], and hence has a Z topological invari-ant. Our model has a chiral symmetry operator for theeffective chain Hamiltonian given by S = τ y σ x . Note,however that this is not a symmetry of the full prob-lem, H + H imp , as the term λ sin k x σ y τ z breaks theTR symmetry. We note here that the two-dimensionallattice described by H is always topologically trivial.Therefore, we can use the topological invariant [54], ν = 14 π i (cid:90) π − π dk tr SH k ∂ k H − k , (7)suitable for a one-dimensional system with chiral sym-metry.The result for the chiral invariant as a function of thechemical potential and the magnetic impurity strengthfor an exemplary set of parameters (∆ = 0 . t and λ = 0 . t ) is given in Fig. 2(a). One can see that a largetopologically non-trivial phase exists, in particular forlarge values of V . II. Testing the formation of Majorana bound states via a probe impurity.
Our second approach is also basedon the effective Green’s function obtained in Eq. (6). Totest the formation of Majorana states we add an infinite-amplitude scalar probe impurity which breaks the wireinto two semi-infinite parts, see Fig. 1(b) [60, 61]. Theperturbed Green’s function of the wire in the presenceof the probe impurity can be obtained via the T -matrixformalism using the effective Green’s function of the in-finite wire in Eq. (6) as an unperturbed Green’s func-tion describing the infinite chain. Assuming that theprobe impurity is localized at y = 0 and is describedby U diag { , , , } δ ( y ), we can find the T -matrix andthe correction to the Green’s function as follows:˜ T = (cid:20) I − U (cid:90) π − π dk π G (0 , k ) (cid:21) − U (8)The local density of states is given by δ ˜ ρ (0 , y ) = − π Im tr G (0 , y ) ˜ T G (0 , − y ) , (9)where G (0 , y ) = (cid:82) π − π dk π G (0 , k ) e i ky . To verify the forma-tion of Majorana states, in Fig. 2(b) we plot the zero-energy local density of states in the presence of the infi-nite impurity, as a function of the chemical potential andthe magnetic impurity strength in the chain. To checkthat these states are indeed forming inside a bulk gapand are not just regular bulk states, in the inset we plotthe zero-energy total DOS in the bulk of the wire: δN (0) = (cid:90) dy δ ˜ ρ (0 , y ) . (10)We confirm that this is indeed negligible except along thebulk-gap closing lines. Thus the states that we observeat zero-energy do not correspond to any bulk states, butcome from adding the impurity and are thus Majoranabound states. This is in very good agreement with thechiral invariant in Fig. 2(a).Our results are checked against standard tight-bindingcalculations, thus in Figs. 2(c,d) we consider a finite-sizesystem of dimensions 21 ×
80, with a magnetic impu-rity chain embedded in it, and we solve the tight-bindingmodel in real space numerically. We focus on the lowest-energy state, and we calculate the corresponding Majo-rana polarization C [55, 62–66], i.e., the normalized ex-pectation value of the particle-hole operator integratedover half the wire. For a localized Majorana bound state,i.e. an eigenstate of P , we should have C = 1, and thisis indeed what we observe in Fig. 2(c). We note thatthere is a perfect agreement between our two new meth-ods introduced above, and the numerical tight-bindinganalysis: the topological regions in the phase diagram aspredicted by our methods are exactly the same as thosethat exhibit a perfect Majorana polarization ( C = 1) anda zero-energy state in the numerical tight-binding calcu-lation. III. Invariant from the quasi-one-dimensional strip.
Wenow consider that the system described by H + H imp has a finite number of sites along the x -direction, N x .The resulting system is a quasi-1D strip described bythe Hamiltonian H Q1D , with the finite-size segment inthe x -direction constituting the unit cell. This configu-ration breaks TR symmetry and so the system belongsto the class D with a Z invariant. This invariant canbe calculated by considering the parity of the bandsat high symmetry points, which encodes the band in-version [67]. At the time reversal invariant momenta,following a rotation the Hamiltonian can be written as H Q1D0 ,π → diag ( ¯ H ,π , − ¯ H ,π ), and calculating the parityof the negative energy bands becomes equivalent to cal-culating δ = ( − ν Q1D = sgn (cid:2) det ¯ H det ¯ H π (cid:3) . FollowingRef. [55] a closed form expression can be found for δ . As ν Q1D is a Z invariant we can unambiguously write that ν Q1D = 1 − sgn (cid:2) det ¯ H (0) det ¯ H ( π ) (cid:3) . (11)More details can be found in the Supplemental Material.The quasi-one-dimensional invariant works very well,and its predictions agree perfectly with the previous cal-culations for the topological phase diagram presented inFig. 2, provided N x is large enough. In Fig. 3(a) weplot the quasi-one-dimensional invariant as a function ofchemical potential and magnetic impurity strength for aa strip of width of N x = 51, while in the right panel weconsider N x = 11. Note that the oscillations present atsmall widths go away for N x = 51. ν Q1D
0. 1. 2. 3. 4.0.1.2.3.4.V / t μ / t ( a ) ν Q1D
0. 1. 2. 3. 4.0.1.2.3.4.V / t μ / t ( b ) Figure 3. The quasi-one-dimensional invariant, ν Q1D , for (a) N x = 51 and (b) N x = 11. Comparison with experiments.
To facilitate the numeri-cal computations, so far we have used unrealistically largevalues of the parameters. In what follows we consider aset of parameters expected to correspond more closelyto realistic experimental parameters. In the experimentsthe value of the superconducting gap varies from 1 . . . λ = ∆ = 0 . t ∼ ν Q1D
0. 1. 2. 3. 4.0.1.2.3.4.V / t μ / t ( a )
0. 0.25 0.5 0.75 1.0.51.52.V / t μ / t ( b )
3. 3.375 3.75 4.125 4.50.0.51.1.52.V / t μ / t ( c ) Figure 4. The quasi-one-dimensional invariant for λ = ∆ =0 . t for N x = 1001. Small panels (b,c) demonstrate theparts of the phase diagram close to µ/t = 0 and µ/t = 4 in ahigher resolution. but the results are comparable to those for larger ∆ and λ . The main differences come from the points close to µ = 0 and µ = 4 and for these points we performed de-tailed analyses in the panels (b,c) of Fig. 4. We notein particular that smaller values of the magnetic impu-rity potential are required close to µ = 4 to give riseto a topological phase. It appears that quite large mag-netic impurity potentials are necessary in the other re-gions of the phase diagram, which may not be alwayseasy to achieve experimentally. However, we note thatthis could be achieved in an alternative configuration, inwhich impurities are deposited above the surface atomsand connected with the substrate via a tunneling cou-pling t (cid:48) . We found that this configuration of suspendedatoms with an impurity potential U is exactly equivalentwith that of embedded impurities with an effective im-purity potential V = t (cid:48) /U . Thus larger values of theeffective impurity potential V can be reached, even fornot so large values of the bare U , see the SupplementalMaterial for more details. Conclusions.
We proposed three new analytical andsemi-analytical techniques to calculate the topologicalphase diagrams of Shiba chains. The first relies on calcu-lating the the effective Green’s function of the impuritychain using the T-matrix formalism and evaluating thecorresponding chiral invariant. The second technique isalso based on the effective GF of the chain, and consistsin calculating the zero-energy DOS in the presence of aprobe impurity. The third method is based on calculat-ing the quasi-one-dimensional invariant of a ribbon withan embedded impurity chain in the limit of infinite rib-bon width. We applied our techniques to the example ofa chain of magnetic impurities on an s -wave supercon-ducting substrate with Rashba spin-orbit coupling. Wecompared the results obtained with the newly introducedtools to those obtained via tight-binding numerical calcu-lations, and we found a perfect agreement. We discussedthe experimental relevance of our results. It will be in-teresting to apply our techniques to problems that aredifficult to solve correctly numerically due to size effects[19, 45].This work was supported by the National Sci-ence Centre (NCN, Poland) under the grant2019/35/B/ST3/03625 (NS). ∗ e-mail: [email protected][1] L. Yu, Acta Physica Sinica , 75 (1965).[2] H. Shiba, Progress of Theoretical Physics , 435 (1968).[3] A. I. Rusinov, Soviet Journal of Experimental and The-oretical Physics Letters , 85 (1969).[4] A. Yazdani, B. A. Jones, C. P. Lutz, M. F. Crommie, andD. M. Eigler, Science , 1767 (1997).[5] G. C. M´enard, S. Guissart, C. Brun, S. Pons, V. S.Stolyarov, F. Debontridder, M. V. Leclerc, E. Janod,L. Cario, D. Roditchev, P. Simon, and T. Cren, NaturePhysics , 1013 (2015).[6] S. Nadj-Perge, I. K. Drozdov, B. A. Bernevig, andA. Yazdani, Phys. Rev. B , 020407 (2013).[7] F. Pientka, L. I. Glazman, and F. Von Oppen, PhysicalReview B , 155420 (2013).[8] B. Braunecker and P. Simon, Phys. Rev. Lett. ,147202 (2013).[9] J. Klinovaja, P. Stano, A. Yazdani, and D. Loss, Phys.Rev. Lett. , 186805 (2013).[10] M. M. Vazifeh and M. Franz, Phys. Rev. Lett. ,206802 (2013).[11] F. Pientka, L. I. Glazman, and F. von Oppen, Phys. Rev.B , 180505 (2014).[12] K. P¨oyh¨onen, A. Weststr¨om, J. R¨ontynen, and T. Oja-nen, Phys. Rev. B , 115109 (2014).[13] I. Reis, D. J. J. Marchand, and M. Franz, Phys. Rev. B , 085124 (2014).[14] Y. Kim, M. Cheng, B. Bauer, R. M. Lutchyn, andS. Das Sarma, Phys. Rev. B , 060401 (2014).[15] J. Li, H. Chen, I. K. Drozdov, A. Yazdani, B. A.Bernevig, and A. H. MacDonald, Phys. Rev. B ,235433 (2014).[16] A. Heimes, P. Kotetes, and G. Sch¨on, Phys. Rev. B ,060507 (2014).[17] P. M. R. Brydon, S. Das Sarma, H.-Y. Hui, and J. D.Sau, Phys. Rev. B , 064505 (2015).[18] A. Weststr¨om, K. P¨oyh¨onen, and T. Ojanen, Phys. Rev.B , 064502 (2015).[19] Y. Peng, F. Pientka, L. I. Glazman, and F. Von Oppen,Physical Review Letters , 106801 (2015).[20] H.-Y. Hui, P. M. R. Brydon, J. D. Sau, S. Tewari, andS. D. Sarma, Scientific Reports , 10.1038/srep08880(2015).[21] J. R¨ontynen and T. Ojanen, Phys. Rev. Lett. , 236803(2015).[22] B. Braunecker and P. Simon, Phys. Rev. B , 241410(2015).[23] K. P¨oyh¨onen, A. Weststr¨om, and T. Ojanen, Phys. Rev.B , 014517 (2016).[24] J. Zhang, Y. Kim, E. Rossi, and R. M. Lutchyn, Phys.Rev. B , 024507 (2016).[25] J. Li, T. Neupert, B. A. Bernevig, and A. Yazdani, Na-ture Communications , 10.1038/ncomms10395 (2016). [26] J. R¨ontynen and T. Ojanen, Phys. Rev. B , 094521(2016).[27] S. Hoffman, J. Klinovaja, and D. Loss, Phys. Rev. B ,165418 (2016).[28] J. Li, T. Neupert, Z. Wang, A. H. MacDonald, A. Yaz-dani, and B. A. Bernevig, Nature Communications ,10.1038/ncomms12297 (2016).[29] M. Schecter, K. Flensberg, M. H. Christensen, B. M. An-dersen, and J. Paaske, Phys. Rev. B , 140503 (2016).[30] M. H. Christensen, M. Schecter, K. Flensberg, B. M. An-dersen, and J. Paaske, Phys. Rev. B , 144509 (2016).[31] V. Kaladzhyan, P. Simon, and M. Trif, Phys. Rev. B ,020507 (2017).[32] G. M. Andolina and P. Simon, Phys. Rev. B , 235411(2017).[33] A. Kobia(cid:32)lka, P. Piekarz, A. M. Ole´s, and A. Ptok, Phys.Rev. B , 205143 (2020).[34] G. C. M´enard, S. Guissart, C. Brun, R. T. Leriche,M. Trif, F. Debontridder, D. Demaille, D. Roditchev,P. Simon, and T. Cren, Nature Communications ,10.1038/s41467-017-02192-x (2017).[35] G. C. M´enard, C. Brun, R. Leriche, M. Trif, F. Debon-tridder, D. Demaille, D. Roditchev, P. Simon, andT. Cren, The European Physical Journal Special Topics , 2303 (2019).[36] S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon,J. Seo, A. H. MacDonald, B. A. Bernevig, and A. Yaz-dani, Science , 602 (2014).[37] R. Pawlak, M. Kisiel, J. Klinovaja, T. Meier, S. Kawai,T. Glatzel, D. Loss, and E. Meyer, npj Quantum Infor-mation , 10.1038/npjqi.2016.35 (2016).[38] M. Ruby, F. Pientka, Y. Peng, F. von Oppen, B. W.Heinrich, and K. J. Franke, Phys. Rev. Lett. , 197204(2015).[39] B. E. Feldman, M. T. Randeria, J. Li, S. Jeon, Y. Xie,Z. Wang, I. K. Drozdov, B. A. Bernevig, and A. Yazdani,Nature Physics , 286 (2016).[40] M. Ruby, B. W. Heinrich, Y. Peng, F. von Oppen, andK. J. Franke, Nano Letters , 4473 (2017).[41] H. Kim, A. Palacio-Morales, T. Posske, L. R´ozsa,K. Palot´as, L. Szunyogh, M. Thorwart, and R. Wiesen-danger, Science Advances , eaar5251 (2018).[42] R. Pawlak, S. Hoffman, J. Klinovaja, D. Loss, andE. Meyer, Progress in Particle and Nuclear Physics ,1 (2019).[43] A. Heimes, D. Mendler, and P. Kotetes, New Journal ofPhysics , 023051 (2015).[44] R.-J. Slager, L. Rademaker, J. Zaanen, and L. Balents,Phys. Rev. B , 085126 (2015).[45] T. D. Stanescu, R. M. Lutchyn, and S. Das Sarma, Phys.Rev. B , 144522 (2011).[46] J. M. Byers, M. E. Flatt´e, and D. J. Scalapino, Phys.Rev. Lett. , 3363 (1993).[47] M. I. Salkola, A. V. Balatsky, and D. J. Scalapino, Phys.Rev. Lett. , 1841 (1996).[48] W. Ziegler, D. Poilblanc, R. Preuss, W. Hanke, and D. J.Scalapino, Phys. Rev. B , 8704 (1996).[49] G. D. Mahan, Many-Particle Physics (Springer US,2000).[50] A. V. Balatsky, I. Vekhter, and J.-X. Zhu, Rev. Mod.Phys. , 373 (2006).[51] A. Zazunov, R. Egger, and A. Levy Yeyati, Phys. Rev. B , 014502 (2016).[52] C. Bena, Comptes Rendus Physique , 302 (2016). [53] M. Alvarado, A. Iks, A. Zazunov, R. Egger, and A. L.Yeyati, Phys. Rev. B , 094511 (2020).[54] V. Gurarie, Physical Review B , 085426 (2011).[55] N. Sedlmayr, J. M. Aguiar-Hualde, and C. Bena, PhysicalReview B , 155425 (2016).[56] S. Pinon, V. Kaladzhyan, and C. Bena, Phys. Rev. B , 115405 (2020).[57] C. K. Chiu, J. C. Teo, A. P. Schnyder, and S. Ryu, Re-views of Modern Physics , 35005 (2016).[58] M. Sato and Y. Ando, Reports on Progress in Physics ,aa6ac7 (2017).[59] S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Lud-wig, New Journal of Physics , 65010 (2010).[60] V. Kaladzhyan and C. Bena, Phys. Rev. B , 081106(2019).[61] A. Zazunov, R. Egger, and A. Levy Yeyati, Phys. Rev. B , 014502 (2016).[62] N. Sedlmayr and C. Bena, Phys. Rev. B , 115115(2015).[63] N. Sedlmayr, V. Kaladzhyan, C. Dutreix, and C. Bena,Phys. Rev. B , 184516 (2017).[64] V. Kaladzhyan, J. Despres, I. Mandal, and C. Bena, TheEuropean Physical Journal B , 10.1140/epjb/e2017-80103-y (2017).[65] V. Kaladzhyan, C. Bena, and P. Simon, Phys. Rev. B ,104512 (2018).[66] S. G(cid:32)lodzik, N. Sedlmayr, and T. Doma´nski, Phys. Rev.B , 085411 (2020).[67] M. Sato, Phys. Rev. B , 214526 (2009).[68] The value provided in the aforementioned reference is0 .
11 eV · ˚A. To convert it into the units of energy, wedivide this value by the lattice constant of lead ( a =4 .
95 ˚A), and we get λ ≈
22 meV. Supplemental Material
Quasi-one-dimensional invariant
We can write the Hamiltonian for H + H Ch as if the x direction is a unit cell, and making a Fourier transformalong y : H Q1D k = f k + V D C y L k − ∆ 0 L † k f k − V D C y − ∆ − ∆ 0 V D C y − f †− k L T − k − ∆ L ∗− k − f †− k − V D C y . (S1)The matrix D C y is a square matrix with a non-zero entry only for the site of the chain C y , and ∆ is a diagonal matrix.Furthermore we have f k = − [2 t cos[ k ] + µ ] I N x − t I + N x − t I − N x , (S2)and L k = − α sin[ k ] I N x − i α I + N x + i α I − N x , (S3)where I N x is an N x × N x identity matrix and [ I ± N x ] ij = δ i,j ± are also N x × N x matrices.We write ˜ H ( k ) = U † H ( k ) U where U = 1 − σ y τ y I N x + σ z τ z − σ x τ x I N x , (S4)we have introduced an N x × N x matrix whose entries are (cid:2) ¯ I N x (cid:3) nn (cid:48) = δ n,N y x +1 − n (cid:48) . Then˜ H (ˆΓ i ) = (cid:18) ¯ H (ˆΓ i ) 00 − ¯ H (ˆΓ i ) (cid:19) , (S5)from which we can calculate the invariant same as in the main text [55, 67]: ν Q1D = 1 − sgn (cid:2) det ¯ H (0) det ¯ H ( π ) (cid:3) . (S6)For example for N x = 3 we find, with t = 1,det ¯ H ,π = − (cid:0) ∆ + ( µ ± (cid:1) (cid:104) ∆ − V (cid:0) ∆ + ( µ ± (cid:1) + (cid:0) − λ + µ ( µ ± (cid:1) + 2∆ (cid:0) λ + µ ( µ ± (cid:1)(cid:105) . (S7)For larger N x the expressions become rapidly rather complicated, but can be calculated numerically relatively quickly.It is worth noting that the effective 1D chain calculated from the Green’s function and the quasi-1D systemconsidered here belong to different symmetry classes with different invariants, Z and Z respectively. In the case inwhich the chiral invariant | ν | ≥ Impurities as a suspended chain
If the magnetic impurity sites for the chain are above the superconductor’s surface, see Fig. S1, then the system isdescribed by the Hamiltonian H + H Ad where H Ad = t (cid:48) (cid:88) r ∈ C Ψ † r τ z Φ r + H.c. + U (cid:88) r ∈ C Φ † r σ z Φ r . (S8)2 Figure S1. A schematic of the magnetic impurity sites (green atoms) above the superconducting surface (dark blue atoms). Φ r = ( a r , ↑ , a r , ↓ , a † r , ↓ , − a † r , ↑ ) T , where a ( † ) r ,σ annihilates (creates) an electron with spin σ = ↑ , ↓ at site r ∈ C . H isgiven by (1), same as in the main text.The effect of the extra impurity sites is to generate an effective Zeeman term on the sites directly underneath givenby V = t (cid:48) /U . One can therefore approximately map this system to H + H imp with V = t (cid:48) /U . To demonstratethis we numerically diagonalized this tight-binding model and calculated the Majorana polarisation and the energycorresponding to the lowest-energy eigenvalue, see the main text for more details. The results can be seen in Fig. S2,where these quantities are plotted as a function of µ and V = t (cid:48) /U . Note that we obtain a perfect agreement toFig. S2. = t ' Ut μ / t ϵ / t t ' Ut μ / t t' = = t ' Ut μ / t ϵ / t t ' Ut μ / t U = Figure S2. A comparison of the lowest energy state energy (cid:15) and its Majorana polarization C for the impurities defined assites above the two-dimensional lattice. Here we take λ = 0 . t , and ∆ = 0 . t . Either t (cid:48) or U are kept constant and the otherterm is varied, as labelled in the figures. Analytical result for the chain Green’s function
We start with the momentum-space Hamiltonian for the two-dimensional superconducting substrate with Rashbaspin-orbit coupling: H k = − ξτ z − λ [sin k x σ y − sin k y σ x ] τ z − ∆ τ x , (S9)where ξ ≡ [ µ + 2 t (cos k x + cos k y )]. In what follows we first calculate the unperturbed Green’s function of the bulkusing G ( ω, k ) = (cid:0) i ω + H k (cid:1) − . For the sake of simplicity we perform this calculation only at ω = 0. Thus we get: G (0 , k ) = − { [ µ + 2 t (cos k x + cos k y )] τ z + 2 λ [sin k x σ y − sin k y σ x ] τ z + ∆ τ x } − (S10)This inversion can be performed using the following identity: (cid:18) M ∆ I ∆ I − M (cid:19) − ≡ (cid:18) M ( M + ∆ I ) − ∆( M + ∆ I ) − ∆( M + ∆ I ) − − M ( M + ∆ I ) − (cid:19) , (S11)3where we denoted M ≡ (cid:18) ξ − λ (sin k x − i sin k y )2i λ (sin k x + i sin k y ) ξ (cid:19) . (S12)Finally, we get: G (0 , k ) = − (cid:20) − i ξ ) + 4 λ (sin k x + sin k y ) (cid:18) M − ∆ − ∆ − − M − (cid:19) + 1(∆ + i ξ ) + 4 λ (sin k x + sin k y ) (cid:18) M + ∆ + ∆ + − M + (cid:19)(cid:21) , (S13)where we defined: M ± = (cid:18) ξ (∆ ± i ξ ) ∓ λ (sin k x + sin k y ) − λ (sin k x − i sin k y )2i∆ λ (sin k x + i sin k y ) ξ (∆ ± i ξ ) ∓ λ (sin k x + sin k y ) (cid:19) (S14)∆ ± = (cid:18) ∆(∆ ± i ξ ) ∓ λ (sin k x − i sin k y ) ± λ (sin k x + i sin k y ) ∆(∆ ± i ξ ) (cid:19) (S15)To calculate the line Green’s function G (0 , k y ) at zero energy we need to integrate the bulk Green’s function G (0 , k ) found above over k x from − π to π . All the odd integrands in k x yield zero, and hence we are left with thethree following integrals to compute: I ± = (cid:90) π − π dk x π ± i ξ ) + 4 λ (sin k x + sin k y ) (S16) I ± = (cid:90) π − π dk x π cos k x (∆ ± i ξ ) + 4 λ (sin k x + sin k y ) (S17) I ± = (cid:90) π − π dk x π cos 2 k x (∆ ± i ξ ) + 4 λ (sin k x + sin k y ) (S18)We note that for every integral above I + i = I − i (∆ → − ∆), therefore, we can compute only I − i , where i = 0 , , z = e i k x . We rewrite thedenominators of the integrands as follows:(∆ − i ξ ) + 4 λ (sin k x + sin k y ) = − ( t + λ ) (cid:18) z + 1 z − w + (cid:19) (cid:18) z + 1 z − w − (cid:19) = (cid:63) (S19)where w ± ≡ − A ± (cid:112) A + B, A ≡ t ( µ + 2 t cos k y + i∆) t + λ , B ≡ λ (1 + sin k y ) − ( µ + 2 t cos k y + i∆) t + λ . (S20)Solving the remaining two quadratic equations, we get: (cid:63) = − ( t + λ ) · z ( z − z )( z − z )( z − z )( z − z ) , z , = 12 (cid:18) w − ∓ (cid:113) w − − (cid:19) , z , = 12 (cid:18) w + ∓ (cid:113) w − (cid:19) . (S21)A simple analysis of the roots shows that:1) | z , | < | z , | > , if 2 λ (cid:113) k y + µ + 2 t cos k y (cid:62) , λ (cid:113) k y − ( µ + 2 t cos k y ) (cid:62) | z , | < | z , | > , if 2 λ (cid:113) k y + µ + 2 t cos k y > , λ (cid:113) k y − ( µ + 2 t cos k y ) < | z , | < | z , | > , if 2 λ (cid:113) k y + µ + 2 t cos k y < , λ (cid:113) k y − ( µ + 2 t cos k y ) > I − = − t + λ π i (cid:73) | z | =1 dz z Π i =1 ( z − z i ) = (S25)1) = − t + λ (cid:20) z ( z − z )( z − z )( z − z ) + z ( z − z )( z − z )( z − z ) (cid:21) (S26)2) = − t + λ (cid:20) z ( z − z )( z − z )( z − z ) + z ( z − z )( z − z )( z − z ) (cid:21) (S27)3) = − t + λ (cid:20) z ( z − z )( z − z )( z − z ) + z ( z − z )( z − z )( z − z ) (cid:21) (S28) I − = −
12 1 t + λ π i (cid:73) | z | =1 dz z + 1Π i =1 ( z − z i ) = (S29)1) = −
12 1 t + λ (cid:20) z + 1( z − z )( z − z )( z − z ) + z + 1( z − z )( z − z )( z − z ) (cid:21) (S30)2) = −
12 1 t + λ (cid:20) z + 1( z − z )( z − z )( z − z ) + z + 1( z − z )( z − z )( z − z ) (cid:21) (S31)3) = −
12 1 t + λ (cid:20) z + 1( z − z )( z − z )( z − z ) + z + 1( z − z )( z − z )( z − z ) (cid:21) (S32) I − = −
12 1 t + λ π i (cid:73) | z | =1 dz z + 1 z Π i =1 ( z − z i ) = (S33)1) = −
12 1 t + λ (cid:20) z + 1 z ( z − z )( z − z )( z − z ) + z + 1 z ( z − z )( z − z )( z − z ) (cid:21) (S34)2) = −
12 1 t + λ (cid:20) z + 1 z ( z − z )( z − z )( z − z ) + z + 1 z ( z − z )( z − z )( z − z ) (cid:21) (S35)3) = −
12 1 t + λ (cid:20) z + 1 z ( z − z )( z − z )( z − z ) + z + 1 z ( z − z )( z − z )( z − z ) (cid:21) (S36)The residue at z = 0 is calculated using the fact that Π i =1 z i = 1 (consequence of the Vieta’s formula). The remainingintegrals I +0 , I +1 , I +2 can be computed by substituting ∆ → − ∆ in the expressions above. Changing the sign of ∆ isequivalent to taking complex conjugates of the roots, i.e., z i → z ∗ i for i = 1 ,
4. Therefore, I + i = (cid:0) I − i (cid:1) ∗ for i = 0 , M ± matrices we get M ± = M ± = µ ∆ ± i( µ + 4 t − λ ) + 2 t (∆ ± µ ) cos k y ± t + λ ) cos 2 k y + [2 t ∆ ± t ( µ + 2 t cos k y )] cos k x ± (cid:0) t + λ (cid:1) cos 2 k x (S37) M ± = − λ (sin k x − i sin k y ) (S38) M ± = +2i∆ λ (sin k x + i sin k y ) (S39)Integrating the expressions above over k x , we get:˜ M ± = ˜ M ± = (cid:2) µ ∆ ± i( µ + 4 t − λ ) + 2 t (∆ ± µ ) cos k y ± t + λ ) cos 2 k y (cid:3) I ± + [2 t ∆ ± t ( µ + 2 t cos k y )] I ± ± (cid:0) t + λ (cid:1) I ± (S40)˜ M ± = − λ sin k y · I ± (S41)˜ M ± = − λ sin k y · I ± (S42)5For the ∆ ± matrices we get: ∆ ± = ∆ ± = ∆ ± i∆( µ + 2 t cos k y ) ± t ∆ cos k x (S43)∆ ± = ∓ λ (sin k x − i sin k y ) (S44)∆ ± = ± λ (sin k x + i sin k y ) (S45)Integrating the expressions above over k x , we get:˜∆ ± = ˜∆ ± = (cid:2) ∆ ± i∆( µ + 2 t cos k y ) (cid:3) I ± ± t ∆ · I ± (S46)˜∆ ± = ± λ sin k y · I ± (S47)˜∆ ± = ± λ sin k y · I ± (S48)Finally, we can write the integral of the Green’s function as follows: G (0 , k y ) = π (cid:90) − π dk x π G (0 , k ) = 12∆ (cid:88) σ = ± (cid:18) ˜ M σ ˜∆ σ ˜∆ σ − ˜ M σ (cid:19) ..