Two-dimensional chiral topological superconductivity in Shiba lattices

The chiral p-wave superconductor is the archetypal example of a state of matter that supports non-Abelian anyons, a highly desired type of exotic quasiparticle. With this, it is foundational for the distant goal of building a topological quantum computer. While some candidate materials for bulk chiral superconductors exist, they are subject of an ongoing debate about their actual paring state. Here we propose an alternative route to chiral superconductivity, consisting of the surface of an ordinary superconductor decorated with a two-dimensional lattice of magnetic impurities. We furthermore identify a promising experimental platform to realize this proposal.

T he chiral p-wave superconductor in two dimensions (2D) and the closely related fractional quantum Hall Pfaffian state at filling fraction n ¼ 5/2 are the archetypal examples of topologically ordered states of matter that support non-Abelian anyonic excitations [1][2][3] . The theoretical exploration of these states has shaped our understanding of topological order and is foundational for the distant goal of building a topological quantum computer [4][5][6] . In contrast to the n ¼ 5/2 fractional quantum Hall state, however, to date chiral p-wave superconductivity has not been confirmed in any experimental system. The most prominent candidate system, superconducting Sr 2 RuO 4 (ref. 7), has been the subject of an ongoing debate about its actual paring state [8][9][10][11] . Chiral superconductors break time-reversal symmetry. This hinders the formation of Cooper pairs, since orbital (and possibly paramagnetic) pair-breaking effects can come into play. Depairing is also the main hurdle for realizing a line of proposals in which layered heterostructures involving ferromagnets and s-wave superconductors are used to build an artificial 2D p-wave superconductor [12][13][14][15] . The guiding principle for these proposals is to design a band structure with a single normal-state Fermi surface with no spin-degeneracy. If Rashba spin-orbit coupling is present so that the states on this Fermi surface are not fully spin-polarized, and if s-wave superconductivity is proximityinduced in such a system, the effective pairing near the single Fermi surface is equivalent to that of a chiral p-wave superconductor.
In one dimension (1D), based on the principle of combining spin-orbit coupling and externally applied magnetic fields, various groups have proposed engineering artificial realizations of p-wave superconductors [16][17][18][19][20][21][22] . An experimental realization of these proposals employing semiconductor nanowires with strong spin-orbit coupling 23 has reported Majorana fermion signatures. In this setup, the externally applied magnetic field has to be rather small (to avoid suppressing superconductivity). As the phasespace for the existence of a topological superconductor is controlled by the Zeeman gap 24 , these systems require a delicate balance of the parameters involved (spin-orbit coupling, magnetic field and chemical potential) in order to create the topological superconductor.
Recently, a 1D topological superconductor was realized in a system that is quite distinct but employs similar microscopic ingredients-spin-orbit coupling, ferromagnetism and s-wave superconductivity 25 . A chain of magnetic Fe atoms is deposited on the surface of an s-wave superconductor with strong spin-orbit interactions. The Fe chain is ferromagnetically ordered 25 with a large magnetic moment, which takes the role of the magnetic field in the nanowire experiments. Unlike previous proposals, this magnetic field is mostly localized on the Fe chain, with small leakage outside. Superconductivity is not destroyed along the chain. In this setup, the energy scale of the exchange coupling of the Fe atoms is typically much larger than that of the Rashba spin-orbit coupling and the superconducting pairing. The ferromagnetically ordered Fe atoms induce localized Shiba states within the gap of the superconductor [26][27][28] . The hybridization of these states forms the band structure of a 1D p-wave superconductor that supports Majorana end states 29,30 . Because the Fe bands are fully spin split, no additional control over the chemical potential is necessary. A similar scenario applies when the Fe orbitals are magnetic but itinerant 24 .
In this article, we point out that this strategy can also be successful in 2D. Magnetic adatoms on the surface of a superconductor with strong spin-orbit coupling, when arranged in a 2D lattice, can yield a 2D topological chiral p-wave superconductor whose chiral Majorana edge modes can be observed in scanning tunneling microscope measurements.
To shed light on the rich range of possibilities, we analyse the topological properties of a system with dense local moments that are exchange coupled to a model 2D superconductor, demonstrating that topological superconductors with higher Chern numbers, and consequently multiple chiral Majorana edge channels, can easily occur. We are also able to analyse the model's dilute magnetic impurity limit analytically and obtain numerical topological phase diagrams for intermediate impurity concentrations. Based on density-functional-theory (DFT) calculations, we further propose realizing 2D topological chiral p-wave superconductors experimentally by depositing transition metal adatoms on superconducting Pb. The type of magnetic ion can be varied to access different strengths of the magnetic moment. In the case of Fe adatoms on a Pb (111) surface, we show that strong magnetic order in general leads to an odd number of 2D Fermi surface segments. As a consequence the proximity-induced superconducting phases can have nonzero Chern numbers and chiral Majorana edge modes.

