Topological superconductivity with deformable magnetic skyrmions

Magnetic skyrmions are nanoscale spin configurations that can be efficiently created and manipulated. They hold great promises for next-generation spintronics applications. In parallel to these developments, 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 magnetic 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 regime of parameters, 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 a deformed skyrmion configuration, while being robust to local disorder.


Introduction
Magnetic skyrmions are nano-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-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 spin-orbit 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. 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 22 or iron adatom clusters on a rhenium surface 23 . 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,24,25 .
Additionally, 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 25 . 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 shape of the skyrmion. Such deformations preserve the number of edge states proportional to the perimeter length of the edge.
For systems with translational symmetry in real space there is a theoretical classification of topological superconducting phases, and predictions for a corresponding CMEM along a given edge of the system [26][27][28] . Furthermore, the existence of an MFB along an edge can be deduced from an appropriate discrete chiral symmetry and topological indices in lower spatial dimension [29][30][31] . 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 skyrmion's shape. 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.

I. RESULTS
A. 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)) , 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 Ψ † (r) = ψ † ↑ (r) , ψ † ↓ (r) , ψ ↓ (r) , −ψ ↑ (r) , where ψ σ (r) annihilates an electron with spin σ at position r in 2D, the total Hamiltonian H can be where the s-wave superconducting order parameter ∆ 0 is taken real without loss of generality, the electron effective mass is m, the chemical potential is µ, and we set = 1. 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 spin-orbit length in the superconductor l so = 2 /(mα) with α the Rashba spin-orbit amplitude, is much larger than the typical lengthscale of the texture so that the magnetoelectric coupling can be neglected [32][33][34] .
Note that the definition in Eq. (1) may describe both Bloch and Néel skyrmions by adding an overall constant angular shift, the helicity 1 , which does not alter our conclusions. 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 introduce the useful lengthscale λ over which a single radial spin flip occurs, i.e. , R sk = pλ. For simplicity, the function f (r) is chosen to be linear between the core and the edge of the skyrmion but its exact shape has no influence on our conclusions.

B. 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 real-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 Sec. A 1). Considering a hard-wall boundary condition at the edge of the skyrmion i.e. at r = R sk , and the regime of J/∆ 0 large enough (estimated as J > ∆ 2 0 + µ 2 below), Fig. 1b shows the resulting strong peak in the local density of states (LDoS), near zero energy and at the skyrmion's edge. Further, the LDoS clearly displays a reduced gap ∆ eff ≈ 5% ∆ 0 consistent with an effective (topological) p-wave superconducting gap 35 . As J/∆ 0 is reduced the effective gap closes (at J = ∆ 2 0 + µ 2 , see estimate below), and a full gap ∆ 0 develops without any edge states. This is expected in a transition from a topological to trivial superconductivity. The regime of topological p-wave superconductivity is also consistent with our finding of other in-gap 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-pairing-dominated trivial phase (J/∆ 0 small enough).
We further clarify the edge states and the 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 and 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 25 , we showcase p = 6).
Note that in contrast, we find a single Majorana zero mode at the core of the skyrmion 25 only if the skyrmion's azimuthal winding number q is even. 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 singlevaluedness of the wavefunction. The existence of edge states indicates that skyrmions of any q induce topological superconductivity. 36

