Topological superconductivity with deformable magnetic skyrmions

Magnetic skyrmions are nanoscale spin configurations that are efficiently created and manipulated. They hold great promises for next-generation spintronics applications. In parallel, the interplay of magnetism, superconductivity and spin-orbit coupling has proved to be a versatile platform for engineering topological superconductivity predicted to host non-abelian excitations, Majorana zero modes. We show that topological superconductivity can be induced by proximitizing skyrmions and conventional superconductors, without need for additional ingredients. Apart from a previously reported Majorana zero mode in the core of the skyrmion, we find a more universal chiral band of Majorana modes on the edge of the skyrmion. We show that the chiral Majorana band is effectively flat in the physically relevant parameter regime, leading to interesting robustness and scaling properties. In particular, the number of Majorana modes in the (nearly-)flat band scales with the perimeter length of the system, while being robust to local disorder. Skyrmions, Majorana fermions and topological superconductivity are active areas of research within the condensed matter physics community and each offer the promise of new physics and potential applications. Here, the authors theoretically demonstrate that by bringing skyrmions into close proximity with a superconductor topological superconductivity can be induced, as well as Majorana states localized at the edges of the skyrmions.

M agnetic skyrmions are nano-scale or meso-scale whirling spin configurations of topological nature which gives them some stability and long lifetime. Magnetic skyrmions have been found in a variety of non-centrosymmetric magnets 1 , in ultrathin magnetic films [2][3][4] , as well as in multiferroic insulators [5][6][7] . Quite remarkably, magnetic skyrmions can be stabilized over a wide temperature domain ranging from room temperature 8,9 to cryogenic temperature 2,3,10 . Evidence that magnetic skyrmions can be driven by ultralow electric current densities 11,12 make them promising candidates for future spintronic applications 13 .
In parallel to these developments, the search for Majorana modes in condensed-matter systems has been the focus of great attention, motivated by their potential application in quantum computation. Various systems have been considered as hosts for topological superconductivity and Majorana modes, based on the paradigm of combining ferromagnetic order with strong spinorbit coupling and conventional superconductivity. The paradigm led to successes in predicting [14][15][16][17][18] and experimentally indicating Majorana zero-energy modes at endpoints of one-dimensional systems, such as iron atomic chains [19][20][21] and semiconducting wires 22 . Recent experiments have extended the paradigm to two dimensions, reporting some evidence for dispersive Majorana edge states around two-dimensional magnetic domains using cobalt atom clusters under monolayer lead 23 or iron adatom clusters on a rhenium surface 24 . Since the long-term goal is a flexible platform for manipulation of Majorana modes, two challenges for the paradigm are that the preformed structures (clusters, wires) are hard to manipulate, and that the systems are constrained by requirement of strong spin-orbit coupling.
An alternative approach to engineering topological superconductors while circumventing these two challenges could be to remove the spin-orbit coupling ingredient, and instead consider a non-collinear magnetic texture proximitized by a conventional superconductor 14,[25][26][27][28][29] . In fact, our results are relevant for a broader class of skyrmion-like textures, such as magnetic bubbles [30][31][32][33] . In addition, a texture such as a skyrmion can be manipulated by external fields, potentially facilitating the manipulation of Majorana states. Yang et al. recently found that skyrmions having an even azimuthal number can indeed bind a single Majorana zero mode in their core 26 . Moreover, very elongated magnetic skyrmions were shown to host Majorana zero modes at their endpoints 34 . In contrast, we find here that a magnetic skyrmion of any azimuthal winding and sufficient radial winding gives rise to a single band of states at the edge of the skyrmion, i.e., a chiral Majorana edge mode (CMEM).
Surprisingly, for the physically relevant range of parameters (skyrmion size, winding numbers, magnetic coupling strength) the CMEM has negligible velocity, i.e., it is nearly a Majorana flat band (MFB). Furthermore, we find that the CMEM is robust to local perturbations, as well as to smooth deformations of the edge geometry. Such deformations preserve the number of edge states proportional to the perimeter length of the edge.
For systems with translational symmetry there is a theoretical classification of topological superconducting phases, and predictions for a corresponding CMEM along a given edge of the system [35][36][37] . Furthermore, the existence of a MFB along an edge can be deduced from an appropriate discrete chiral symmetry and topological indices in lower spatial dimension [38][39][40] . In our case, the skyrmion is an inhomogeneous texture so these methods cannot be directly used to explain the observed robustness and near-flatness of the CMEM. We however deduce the underlying topological protection of the skyrmion's CMEM by a mapping to a cylinder geometry. Although this construction requires rotational symmetry of the skyrmion, the CMEM by its nature provides robustness against small deformations of the shape of the system. Further, we identify the chiral symmetry that would protect a strict MFB (instead of a CMEM), and show that this symmetry is only weakly broken by the skyrmion texture, leading to a nearly flat CMEM and providing further protection against low-energy perturbations. Finally, we will discuss potential material realizations, and possibilities for manipulation of Majorana states within the nearly-flat CMEM.