Results
Model Hamiltonian. We first present a model system that bears a number of generic features of superconductor surfaces with ferromagnetically ordered magnetic adatoms (see Fig. 1). To render the analytical calculations tractable, we consider a Hamiltonian that models only the surface layer of a bulk 3D s-wave superconductor on which the s-wave superconducting order parameter D is induced from the bulk. On this superconducting layer, we model the magnetic impurities as classical spins whose only interaction with the electrons in the superconductor is through Zeeman-like couplings [26][27][28] . Employing a tight-binding description on a 2D square lattice L that is spanned by the primitive lattice vectorsê 1 andê 2 , we consider the mean-field Hamiltonian Here, c w r ¼ c w r;" ; c w r;# is a spinor of the creation operators for electrons at site r with spin m, k, and s i ,i ¼ 1,2,3, are the three  Pauli matrices. We denote by t the nearest neighbor hopping integral in the superconductor, m the chemical potential, and a the strength of the Rashba spin-orbit coupling. Classical magnetic moments ferromagnetically aligned normal to the plane in the s 3 direction are positioned on a sublattice L* of L and are exchangecoupled to electrons via the term proportional to J. Drawing experience from the 1D situation and the ab-initio calculations presented below, the physically relevant hierarchy of energy scales that we consider here is given by Dense impurity limit. As a warmup, it is instructive to consider the simplest situation in which each lattice site is coupled to a magnetic moment, i.e., L ¼ L*. This limit is representative of system with self-assembled islands of magnetic adatoms. Its consideration allows us to highlight the difference between the regime in equation (2), and that of small J that has been studied previously 12,14,15 . In particular, it is possible to access a phase with Chern number 2 in the large J limit. For the case L ¼ L* and a ¼ 0, the Hamiltonian in equation (1) is gapless at zero energy with nodal lines in momentum space defined by (see Supplementary Note 1) where e k ¼ 2t(cos k 1 þ cos k 2 ) À m. Adding spin-orbit coupling aa0 will generically lift these degeneracy lines and gap the spectrum, except if the degeneracy occurs at the four inversion- where the spin-orbit coupling vanishes. The occurrence of these nodal points is a signature of transitions between phases characterized by different Chern numbers. The critical chemical potentials are therefore determined by equation (3) with k ¼ k ij , which yields Around each of the nodal points in the three-dimensional k-m parameter space, we can reduce the Hamiltonian to an effective two-band model and expand it to linear order in the deviations dk from k ij and the deviations dm from m ijl yielding Here the Pauli matrices act on a subspace defined by the two bands that satisfy condition in equation (3). In the k-m space, the Hamiltonian in equation (5) is a Weyl Hamiltonian that is characterized by a unit topological charge i Â j Â l ¼ ± 1. As the chemical potential ramps through a critical point at m ijl , the Chern number associated with the original Hamiltonian in equation (1) changes precisely by the value of the topological charge, or the total topological charges when multiple nodal points occur at the same m. Therefore, by increasing m and assuming 2t4 ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi J 2 À D 2 p , the system exhibits phases with the Chern numbers (see Supplementary Note 1) where we have neglected the trivial phases with m falling outside the band width. Hence, with homogeneous magnetization, the superconductor may already exhibit Chern number equal to 2. In cases where magnetic impurities are spaced more sparsely, that is, if L* is a sublattice of L, even higher Chern numbers can be obtained. We solved the Hamiltonian in equation (1) numerically for the case of one magnetic impurity every 2 Â 2 and 3 Â 3 plaquettes of the square lattice, and show the phase diagrams in Fig. 2. From the phase diagrams, we can read off three general features. First, at small chemical potentials around the band bottom, where the Fermi wavelengths are larger than or comparable to the lattice spacing of L*, the sequence of Chern numbers ( À 1, 0, þ 2) always occurs when J/t is not too large-this universal feature corresponds to the dense limit that we have discussed above; Second, at larger chemical potentials, more than two Fermi surfaces can exist in the reduced Brillouin zone defined by L* as the Fermi wavelengths are significantly smaller than the lattice spacing of L*, as a consequence higher Chern numbers can occur (for example, 8 in Fig. 2a and 15 in Fig. 2b) but the trade-off is an overall smaller induced gap; Third, the phases with different Chern numbers are generally separated by lines (in the 2D parameter space) defined by conditions similar to equation (4), and across each specific line the change of Chern numbers is a constant determined essentially the same way as in our preceding analysis.
Dilute impurity limit. To better understand the dilute impurity limit, we complement our results on the lattice by a calculation in which we treat the underlying 2D superconductor in the continuum limit and consider sparsely distributed Shiba impurities that are arranged in a square lattice on top of it. This allows us to derive an effective two-band model for the hybridizing Shiba states. This effective Hamiltonian represents a chiral p-wave superconductor in the appropriate parameter regime.
The strategy of our derivation is inspired by the calculation for 1D Shiba chains of Pientka et al. 29 (for details, see Methods and Supplementary Note 2). We start from a 4 Â 4 Bogoliubov-de Gennes (BdG) Hamiltonian that acts on 4-spinor valued wave functions cðrÞ ¼ ðc " ; c # ; c w # ; À c w " Þ T ðrÞ, r 2 R 2 , with x k ¼ k 2 /(2m) À m þ a(k 1 s 2 À k 2 s 1 ). In addition, the magnetic impurities are represented by the Hamiltonian where S 3 ¼ s 3 #t 0 , with t 0 the identity matrix acting on particlehole space. If we restrict the wave-function C(r) to the locations r*AL* of the impurities, we obtain the self-consistency equation for the Fourier transformsCðqÞ ¼ P r Ã 2L Ã e À iqÁr Ã Cðr Ã Þ, where the momentum qA[0,2p) 2 now belongs to the L* Brillouin zone and we have set the impurity spacing to unity. We need two more steps to reduce equation (9) to an effective two-band model for the Shiba states, assuming they are deep in the superconducting gap and dilute compared with the Fermi wavelength. First, the left-hand side is expanded to linear order in the energy E to cast the equation in the form of the timeindependent Schroedinger equation. Second, we project the effective Hamiltonian into the eigenstates of an isolated Shiba impurity on every site, given by C þ ¼ 1; 0; 1; 0 ð Þ T = ffiffi ffi 2 p and C À ¼ 0; 1; 0; À 1 ð Þ T = ffiffi ffi 2 p (ref. 29). We obtain the effective two-band Hamiltonian are defined in terms of the functions In the equation above, x is the coherence length of the 2D superconductor (without the magnetic impurities), Ç ma are the Fermi wave-vectors for its two spin-split bands, and Z ¼ J½k is a dimensionless parameter. In the deep Shiba limit in which the projection in the states C þ and C À is justified, we have ZB1.
The Hamiltonian in equation (10) represents the effective superconductor formed by the Shiba bound states within the gap of the underlying s-wave superconductor. Similar to the case of the effective two-band model in equation (5), this Hamiltonian can have nodal points in k-m-space at which the Chern number changes. However, unlike in equation (5), the nodes can occur at any point in the Brillouin zone, making an analytic treatment intractable. In addition, the validity of equation (10) requires a self-consistency that permits the low-energy expansion of equation (9). Therefore, instead of computing the Chern numbers in an extended parameter space (for examples of Chern numbers along several linecuts of the parameter space, see Supplementary  Fig. 1), we focus on information that can be obtained at special points of the Brillouin zone at infinitesimal energy. To that end, observe that the Hamiltonian in equation (10) has C 4 rotational symmetry. Thus, any gap closing at points other than the C 4symmetric momenta k ¼ (0, 0) and k ¼ (p, p) changes the Chern number by an even integer due to the symmetry-imposed multiplicity of the nodal points. By expanding the Hamiltonian into Dirac form around the C 4 -symmetric momenta, we obtain the expression (see Supplementary Note 2) Figure 3 | Phase diagram in the dilute impurity limit. Parity of the Chern number (white for even and black for odd) for the Bogoliubov-de Gennes band structure of an s-wave superconductor decorated with dilute Shiba impurities, following equation (14). Here the spin-orbit coupling strength a is scaled by the Fermi velocity v F , and the lattice spacing of the Shiba impurities a L* is scaled by the Fermi wavelength l F . Both v F and l F are taken in the limit a-0 and we have assumed the superconducting coherence for the parity of the Chern number. The numerical evaluation of this equation is shown in Fig. 3 in the form of a phase diagram.
Helical magnetic order. We have also performed calculations for magnetic orders other than simple ferromagnetism. In particular, the case where the magnetic configurations corresponds to 2D helices is related to previous studies on 1D helices [31][32][33][34][35][36][37][38][39][40] . We obtained criteria for such a system to be fully gapped by proximity effect, and found that the fully-gapped superconducting phases can be generically topologically nontrivial. The results and phase diagrams are presented in Supplementary Note 3 and Supplementary Figs 2-4.
Material proposal. We complement our simple model considerations with a specific material proposal to realize a superconductor with nonzero Chern numbers. For that, we consider transition metal atoms, in particular Fe, deposited on the (111) surface of Pb, a strong type-II superconductor. The same combination of materials, but a different surface of Pb, was used in the experimental realization for the 1D p-wave superconductor 25 . The Pb atoms on the (111) surface form a triangular lattice.
Through ab-initio calculations that assume a ferromagnetic alignment of the Fe magnetic moments and include spin-orbit coupling, we compared the relaxation energies for various densities and arrangements of Fe adatoms on the Pb (111) surface, and found that a deposition with one Fe atom in each triangular plaquette is particularly favorable (see Supplementary  Figs 5 and 6). In this case, the Fe atoms form a honeycomb lattice in which the atoms sit at different heights in each sublattice (see Fig. 4a). We further performed DFT calculations of the electronic structure. The resulting (not spin-degenerate) Fermi surface and band structure restricted to the Fe d-orbitals is shown in Fig. 4b,c. Critically, we find an odd number of Fermi surfaces (for examples of the Fermi surfaces at several Fermi energies, see Supplementary  Fig. 7). The Chern number of the corresponding BdG Hamiltonian is indeed nonzero over a large range of the chemical potential (see Fig. 4d).

