Topological Superconductivity in Skyrmion Lattices

Atomic manipulation and interface engineering techniques have provided a novel approach to custom-designing topological superconductors and the ensuing Majorana zero modes, representing a new paradigm for the realization of topological quantum computing and topology-based devices. Magnet-superconductor hybrid (MSH) systems have proven to be experimentally suitable to engineer topological superconductivity through the control of both the complex structure of its magnetic layer and the interface properties of the superconducting surface. Here, we demonstrate that two-dimensional MSH systems containing a magnetic skyrmion lattice provide an unprecedented ability to control the emergence of topological phases. By changing the skyrmion radius, which can be achieved experimentally through an external magnetic field, one can tune between different topological superconducting phases, allowing one to explore their unique properties and the transitions between them. In these MSH systems, Josephson scanning tunneling spectroscopy spatially visualizes one of the most crucial aspects underlying the emergence of topological superconductivity, the spatial structure of the induced spin-triplet correlations.

Atomic manipulation and interface engineering techniques have provided a novel approach to custom-designing topological superconductors and the ensuing Majorana zero modes, representing a new paradigm for the realization of topological quantum computing and topologybased devices. Magnet-superconductor hybrid (MSH) systems have proven to be experimentally suitable to engineer topological superconductivity through the control of both the complex structure of its magnetic layer and the interface properties of the superconducting surface. Here, we demonstrate that two-dimensional MSH systems containing a magnetic skyrmion lattice provide an unprecedented ability to control the emergence of topological phases. By changing the skyrmion radius, which can be achieved experimentally through an external magnetic field, one can tune between different topological superconducting phases, allowing one to explore their unique properties and the transitions between them. In these MSH systems, Josephson scanning tunneling spectroscopy spatially visualizes one of the most crucial aspects underlying the emergence of topological superconductivity, the spatial structure of the induced spin-triplet correlations.
The ability to create, control and manipulate topological superconducting phases is quintessential for the realization of topological quantum computing using the non-Abelian braiding statistics of Majorana zero modes [1]. These modes have been observed in one- [2][3][4][5][6] and two-dimensional (2D) topological superconductors [7][8][9], with the latter also exhibiting chiral Majorana edge modes [10][11][12]. Magnet-superconductor hybrid (MSH) systems consisting of chains, islands or layers of magnetic adatoms deposited on the surface of conventional s-wave superconductors, have proven to be suitable experimental systems for (a) the creation of topological superconductivity using atomic manipulation [6] or interface engineering techniques [11], and (b) the study of Majorana modes using scanning tunneling spectroscopy (STS). In particular, 2D MSH systems, with their topological invariant given by the Chern number, are predicted to exhibit a rich topological phase diagram [13][14][15]. However, the experimental ability to tune between different topological phases in 2D MSH systems, essential for exploring the nature of topological superconductivity, has not yet been realized.
In the following, we demonstrate that the ability to tune between topological phases can be achieved in 2D MSH systems containing a magnetic skyrmion lattice by varying the skyrmion radius. As the latter can be experimentally controlled through the application of an external magnetic field [16], the skyrmion MSH system presents an unprecedented opportunity to explore a rich phase diagram of topological superconducting phases, and the transitions between them. The underlying origin for the ability to control the topological phases lies in a spatially inhomogeneous Rashba spin-orbit (RSO) interaction that is induced by the magnetic skyrmion lattice. The induced RSO interaction images the local topological skyrmion charge -the skyrmion number density -and possesses a characteristic spatial signature in the zero-energy local density of states (LDOS) which can be observed at a topological phase transition as well as in the LDOS of chiral Majorana edge modes. Finally, we demonstrate that Josephson scanning tunneling spectroscopy can be employed to visualize one of the most fundamental aspects underlying the emergence of topological superconductivity, the existence of induced spin-triplet superconducting correlations. As 2D skyrmion MSH systems can be built with currently available experimental techniques, our results open new venues for the exploration and manipulation of topological superconductivity and Majorana zero modes.