C. Topological origin and the near-flatness of edge mode
In order to explain the origin of the edge states we use a procedure introduced in Ref. 37 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 Sec. B) 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, 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 radial distance r, and finally we apply the mapping to cylinder. The resulting Hamiltonian H cyl m J (r) can be written as the sum of three parts, introduces the mapping: ϕ = π/2 realizes the disk geometry while the limit ϕ → 0 with r → ∞ and r sin ϕ = R sk realizes the cylinder geometry, where r is the distance with respect to the tip of the cone. The core is covered by a white disc for clarity. The excitation energy spectrum ε 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 in a trivial state (J < ∆ 2 0 + µ(m J ) 2 ) or in a topological state 38,39 (J > ∆ 2 0 + µ(m J ) 2 ). For each topological wire there is a single Majorana zero mode localized at the end of the wire, i.e. a single zero mode at the edge of the skyrmion. Due to the variation of µ(m J ), there is generically a flat zero-energy band of edge modes, i.e. an MFB, for a range of |m J | < |m * J |, where where all energies are in units of the bandwidth t, all distances are in units of the lattice spacing a (see Supplementary Materials (SM) Sec. A for details). For precisely |m 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 30,31,40 . Note that the wire Hamiltonian and its MFB become a correct model for a skyrmion if we choose q = 0 and thereby nullify the H slope m J (r). Physically, this is a special case where the skyrmion texture becomes coplanar (in x − z plane, see Eq. (1)), and the orthogonal direction provides a chiral operator Ξ = τ y σ y 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 40 , namely, a perturbation must have energy larger than the effective gap; or, it should hybridize the MFB with low energy bulk states at |m J | ≈ |m * J |, which are few; or, chirality symmetry must be broken (out-of-plane exchange field). Consider now deformations of the shape of the skyrmion edge. The proof of existence of the MFB rests on the rotational symmetry of the (q = 0) skyrmion, providing the m J quantum numbers, which would generally mix under such deformations. Yet, the described stability of the MFB means that a geometric perturbation would be inefficient in removing the MFB states.
We can now proceed to the relevant model for a skyrmion with arbitrary q = 0. 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 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. The CMEM therefore behaves like the MFB it derives from, with weak breaking of chiral symmetry. 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. , ε * ≡ |ε edgestate (|m * J |)|. Treating H slope m J (r) as a first order perturbation to the MFB (see SM Sec. B), the estimate ε * = q R sk µ + J 2 − ∆ 2 0 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) 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, ∂ε edgestate (m J )/∂(m J /R sk ) ∼ q 2mR sk . 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

D. Edge states of deformed skyrmions
The number of states in the single CMEM of a perfectly rotationally symmetric skyrmion is given simply by the highest angular momentum which is reached by the edge states, i.e. m * J . By mapping the skyrmion to Rashba wires as introduced above, we have seen that m * J is estimated by the relation J 2 = ∆ 2 0 + µ (m * J ) 2 , yielding the value for m * J in Eq. (6) up to negligible corrections of order 1. The number of states in CMEM therefore scales linearly with the perimeter of the perfect skyrmion. If the nearly-flat CMEM is indeed robust, we hypothesize that geometric deformations of the edge of the skyrmion texture would preserve the scaling of number of states in CMEM with the perimeter of the deformed texture. We deduced the robustness of CMEM based on the mapping to wires, but the mapping itself concerns low energies and small momenta, so our perimeter hypothesis is far from obviously true.
We substantiate the perimeter hypothesis with an extensive analysis of a 2D squarelattice tight-binding discretization of the skyrmion model Eq. (2) (see Methods Sec. A 2).
Setting the skyrmion exchange strength J to zero outside the skyrmion radius R sk reveals that in the regime of J/∆ 0 large enough (estimated above as J > ∆ 2 0 + µ 2 ), we indeed find weakly-dispersing states extended around the skyrmion's edge and localized near R sk in the radial direction. They occur at low energies whatever the value of q (even or odd), while a Majorana zero more occurs at the core for q even (see SM Sec. D 2 and particularly Fig. D.2).
We firstly investigate the original skyrmion texture that is close to perfectly rotationally symmetric (except for the breaking of the spatial rotation symmetry down to the square lattice's discrete one), dubbed "geometry 1". Then we consider other geometries where the edge of texture is far from a circle. In particular, we consider geometries composed of the skyrmion with an edge shaped by two overlapping circles, one of which has a varying radius.
We consider two cases: "geometry 2", where the circles' radii are equal; and "geometry 3", where the smaller radius is half of the other. The textures are practically obtained by defining a skyrmion texture on a very large area, defining the desired geometric shape of the edge, and then setting the exchange strength to zero outside the edge (see SM Sec. D 1). In the original skyrmion case ("geometry 1"), a state is defined as an edge state if the maximum of its wavefunction lies in some fixed range (a few lattice constants) surrounding the edge of the skyrmion. The same definition is appropriately modified for the other geometries we investigate (see SM Sec. D). We count the edge states as the overall size of the given texture is varied, and 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, which 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. The black horizontal line is the average slope estimated to be 0.12(13) a −1 .
To further investigate the robustness of the states forming the single nearly-flat CMEM, we notice that the states in 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 SM Sec. D 3).
Sharp corners in the shape of the edge seem a stronger perturbations than scalar disorder, to which the nearly-flat CMEM should be very robust. We numerically confirm using the 2D tight-binding model that the nearly-flat CMEM is indeed robust to uncorrelated scalar disorder (see SM Sec. E).