Results
Setup and model. Consider a two-dimensional (2D) magnetic thin film hosting a skyrmion, which is represented by a classical magnetization texture n r ð Þ ¼ sin f ðrÞcosðqθÞ; sin f ðrÞsinðqθÞ; cos f ðrÞ ð Þ ; ð1Þ written in polar coordinates r = (r, θ), where f(r) is a radial profile that we will specify shortly. We study such a thin film proximitized by a conventional s-wave superconductor (Fig. 1a). The electrons interact with the skyrmion texture via a direct exchange interaction of strength J. In the Nambu basis Ψ y r ð Þ ¼ ψ y " r ð Þ; ψ y # r ð Þ; ψ # r ð Þ; Àψ " r ð Þ , where ψ σ (r) annihilates an electron with spin σ at position r in 2D, the total Hamiltonian H can be written as Bogoliubov-de-Gennes (BdG) Hamiltonian H r ð Þ defined as where the s-wave superconducting order parameter Δ 0 is taken real without loss of generality, the electron effective mass is m and the chemical potential is μ. We set ħ = 1 unless explicitly written otherwise. The σ α and τ α (α = x, y, z) are Pauli matrices acting in spin and particle-hole space, respectively. We thus assume that the skyrmion affects the electrons only through the exchange field, which is justified in the limit of strong local exchange interaction. We further consider the limit where the effective spin-orbit coupling is dominated by the one induced by the skyrmion exchange field. This is justified if both intrinsic and Rashba spin-orbit couplings are relatively weak, which we argue in the discussion section to be the case in a typical superconducting material such as aluminum. The inclusion of intrinsic and Rashba spin-orbit couplings would require another in-depth study due to the loss of rotational symmetry, and due to the non-trivial interplay of different effective spin-triplet pairings introduced by these couplings. In the limit that we consider in this work, in absence of skyrmion the spin-orbit length in the superconductor l so = 1/(mα), with α a spin-orbit amplitude, is much larger than the typical lengthscale of the skyrmion, so that the magnetoelectric coupling and appearance of vortices can also be neglected [41][42][43][44][45] .
The skyrmion is parametrized by three numbers: the radial winding number p, which counts the number of spin flips as one moves radially away from the core of the skyrmion; the azimuthal winding number q, which counts the number of spin flips as one winds around the origin; and finally the skyrmion radius R sk , which determines its size. We consider a hard-wall boundary condition at the edge i.e., at r = R sk . Formally this can be realized by having the exchange J = 0 in the magnetic insulator outside the edge, which might be experimentally unattainable. However a simple alternative is to deposit the superconductor in form of an island, whose edge would naturally become the edge in our model. In such a setup it is natural to consider various geometrical shapes of the edge given a fixed underlying magnetic texture.
For simplicity, the function f(r) is chosen to be linear, defining a straightforward skyrmion texture as in Fig. 1a. As we show in the discussion section, the exact shape of f(r) has weak influence on our conclusions, allowing us to extend our results to a broader class of textures including magnetic bubbles.
Skyrmion edge states and topological superconductivity. We first solve the model in Eq. (2) by using rotational symmetry. Our model in Eq. (2) has a rotational symmetry in the combined realspace and spin-space, given by the conserved total angular momentum J z ¼ L z þ q 2 σ z where L z = −i∂ θ is the orbital angular momentum. Using the eigenvalues of J z , denoted as m J , the Hamiltonian in Eq. (2) becomes an effective one-dimensional radial model that can further be discretized and diagonalized numerically (technical details are given in "Methods" section "Radial tight-binding model"). In the regime of J/Δ 0 large enough Fig. 1b shows the resulting strong peak in the local density of states (LDoS), near zero energy and at the skyrmion's edge. Such a spectral feature was observed before in related models 14,25,26,46 . Further, the LDoS clearly displays a reduced gap Δ eff ≈ 5%Δ 0 consistent with an effective (topological) p-wave superconducting gap 47 , see estimate below), and a full gap Δ 0 develops without any edge states. This is expected in a transition from a topological superconducting phase to trivial superconductivity. The regime of topological p-wave superconductivity is also consistent with our finding of other ingap states localized near the skyrmion core that only appear when the edge modes appear. We therefore interpret these states near the core as analogs of states bound to magnetic impurities (here, inhomogeneities of the skyrmion texture), which are only expected for p-wave pairing, but are absent in the s-wave-pairingdominated trivial phase (J/Δ 0 small enough).
We further clarify the edge states and topological superconductivity by looking at the spectrum ε(m J ), in which edge states form a seemingly flat band in a range of m J values, see Fig. 2a, b. Importantly, the edge states appear for any value of the azimuthal winding number q (on the other hand, p always needs to be high enough 26 , we showcase p = 6). Note that in contrast, we find a single Majorana zero mode at the core of the skyrmion only if the skyrmion's azimuthal winding number q is even 26 . This is easily understood since the zero mode must appear in the self-conjugate angular momentum m J = 0 sector, while m J is quantized to be integer (resp. half-odd-integer) when q is even (odd) due to the single-valuedness of the wavefunction. The existence of edge states indicates that skyrmions of any q induce topological superconductivity.
Topological origin and the near-flatness of edge mode. In order to explain the origin of the edge states we use a procedure introduced by Wu et al. 48 to smoothly deform the model in Eq.
(2) defined on the disk to another model defined on the cylinder via the cone geometry as represented in Fig. 2a. Taking the cylinder limit (see "Methods" section "Gradient and Laplace operators in the cone geometry") effectively focuses on the edge of the skyrmion at the price of disregarding the skyrmion core area, which is replaced by an artificial edge.
Explicitly, we use the rotation symmetry, i.e., the total angular momentum m J basis (see "Methods" section "Radial tight-binding model"), then we apply the unitary transformation U(r) = exp (iσ y f(r)/2) to align the exchange field with the z-axis at each point of radial distance r, and finally we apply the mapping to the cylinder. The resulting HamiltonianH cyl m J ðrÞ can be written as the sum of three parts,H cyl For our purpose it is sufficient to show that the edge modes and the effective gap (the energy gap in the m J = 0 sector) of the original model Eq. (2) are connected to such features of a model derived from the cylinder mapping. Therefore, in what follows we safely neglect the part in Eq. (5) since these are a small overall chemical potential renormalization and small overall boundary term. For a given angular momentum m J , the Hamiltonian H wire m J ðrÞ can be interpreted as the extensively studied Hamiltonian of a Rashba wire 49,50 upon introducing a momentumdependent chemical potential μ m (Note: however that the skyrmion-induced effective spin-orbit coupling in 2D is not of a simple Rashba type.) At each m J the superconducting wire Hamiltonian H wire m J ðrÞ is well known to be where all energies are in units of the bandwidth t, all distances are in units of the lattice spacing a (see Supplementary Note 1 and Supplementary Fig. 1 for details). For precisely jm J j ¼ m Ã J the wire is at the topological transition and has a gapless spectrum, giving our model a bulk-gap-closing point as shown in Fig. 2c. The MFB found here has a protection by a chiral symmetry, as MFB's were found to have in models with translational symmetries 39,40,51 . Note that the wire Hamiltonian and its MFB become a correct model for our texture Eq. (1) if we choose q = 0 and thereby nullify the H slope m J ðrÞ term. Physically, this is a special case where instead of the skyrmion shape the texture becomes coplanar (in the xz-plane, see Eq. (1)), and the orthogonal direction provides a chiral operator that anticommutes with the Hamiltonian (see Eq. (2)). Since all the MFB states have the same chirality, they cannot hybridize among themselves. It is difficult to remove the MFB states 51 , namely, a perturbation must have energy larger than the effective gap; or, it should hybridize the MFB with low energy bulk states at jm J j ¼ m Ã J , which are few; or, chirality symmetry must be broken (out-of-xz-plane exchange field). We note that the proof of existence of the MFB rests on the rotational symmetry of the q = 0 coplanar texture, since this symmetry provides the m J quantum number. Consider now deformations of the shape of the edge imposed on our q = 0 coplanar texture. These geometric deformations would generally mix the m J sectors, yet the described stability of the MFB implies that the deformations would be inefficient in removing the MFB states.
We can now proceed to the relevant model for a skyrmion with arbitrary q ≠ 0:H cyl;eff The single term H slope m J ðrÞ breaks the chiral symmetry Ξ, and there are no other chiral operators. The term H slope m J ðrÞ exactly contributes an energy ε edgestate (m J )~m J to an MFB state at m J , making the flat MFB into a linearly dispersing chiral Majorana edge mode (CMEM) of the q ≠ 0 skyrmion (Fig. 2c). The single CMEM itself has general robustness to perturbations, however, we additionally find that the velocity of the CMEM is very small in the relevant physical regime, i.e., the breaking of chiral symmetry is very weak. Qualitatively, we can estimate the upper limit on energy ε * that the CMEM can have, which occurs at the maximal m J of the CMEM, i.e., ε Ã jε edgestate ðjm Ã J jÞj. Treating H slope m J ðrÞ as a first order perturbation to the MFB (see Supplementary Note 2 and Supplementary Fig. 2), the estimate scales the same way with skyrmion size as the estimate of the effective gap Δ eff~p /R sk . For the relevant regime of J, Δ 0 , μ (see "Discussion" section) the quantitative ratio is at most ε * /Δ eff~0 .1. The corresponding Fermi velocity of the CMEM is therefore small and suppressed by the skyrmion size, We thus demonstrated that at low energy the single edge mode of the skyrmion can be connected to the single CMEM of a cylinder made of Rashba wires, and the CMEM is nearly a MFB. The phase diagram of both models (skyrmion model vs. wires on cylinder) obtained by varying J/t in the radial tight-binding setup are compared in Supplementary Note 3 and Supplementary Fig. 3 and show excellent agreement. Importantly, in both the original skyrmion and the cylinder model and for small enough systems as shown in Fig. 2b, we observe the angular momentum value m Ã J in accordance to predictions in Eq. (6), and we observe the nearflatness of the edge mode.
Edge states on deformed edges. The number of states in the single CMEM of a perfectly rotationally symmetric system is given simply by the highest angular momentum that is reached by the edge states, i.e., m Ã J , and therefore scales linearly with the perimeter of a disk-shaped system centered on the rotationally symmetric skyrmion (see Eq. (6), neglecting corrections of order 1 due to μ depending on m Ã J ). We remind that the edge of the system is defined by setting the exchange J = 0 outside it, or equivalently, by depositing a superconducting island with that edge shape on top of the underlying magnetic material. If the nearly-flat CMEM is indeed robust, we hypothesize that geometric deformations of the edge would preserve the scaling of number of states in the CMEM with the perimeter of the deformed edge.
We substantiate the perimeter hypothesis with an extensive analysis of a 2D square-lattice tight-binding discretization of the skyrmion model Eq. (2), which upon setting the skyrmion exchange strength J to zero outside the skyrmion edge, i.e., radius R sk , gives consistent results with the radial model (see "Methods" section "2D tight-binding Hamiltonian", Supplementary Note 4 Counting the number of edge states of circular skyrmion and Supplementary Fig. 4).
Next we consider two more geometries where the edge of the system is far from a circle and count their edge states as the overall system size is varied (see Supplementary Note 4 Defining the geometries and edge state counting and Supplementary  Fig. 5).
The results for the number of edge states vs. the perimeter of the edge are displayed in Fig. 3. It clearly shows that the number of edge states scales linearly with the perimeter of the edge for all three geometries considered, with a mean slope of 0.12(13)a −1 . The inverse slope is a lengthscale ξ associated to the edge state. We find that ξ ≈ 0.5λ for the parameters considered, where the lengthscale λ measures the distance for a single radial spin flip, i.e., R sk = pλ. This is consistent with the observed localization length of edge states in the radial direction. This typical radial width of the edge states thus ranges from a few nanometers for the skyrmions depicted in Fig. 3 to a few tens of nanometers for the skyrmion depicted in Fig. 1.
To further investigate the robustness of the states forming the single nearly-flat CMEM, we notice that the states in the CMEM seem to locally hybridize where the shape of the edge has sharp features. Sharp features in the edge shape allow the edge-state wavefunctions to overlap as they decay perpendicularly to the edge. Therefore the "elastic perimeter law" demonstrated in Fig. 3 is best exhibited when the curvature of the edge is constant on lengthscales comparable to the extent of a single edge-state wavefunction ξ, as we additionally confirm through an investigation of elliptical skyrmion geometries (see Supplementary Note 4 Elliptic geometry and Supplementary Fig. 6).
Sharp corners in the shape of the edge seem a stronger perturbation than uncorrelated scalar disorder, since we numerically show using the 2D tight-binding model that the nearly-flat CMEM is indeed robust to uncorrelated scalar disorder (see Supplementary Note 5 and Supplementary Fig. 7).