Results
Theoretical Model. We investigate the emergence of topological superconductivity in a 2D MSH system, in which a magnetic skyrmion lattice (see Fig. 1a) is placed on the surface of a conventional s-wave superconductor, as described by the Hamiltonian where c † rα creates an electron at lattice site r with spin α, and σ is the vector of spin Pauli matrices. We consider a triangular lattice with lattice constant a 0 , chemical potential µ, and hopping amplitude −t rr = −t between nearest-neighbor sites only. ∆ is the superconducting swave order parameter. The spatial spin structure of the skyrmion lattice is encoded in S r [see Supplementary Information (SI) section 1], which represents the direction of an adatom's spin located at site r, and J is its exchange coupling with the conduction electron spin. Note that the creation of Majorana modes in single skyrmions has previously been discussed in Refs. [17,18]. As Kondo screening is suppressed by the full superconducting gap, the spins S r are taken to be classical vectors of length S. We assume that the triangular lattice of skyrmions is commensurate with the underlying triangular surface lattice, thus allowing the skyrmion radius R to take integer and half-integer values of a 0 . Note that in contrast to earlier studies of 2D MSH systems [13][14][15], the above Hamiltonian does not contain an intrinsic Rashba spin-orbit (RSO) interaction. Moreover, due to the broken time-reversal symmetry arising from the presence of magnetic moments, and the particle-hole symmetry of the superconducting state, the topological superconductor belongs to class D [19].
To characterize the topological superconducting phases of the system, we compute the topological invariant -the Chern number C -given by where E n (k) and |Ψ n (k) are the eigenenergies and the eigenvectors of the Hamiltonian in Eq.(1), with n being a band index, and the trace is taken over Nambu and spin-space. Further insight into the origin underlying the emergence of topological superconductivity in skyrmion MSH systems can be gained by considering the spatial structure of the skyrmion and Chern number densities, n s (r) and C(r), respectively. The former is given by yielding a skyrmion number n s = r n s (r). The latter, C(r) [20,21] (see SI section 2), represents the real-space analog of the Berry curvature, and allows a real space calculation of the Chern number [22,23] C = 1/N 2 r C(r) that coincides with that obtained from Eq.(2).
Topological phase diagram. A crucial aspect for the emergence of topological superconductivity in 2D skyrmion MSH systems is that the magnetic skyrmion lattice induces an effective, spatially varying Rashba spin-orbit interaction. To demonstrate this, we apply a unitary transformation [24] to the Hamiltonian in Eq.(1) (see SI section 2) that rotates the local spin S r to thê z axis, yielding an out-of-plane ferromagnetic order and a spatially inhomogeneous RSO interaction, α(r) (see Fig.1b). α(r) possesses the same spatial structure as the skyrmion number density, n s (r), (see Fig.1c) -reflecting its origin in the local topological charge of the skyrmon lattice -with its largest value, α max = πa 0 t/(2R), in the center of the skyrmion, and a vanishing α(r) at the corners of the skyrmion lattice Wigner-Seitz unit cell. The existence of a non-zero α(r), of an out-of-plane ferromagnetic order in the rotated basis, and of a hard s-wave gap, satisfies all necessary requirements for the emergence of topological superconductivity [13][14][15], resulting in the rich topological phase diagram shown in Fig. 2.
The phase diagram in the (µ, R) plane (see Fig. 2a) reveals the intriguing result that it is possible to tune a skyrmion MSH system between different topological phases by changing the skyrmion radius R, which can be experimentally achieved through the application of an external magnetic field [16]. This unprecedented ability arises from the facts that (a) varying the skyrmion radius leads to changes in the induced α(r), and (b) in contrast to MSH systems with a homogeneous ferromagnetic structure, topological phase transitions in magnetically inhomogeneous MSH systems (as given here) are controlled not only by µ and J, but also by α. Indeed, the results in Fig. 2a reveal that the phase transition lines in the (µ, R) plane are determined by µ = A i +B i /R 2 (see dashed lines) with constants A i , B i . Since α max ∼ 1/R, our result suggests that the induced RSO interaction leads to an effective renormalization of the chemical potential [17], thus facilitating the ability to tune between topological phases. This dependence of the phase transition lines on R is also revealed when considering the phase diagrams in the (µ, JS) plane for different skyrmion radii (see Fig. 2b). These phase diagrams show a very similar structure of topological phases for different R, with the phases moving to lower values of µ with increasing R. We note that the topological phases that are accessible through tuning of R strongly depend on JS (see Fig. 2c): for sufficiently large JS, every change in the skyrmion radius by a half-integer leads to a change in the system's Chern number. Thus, a rich topological phase diagram can be accessed and explored through changes in the skyrmion radius R.
The inhomogeneous magnetic structure of the MSH system also allows us to reveal an intriguing connection between the Chern number density, C(r), which is a local marker for the topological nature of the system, and the Berry curvature in momentum space. In particular, the spatial structure of C(r) (see Fig. 1d) reflects that of the skyrmion lattice, but is complementary to that of the induced α(r) (see Fig. 1b), with the maximum in C(r) occurring at the corners of the skyrmion lattice unit cell. However, it is in these regions that the lowest energy states possess their largest spectral weight (see discussion below), establishing a real space analogue of the observation that the lowest energy states in momentum space in general possess the largest Berry curvature. Electronic structure at a topological phase transition The real space structure of the induced RSO interaction, and hence that of the local topological skyrmion charge, is reflected in the electronic structure of the MSH system, and becomes particularly evident at a topological phase transition. To demonstrate this, we consider the transition between a C = 8 and C = 6 phase, as indicated by the solid black dot in Fig. 2c. While the system possesses a topological gap on either side of the transition, the gap at the transition closes at the K, K -points (see Fig. 3a of the lowest energy band (see Fig. 3b) in the reduced Brillouin zone (RBZ). This gap closing is reflected in a unique spatial and energy structure of the zero-energy local density of states (LDOS) [see (xy)-plane in Fig. 3c]. In particular, the spatial structure of the LDOS reveals that the largest (smallest) spectral weight of the zeroenergy state, associated with the phase transition, is located where the induced RSO interaction is the smallest (largest), at the corners of the Wigner-Seitz unit cell (the skyrmion center). Thus, the spatial pattern of the zero-energy LDOS is complementary to that of the local topological skyrmion charge, n s (r). Moreover, as the topological gap in general increases with increasing RSO interaction, we find that the large induced RSO interaction in the skyrmion center leads to a dome-like region in energy in which the LDOS is suppressed [see (x, E)-and (y, E)-planes in Fig. 3c]. The electronic structure of the skyrmion MSH systems also provides a unique example to demonstrate that the multiplicity m of the momenta in the Brillouin zone, at which the gap closing occurs, determines and is equal to the change in the Chern number at the transition. For the time-reversal invariant (TRI) Γ, M, (K, K ) points, the multiplicity is m = 1, 3, 2 (note that by symmetry, a gap closing at the K point implies a gap closing at K as well), respectively, as each M (K, K ) point is shared by 2 (3) BZs, leading to a change in the Chern number by ∆C = 1, 3, 2 at the transition.
Gap closings can also occur at non-TRI points, e.g., at points along the Γ − M line (see SI section 3), which possess a multiplicity of m = 6, resulting in a change of the Chern number by ∆C = 6. While all of the above gap closings exhibit a Dirac cone (see Fig. 3a), there also exist gap closings that exhibit lines of zero-energy (see SI section 3), rather than discrete zero energy Dirac points. These gap closings, however, are not accompanied by a change in the Chern number.
MSH system with a skyrmion ribbon To study the emergence of chiral Majorana edge modes in a skyrmion MSH system, we next consider a skyrmion ribbon placed on the surface of an s-wave superconductor (see Fig. 4a). In a topological phase with Chern number C, the bulkboundary requires that such a MSH system possess |C| chiral Majorana edge modes per edge. These modes traverse the superconducting gap and disperse linearly near the chemical potential as a function of the momentum along the ribbon edge, as shown in the inset of Fig. 4b for the C = 3 phase. A spatial plot of the zero-energy LDOS (see Fig. 4b) demonstrates that the chiral Majorana mode is as expected localized along the edges of the ribbon, and that its spatial structure is complementary to that of the local skyrmion topological charge. The spatial structure of the skyrmion lattice, and hence of the induced α(r), is also reflected in the combined energy and spatial dependence of the LDOS (see Fig. 4c) as revealed by a line-cut of the LDOS from the bottom to the top of the ribbon along x = 0 (left edge of Fig. 4b).
In particular, in the center of the skyrmions, where α is the largest, the spectral weight in the LDOS is pushed to higher energies. The spatial structure of the LDOS is therefore similar to that exhibited by the MSH system at a phase transition (see Fig. 3c).
In addition to the chiral Majorana edge modes, the magnetic structure of the skyrmion ribbon leads to two unique physical features. The first one is the spatial form of persistent supercurrents that are induced by the broken time-reversal symmetry. While these supercurrents are generally confined to the edges of an MSH system [15], the inhomogeneous magnetic structure of the skyrmion lattice leads to supercurrents (see SI section 4) that circulate each skyrmion, not only along the ribbon's edge, but also in its interior (see Fig. 4d). These supercurrents screen the out-of-plane component of the local magnetic moments, similar to the case of a vortex lattice, and are carried by both the in-gap and bulk states. The second unique feature is the presence of spintriplet superconducting correlations which are a crucial requirement for the emergence of topological superconductivity [15]. The development of Josephson scanning tunneling spectroscopy (JSTS) [25][26][27][28][29] has provided a unique opportunity to visualize not only these correlations in real space at the atomic level, but also to investigate the effects of the inhomogeneous magnetic structure of the skyrmion lattice on the superconducting swave order parameter, ∆(r) [30]. Specifically, pair breaking effects of the magnetic moments lead to a spatially non-uniform suppression of ∆(r) inside the skyrmion ribbon (see Fig. 4e), with the largest suppression occurring where the induced RSO interaction is the weakest. This spatial structure of ∆(r) is well imaged by that of the critical Josephson current, I c (r) (see Fig. 4f), measured via JSTS using a tip with an s-wave superconducting order parameter [30], thus providing direct insight into the strength of local pair breaking effects. Moreover, the inhomogeneous magnetic structure of the skyrmion lattice enables the emergence of superconducting spintriplet correlations not only in the equal-spin channels |↑↑ and |↓↓ (corresponding to Cooper pair spin states S z = ±1), but also in the mixed-spin (S z = 0) channel, |↑↓ + |↓↑ (see SI section 5). The spatial structure of the real and imaginary parts of these correlations in the |↑↑ channel are shown in Figs. 4g and h, respectively (the correlations in the S z = 0, −1 channels are shown in SI section 5). These correlations are a direct consequence of the magnetic structure of the skyrmions, and thus vanish outside the ribbon. To image the spatial structure of these non-local triplet correlations, we compute I c (r) using an extended (2 × 1) JSTS tip with a superconducting triplet order parameter (see SI section 5). If the tip's order parameter is chosen to be either purely real (see Fig. 4i) or purely imaginary (see Fig. 4j), the resulting Josephson current very well images the spatial structure of the real or imaginary parts, respectively, of the superconducting triplet correlations. We note that these triplet correlations can be imaged despite the fact that the MSH system possesses neither a long-range nor a local triplet superconducting order parameter. Thus JSTS can provide unprecedented insight into the existence of one of the most crucial aspects of topological superconductivity, the existence of spin-triplet correlations.
Discussion MSH systems containing a magnetic skyrmion layer are suitable candidate systems to explore a rich topological phase diagram. By varying the skyrmion radius, which can be achieved through the application of an external magnetic field, it is possible to tune these systems between different topological phases, and explore not only their unique properties, but also the transitions between them. The origin of this tunability lies in the spatially inhomogeneous Rashba spin-orbit interaction, which is induced by the magnetic skyrmion lattice and which carries a spatial signature in the zero-energy LDOS that can be observed at a topological phase transition. Skyrmion MSH systems also provide a unique opportunity to employ Josephson scanning tunneling spectroscopy to visualize one of the most fundamental aspects underlying the emergence of topological superconductivity, the existence of induced spin-triplet superconducting correlations.
This, in turn, yields a new experimental approach to identifying topological superconducting phases, a possibility that needs to be further explored in future studies. These results demonstrate that the tunability of the magnetic structure in MSH systems opens a new venue for the quantum engineering and exploration of topological superconductivity, and the ability to engineer Majorana zero modes and chiral Majorana edge modes. This raises the intriguing question of whether a tuning of topological phases, similar to the one discussed here, could also be achieved using other non-collinear magnetic structures [24,31,32].

Section 1: Spin structure of the skyrmion lattice and the induced Rashba spin-orbit interaction
The local spin at site r in the skyrmion lattice is given by S r = S(sin θ r cos φ r , sin θ r sin φ r , cos θ r ) (S1) with θ r and φ r being the polar and azimuthal angles, defined by where r = (x, y), R i = (x i , y i ) is the position of the skyrmion center closest to r, and R is the skyrmion radius.
To demonstrate that the skyrmion lattice induces an effective Rashba spin-orbit interaction, we perform a local unitary transformation that rotates the local spin S r into theẑ-axis, perpendicular to the plane of the MSH system. The corresponding unitary transformation is defined via where the unitary matrixÛ r is given byÛ The Hamiltonian in this basis is then given by The effective hopping is with α rr being the induced Rashba spin-orbit interaction.

Section 2: Chern number density in real space
Insight into the interplay between the inhomogeneous magnetic structure of the skyrmion lattice and the local topological nature of the system can be gained by considering the spatial Chern number density [1][2][3][4][5] given by Tr τ,σ [P [δ x P, δ y P ]] r,r .
The Chern number is then computed via C = 1/N 2 r C(r), and coincides with that obtained from Eq.(2) in the arXiv:2005.00027v1 [cond-mat.supr-con] 30 Apr 2020 main text. Here, Tr τ,σ denotes the partial trace over spin σ and Nambu space τ , where the projector P onto the occupied bands is given by P = α=occ. |ψ α ψ α | for a real-space N × N lattice. The c m 's are central finite difference coefficients for approximating the partial derivatives. The coefficients for positive m can be calculated by solving the following linear set of equations for c = (c 1 , . . . , c Q ): For negative m, we have c −m = −c m . We have taken the largest possible value of Q = N/2. C(r) as defined above thus represents the real-space analog of the Berry curvature F(k), and was previously introduced to discuss the topological phases of Chern insulators [1] and of disordered topological superconductors [5].
An alternative form of the Chern number density can be derived by realizing that the periodicity of the skyrmion lattice can be employed to rewrite the Hamiltonian of Eq.(1) of the main text in the following form (we consider a hopping −t between nearest-neighbor sites only, with δ being the lattice vector between nearest neighbor sites) Here, the primed sum runs over all locations of the unit cell of the skyrmion lattice, and Ψ † k , Ψ k are spinors in Nambu space, spin space, and positions r in the unit cell. Moreover, RBZ is the reduced Brillouin zone whose wave-vectors are defined via where M i (N i ) is the number of skyrmion unit cells (number of lattice sites) along the lattice vector a i defined via and b i are the reciprocal lattice vectors given by Finally, the fermionic operators are defined via with N = N 1 N 2 , and L i is the length of the skyrmion unit cell (in units of a 0 ) along a i .
In analogy to the calculation of the Chern number in Eq.(2) of the main text, we can now define the Chern number density (i.e., a "local" Chern number), by diagonalizingĤ(k) in Eq.(S9), yielding r, τ, σ|P k [∂ kx P k , ∂ ky P k ]|r, τ, σ (S15a) where E n (k) and |Ψ n (k) are the eigenenergies and eigenvectors ofĤ(k), respectively, with n being a band index, and the summation is taken over Nambu and spin-space.

Section 3: Topological Phase Transitions and Gap Closings
In the topological phase diagram of the skyrmion MSH system, we observe two different types of gap closings. The first one occurs via a Dirac cone, with the Dirac point being located at (a) time-reversal invariant (TRI) Γ, M, K, K points in the Brillouin zone (an example of this is shown in Fig.3 of the main text), or (b) at a non-time-reversal invariant point along the Γ − M line. An example of the latter is shown in Fig. S1a, where we present the dispersion E k of the lowest energy band at the transition between the C = 4 and C = −2 phases (see black line labelled a in Fig. S1e). Here, the multiplicity of the points in the BZ where the gap closes is m = 6 (as none of the points lie on the boundary with another BZ), implying that the Chern number changes by ∆C = 6 at the transition. To gain further insight into the nature of the topological phase transition, we rewrite the expression for the Chern number in Eq.(2) of the main text as and present the momentum dependence of C(k) in Fig. S1b on either side of the transition. This plot reveals that the change in the Chern number by ∆C = 6 arises from the same change in C(k) near each of the gap closing points. The second type of gap closing in the phase diagram (see black line labelled c in Fig. S1e) involves a line of momenta along which the gap vanishes (and thus not discrete Dirac points) in the Brillouin zone, as shown in Fig. S1c. Such a gap closing does not lead to a change in the Chern number between the topological phases adjacent to it, which is also reflected in the momentum dependence of C(k) (see Fig. S1d), which does not show any qualitative change between the two phases.

Section 4: Supercurrents in the skyrmion lattice
The persistent supercurrent associated with the hopping of an electron with spin σ from site r s to a nearest neighbor site r s + δ can be computed via where −t is the spin-independent hopping parameter between sites r and r s + δ, and g r,a,< mn (r, r + δ, ω) are the (m, n) elements in Nambu space of the retarded, advanced, or lesser Green's function matrices. The Green's function matrix in Matsubara time is defined viaĝ where the spinors are defined via To obtain the above Greens functions for the system, we diagonalize the real space Hamiltonian in Eq.(1) of the main text, yielding energy eigenvalues E k and eigenvectors u mk (r). The Greens functions can then be computed using Note that the energy eigenvalues come in pairs ±E k , and that the summation in the above equation runs over all eigenvalues. The elements (m, n) of the Greens function matrix used in Eq.(S17) is determined by the spin index σ with n = m = 1 for σ =↑, and n = m = 2 for σ =↓. The total supercurrent between neighboring sites is then given by and is presented in Fig. 4d of the main text.
The spatial form of the local superconducting correlations in the s-wave channel are shown in Fig. 4e, and for the real and imaginary parts of c r+δ,↑ c r,↑ in Figs. 4g and h, respectively, of the main text. The real and imaginary parts of c r+δ,↓ c r,↓ and of c r+δ,↑ c r,↓ + c r+δ,↓ c r,↑ are shown in Figs. S2a, b and Figs. S2c, d, respectively. To spatially image these correlations via JSTS, the superconducting order parameter in the JSTS tip needs to possess the same symmetry and spin-structure as the correlations we intend to probe. For the spin-singlet, s-wave channel, this has previously been described in Ref. [7]. In this case, a JSTS tip ending in a single site is sufficient to probe the s-wave correlations, and the resulting spatial structure of the critical Josephson current, I c , is shown in Fig. 4f of the main text. In contrast, the triplet superconducting correlations are by definition odd in real space, implying that they are non-local. To probe them via JSTS thus requires a tip that also exhibits a non-local superconducting order parameter (similar to the case of a d x 2 −y 2 -wave order parameter, discussed in Ref. [6]). This implies that the JSTS tip has to end in at least 2 sites from which electrons can tunnel into the system, a case we will consider in the following (we refer to such a JSTS tip as a 2-site tip). As we have shown before [6], using tips with a larger number