{'fraction_non_alphanumeric': 0.04594615931159706, 'fraction_numerical': 0.013447014552248625, 'mean_word_length': 4.275857281688186, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Topological superconductivity in a multichannel Yu-Shiba-Rusinov chain

Chains of magnetic atoms placed on the surface of an s-wave superconductor with large spin-orbit coupling provide a promising platform for the realization of topological superconducting states characterized by the presence of Majorana zero-energy modes. In this work we study the properties of the one-dimensional chain of Yu-Shiba-Rusinov states induced by magnetic impurities using a realistic model for the magnetic atoms that include the presence of multiple scattering channels. These channels are mixed by the spin-orbit coupling and, via the hybridization of the Yu-Shiba-Rusinov states at different sites of the chain, result in a multi-band structure for the chain. We obtain the topological phase diagram for such band structure. We identify the parameter regimes for which the different bands lead to a topological phase and show that the inclusion of higher bands can greatly enlarge the phase space for the realization of topological states.


I. INTRODUCTION
The search for Majorana zero modes (Majoranas) in condensed matter systems has been an active and exciting pursuit [1][2][3][4]. The reason for much excitement is due to the theoretical prediction that these modes manifest non-Abelian quantum statistics [5][6][7][8], and, as such, would open the possibility to realize topological quantum computing [9][10][11][12]. Currently, errors caused by the decoherence of the quantum states used to encode the data constitute the biggest fundamental obstacle for the realization of a scalable quantum computer. In topological quantum computing this obstacle is overcome by the topological protection of the quantum states used to encode the information. Most platforms for realizing topological phases of matter supporting Ising anyons (i.e. exotic defects binding Majorana zero modes) involve superconducting heterostructures . Recently, several works [34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51] have proposed that a chain of magnetic atoms placed on the surface of a superconductor can be in a robust topological phase characterized by the presence of Majorana modes located at its ends. In addition, recent experimental results [52] have shown that a chain of Fe atoms placed on the surface of superconducting Pb exhibits a zero-bias peak localized at its ends, consistent with the presence of Majorana modes. A very recent preprint also presents experimental results for a similar system [53]. A full understanding of the experiment presented in Refs. 52 and 53 is still being developed, for example, the height of the zero-bias peak is a small fraction of the predicted universal value of 2e 2 /h [54,55], a fact that could be attributed, for example, to finite temperature broadening or disorder.
The potential of a chain of magnetic impurities, placed on the surface of a superconductor, to be in a robust topological phase, calls for a thorough theoretical understanding of this system. So far most of the works have assumed one bound state per magnetic atom, corresponding to the zeroth angular momentum channel (l = 0) of Yu-Shiba-Rusinov(YSR) states [56][57][58]. However, realistic adatoms are expected to induce several bound states corresponding to different angular momentum scattering channels (i.e. l = 0, ±1, etc.). It has been shown experimentally that partial waves beyond s-wave are essential to explain the energy spectrum of magnetic atoms such as Mn, Cr and Fe [59][60][61]. Furthermore, the YSR states originating from the l = 0 channel (s-wave) are not always the lowest energy eigenstates. It is therefore necessary to understand the interplay of different angular momentum channels, in particular in the presence of significant Rashba spin-orbit coupling(SOC) [62]. The Rashba SOC term in the effective low-energy Hamiltonian describing the fermionic degrees of freedom is expected to be present at the surface due to broken inversion symmetry [63]. Moreover, the presence of SOC is required in order to have a stable topological phase for a ferromagnetically ordered chain [42,52] placed on the surface of an s-wave superconductor.
In this work we study a realistic model for a chain of magnetic impurities placed on the surface of an s-wave superconductor, see Fig. 1, that takes into account multiple scattering channels for the adatoms and the presence of SOC. Treating the magnetic impurities classically we study the properties of the bands formed by the hybridization of the YSR states bound to the different adatoms forming the chain. The multichannel treatment of the scattering potential of a single impurity implies that, for the chain, we obtain a multiband model. In the remainder, to simplify the presentation and to be able to clearly point out the main qualitative features of the multiband structure resulting from the hybridization of multichannel YSR states at different impurity sites, we assume that the YSR states corresponding to different values of |l| are well separated in energy. In this limit we immediately notice a fundamental difference between the band, s-band, formed from the |l| = 0 states and the bands formed from |l| > 0 states: the s-band can be assumed to be well separated from all the other bands and therefore can be treated effectively as a single isolated band; however, given the degeneracy, in the limit of no SOC, of the +|l|, −|l| states, we have that the two bands formed from +|l|, and −|l| states are always very close in energy and therefore that for bands formed from |l| > 0 states, a multiband treatment is necessary. Henceforth, we limit ourselves to the case in which max(|l|) = 1, and consider the minimal model that captures the aforementioned features. arXiv:1505.05862v2 [cond-mat.supr-con] 11 Nov 2015 The topological properties of the chain are determined by the band (bands) that is (are) closest to the midgap energy of the superconductor. To exemplify the main properties of the chain of multichannel magnetic impurities we consider two limits: (i) the limit in which the l = 0 YSR states are closest to the midgap energy of the superconductor so that for the chain the resulting s-band is also the closest to the midgap region, see Fig. 1 (b); we call this case "deep s-band" limit; (ii) the limit in which the |l| = 1 YSR states are closest to the midgap energy of the superconductor so that for the chain the resulting bands, "p-bands", are also the closest to the midgap region see Fig. 1 (c); we call this case "deep p-band" limit. We obtain the topological phase diagram for both the deep s-band and the deep p-band limit. The presence in the deep p-band limit of two bands close in energy would suggest that in this limit the topological phase could be strongly suppressed. Contrary to this naïve expectation we find that in the deep p-band limit the phase space in which the topological phase is present can be even larger than in the deep s-band limit. We also find that in the multichannel case the presence of SOC leads to the dependence of the chemical potential on the direction of the chain magnetization, which, in principle, allows one to tune between topological and non-topological phases. This is an important feature for braiding Majoranas [64].
The paper is organized as follows. In Sec. II we introduce our model and explain the general framework for the calculation. In Secs. III and IV, we derive an effective Hamiltonian and calculate the topological phase diagram as well as quasiparticle gap for the deep sand p-band limits, respectively. In Sec.V we discuss the qualitative difference between those two limits and compare their topological phase diagrams. The technical details are presented in the Appendices.

II. THEORETICAL MODEL FOR MULTICHANNEL YSR CHAIN
We consider a chain of magnetic impurities separated by a distance a and placed on top of an s-wave superconductor with Rashba spin-orbit coupling (SOC). The corresponding Hamiltonian describing an effectively two-dimensional superconducting film with Rashba SOC reads ( = 1): where k = (k x , k y ) is the electron momentum, ε k = k 2 /2m − µ with µ the chemical potential and m the effective electron mass,α is the strength of the Rashba SOC,ẑ is the unit vector normal to the plane, see Fig.1, and ∆ is the superconducting gap. The Pauli matrices σ i and τ i operate in spin and particle-hole space, respectively. The presence of the magnetic impurities is taken into account via the Hamiltonian where V j are the individual impurity magnetic potentials with R j , S j andJ being the impurity position, its classical spin (Color online) Schematic setup of the proposed structure supporting multichannel Yu-Shiba-Rusinov chain. Atoms with high magnetic moments form a ferromagnetic chain on the surface of thin film SC [42,52]. Each magnetic atom creates multiple Yu-Shiba-Rusinov states with distinct angular momentum quantum number l.
Here we consider l = −1, 0, and 1 states which forms three bands in a quasi one-dimensional system. In the topological phase, there is an odd number of Majorana zero-energy modes at the opposite ends of the chain. and its exchange coupling to the host quasiparticles, respectively. In order to find the band structure of the chain of YSR states induced by the magnetic impurities, we need to solve the Schrodinger equation: Here The spectrum of the subgap states is determined by the pole of the T -matrix [65]. In terms of the Green's function for the superconductor G = [E − H SC ] −1 , the eigenvalue problem for the subgap states is given by (1 − GH imp )ψ(r) = 0. We note that the Green's function of an s-wave superconductor in the presence of SOC has both even-and odd-parity components [66], and can be written as G(k; E) = l G l (k; E)e ilθ k where k = |k| and θ k = arctan k x /k y . In the presence of Rashba SOC, the G l are non-zero only for l = −1, 0, 1 and are given by: where σ ± = (σ x ± iσ y ) /2, and ε λ (k) = k 2 /2m − λαk − µ are the dispersions of the helical bands with λ = ±.
Transforming to momentum space, Eq. (3) can then be rewritten as: where ψ i (k) is the Fourier transform of the wave function centered at R i . Assuming that at the Fermi surface the scattering potential is weakly dependent on p ≡ |p| we have The magnetic potential can be decomposed into independent angular momentum channels: withJ l the angular momentum components ofJ(θ k − θ p ).
Since V j (θ k − θ p ) is Hermitian and an even function of θ k −θ p , we haveJ l =J −l and V j,l = −J l (S j · σ) τ 0 . We assume that l = 0, the s channel, and l = ±1, the p channels, are the dominant scattering channels: V j = l=0,±1 V j,l . For the impurity chain arranged along thex direction in a ferromagnetic pattern, i.e., S j = (cos(ϕ) sin θ, sin(ϕ) sin θ, cos(θ)), (illustrated in Fig. 1(a)), Eq. (7) in the angular mode representation takes the form where an overline indicates an integration over momentum amplitude: ψ i,l ≡ kdk 2π ψ i,l (k), with ψ i,l the angular momentum components of Using Eq. (4)-(6) one can rewrite G ij l−l (E) in terms of the dimensionless integrals: 1+α 2 is the density of states of the λ helical band at the Fermi level in the normal state, and N F = (N + + N − )/2. The analytic results for the above integrals in the limit D → ∞ and for G ij l−l (E) are presented in Appendix A and B. For convenience, we define the dimensionless exchange couplings J l ≡J l |S|πN F which will be used henceforth.
It is convenient to rewrite Eq. (10) in the following form: Here the local part of the matrix M ii l,l determines the YSR spectrum of a single magnetic atom [62] whereas the non-local part M ij l,l describes the hybridization between YSR states induced by the magnetic atoms at i and j sites. For an equally spaced magnetic atom chain with distance a between the two nearest atoms, this hybridization leads to the formation of the YSR bands. In the limit of k F a 1, which we consider henceforth, the hopping energy scale is proportional to 1/ √ k F a and, thus, the bandwidth W is small, i.e. W ∆. In this limit, the bands maintain the character of the single impurity YSR states and, thus, we refer to them as s or p-bands. Strictly speaking, SOC mixes different angular momentum states but, since we assume that α 1, this terminology is justified. When s and p bands are well-separated by a gap that is much larger than the temperature, see Fig. 1(b) and (c), the problem can be considerably simplified by integrating out the higher-energy bands (i.e. by taking into account virtual scattering processes to the bands higher in energy). In the following, we consider two limiting cases corresponding to the deep sand p-band limits and discuss the corresponding topological phase diagrams. We show that these two cases are qualitatively different since deep p band limit consists of two bands originating from the l = ±1 YSR states.

A. Derivation of the Effective Hamiltonian
We first consider the deep s-band limit such that the energy of the l = 0 state is close to the midgap, i.e., J 0 ∼ 1 with the on-site energy 0 ≈ ∆(1 − J 0 ) + O(α 2 ) → 0. Provided J 1 J 0 the l = 0 band is well separated from the p bands, i.e., the bandwidth W ∆ |J 0 − J 1 |. After integrating out the l = ±1 states, we obtain a tight-binding description for the single s-band with the p channels taken into account perturbatively, by allowing for virtual transitions through the p channels. This can be done by first solving for ψ i,±1 using Eq. (13) and substituting above expressions into the equation for l = 0 component and keeping terms up to the linear order in intersite coupling, we obtain where the matrix M ij s (E) is given in the Appendix C. In order to solve Eq. (15) analytically, we expand the local on-site matrix to the linear order in E around E = 0, assuming that and set E = 0 in the inter-site matrix: In doing so we ignore terms With these approximations, Eq.(15) can be written as where the local and non-local contributions are given by H ii The tight-binding Hamiltonian H s (i, j) is obtained by projecting Eq. (16) onto the local YSR states: the eigen spinors of the single-impurity bound states with energy ± 0 . The local basis can be found by solving the single-site equation M ii s (E)ϕ ± = 0 as a special case of Eq. (15), where the bound state energies are determined from Det M ii s (E) = 0, see Appendix D. We first assume that the magnetic atom moments are aligned ferromagnetically along a)ẑ-(out-of-plane), b)x-(along the chain direction), and c)ŷ-(in-plane but normal to the chain) axis, and present explicit expressions for the corresponding effective Hamiltonian. We note that due to the presence of SOC, the effective Hamiltonian is anisotropic which can be readily seen already at the single-impurity level [62].
In the a) and b) cases, by transforming the effective tightbinding Hamiltonian H s (i, j) to momentum space, we find that the corresponding Bogoliubov-de Gennes (BdG) Hamiltonian reads To order α 2 and α/ √ k F a, the effective hopping energy is given by where the functions I 0,± (k) and I 2,± (k) are obtained by taking the E → 0 limit in Eqs.(A17). The on-site energy is The effective p-wave pairing take the form The functions I n,± (k) ≡ I n,± (k, E = 0) and K n,± (k) ≡ K n,± (k, E = 0) are defined in the Appendix A, see Eqs.(A17-A20). The functions h(k) and∆(k) have the following properties h(k) = h(−k) and∆(−k) = −∆(k), and, thus, the gap is generically vanishing at k = 0, π/a. Moreover, the effective pairing∆(k) is vanishing for α → 0 analogous to the semiconductor nanowire proposal [17,18] in that the SOC controls the excitation gap.
To have a better understanding of the Hamiltonian structure in Eq. (17), it is instructive to perform a perturbative expansion of h(k) and∆(k), for example, around k = 0, see Eqs. (18) and (21) where the expressions for h (0,2) z,x and ∆ (1) are given in the Appendix E. In the limit a/ξ 0 → 0, where ξ 0 = v F /∆ is the superconducting coherence length, these functions have singular points for some values of k F a which is a consequence of the long-range nature of the hopping matrix element in the effective Hamiltonian. The presence of a finite coherence length ξ 0 , however, regularizes the singularities. Nevertheless, such a strong dependence on k F a leads to significant variations of the effective mass and Fermi velocity. Another important feature that we find is that the effective Hamiltonian is anisotropic due to the SOC (cf. Eq. (20)) which might be helpful to drive the topological transition by changing the direction of the magnetization of the impurities forming the chain. We note that this effect is absent for J 1 = 0, in which case we recover the results of Ref. [43]. Thus, the dependence of the effective chemical potential on the angle θ, which is the only tuning parameter in the Hamiltonian (17), is a feature of the multichannel magnetic impurity model. Finally, we discuss the case in which S ŷ. In this case the magnetic exchange energy term and the SOC term commute, and, as a result, the system is gapless. Thus, this case is not interesting from the point of view of Majorana physics.  Having derived the effective Hamiltonian (17), we can now calculate the topological phase diagram. The Hamiltonian (17) for a generic direction of magnetization is in the symmetry class D [67][68][69], and, thus, is characterized by the Z 2 topological invariant, the so-called Majorana number M [70]: The system is in the topological superconducting phase when M = −1, whereas M = +1 indicates a non-topological phase. We obtain the topological phase diagram by calculating M. Figure 2 (a) shows the topological phase diagram in the (k F a, J 0 ) plane for the deep s-band limit for the case in which the magnetic moments of the impurities forming the chain are aligned along the z direction and α = 0.3. The range of values of k F a has been chosen so that the inequality 1/ √ k F a 1, on which the treatment of the previous section relies, is well satisfied. From Fig. 2 (a) we see that, for α = 0.3 there is a large fraction of the (k F a, J 0 ) in which the chain is expected to be in a topological phase characterized by odd number of Majoranas at its ends. It is interesting to ask how the topological phase diagram is affected by the strength of the SOC. This question is addressed by the results presented in Fig. 2 (b) that shows the dependence of M on α and k F a for a fixed value of J 0 . Again we notice that there is a large fraction of the (α, k F a) space in which the chain is expected to be in a topological phase.
Experimentally it can be challenging to vary in a controlled way parameters such as α, J 0 , and k F a and therefore to verify the theoretical predictions shown in Fig. 2 (a), (b). However, our multichannel treatment, contrary to the single channel treatment [39,43,48], reveals that, the topological character of the chain state also depends on the direction of the magnetization. This is shown in Fig. 2 (c) in which we can observe that the boundaries of the topological phase in the (k F a, J 0 ) plane are different depending on the direction, z or x, of the magnetic moment of the impurities forming the chain. This result can be easily understood by considering that the onsite energy, Eqs. (20), depends on the direction of the impurity magnetization. The dependence of the topological index on the direction of the chain magnetization is very important because it allows, in principle, to tune the chain in and out of a topologically nontrivial phase by varying a quantity that can be tuned and controlled experimentally by using an external magnetic field. The appearance, disappearance, of a zero bias peak as a function of the direction of the magnetization of the chain according to theoretical predictions like the ones presented in Fig. 2 (c) would provide compelling evidence of the Majorana character of the observed zero energy states.
In addition to the topological index (Majorana number), we have also calculated quasiparticle excitation gap as a function of J 0 , α, and k F a, see Fig. 2(d) and Fig. 3(a)-(c). One can notice that the closing of the gap is consistent with the phase diagram shown in Fig. 2(a)-(b). Additionally, one has information regarding the magnitude of the gap deep in the topological phase which is crucial for understanding the stability of the topological phase. Fig. 3(c) shows there also exist gap closing points inside the topological phase, not determined by calculating the topological index. To identify the location of these gap closing points, we plot the hopping energy spectrum h(k), the effective pairing energy spectrum Im∆(k), and the energy spectrum of the particle band as a function of k and α in Fig. 3(d)-(f), respectively. Clearly, for certain values of α both hopping and pairing energies vanish at certain momentum points between 0 and π resulting in the gap closing. In the limit J 1 → 0 it can be shown analytically that these accidental zeros occur when the k F,± are commensurate with each other, i.e. when α = nπ/(k F a), with n ∈ N. We will discuss the origin of these accidental gap closing points in more details in the next section.  .