Discussion
In conclusion, we have proposed a versatile platform for realizing chiral superconductors in 2D. We have obtained analytically the topological phase diagram (Chern number and gaps of the superconductor) of the dilute and dense limit, and numerically evaluated the phase diagram in the intermediate regime. We then showed through a more realistic ab-initio calculation that ferromagnetically ordered Fe atoms on the (111) surface of Pb in the dense limit could give rise to a chiral superconductor.
To find flat islands of magnetic adatoms on the Pb substrate, however, is currently an experimental challenge because under standard growth conditions the magnetic adatoms tend to ball up instead of making flat islands on the (111) surface of Pb (see Supplementary Fig. 5).
The presence of a 2D chiral superconductor could be established experimentally by tunneling into the chiral Majorana modes, whose number is equal to the Chern number of the phase, and which would take place only on the edge of a 2D thin island of Fe on the surface of Pb. Such a technique has been used to image 1D topological edge states of bismuth bilayers in the absence of superconductivity 41 . Similar observations of quasiparticle edge modes inside the superconducting gap will be a strong signature of topological superconductivity proposed in this paper.

Methods
Effective Hamiltonian in the dilute impurity limit. We outline the general formalism by which an effective Hamiltonian for bound states of a Shiba lattice can be derived in the limit where the impurity spacing is large compared with the spacial extend of the bound states of an isolated impurity. We want to derive the low-energy effective theory for a Hamiltonian of the general form where H 0 is the original superconducting BdG Hamiltonian which is gapped (ED) around zero energy, and H 1 is the Hamiltonian for a collection of delta-function impurities at positions r n , where n takes values in a set c, for example, a lattice.
Here, V n are matrices associated with the local degrees of freedom (such as spin and particle-hole) which can induce in-gap states. We implicitly keep r as a d-dimensional vector, so that the formalism is applicable for systems in any dimension d (the same applies to r n , n, m, k, q, R and n b below). We start with the Schrödinger equation for bound state wave functions C It follows that where G 0 (E) ¼ (E À H 0 ) À 1 is the Green function.
Because H 1 is composed of delta-functions for a small set c, G 0 H 1 is nonzero only in the columns corresponding to c, thus equation 17 is equivalent to X n2c G 0 ðE; r; r n ÞV n cðr n Þ ¼ cðrÞ: In the simplest case, if c contains only one single point, labelled by 0, then equation 18 implies The energy of the excitations can be obtained by solving In more complicated cases, the following equation, again implied by equation 18, can serve as the starting point to extract an effective Hamiltonian 8m 2 c : X n2c G 0 ðE; r m ; r n ÞV n cðr n Þ ¼ cðr m Þ; ð22Þ G 0 ðE; r m ; r n Þ ¼G 0 ðE; r m À r n Þ In addition, if 8n: V n ¼ V 0 and r n ¼ nR (nAZ d with d the dimension), equation (22) can be transformed to k-space: where G 0 ðE; qÞ :¼ X nb G 0 ðE; q þ 2n b p=RÞ: ð29Þ n b can be interpreted as the band index when the k-space is folded into the Brillouin zone defined by [ À p/R, p/R). The application of this formalism to the Hamiltonian in equation (7) is detailed in Supplementary Note 2 and results in the Hamiltonian in equation (10).
First principle calculations. We performed electronic structure calculations within the DFT formalism as implemented in the Vienna ab initio simulation package 42 . We used the all-electron projector augmented wave 43,44 basis sets with the generalized gradient approximation of Perdew, Burke and Ernzerhof 45 to the exchange correlation potential. The Hamiltonian contains scalar relativistic corrections, and the spin-orbit coupling was taken into account by the second variation method 46 .
In this work, we chose the host superconductor to be Pb thin film with a (111) surface, and consider different transition metal adatoms. We started by finding the stable configurations of the adatoms on top of the Pb surface. To this end we have compared the relaxation energy (per atom) for an extensive collection of possible configurations, four of which are shown in Supplementary Fig. 5. Based on this energetic consideration and for simplicity, we adopted the configuration with two adatoms per Pb unit cell (highlighted in Supplementary Fig. 5) in our following simulations. When the transition metal element is chosen to be Fe, the details of the configuration are shown in Supplementary Fig. 6. In this configuration, two species of Fe atoms form a buckled honeycomb structure which gains bonding energy due to the short nearest-neighbor distance.
To simulate the composite system and consider the effect of Pb, we used six layers of Pb atoms as substrate in the relaxation calculations with roughly 15 Å vacuum space, taking into account spin-orbital coupling and the ferromagnetic alignment of the Fe moments. We performed DFT calculations with the abovementioned stable configuration. Based on the DFT calculations, we constructed the maximally localized Wannier functions 47,48 for Fe, and obtained a tight-binding model with a band structure that agrees well with the DFT result. We then used the tight-binding model and added a small s-wave superconducting pairing term to it. We computed the Chern numbers of the thus-obtained Bogoliubov-deGennes Hamiltonian to be nonvanishing, as shown in Fig. 4. For completeness, we present a few more Fermi surfaces with different values of E F in Supplementary Fig. 7 to complement Fig. 4.