Discussion
In summary, we have shown that a system composed of a magnetic skyrmion coupled to a conventional s-wave superconductor realizes a topological superconducting phase with a nearly dispersionless chiral Majorana mode at its edge. Deforming the edge of the skyrmion away from a circular shape shows that the number of edge states can be tuned and scales linearly with the perimeter of the edge.
As skyrmions usually appear in ferromagnetic thin films, we also considered the effect of a ferromagnetic background on the edge states. For this purpose, in the radial tight-binding model we move the boundary of the system farther than the edge of the skyrmion, filling the added space with a ferromagnetic exchange field without changing the strength of the interaction J. We find that the edge states initially localized at the edge of the skyrmion delocalize in the background, as seen in Supplementary Note 6 and Supplementary Fig. 8. This can be understood rather simply because the superconductor is gapless in that region. The delocalization of the edge states is consistent with the analytical treatment of Yang et al. 26 .
Our analysis is carried out in the Bogoliubov-de Gennes formalism without self-consistency. We believe that a selfconsistent calculation would not change our main conclusions since self-consistent calculations 52,53 on similar systems related to one-dimensional wire Hamiltonians, to which we map the skyrmion, didn't show any qualitative change of the physics. The only effects would then be expected near the topological phase transition where the gap is small. In our case, this may for Perimeter/a Number of edge states Fig. 3 Edge states of deformed structures. Number of low-energy edge states as a function of the perimeter for each of the three geometries of the edge geometry 1, 2, and 3. Geometry 1 consists in a circular skyrmion while geometry 2 and 3 are made of two overlapping disks obtained from a larger skyrmion by setting the exchange field to 0 outside the desired region. The radius of the additional disk is equal to that of the central disk for geometry 2 and half this for geometry 3 (see Supplementary Note 4 Defining the geometries and edge state counting). In all geometries, the underlying skyrmion texture has azimuthal winding number q = 1. The other model parameters are the spin-flip length λ/a = 16, the s-wave order parameter Δ 0 /t = 0.1, the chemical potential μ/t = 0 and the exchange interaction strength J/t = 0.2. For geometry 1, the radial winding number is p = 9 while for geometry 2 and 3, the skyrmion has p = 6 and p = 9, respectively. The graphic for black disks (geometry 1) shows the real-space image of the local density of states of one typical low-energy state. Graphics for red triangles and blue squares (geometry 2 and 3, respectively) show the local density of states averaged over the 30 lowest-energy states. The colorbar represents the local density of states from low (black) to high (red) values with different scales for the three geometries. Top inset: linear slope extracted for each geometry. The gray shading indicates error bars from the fitting procedure. The black horizontal line is the average slope estimated to be 0.12(1)a −1 where a is the lattice spacing example slightly shift the value of m J * defined in Eq. (6) at which the gap closes. Furthermore, we assumed that our skyrmion arises in a magnetic insulator which guarantees that the mutual interplay between the magnetic insulator and the superconductor is weak.
The chirality of our CMEM is determined by the azimuthal winding number q, so the questions arise whether the CMEM appears in other textures that have azimuthal winding, and whether chiral materials are necessary. First, note that our texture definition in Eq. (1) may describe both Bloch and Néel skyrmions by adding a constant phase shift, named "helicity" 1 that can be unitarily removed from our model and does not affect the discovered spectrum nor the wavefunction localization. Second, we have so far focused on skyrmions but our findings also apply to magnetic bubbles that have a different microscopic stabilizing mechanism but have the same topology [30][31][32][33] . For our purpose, the key distinguishing aspect of bubbles is their spatial profile: bubbles are essentially annulus-shaped domains of uniform polarization separated by ring-shaped domain walls [30][31][32][33] . We include this spatial feature of a bubble directly in our exchange field model Eq. (1) by tuning the function f(r), and we show in Supplementary Note 7 and Supplementary Figs. 9 and 10 that our main results, i.e., the gapped spectrum with a CMEM and wavefunction localization, appear in the bubble model too. This indicates that a wider spectrum of materials and textures could be experimentally explored for realization of our predictions.
The realization of topological superconductivity and the edge states in our system puts constraints on the parameter values. We consider three requirements for successful realization: (i) A substantial effective p-wave gap in the m J = 0 sector, e.g., Δ eff /Δ 0 > 5%. An estimate of the effective gap 47 based on the skyrmioninduced spin-orbit coupling and chemical potential μ m J À Á ¼ μ À m 2 J þ q 2 4 =R 2 sk % μ when R sk large, is given by where λ = R sk /p is the spin-flip length. In this formula, all energy scales are in units of the hopping energy t which determines the bandwidth and we take it to be t~1 eV; the unit lengthscale in the formula is the lattice spacing a whose dimensionful value should correspond to the microscopic electron lengthscale, so we take ã 0.1 nm. The requirement (i) now says that the exchange strength cannot be too large, i.e., J/t < 600 a/R sk , assuming p ≲ 10. Since in materials generically J~1 -10 meV, the skyrmion size is allowed to reach micrometers. The second requirement is that: (ii) The topological regime is reached, so that the exchange scale J surpasses the superconducting pairing Δ 0 . This means Δ 0 is below the 1-10 meV range, or the coherence length is in the range 10-100 nm, which is generally realistic. The final requirement is that: (iii) The CMEM is localized at the edge, i.e., the localization length of the edge-state wavefunctions (estimated to be a ⋅ t/Δ eff ) has to be an order of magnitude smaller than the skyrmion radius R sk . From Eq. (9) using J ≈ Δ 0 we get the constraint that radial winding p~10, consistent with Yang et al. 26 . One may try to relax this constraint by increasing the exchange strength.
For the superconducting part of our setup, we propose aluminum which is a known superconductor and has negligible atomic spin-orbit coupling, in accord with our general assumption that Rashba and intrinsic spin-orbit couplings are sufficiently weak. First of all, disordered thin films of aluminum have a critical temperature of the order of 3 K with a coherence length of around 50 nm 54,55 or less, which is within our theoretically relevant range. Second, direct measurements of the Rashba spinorbit coupling in such thin films are hard to come by, but we find an estimate of E so /Δ 0 = 5% for the ratio of energy scale E so of spin-orbit scattering to the energy scale Δ 0 of s-wave pairing in thin-film aluminum 56 . The skyrmion-induced spin-orbit energy scale E so = (π 2 /2m) ⋅ (p 2 /R 2 ) [cf. derivation of the wire model, Eq. (3)] can without problem reach 5% Δ 0 or several times higher for theoretically relevant values of the parameters p ≲ 10, R sk~1 0-100 nm and Δ 0~1 meV.
For the experimental realization of our findings we propose that the magnetic material be insulating so as to protect the CMEM. From the materials perspective, there are currently two known insulators hosting skyrmions, Cu 2 OSeO 3 and BaFe 12−x −0.05 Sc x Mg 0.05 O 19 (x = 1.6) [5][6][7] . In terms of their parameters, 100 nm-thick Cu 2 OSeO 3 films host skyrmions of radius 25 nm at temperatures ranging from a few Kelvins up to 57 K 57 . There is a sizeable electronic gap of 2.5 eV at 15 K 58 , while the lattice constant is 8.925 Å 59 . All these parameters are within the ranges for which our results are relevant, as detailed in the previous paragraphs. We note that in BaFe 12−x−0.05 Sc x Mg 0.05 O 19 (x = 1.6) the skyrmions are larger, but could be within the upper limit of the tens-ofnanometers range we consider for this parameter. If these magnetic insulators could be grown on a metallic substrate, then one may consider a finite superconducting island deposited on top of the system making the system suitable for Scanning Tunneling Microscopy/Spectroscopy (STM/STS) experiments. In addition, our model might also apply to the case of a metallic magnet, although feedback effects between the texture and the electrons (not considered here) can be important 60,61 . In that regard, we note that skyrmions displaying a three-ring structure where observed experimentally, albeit with a change in the helicity 6 . Further, magnetic skyrmions with q = 2 have also been predicted in frustrated 62 and itinerant 63 magnets. An alternative platform to consider would be thick permalloy (Ni 81 Fe 19 ) disks 64 , since the existence of skyrmions with p up to 3 was recently shown in them, although this would require a different setup. High-p skyrmions were also recently observed in Pd/Fe/Ir(111) magnetic islands 65 . These systems, albeit metallic, naturally provide an edge to localize the CMEM and remove the need to grow a superconducting island. These results are important developments since the larger p also ensures the localization of the CMEM.
The biggest challenge in the experimental verification of our findings lies in the choice of the materials. Indeed, both ingredients (skyrmions and superconductivity) are separately wellcontrolled and well-understood, but little is known about their combination. In particular, we expect that the strength of the exchange field will depend on the achieved interfacing between the magnetic and superconducting materials, which is hard to predict. Recent works aiming at engineering topological superconductivity by using magnetic adatoms or external magnetic fields have shown interesting possibilities, which means that bringing together the magnetism/spintronics and topological superconductivity communities holds great promises.