II. DISCUSSION
In conclusion, 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 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 /R 2 sk ≈ µ when R sk large, is given by where λ = R sk /p is the spin-flip length. Requiring a substantial gap, e.g. ∆ eff /∆ 0 > 5%, leads to the requirement that the exchange strength cannot be too large, i.e. J/t < 600 a/R sk , with t the bandwidth and a the microscopic electronic lengthscale. In materials generically, J ∼ 10 meV, t ∼ eV, a ∼ 0.1 nm so the skyrmion size is allowed to reach micrometers. Next, for topological superconductivity the exchange needs to surpass the superconducting pairing, which constrains ∆ 0 into the meV range, which is not unrealistic. Finally, to have a CMEM localized at the edge, the localization length of the edge-state wavefunctions (estimated to be 1/∆ eff ) has to be an order of magnitude smaller than the skyrmion radius R sk , which leads to the requirement of the radial winding p ∼ 10, consistent with Ref. 25. One may try to lower this bound by increasing the exchange strength.
As far as the experimental realization of our findings is concerned, we propose that the magnetic material be insulating so as to protect the CMEM. Additionally, our model might also apply to the case of a metallic magnet, although feedback effects between the 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 We conveniently rescale the spinor by √ r 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 H m J (r) reads We discretize the remaining polar r variable by introducing a lattice spacing a so that In what follows, we set a = 1. To do so, we write a nearest-neighbor tight-binding Hamiltonian in the Nambu basis C † j = c † ↑ (ja) , c † ↓ (ja) , c ↓ (ja) , −c ↑ (ja) . We parametrize the tight-binding Hamiltonian as We exactly diagonalize the Hamiltonian in the form of Eq. (A5).

2D tight-binding Hamiltonian
On the square lattice r = (xa, ya) where a ≡ 1 is the lattice spacing, and x, y are integers labelling the sites of the lattice, the two-dimensional tight-binding Hamiltonian is 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 s-wave gap and J the exchange coupling with the texture. The unit vector in the x (resp. y) direction is denoted asx (resp.ŷ).
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θ.  where m * J is given by Eq. (6) in the main text, assuming the case |µ| < µ * . We have also assumed τ z σ z r −2 MFBstate = cR −2 sk , taking into account that the edge states are very localized around R sk for the relevant range of model parameters, and expecting that c is a constant of order unity.
The effective p-wave gap ∆ eff being estimated by Eq. (9) in the main text, we obtain for the ratio between the maximal CMEM energy and the effective gap: where all energies are in units of the bandwidth t. We next compare this estimate to the numerical results in the radial tight-binding skyrmion model, see The vertical dashed gray lines mark the |µ| = µ * points.