A. Derivation of the Effective Hamiltonian
We now discuss the deep p-band limit assuming that the energy of the l = ±1 states is lower than that of l = 0, i.e., J 0 J 1 . In the limit J 1 ∼ 1 and α 2 J 0 1, the on-site energy ≈ ∆(1−J 1 )+O(α 2 ) is close to the midgap [62]. Once again, we assume that the l = ±1 states are well separated from the l = 0 state, and, therefore, l = 0 states can be integrated out. As a result, we obtain a tight-binding description for the p-bands with the s channel taken into account perturbatively by allowing for the transitions through intermediate virtual l = 0 states. Note that there is a significant difference with respect to the calculation in Sec. III since there are now two low-energy p-bands. Following the same procedure as in the previous section, we obtain the 8 dimensional matrix equation for the deep p-band limit: where . We then project H ij p onto the local basis of YSR states: φ 1,+ , φ 2,+ , φ 1,− , φ 2,− T where φ 1(2),± are the eigen spinors of the single-impurity bound states with energy ± 1(2) in channel 1 (l = −1) and channel 2 (l = 1). The local basis can be found by solving the single-site equation M ii p (E)φ = 0 as a special case of Eq. (24), where the bound state energies are determined from Det M ii p (E) = 0, see Appendix G. After the projection onto the local basis we obtain the effective Hamiltonian H p (k) describing the two coupled bands of the YSR chain in the deep-p band limit.
Simple analytical expressions for H p (k) can be obtained when the impurity spins are polarized normal to the plane, i.e. alongẑ direction. In the case that the impurities are polarized inẑ direction, transforming the effective tight-binding Hamiltonian H p (i, j) to momentum space, we obtain a two-band BdG Hamiltonian, The coefficients here satisfy the following properties: h ij (k) = h ij (−k),∆ ij (−k) = −∆ ij (k) and, therefore, ∆ ij (p) = 0 at p = 0, π/a. The effective hopping energies include both intra-channel hopping and inter-channel hopping with the on-site energies The effective p-wave pairing also contains both intra-channel pairing∆ 11 k.
22 k. and inter-channel pairing The coefficients of the small k expansions are explained in the Appendix E.
In order to understand the physics described by the Hamiltonian (26), we first discuss the effect of SOC on the normalstate band structure (i.e. ∆ ij = 0). The spectrum for the two bands reads As shown in the Appendix A, to leading order in 1/ √ k F a, I 0,λ (k) ≈ I 2,λ (k). Hence h 12 is approximately the same as h 11 and h 22 . In the absence of SOC, δ 12 vanishes, and the band structure is characterized by a heavy- with the bandwidth doubled compared to s-band, as shown in Fig. 4(a). The on-site orbital structure of these two bands are symmetric(light) and anti-symmetric(heavy) combinations of l = ±1 states. The physical origin of these orbital structures reflects the degeneracy due to the isotropic magnetic potential and the asymptotically equal hopping amplitudes.
In the presence of SOC and a finite s channel coupling (i.e. J 0 = 0 and α = 0) δ 12 becomes finite, which induces a hybridization gap between the two bands leading to an avoided level crossing, as shown in Fig. 4(b). The induced hybridization gap gives rise to an interesting feature in the topological phase diagram as discussed below.