Methods
Radial tight-binding model. The rotational symmetry of the problem can be exploited by defining the total angular momentum operator J z around the z axis perpendicular to the plane of motion of the electrons. In polar coordinates (r, θ), it is defined as J z ¼ L z þ q 2 σ z where L z = −i∂ θ is the orbital angular momentum. Denoting the eigenvalues of J z as m J , we can expand the electronic field operators as The Nambu spinor Ψ(r) can thus be expanded as whereΨ m J r ð Þ ¼ψ m J ;" r ð Þ;ψ m J ;# r ð Þ;ψ y Àm J ;# r ð Þ; Àψ y Àm J ;" r ð Þ T . We conveniently rescale the spinor by ffiffi r p so that the rdrdθ measure simplifies to drdθ. After all these transformations, the BdG Hamiltonian is block-diagonal in angular momentum space and a single blockĤ m J r ð Þ readŝ We discretize the remaining polar r variable by introducing a lattice spacing a so that r → r j = ja, and in numerical calculations we set a = 1. The nearest-neighbor tightbinding Hamiltonian uses the Nambu basis C y j ¼ c y " ja ð Þ; c y # ja ð Þ; c # ja ð Þ; Àc " ja ð Þ .
We parametrize the tight-binding Hamiltonian aŝ Now, we Taylor expand Eq. (13) to second order, integrate by parts and identify the matrices M and C from Eq. (12). This leads tô ÀtC y jþ1 τ z C j þ h:c: þ C y j 2t À μ À t 4j 2 The tight-binding hopping energy is t = 1/(2ma 2 ) in terms of the effective electron mass m. We exactly diagonalize the Hamiltonian in the form of Eq. (14) without implementing the self-consistency inherent to the Bogoliubov-de-Gennes formalism.
2D tight-binding Hamiltonian. On the square lattice r = (xa, ya) where a ≡ 1 is the lattice spacing, and x, y are integers labeling the sites of the lattice, the twodimensional tight-binding Hamiltonian is H 2D TB ¼ P r¼x;y P σ¼";# Àtc y rþxσ c rσ À tc y rþŷσ c rσ þ ð4t À μÞc y rσ c rσ " þ Δ 0 c y r" c y r# þ h:c: where the parameters are the same as in the main text, and t is the hopping amplitude, μ the chemical potential measured from the bottom of the band, Δ 0 the swave gap and J the exchange coupling with the texture. The unit vector in the x (resp. y) direction is denoted asx (resp.ŷ). Exact diagonalization is then performed without implementing the self-consistency inherent to the Bogoliubov-de-Gennes formalism.
Consistently with the radial model, in the regime of J/Δ 0 large enough (estimated as J > ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Δ 2 0 þ μ 2 q ) and for any q we find weakly-dispersing states extended around the edge and localized in the radial direction near R sk , while only for q even there is a zero energy state localized at the skyrmion center.
Gradient and Laplace operators in the cone geometry. As in Wu et al. 48 , consider a cone of half-opening angle φ and base radius R sk where the coordinates r and θ, respectively denote the distance measured from the tip of the cone and the usual polar angle. Denoting byr andθ the unit vectors on the cone, the gradient and Laplace operators read The cylinder limit is φ = 0, r → ∞ while keeping r sin φ = const = R sk . Under the transformation from the disk to the cylinder via the cone, the surface element varies like r dr dθ ! r sin φ dr dθ ! R sk dr dθ: ð17Þ

Data availability
Codes and datasets used in this study are available from the corresponding author upon reasonable request.