Supplementary Materials C: Phase diagrams of the skyrmion model and the cylinder model
To support the use of our mapping from the disk to the cylinder, we here plot the two phase diagrams at a fixed value of the chemical potential µ/t = 0. To do so, we tune the exchange coupling J and measure the effective gap, i.e. the gap in the m J = 0 sector. Fig. C.1 clearly shows that the effective gap behaves the same way in both models, including the topological phase transition at the analytically predicted value J = ∆ 2 0 + µ 2 . Moreover, the agreement in the topological regime is quantitative, while the small discrepancy can be accounted for by the fact that in the cylinder model we neglected the small chemical potential renormalization and the boundary term, Eq. (5). L e 6 t L + b j + j u R c 7 s 6 r W e 2 Y E X e K j T 9 D I 2 s o n M s Z F N K V k k 4 T x j 4 x Z y L I 7 V F z V 4 6 5 W q P H + x b / f 7 t u g + 3 1 Q R w N 4 i 9 v + 5 v v / c 2 b x Q p e 4 Q 1 v 1 z t s 0 n E L Q / q X + I G f + N X 9 E 0 w F c 8 H 8 V W q 3 4 z U v 8 M 8 I X v 4 F 3 / + r v Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " L X D F h Q x k A h E Q C K / Y j A B D D V Q p 9 s Q = " > A A A D X H i c j V J d S 9 x A F D 2 7 s W r X + t V C X w Q J X Q p 9 W h M p 1 E d p + 1 D 6 Z K G r U h V J s q M G J x 9 M J s q y 9 F / 0 t f 1 f f e m z / o u e u Y 5 g K 3 5 M S H L m 3 H v O z L 0 z a a 3 z x k b R 7 0 4 3 m H o y P T P 7 t D f 3 b H 5 h c W n 5 + X Z T t S Z T w 6 z S l d l N k 0 b p v F R D m 1 u t d m u j k i L V a i c 9 / e D i O 2 f K N H l V f r X j W h 0 U y X G Z H + V Z Y k l 9 + 7 y 2 / 1 F p m x x G h 0 v 9 a B D J C G + D 2 I M + / N i q l j u r 2 M c I F T K 0 K K B Q w h J r J G j 4 7 C F G h J r c A S b k D F E u c Y X v 6 F H b M k s x I y F 7 y u 8 x Z 3 u e L T l 3 n o 2 o M 6 6 i + R o q Q 7 y m p m K e I X a r h R J v x d m x d 3 l P x N P t b c x / 6 r 0 K s h Y n Z B / S X W c + V u d q s T j C h t S Q s 6 Z a G F d d 5 l 1 a 6 Y r b e X i j K k u H m p z D I 8 Y N c S b K 6 z 6 H o m m k d t f b R O I X k u l Y N 8 9 8 b o v L e 6 t L + b j + j u R c 7 s 6 r W e 2 Y E X e K j T 9 D I 2 s o n M s Z F N K V k k 4 T x j 4 x Z y L I 7 V F z V 4 6 5 W q P H + x b / f 7 t u g + 3 1 Q R w N 4 i 9 v + 5 v v / c 2 b x Q p e 4 Q 1 v 1 z t s 0 n E L Q / q X + I G f + N X 9 E 0 w F c 8 H 8 V W q 3 4 z U v 8 M 8 I X v 4 F 3 / + r v Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " L X D F h Q x k A h E Q C K / Y j A B D D V Q p 9 s Q = " > A A A D X H i c j V J d S 9 x A F D 2 7 s W r X + t V C X w Q J X Q p 9 W h M p 1 E d p + 1 D 6 Z K G r U h V J s q M G J x 9 M J s q y 9 F / 0 t f 1 f f e m z / o u e u Y 5 g K 3 5 M S H L m 3 H v O z L 0 z a a 3 z x k b R 7 0 4 3 m H o y P T P 7 t D f 3 b H 5 h c W n 5 + X Z T t S Z T w 6 z S l d l N k 0 b p v F R D m 1 u t d m u j k i L V a i c 9 / e D i O 2 f K N H l V f r X j W h 0 U y X G Z H + V Z Y k l 9 + 7 y 2 / 1 F p m x x G h 0 v 9 a B D J C G + D 2 I M + / N i q l j u r 2 M c I F T K 0 K K B Q w h J r J G j 4 7 C F G h J r c A S b k D F E u c Y X v 6 F H b M k s x I y F 7 y u 8 x Z 3 u e L T l 3 n o 2 o M 6 6 i + R o q Q 7 y m p m K e I X a r h R J v x d m x d 3 l P x N P t b c x / 6 r 0 K s h Y n Z B / S X W c + V u d q s T j C h t S Q s 6 Z a G F d d 5 l 1 a 6 Y r b e X i j K k u H m p z D I 8 Y N c S b K 6 z 6 H o m m k d t f b R O I X k u l Y N 8 9 8 b o v L e 6 t L + b j + j u R c 7 s 6 r W e 2 Y E X e K j T 9 D I 2 s o n M s Z F N K V k k 4 T x j 4 x Z y L I 7 V F z V 4 6 5 W q P H + x b / f 7 t u g + 3 1 Q R w N 4 i 9 v + 5 v v / c 2 b x Q p e 4 Q 1 v 1 z t s 0 n E L Q / q X + I G f + N X 9 E 0 w F c 8 H 8 V W q 3 4 z U v 8 M 8 I X v 4 F 3 / + r v Q = = < / l a t e x i t > J/ 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " L X D F h Q x k A h E Q C K / Y j A B D D V Q p 9 s Q = " > A A A D X H i c j V J d S 9 x A F D 2 7 s W r X + t V C X w Q J X Q p 9 W h M p 1 E d p + 1 D 6 Z K G r U h V J s q M G J x 9 M J s q y 9 F / 0 t f 1 f f e m z / o u e u Y 5 g K 3 5 M S H L m 3 H v O z L 0 z a a 3 z x k b R 7 0 4 3 m H o y P T P 7 t D f 3 b H 5 h c W n 5 + X Z T t S Z T w 6 z S l d l N k 0 b p v F R D m 1 u t d m u j k i L V a i c 9 / e D i O 2 f K N H l V f r X j W h 0 U y X G Z H + V Z Y k l 9 + 7 y 2 / 1 F p m x x G h 0 v 9 a B D J C G + D 2 I M + / N i q l j u r 2 M c I F T K 0 K K B Q w h J r J G j 4 7 C F G h J r c A S b k D F E u c Y X v 6 F H b M k s x I y F 7 y u 8 x Z 3 u e L T l 3 n o 2 o M 6 6 i + R o q Q 7 y m p m K e I X a r h R J v x d m x d 3 l P x N P t b c x / 6 r 0 K s h Y n Z B / S X W c + V u d q s T j C h t S Q s 6 Z a G F d d 5 l 1 a 6 Y r b e X i j K k u H m p z D I 8 Y N c S b K 6 z 6 H o m m k d t f b R O I X k u l Y N 8 9 8 b o v L e 6 t L + b j + j u R c 7 s 6 r W e 2 Y E X e K j T 9 D I 2 s o n M s Z F N K V k k 4 T x j 4 x Z y L I 7 V F z V 4 6 5 W q P H + x b / f 7 t u g + 3 1 Q R w N 4 i 9 v + 5 v v / c 2 b x Q p e 4 Q 1 v 1 z t s 0 n E L Q / q X + I G f + N X 9 E 0 w F c 8 H 8 V W q 3 4 z U v 8 M 8 I X v 4 F 3 / + r v Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " L X D F h Q x k A h E Q C K / Y j A B D D V Q p 9 s Q = " > A A A D X H i c j V J d S 9 x A F D 2 7 s W r X + t V C X w Q J X Q p 9 W h M p 1 E d p + 1 D 6 Z K G r U h V J s q M G J x 9 M J s q y 9 F / 0 We start by defining a skyrmion in a large simulation window, and then defining the desired texture's edge by setting the exchange J to zero outside the edge curve, as shown in disks. The second disk, of radius r is positioned so that its center lies on the perimeter of the first disk, which has radius R (Fig. D.1). The exact perimeter of the resulting texture's edge (i.e. the outside perimeter of the overlapping disks), P r R , reads Alongside geometry 1, defined by a single circular edge of radius R sk (i.e. the original skyrmion), the other two geometries presented in the main text are defined by disk radii r = R (geometry 2) and r = R/2 (geometry 3). For comparisons, we set R ≡ R sk .