B. Topological Properties
The topological phase diagram in the deep p-band limit (26) involves two bands which are hybridized by the SOC. As a result, it exhibits a more intricate dependence on the parameters compared with the deep s-band limit. In order to calculate the Z 2 topological invariant M, we use the method developed for the multiband semiconductor nanowire system [71]: M = sgn [PfB(0)PfB(π/a)] = ±1, where the antisymmetric matrix B(p) = Hẑ p (p)τ x is calculated at the particle-hole invariant points: p = 0, π/a. For the twoband system, the corresponding expression for the Pfaffian is The topological phase diagram for the effective Hamiltonian (26) as a function of J 1 and k F a is shown in Fig. 5 (a). We immediately notice that in the deep p-band limit there is a large fraction of the (J 1 , k F a) plane in which the chain is in a topological phase. This is somewhat surprising given that in the deep p-band we have two bands close in energy and one could expect that the hybridization of the two bands could (a)  lead for most values of J 1 and k F a to a situation in which the Fermi energy intersects the bands an even number of times. The reason why this is not the case is evident from the plots of the two bands for a typical situation shown in Fig. 4: one can see that the result of the hybridization of the bands leads to the formation of the light-and heavy-fermion bands. As a consequence, for most values of J 1 and k F a the Fermi energy crosses only once the light-fermion band. On the other hand, we expect that when J 1 and k F a are such that the Fermi energy is very close to the heavy band the chain should be in a topologically trivial phase. Given the flatness of the heavy band we expect that this will happen only for a very small range of values of J 1 . Indeed, a close inspection of the topological phase diagram of Fig. 5 (a) shows that for J 1 very close to 1 there is a narrow region, highlighted in the figure by a rectangular box, in which the chain is in a topologically trivial phase.
To understand the re-entrance region from the topological phase to the trivial phase when J 1 ∼ 1, in Fig. 5 (b) we show a zoom-in of the topological phase diagram in the region surrounded by the rectangular box in Fig. 5 (a). In the zoom-in plot the presence of a topologically trivial phase is clearly visible. In this region, due the fact that J 1 → 1 the effective chemical potential is close to zero so that the heavy-fermion band becomes important: in this situation the chemical potential lies inside the hybridization gap, as shown in Fig. 4 (b), which results in the re-entrance to the trivial phase. Outside this narrow range of values of J 1 the heavy-fermion band is either completely empty or filled and the chemical potential crosses the light-fermion band an odd number of times, thus, there is an odd number of Majorana zero modes per end originating from the light-fermion band. Due to this interplay, the two YSR bands do not annihilate each other (except for the small region close to J 1 → 1) which is a peculiar feature of this model.
In order to understand the stability of the topological phase and to corroborate the results discussed above, we also compute the quasiparticle excitation gap E g , see Fig. 5(c). Figure 5(d)-(e) plot the value of E g along the two line-cuts on the phase diagram near the re-entrance region. One can see that the quasiparticle gap closing is consistent with the topological phase diagram in Fig. 5(b). We once again find that SOC controls the magnitude of the quasiparticle gap and, as such, is crucial for the stability of the topological phase.
In addition to the quasiparticle gap closing at the topological phase boundary, one can notice that there are also points where the gap vanishes in the topological phase, see Fig.5(d) at J 1 ∼ 1.005 and Fig.5(e) at J 1 ∼ 1.017. We now investigate in detail the two regions across these points using an additional symmetry of our effective Hamiltonian (26). In addition to the particle-hole symmetry P = τ x K where K refers to complex conjugation, our effective spinless Hamiltonian also has a pseudo-time reversal symmetry T = K. Using these two symmetries, one can construct another symmetry -chiral symmetry S = T P = τ x which anticommutes with the Hamiltonian (26). Thus, the Hamiltonian (26) belongs to the BDI symmetry class [67][68][69] which is characterized by the integer invariant W and supports multiple spatiallyoverlapping Majorana zero modes [72]. In order to calculate the topological index W, it is convenient to transform the Hamiltonian (26) into a chirality basis using a unitary transformation U = e −i π 4 τy which converts the Hamiltonian to the off-diagonal form: Then, the winding number (i.e. the number of Majorana modes per each end) can be calculated by introducing a complex variable z(k) = det[A(k)]/| det[A(k)]|, and calculating the integral Using this analysis we find that, for example, the phases at J 1 = 1.001 and J 1 = 1.007 in Fig. 5(d) have different winding numbers W(J 1 = 1.001) = 1 and W(J 1 = 1.007) = −1.
Thus, gap closing between these two regions corresponds to the transition between W = ±1. The same argument holds for J 1 = 1.015. Thus, accidental gap closing points inside of the topological or non-topological phases are not really accidental but represent the change of the winding number by an even integer. The analysis above relies on the chiral symmetry. However, in realistic systems the chiral symmetry can be easily broken by allowing, for example, for a generic direction of magnetic chain polarization (i.e. along y-axis). The precise magnitude for the Majorana splitting energy, which is important for tunneling transport measurements, depends on the details of the chiral-symmetry-breaking perturbations, and we refer a reader to Refs. [45,48,[73][74][75] for more details. As a consequence, the topological phases identified by the parity of the topological index are expected to be much more robust and this is the reason that our analysis has been focused mostly on characterizing the dependence of such index on the parameters of the system.

V. CONCLUSIONS
We have studied the topological properties of a chain of magnetic impurities placed on the surface of an s-wave superconductor with Rashba spin-orbit coupling taking into account the presence of multiple scattering channels, in the limit in which the states induced by isolated impurities are well described as Yu-Shiba-Rusinov states and the distance a between the impurities forming the chain is such that k F a 1. The inclusion of multiple angular momentum scattering channels, l, implies that for the chain of YSR states we have multiple bands. We have shown that the multiband character of the bands strongly affects the topological properties of the chain for the case when the lowest energy bands are the ones arising from the hybridization of YSR states with |l| > 0. Considering the lowest l = 0 and |l| = 1 channels we have obtained the topological phase diagram in the deep s-band and deep p-band limits. Our results show, somehow unexpectedly, that even in the deep p-band limit there is a large region of parameters space in which the chain is in a topological phase. Moreover, we find that even though the deep p-band case involves two bands, and it is not obvious a priori that these two bands do not "annihilate" each other, this limit seems to be more favorable for the observation of Majorana zero modes. This can be seen from the comparison, Fig. 6, of the topological phase diagram obtained in the deep s-band and deep p-band limit: we see that, for the same range of values of k F a and relevant coupling constants (J 0 in the deep s-band, J 1 in the deep p-band limit) the phase space where the chain is in a topologically non-trivial state is larger in the deep p-band limit than in the deep s-band limit. We have also characterized the stability of the topological states by computing the quasiparticle excitation gap. Our results have important implications for the ongoing experimental search for Majorana zero modes in this system. The interatomic spacing in Ref. [52] is of the order of the Fermi wave length (k F a ∼ 1), in which case direct tunneling between iron atoms needs to be included, whereas our calculation assumes k F a 1. However, at the qualitative level one can already draw a number of conclusions. As k F a decreases the bandwidth of subgap YSR bands becomes larger and, thus, the mixing between them becomes even more important. Therefore, it is conceivable that the subgap features reported recently by Yazdani et al. [76] might correspond to l > 0 YSR bands. In order to understand these additional subgap features within the framework of our model, one would need to systematically measure the spectrum of the YSR states of single-and two-atom structures [62] to identify the angular momentum channels that are the most relevant for the onedimensional chains. Another important feature that we predict is that SOC allows one to tune the effective chemical potential for the YSR bands and drive the topological phase transition using an external magnetic field. As such, large SOC might be helpful for the manipulation of the Majorana zero modes in this system. At the theoretical level, it would be interesting to establish the correspondence between our results involving classical magnetic impurities with the more microscopic multi-orbital Anderson model for magnetic impurities.
Appendix A: Analytic Expressions of Integrals I l,λ and K l,λ In this Appendix, we provide analytic expressions of the integral functions defined in Eq. (11) and Eq. (12): Here J n (z) and H n (z) are Bessel and Struve functions of order n, respectively; Note that the expressions for K l,λ (x; E) given above are valid for x = 0, and the integral K l,λ (0; E) = 0 for x = 0. Assuming k F |x| 1 and ζ −1 k F,λ , we can use the asymptotic forms of the Bessel and Struve functions [77]. In the limit k F x 1, one can find approximate expressions up to the order 1/(k F x) 2 : The corresponding Fourier transforms of the above asymptotic forms to the leading order of 1 √ k F a are given by where Li s (z) is the polylogarithm function The non-local Green's function (i = j) is defined as where m = l − l = 0, ±1, ±2, and the functions I n,λ (x ij ; E) and K n,λ (x ij ; E) are given in Appendix A.
Appendix C: Derivation of effective Hamiltonian in deep s-band limit In this Appendix we provide the details of the derivation of effective Hamiltonian in deep s-band limit. Using the method outlined in Sec. II, we find that the effective eigenvalue equation for the deep s-band limit reads After substituting Eq. (14) back into the above equation, one finds The assumption k F a 1 allows one to neglect the terms O((M i =j ) 2 ). Using the tight-binding approximation, one finally arrives at Appendix D: Local basis in the deep s-band limit In this Appendix, we provide details for the projection procedure used to the derive an effective Hamiltonian for the deep s-band limit. Assuming that the magnetic-atom polarization is alongẑ axis, the unnormalized local basis for magnetic impurity reads: When the polarization is alongx-axis, the unnormalized local basis for magnetic impurity is given by For the sake of completeness, we also derive effective Hamiltonian when the magnetization is along y-axis. In this case, the unnormalized local basis for magnetic impurity becomes After the projection onto the local basis, the effective Hamiltonian is given by with the functions h y (k) and d y (k) being Here the on-site energy y = x . One can see that quasiparticle spectrum in this case is indeed gapless.
Appendix E: Effective Hamiltonian in the long wavelength limit k → 0 It is instructive to expand the functions I(n, k) and K(n, k) appearing in our effective Hamiltonian close to k = 0 in order to understand the spectrum qualitatively. After some algebra, one finds The dependence of the functions A 0 (z), A 2 (z), B 1 (z) and C 1 (z) on the external parameters is shown in Fig. 7. One can notice that when k F,λ a = 2πn with n being an integer, these functions have singularities which follows from the definition of polylogarithm function. These singularities are cutoff by the finite coherence length. In realistic systems, however, the superconducting coherence length is much larger than the interatomic spacing, and, thus, the parameters such as effective mass and Fermi velocity are strongly dependent on k F a, see Fig. 7. Finally, the expansion of the coefficients in the deep s-band Hamiltonian at k → 0 becomes The expansion coefficients in the deep p-band Hamiltonian are Appendix F: Derivation of effective Hamiltonian in deep p-band limit In this Appendix, we provide the details of the derivation of effective Hamiltonian in deep p-band limit. Using the method outlined in Sec. II, we find the effective eigenvalue equation for the deep p-band limit. The corresponding equations for p-wave bands are given by j,l M ij −1,l (E)ψ j,l = 0 (F1) j,l M ij 1,l (E)ψ j,l = 0 (F2) In order to integrate out s-channel, we have to solve for ψ i,0 finding that Substituting Eq. (F3) into Eq. (F1, F2) and following the same procedure as in Appendix C, we eventually obtain two coupled equations for the p-wave bands After some manipulations, one can write eigenvalue equations in the compact form, see Eq. (24).