Counting the number of edge states
As mentioned in the main text, for a circular edge of the texture ("geometry 1") an eigenstate of the 2D tight-binding model is defined as an edge state if the maximum of its wavefunction lies in a certain corona of width 2l around the edge of the skyrmion at r = R sk .
The position R max of the maximum of the wavefunction is found using the angularly-averaged wavefunction. Precisely, a state is an edge state if |R sk − R max | ≤ l. We will comment on the chosen values for l below. Based on this definition, the count of edge states for the circular skyrmion is presented in Fig. D.2, where the cut-off l is 3 lattice sites so that l/λ ≈ 0.19, given the skyrmion's value λ/a = 16. In the case of non-circular edge shape, we straightforwardly generalize the above definition. Note that our geometries 2 and 3 are formed by adding an outward bulge to the circular edge of geometry 1. Therefore for simplicity we define as an edge state any state that satisfies R max ≥ R sk − l. This criterion is less precise and may slightly overestimate the number of edge states. For simplicity we examine only the eigenstates whose energy is below ∆ eff /4, recalling that the highest edge state energy is expected to be ε * ≈ ∆ eff /10. In all studied cases the energy cutoff is ≈ 0.05∆ 0 .
In all geometries presented here (1, 2 and 3), the variation of the perimeter was achieved by changing the overall spatial scale while keeping the number of sites per spin-flip the same.

Elliptic geometry
The elliptic geometry is defined simply by replacing the two-disk construction of Fig. D.1 (see SM Sec. D 1) by a single ellipse at the center. We use the perimeter P (a, b) of an ellipse of semi-major axis a and semi-minor axis b given as where E (m) is the complete elliptic integral of the second kind defined as The selection criterion for edge states is a simple generalization of the circular case (geometry The data is clearly more noisy than in the case of other geometries, and the fit slopes deviate from the results of the main text (even within error bars, see Fig. 3 in the main text). We ascribe the discrepancy to the varying curvature of the texture's edge that causes hybridization between the wavefunctions on the edge. Note that in geometries 2 and 3 the edge has a non-constant curvature only at isolated points (where two circles meet).

Supplementary Materials E: Stability to disorder in the 2D tight-binding model
The robustness of the chiral Majorana edge mode is further confirmed by an analysis of the effect of scalar disorder in the 2D tight-binding model. The disorder is a random, spatially uncorrelated variation of the on-site energy applied throughout the system: µ → µ + δµ where δµ is distributed according to a normal law of mean 0 and standard deviation σ µ .
In Fig. E.1 we show the results for varying disorder strengths σ µ /∆ 0 = 0.2, 0.5, 0.7 in the case of the p = 9 and λ/a = 16 skyrmion. An example of a clean system's density of states is shown in Fig. D.2. While for the particular clean system whose parameters exactly correspond to the disordered systems here, there are 90 edge states in the energy range ±0.05 ∆ 0 . The density of states of the disordered systems, including the count of edge states in the energy range ±0.05 ∆ 0 , namely, 90, 106, and 98 states for Fig. E.1a, b and c, respectively, confirm that the edge mode is stable to relatively high disorder strengths. Ref. 25 predicts that if a constant magnetic background is added outside the skyrmion then the edge states delocalize from the edge into this background. The reason is simply that the superconductor is gapless in the background region. This is readily verified numerically in the radial tight-binding model of the skyrmion, as shown in Fig. F.1 where we plot the spectral weight of the peak of the edge state as the size of the background region is increased.
Defining the wavefunction as (u ↑ (r), u ↓ (r), v ↓ (r), v ↑ (r)), we compute the spectral weight S peak of the state's peak through the u ↑ component, as where W peak is the estimated width (i.e. length in radial direction) of the peak of the edge state. We plot the results for S peak as a function of the background region size (i.e. length in radial direction) in Fig. F