Topological superconductivity with deformable magnetic skyrmions

TL;DR

This paper demonstrates that magnetic skyrmions can induce topological superconductivity and host Majorana modes without extra ingredients, with a robust, scalable chiral Majorana band along skyrmion edges.

Contribution

It reveals a new mechanism for topological superconductivity using deformable magnetic skyrmions proximitized with conventional superconductors.

Findings

Majorana zero mode in skyrmion core

Chiral Majorana band on skyrmion edge

Majorana modes scale with skyrmion perimeter

Abstract

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.

Equations

$$\mathbf{n}\left(\mathbf{r}\right)=\left(\sin f(r)\cos (q\theta ),\sin f(r)\sin (q\theta ),\cos f(r)\right),$$

$$\mathcal{H}\left(\mathbf{r}\right)=\left(-\frac{{\nabla }^2}{2m}-\mu \right){\tau }_z+J\mathbf{\sigma }\cdot \mathbf{n}\left(\mathbf{r}\right)+{\Delta }_0{\tau }_x,$$

$${\mathcal{H}}_{m_{\mathrm{J}}}^{\mathrm{wire}}(r)=\left[-\frac{1}{2m}{\partial }_r^2-\mu \right]{\tau }_z+\frac{1}{2m{R}_{\mathrm{sk}}^2}\left(m_{\mathrm{J}}^2+\frac{q^2}{4}\right){\tau }_z+\frac{f^{\prime }}{2m}{\partial }_{r}i{\sigma }_y{\tau }_z+J{\sigma }_z+{\Delta }_0{\tau }_x$$

$${\mathcal{H}}_{m_{\mathrm{J}}}^{\mathrm{slope}}(r)=-\frac{q{m}_{\mathrm{J}}}{2m{R}_{\mathrm{sk}}^2}(-1{)}^p{\sigma }_z{\tau }_z$$

$${\mathcal{H}}_{m_{\mathrm{J}}}^{\prime }\left(r\right)=\frac{f^{\prime 2}}{8m}{\tau }_z+\frac{f^{\prime \prime }}{4m}i{\sigma }_y{\tau }_z$$

$$\left|m_{\mathrm{J}}^{\ast }\right|=R_{\mathrm{sk}}\sqrt{\mu +\sqrt{J^2-{\Delta }_0^2}},$$

$$\Xi ={\tau }_y{\sigma }_y$$

$${\tilde{\mathcal{H}}}_{m_{\mathrm{J}}}^{\mathrm{cyl},\mathrm{eff}}(r)={\mathcal{H}}_{m_{\mathrm{J}}}^{\mathrm{wire}}(r)+{\mathcal{H}}_{m_{\mathrm{J}}}^{\mathrm{slope}}(r).$$

$${\Delta }_{\mathrm{eff}}=\frac{\pi }{\lambda }\frac{{\Delta }_0}{J}\sqrt{J+\mu }$$

$${\psi }_{\sigma }\left(\mathbf{r}\right)=\sum _{m_{\mathrm{J}}=-\infty }^{+\infty }{e}^{i\left[m_{\mathrm{J}}-\frac{q}{2}{\left({\sigma }_z\right)}_{\sigma \sigma }\right]\theta }{\tilde{\psi }}_{m_{\mathrm{J}},\sigma }\left(r\right)$$

$$\Psi \left(\mathbf{r}\right)=\sum _{m_{\mathrm{J}}=-\infty }^{+\infty }{e}^{i\left[m_{\mathrm{J}}-\frac{q}{2}{\sigma }_z\right]\theta }{\tilde{\Psi }}_{m_{\mathrm{J}}}\left(r\right)$$

$${\hat{\mathcal{H}}}_{m_{\mathrm{J}}}\left(r\right)=-\frac{1}{2m}\left[{\partial }_r^2+\frac{1}{4r^2}\left(1-q^2-4m_{\mathrm{J}}^2+4q{m}_{\mathrm{J}}{\sigma }_z\right)-\mu \right]{\tau }_z+J{\sigma }_z\cos f+J{\sigma }_x\sin f+{\Delta }_0{\tau }_x$$

$$\begin{array}{rl}{\hat{H}}_{m_{\mathrm{J}}}^{\mathrm{TB}}=\sum _{j=1}^L{\mathcal{C}}_{j+1}^{†}M{\mathcal{C}}_j& +h.c.+{\mathcal{C}}_j^{†}C{\mathcal{C}}_j\end{array}$$

$$\begin{array}{rl}{\hat{H}}_{m_{\mathrm{J}}}^{\mathrm{TB}}\approx \sum _{j=1}^L-t{\mathcal{C}}_{j+1}^{†}{\tau }_z{\mathcal{C}}_j& +h.c.+{\mathcal{C}}_j^{†}\left[2t-\mu -\frac{t}{4j^2}\left(1-q^2-4m_{\mathrm{J}}^2+4q{m}_{\mathrm{J}}{\sigma }_z\right)\right]{\tau }_z{\mathcal{C}}_j\\ & +{\mathcal{C}}_j^{†}\left[J{\sigma }_z\cos f+J{\sigma }_x\sin f+{\Delta }_0{\tau }_x\right]{\mathcal{C}}_j.\end{array}$$

$$\begin{array}{rl}{H}^{2D\mathrm{TB}}& =\sum _{\mathbf{r}=x,y}[\sum _{\sigma =\uparrow ,\downarrow }-t{c}_{\mathbf{r}+\hat{\mathbf{x}}\sigma }^{†}{c}_{\mathbf{r}\sigma }-t{c}_{\mathbf{r}+\hat{\mathbf{y}}\sigma }^{†}{c}_{\mathbf{r}\sigma }+(4t-\mu )c_{\mathbf{r}\sigma }^{†}{c}_{\mathbf{r}\sigma }\\ & +{\Delta }_0{c}_{\mathbf{r}\uparrow }^{†}{c}_{\mathbf{r}\downarrow }^{†}+h.c.\\ & +J\sum _{\sigma ,{\sigma }^{\prime }}{c}_{\mathbf{r}\sigma }^{†}{\left(\mathbf{n}\left(\mathbf{r}\right)\cdot \mathbf{\sigma }\right)}_{\sigma {\sigma }^{\prime }}{c}_{\mathbf{r}{\sigma }^{\prime }}],\end{array}$$

$$\nabla ={\partial }_r\hat{\mathbf{r}}+\frac{1}{r\sin \phi }{\partial }_{\theta }\hat{\mathbf{\theta }}$$

$$rdrd\theta \to r\sin \phi drd\theta \to R_{\mathrm{sk}}drd\theta .$$

$$\left|m_{J,\pm }^{\ast }\right|\approx R_{\mathrm{sk}}\sqrt{\left(\mu \pm \sqrt{J^2-{\Delta }_0^2}\right)}+\mathcal{O}\left(1\right).$$

$${\epsilon }^{\ast }=c\frac{q{m}_{\mathrm{J}}^{\ast }}{2m{R}_{\mathrm{sk}}^2}$$

$$\frac{{\epsilon }^{\ast }}{{\Delta }_{\mathrm{eff}}}=c\frac{q}{p}\frac{J}{\pi {\Delta }_0}\sqrt{\frac{\mu +\sqrt{J^2-{\Delta }_0^2}}{J+\mu }},$$

$$\frac{P\left(\frac{r}{R}\right)}{2\pi R}=1+\frac{r}{R}-\frac{1}{\pi }\arccos \left(1-\frac{r^2}{2R^2}\right)-\frac{1}{\pi }\frac{r}{R}\arcsin \left(\sqrt{1-\frac{r^2}{4R^2}}\right).$$

$$P(a,b)=4aE\left(\sqrt{1-\frac{b^2}{a^2}}\right),$$

$$E\left(m\right)={\int }_0^{\frac{\pi }{2}}\sqrt{1-m^2{\sin }^2\theta }d\theta ={\int }_0^1\frac{\sqrt{1-m^2{t}^2}}{\sqrt{1-t^2}}dt.$$

$$S_{\mathrm{peak}}\approx \frac{1}{W_{\mathrm{peak}}}\sum _{j\in \mathrm{peak}}{\left|u_{\uparrow }(j)\right|}^2$$

$$f\left(r\right)=\pi \sum _{i}g\left(r-r_0-r_i,w\right)$$

Full text

I 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 TokNag2013 , in ultrathin magnetic films Heinze2011 ; Romming2013 ; Gross2018 as well as in multiferroic insulators Seki2012 ; Tokura2012 ; White2012 . Quite remarkably, magnetic skyrmions can be stabilized over a wide temperature domain ranging from room temperature Yu2011 ; Parkin2017 to cryogenic temperature Heinze2011 ; Romming2013 ; Romming2015 . Evidence that magnetic skyrmions can be driven by ultralow electric current densities Jonietz2010 ; Yu2012 make them promising candidates for future spintronic applications Fert2013 .

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 Nakosai2013 ; NP2013 ; Braunecker2013 ; Klinovaja2013 ; Pientka2013 and experimentally indicating Majorana zero-energy modes at endpoints of one-dimensional systems, such as iron atomic chains Nadj-Perge_Yazdani_adatom_chain_2014 ; Pawlak2016 ; Wiesendanger_Chain_2018 and semiconducting wiresLutchyn:2018 . 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 leadMenardPbCoSi or iron adatom clusters on a rhenium surfacePalacioMoralesWiesendanger_Island_2018 . 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 superconductorNakosai2013 ; Schnyder2015 ; Loss_Majorana_skyrmion ; Zutic_PRL ; Mohanta2019 ; Kontos2019 . In fact, our results are relevant for a broader class of skyrmion-like textures, such as magnetic bubblesKiselev:2011 ; Tchernyshyov:2012 ; Buttner:2018 ; BernandMantel:2018 . 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 Loss_Majorana_skyrmion . Moreover, very elongated magnetic skyrmions were shown to host Majorana zero modes at their endpointsKovalevPRB2018 . 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 systemWongLee12 ; HaoTing16 ; DaidoYanase17 . Furthermore, the existence of a MFB along an edge can be deduced from an appropriate discrete chiral symmetry and topological indices in lower spatial dimensionVolovik11 ; Sato11 ; WangLee12 . 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.

II Results

II.1 Setup and model

Consider a two-dimensional (2D) magnetic thin film hosting a skyrmion, which is represented by a classical magnetization texture

[TABLE]

written in polar coordinates ${\bf r}=(r,\theta)$, 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 $\Psi^{\dagger}\left({\bf r}\right)=\left(\psi_{\uparrow}^{\dagger}\left({\bf r}\right),\psi_{\downarrow}^{\dagger}\left({\bf r}\right),\psi_{\downarrow}\left({\bf r}\right),-\psi_{\uparrow}\left({\bf r}\right)\right)$, where $\psi_{\sigma}\left({\bf r}\right)$ annihilates an electron with spin $\sigma$ at position ${\bf r}$ in 2D, the total Hamiltonian $H$ can be written as $H=\frac{1}{2}\int d{\bf r}\,\Psi^{\dagger}\left({\bf r}\right){\cal H}\left({\bf r}\right)\Psi\left({\bf r}\right)$, with the Bogoliubov-de-Gennes (BdG) Hamiltonian ${\cal H}\left({\bf r}\right)$ defined as

[TABLE]

where the $s$-wave superconducting order parameter $\Delta_{0}$ is taken real without loss of generality, the electron effective mass is $m$ and the chemical potential is $\mu$. We set $\hbar=1$ unless explicitly written otherwise. The $\sigma_{\alpha}$ and $\tau_{\alpha}$ ($\alpha=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 aluminium. 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_{\rm so}=1/(m\alpha)$, with $\alpha$ 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 Rudner2016 ; Buzdin2018 ; Takashima2016 ; Eremin2018 ; Rex_SK_Vortex_2019 .

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_{\rm sk}$, which determines its size. We consider a hard-wall boundary condition at the edge i.e. at $r=R_{\rm 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.

II.2 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}+\frac{q}{2}\sigma_{z}$ where $L_{z}=-i\partial_{\theta}$ is the orbital angular momentum. Using the eigenvalues of $J_{z}$, denoted as $m_{\rm 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 Radial tight-binding model). In the regime of $J/\Delta_{0}$ large enough (estimated as $J>\sqrt{\Delta_{0}^{2}+\mu^{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. Such a spectral feature was observed before in related models Nakosai2013 ; Schnyder2015 ; Loss_Majorana_skyrmion ; Morr_engineering . Further, the LDoS clearly displays a reduced gap $\Delta_{\rm eff}\approx 5\%\,\Delta_{0}$ consistent with an effective (topological) $p$-wave superconducting gap(AliceaPRB, ). As $J/\Delta_{0}$ is reduced the effective gap closes (at $J=\sqrt{\Delta_{0}^{2}+\mu^{2}}$, see estimate below), and a full gap $\Delta_{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 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/\Delta_{0}$ small enough).

We further clarify the edge states and topological superconductivity by looking at the spectrum $\varepsilon(m_{\rm J})$, in which edge states form a seemingly flat band in a range of $m_{\rm 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 enoughLoss_Majorana_skyrmion , 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 evenLoss_Majorana_skyrmion . This is easily understood since the zero mode must appear in the self-conjugate angular momentum $m_{\rm J}=0$ sector, while $m_{\rm 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.

II.3 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.MappingMartin 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 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_{\rm J}$ basis (see Methods Radial tight-binding model), then we apply the unitary transformation $U(r)=\exp(i\sigma_{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 Hamiltonian $\widetilde{{\cal H}}_{m_{\rm J}}^{\rm cyl}(r)$ can be written as the sum of three parts, $\widetilde{{\cal H}}_{m_{\rm J}}^{\rm cyl}(r)={\cal H}_{m_{\rm J}}^{\rm wire}(r)+{\cal H}_{m_{\rm J}}^{\rm slope}(r)+{\cal H}_{m_{\rm J}}^{\prime}(r)$, where:

[TABLE]

For our purpose it is sufficient to show that the edge modes and the effective gap (the energy gap in the $m_{\rm 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_{\rm J}$, the ${\cal H}_{m_{\rm J}}^{\rm wire}(r)$ can be interpreted as the extensively studied Hamiltonian of a Rashba wire Lutchyn2010 ; Oreg2010 upon introducing a momentum-dependent chemical potential $\mu\left(m_{\rm J}\right)=\mu-\left(m_{\rm J}^{2}+\frac{q^{2}}{4}\right)/\left(2mR_{\rm sk}^{2}\right)$. (Note however that the skyrmion-induced effective spin-orbit coupling in 2D is not of a simple Rashba type.) At each $m_{\rm J}$ the superconducting wire Hamiltonian ${\cal H}_{m_{\rm J}}^{\rm wire}(r)$ is well known to be in a trivial state ($J<\sqrt{\Delta_{0}^{2}+\mu(m_{\rm J})^{2}}$) or in a topological stateLutchyn2010 ; Oreg2010 ($J>\sqrt{\Delta_{0}^{2}+\mu(m_{\rm 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 $\mu(m_{\rm J})$, there is generically a flat zero-energy band of edge modes, i.e. a MFB, for a range of $|m_{\rm J}|<\left|m_{\rm J}^{*}\right|$, where

[TABLE]

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 Figure 1 for details). For precisely $|m_{\rm J}|=\left|m_{\rm J}^{*}\right|$ 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 symmetriesSedlmayerBena11 ; Sato11 ; WangLee12 . 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 ${\cal H}_{m_{\rm J}}^{\rm slope}(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

[TABLE]

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 statesSedlmayerBena11 , namely, a perturbation must have energy larger than the effective gap; or, it should hybridize the MFB with low energy bulk states at $|m_{\rm J}|\approx\left|m_{\rm J}^{*}\right|$, 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_{\rm 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_{\rm 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\neq 0$:

[TABLE]

The single term ${\cal H}_{m_{\rm J}}^{\rm slope}(r)$ breaks the chiral symmetry $\Xi$, and there are no other chiral operators. The term ${\cal H}_{m_{\rm J}}^{\rm slope}(r)$ exactly contributes an energy $\varepsilon^{\rm edgestate}(m_{\rm J})\sim m_{\rm J}$ to an MFB state at $m_{\rm J}$, making the flat MFB into a linearly dispersing chiral Majorana edge mode (CMEM) of the $q\neq 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 $\varepsilon^{*}$ that the CMEM can have, which occurs at the maximal $m_{\rm J}$ of the CMEM, i.e. , $\varepsilon^{*}\equiv|\varepsilon^{\rm edgestate}(\left|m_{\rm J}^{*}\right|)|$. Treating ${\cal H}_{m_{\rm J}}^{\rm slope}(r)$ as a first order perturbation to the MFB (see Supplementary Note 2 and Supplementary Figure 2), the estimate $\varepsilon^{*}=\frac{q}{R_{\rm sk}}\sqrt{\mu+\sqrt{J^{2}-\Delta_{0}^{2}}}$ scales the same way with skyrmion size as the estimate of the effective gap $\Delta_{\rm eff}\sim p/R_{\rm sk}$. For the relevant regime of $J,\,\Delta_{0},\,\mu$ (see Discussion) the quantitative ratio is at most $\varepsilon^{*}/\Delta_{\rm eff}\sim 0.1$. The corresponding Fermi velocity of the CMEM is therefore small and suppressed by the skyrmion size, $\partial\varepsilon^{\rm edgestate}(m_{\rm J})/\partial(m_{\rm J}/R_{\rm sk})\sim\frac{q}{2mR_{\rm 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 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 Figure 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_{\rm J}^{*}$ in accordance to predictions in Eq. (6), and we observe the near-flatness of the edge mode.

II.4 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_{\rm 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 $\mu$ depending on $m_{\rm J}^{*}$). We remind that the edge of system is defined by setting 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_{\rm sk}$, gives consistent results with the radial model (see Methods 2D tight-binding Hamiltonian, Supplementary Note 4 Counting the number of edge states of circular skyrmion and Supplementary Figure 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 Figure 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 $\xi$ associated to the edge state. We find that $\xi\approx 0.5\lambda$ for the parameters considered, where the lengthscale $\lambda$ measures the distance for a single radial spin flip, i.e. , $R_{\rm sk}=p\lambda$. 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 $\xi$, as we additionally confirm through an investigation of elliptical skyrmion geometries (see Supplementary Note 4 Elliptic geometry and Supplementary Figure 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 Figure 7).

III 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 Figure 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 Loss_Majorana_skyrmion .

Our analysis is carried out in the Bogoliubov-de Gennes formalism without self-consistency. We believe that a self-consistent calculation would not change our main conclusions since self-consistent calculationsBjornson_selfconst_2015 ; Christen2016 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 example slightly shift the value of $m_{\rm 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”TokNag2013 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 topologyKiselev:2011 ; Tchernyshyov:2012 ; Buttner:2018 ; BernandMantel:2018 . 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 wallsKiselev:2011 ; Tchernyshyov:2012 ; Buttner:2018 ; BernandMantel:2018 . 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 Figure 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_{\rm J}=0$ sector, e.g. $\Delta_{\rm eff}/\Delta_{0}>5\%$. An estimate of the effective gapAliceaPRB based on the skyrmion-induced spin-orbit coupling and chemical potential $\mu\left(m_{\rm J}\right)=\mu-\left(m_{\rm J}^{2}+\frac{q^{2}}{4}\right)/R_{\rm sk}^{2}\approx\mu$ when $R_{\rm sk}$ large, is given by

[TABLE]

where $\lambda=R_{\rm sk}/p$ is the spin-flip length. In this formula, all energy scales are in units of hopping energy $t$ which determines the bandwidth and we take it to be $t\sim 1\,{\rm 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 $a\sim 0.1\,{\rm nm}$. The requirement (i) now says that the exchange strength cannot be too large, i.e. $J/t<600\,a/R_{\rm sk}$, assuming $p\lesssim 10$. Since in materials generically $J\sim 1-10\,{\rm 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 $\Delta_{0}$. This means $\Delta_{0}$ is below the $1-10\,{\rm meV}$ range, or the coherence length is in the range $10-100\,{\rm 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\cdot t/\Delta_{\rm eff}$) has to be an order of magnitude smaller than the skyrmion radius $R_{\rm sk}$. From Eq. (9) using $J\approx\Delta_{0}$ we get the constraint that radial winding $p\sim 10$, consistent with Yang et al Loss_Majorana_skyrmion . 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 aluminium have a critical temperature of the order of 3K with a coherence length of around 50 nm(Meservey:1971, ; Dressel:2008, ) or less, which is within our theoretically relevant range. Second, direct measurements of the Rashba spin-orbit coupling in such thin films are hard to come by, but we find an estimate of $E_{\rm so}/\Delta_{0}=5\%$ for the ratio of energy scale $E_{\rm so}$ of spin-orbit scattering to the energy scale $\Delta_{0}$ of $s$-wave pairing in thin-film aluminumMeservey:1994 . The skyrmion-induced spin-orbit energy scale $E_{\rm so}=(\pi^{2}/2m)\cdot(p^{2}/R^{2})$ [cf.  derivation of the wire model, Eq. (3)] can without problem reach $5\%\Delta_{0}$ or several times higher for theoretically relevant values of the parameters $p\lesssim 10$, $R_{\rm sk}\sim 10-100\,{\rm nm}$ and $\Delta_{0}\sim 1\,{\rm 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, ${\rm Cu}_{2}{\rm O}{\rm Se}{\rm O}_{3}$ and ${\rm Ba}{\rm Fe}_{12-x-0.05}{\rm Sc}_{x}{\rm Mg}_{0.05}{\rm O}_{19}$ ($x=1.6$) Seki2012 ; Tokura2012 ; White2012 . In terms of their parameters, 100 nm-thick ${\rm Cu}_{2}{\rm O}{\rm Se}{\rm O}_{3}$ films host skyrmions of radius 25 nm at temperatures ranging from a few Kelvins up to 57 KEverschor-Sitte:2018 . There is a sizeable electronic gap of 2.5 eV at 15 KVersteeg:2016 , while the lattice constant is 8.925 $\AA$ Portnichenko:2016 . All these parameters are within the ranges for which our results are relevant, as detailed in the previous paragraphs. We note that in ${\rm Ba}{\rm Fe}_{12-x-0.05}{\rm Sc}_{x}{\rm Mg}_{0.05}{\rm O}_{19}$ ($x=1.6$) the skyrmions are larger, but could be within the upper limit of the tens-of-nanometers 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. Additionally, 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 Motome2018 ; NogueraPRB2018 . In that regard, we note that skyrmions displaying a three-ring structure where observed experimentally, albeit with a change in the helicity Tokura2012 . Further, magnetic skyrmions with $q=2$ have also been predicted in frustratedMostovoyNC2015 and itinerantMotomeLargeQ2017 magnets. An alternative platform to consider would be thick permalloy (${\rm Ni}_{81}{\rm Fe}_{19}$) disksTargetSK-permalloy2018 , 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 islandsWiesendanger_target_exp . 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 well-controlled 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.

Acknowledgement:

We would like to acknowledge useful conversations with Marco Aprili, Freek Massee, Stanislas Rohart, Nicholas Sedlmayr and Silas Hoffman. This work has been partially supported by French Agence Nationale de la Recherche through the contract ANR Mistral.

Author contributions:

M. G. performed the numerical and analytical calculations under the supervision of A. M. and P. S. All authors discussed the results and contributed to the final manuscript.

Data availability:

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

Competing interests:

The authors declare that they have no competing interests.

Methods

I 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 $\left(r,\theta\right)$, it is defined as $J_{z}=L_{z}+\frac{q}{2}\sigma_{z}$ where $L_{z}=-i\partial_{\theta}$ is the orbital angular momentum. Denoting the eigenvalues of $J_{z}$ as $m_{\rm J}$, we can expand the electronic field operators as

[TABLE]

The Nambu spinor $\Psi\left({\bf r}\right)$ can thus be expanded as

[TABLE]

where $\widetilde{\Psi}_{m_{\rm J}}\left(r\right)=\left(\widetilde{\psi}_{m_{\rm J},\uparrow}\left(r\right),\widetilde{\psi}_{m_{\rm J},\downarrow}\left(r\right),\widetilde{\psi}^{\dagger}_{-m_{\rm J},\downarrow}\left(r\right),-\widetilde{\psi}^{\dagger}_{-m_{\rm J},\uparrow}\left(r\right)\right)^{T}$. We conveniently rescale the spinor by $\sqrt{r}$ so that the $rdrd\theta$ measure simplifies to $drd\theta$. After all these transformations, the BdG Hamiltonian is block-diagonal in angular momentum space and a single block $\widehat{{\cal H}}_{m_{\rm J}}\left(r\right)$ reads

[TABLE]

We discretize the remaining polar $r$ variable by introducing a lattice spacing $a$ so that $r\to r_{j}=ja$, and in numerical calculations we set $a=1$. The nearest-neighbor tight-binding Hamiltonian uses the Nambu basis ${\cal C}^{\dagger}_{j}=\left(c^{\dagger}_{\uparrow}\left(ja\right),c^{\dagger}_{\downarrow}\left(ja\right),c_{\downarrow}\left(ja\right),-c_{\uparrow}\left(ja\right)\right)$. We parametrize the tight-binding Hamiltonian as

[TABLE]

Now, we Taylor expand Eq. (13) to second order, integrate by parts and identify the matrices $M$ and $C$ from Eq. (12). This leads to

[TABLE]

The tight-binding hopping energy $t=1/\left(2ma^{2}\right)$ 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.

II 2D tight-binding Hamiltonian

On the square lattice ${\bf r}=\left(xa,ya\right)$ where $a\equiv 1$ is the lattice spacing, and $x,y$ are integers labelling the sites of the lattice, the two-dimensional tight-binding Hamiltonian is

[TABLE]

where the parameters are the same as in the main text, and $t$ is the hopping amplitude, $\mu$ the chemical potential measured from the bottom of the band, $\Delta_{0}$ the $s$-wave gap and $J$ the exchange coupling with the texture. The unit vector in the $x$ (resp. $y$) direction is denoted as ${\bf\hat{x}}$ (resp. ${\bf\hat{y}}$). 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/\Delta_{0}$ large enough (estimated as $J>\sqrt{\Delta_{0}^{2}+\mu^{2}}$) and for any $q$ we find weakly-dispersing states extended around the edge and localized in the radial direction near $R_{\rm sk}$, while only for $q$ even there is a zero energy state localized at the skyrmion center.

III Gradient and Laplace operators in the cone geometry

As in Wu et alMappingMartin , consider a cone of half-opening angle $\varphi$ and base radius $R_{\rm sk}$ where the coordinates $r$ and $\theta$ respectively denote the distance measured from the tip of the cone and the usual polar angle. Denoting by ${\bf\hat{r}}$ and $\hat{\bm{\theta}}$ the unit vectors on the cone, the gradient and Laplace operators read

[TABLE]

The cylinder limit is $\varphi=0$, $r\to\infty$ while keeping $r\sin\varphi={\rm const}=R_{\rm sk}$. Under the transformation from the disk to the cylinder via the cone, the surface element varies like

[TABLE]

Supplementary Note 1 Topological transitions in the wire model

The wire Hamiltonian ${\cal H}_{m_{\rm J}}^{\rm wire}(r)$ (see Eq. (3) in the main text) corresponding to angular momentum sector $m_{\rm J}$ is in a topological phase when $\mu\left(m_{\rm J}\right)^{2}-\left(\mu^{*}\right)^{2}<0$, where we introduced the quantity $\mu^{*}=\sqrt{J^{2}-\Delta_{0}^{2}}>0$. This leads to two cases. Firstly, Eq. (6) of the main text is only valid in the case $\left|\mu\right|<\mu^{*}$, and gives the two solutions presented there. Secondly, if $\mu>\mu^{*}$, solving for the gap-closing momenta yields four solutions $\pm\left|m_{{\rm J},\pm}^{*}\right|$, of the form

[TABLE]

Supplementary Figure 1 shows the different profiles of $\mu\left(m_{\rm J}\right)^{2}-\left(\mu^{*}\right)^{2}$ depending on the ratio of $\mu$ and $\mu^{*}$, where $\mu\left(m_{\rm J}\right)\equiv\mu-\left(m_{\rm J}^{2}+\frac{q^{2}}{4}\right)/\left(2mR_{\rm sk}^{2}\right)$. Supplementary Figure 1c shows that in the second case discussed above (not presented in the main text), the momentum range $m_{\rm J}\in\left[-m_{{\rm J},-}^{*},m_{{\rm J},-}^{*}\right]$ is topologically trivial while the ranges $m_{\rm J}\in\left[m_{{\rm J},-}^{*},m_{{\rm J},+}^{*}\right]$ and $m_{\rm J}\in\left[-m_{\rm{J},+}^{*},-m_{\rm{J},-}^{*}\right]$ are non-trivial as denoted by the gray filling of the curve. Consequently, in this regime of $\mu>\mu^{*}$ there is no Majorana zero-mode in the core of the skyrmion even if the azimuthal winding number $q$ is even. The case $\mu<-\mu^{*}$ corresponds to the case of a fully topologically trivial system. These considerations were confirmed numerically in the radial tight-binding setup of the skyrmion model.

Supplementary Note 2 Velocity of the chiral Majorana edge mode

Treating ${\cal H}_{m_{\rm J}}^{\rm slope}(r)$ in Eq. (4) in the main text as a first order perturbation to the Majorana flat band (MFB) of the $q=0$ model, we can estimate the upper limit on the energy $\varepsilon^{*}$ that the CMEM reaches, which occurs at the maximal $m_{\rm J}$ of the CMEM, i.e. , $\varepsilon^{*}\equiv|\varepsilon^{\rm edgestate}(\left|m_{\rm J}^{*}*\right|)|$. This reads

[TABLE]

where $m_{\rm J}^{*}$ is given by Eq. (6) in the main text, assuming the case $\left|\mu\right|<\mu^{*}$. We have also assumed $\left\langle\tau_{z}\sigma_{z}r^{-2}\right\rangle_{\rm MFBstate}=cR_{\rm sk}^{-2}$, taking into account that the edge states are very localized around $R_{\rm sk}$ for the relevant range of model parameters, and expecting that $c$ is a constant of order unity.

The effective $p$-wave gap $\Delta_{\rm eff}$ being estimated by Eq. (9) in the main text, we obtain for the ratio between the maximal CMEM energy and the effective gap:

[TABLE]

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 Supplementary Figure 2. First, we have calculated numerically the expectation value $\left\langle\sigma_{z}\tau_{z}\right\rangle_{\rm MFBstate}$ and indeed found typical values $\approx\!0.44$ of order unity. Supplementary Figure 2 shows a good agreement between our numerics and our analytical estimate for $c\equiv 0.44$. Due to the large size of the skyrmion used here to minimize finite-size effects, the ratio Supplementary Equation SM 2.3 is of order $5\%$. From the computation on smaller skyrmions and by varying parameters $J,\,\mu$, we find that this value can be increased to at most $\approx\!10\%$ and the CMEM therefore universally appears almost flat. In fact, our analytical estimate confirms that the near-flatness of the CMEM cannot significantly change by varying the parameters $J,\,\Delta_{0}$ in their respective ranges under consideration (see Discussion).

Supplementary Note 3 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 $\mu/t=0$. To do so, we tune the exchange coupling $J$ and measure the effective gap, i.e. the gap in the $m_{\rm J}=0$ sector. Supplementary Figure 3 clearly shows that the effective gap behaves the same way in both models, including the topological phase transition at the analytically predicted value $J=\sqrt{\Delta_{0}^{2}+\mu^{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) in the main text.

Supplementary Note 4 Different edge geometries: counting the states in the chiral Majorana edge mode

In this section we define the different geometries mentioned in the main text, as well as our technique for counting the edge states.

Counting the number of edge states of a circular skyrmion

We firstly investigate the original skyrmion texture that, dubbed “geometry 1”. The circular edge of the original skyrmion texture (“geometry 1”) is close to perfectly rotationally symmetric on the 2D lattice (except for the breaking of the spatial rotation symmetry down to the square lattice’s discrete one). 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_{\rm sk}$. The position $R_{\rm max}$ of the maximum of the wavefunction is found using the angularly-averaged wavefunction. Precisely, a state is an edge state if $\left|R_{\rm sk}-R_{\rm max}\right|\leq 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 Supplementary Figure 4, where the cut-off $l$ is 3 lattice sites so that $l/\lambda\approx 0.19$, given the skyrmion’s value $\lambda/a=16$.

Defining the geometries and edge state counting

Geometries 2 and 3 are obtained 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 Supplementary Figure 5. The shape of the edge curve is defined as the outside perimeter of two overlapping 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$ (Supplementary Figure 4). The exact perimeter of the resulting texture’s edge (i.e. the outside perimeter of the overlapping disks), $P\left(\frac{r}{R}\right)$, reads

[TABLE]

Alongside geometry 1, defined by a single circular edge of radius $R_{\rm 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\equiv R_{\rm sk}$.

In the case of non-circular edge shape, we straightforwardly generalize the edge state counting . 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_{\rm max}\geq R_{\rm 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 $\Delta_{\rm eff}/4$, recalling that the highest edge state energy is expected to be $\varepsilon^{*}\approx\Delta_{\rm eff}/10$. In all studied cases the energy cutoff is $\approx 0.05\Delta_{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 Supplementary Figure 5 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

[TABLE]

where $E\left(m\right)$ is the complete elliptic integral of the second kind defined as

[TABLE]

The selection criterion for edge states is a simple generalization of the circular case (geometry 1) since the location of points on the edge of the ellipse is simply defined in cartesian coordinates. The perimeter is increased by changing the aspect ratio of the ellipse while keeping its surface area constant, starting from $p=16$, $q=2$ and $\lambda/a=12$. The results and linear fits for different edge state selection criteria $l/a=3,\,4,\,5$ ($l/\lambda=0.25,\,0.33,\,0.42$) are shown in Supplementary Figure 6.

The fit results are:

• $\bullet$

$l/\lambda=0.25$: $0.31(36)\pm 0.05(35)\,a^{-1}$

• $\bullet$

$l/\lambda=0.33$: $0.30(95)\pm 0.13(74)\,a^{-1}$

• $\bullet$

$l/\lambda=0.42$: $0.56(93)\pm 0.21(87)\,a^{-1}$.

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 Note 5 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: $\mu\to\mu+\delta\mu$ where $\delta\mu$ is distributed according to a normal law of mean [math] and standard deviation $\sigma_{\mu}$.

In Supplementary Figure 7 we show the results for varying disorder strengths $\sigma_{\mu}/\Delta_{0}=0.2,\,0.5,\,0.7$ in the case of the $p=9$ and $\lambda/a=16$ skyrmion.

An example of a clean system’s density of states is shown in Supplementary Figure 4. While for the particular clean system whose parameters exactly correspond to the disordered systems here, there are 90 edge states in the energy range $\pm 0.05\,\Delta_{0}$. The density of states of the disordered systems, including the count of edge states in the energy range $\pm 0.05\,\Delta_{0}$, namely, 90, 106, and 98 states

for Supplementary Figure 7a, b and c, respectively, confirm that the edge mode is stable to relatively high disorder strengths.

Supplementary Note 6 Delocalization of the edge states in a magnetic background

Yang et al (Supplementary Reference Loss_Majorana_skyrmion, ) predict 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 Supplementary Figure 8 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 $\left(u_{\uparrow}(r),u_{\downarrow}(r),v_{\downarrow}(r),v_{\uparrow}(r)\right)$, we compute the spectral weight $S_{\rm peak}$ of the state’s peak through the $u_{\uparrow}$ component, as

[TABLE]

where $W_{\rm peak}$ is the estimated width (i.e. length in radial direction) of the peak of the edge state. We plot the results for $S_{\rm peak}$ as a function of the background region size (i.e. length in radial direction) in Supplementary Figure 8.

Supplementary Note 7 Magnetic bubbles

For our model the most important spatial feature of bubbles is that they are essentially annulus-shaped domains of uniform polarization separated by ring-shaped domain walls Supplementary Reference Kiselev:2011, ; Tchernyshyov:2012, ; Buttner:2018, ; BernandMantel:2018, , see Supplementary Figure 9a. We therefore contrast such a straightforward model of a bubble with our model of the skyrmion in Supplementary Figure 9b,c. In Supplementary Figure 9b the skyrmion shows smooth oscillation between out-of-plane (blue line is $z$-component) and in-plane magnetization (orange line). In contrast, the bubble, Supplementary Figure 9c, has magnetization roughly out-of-plane and constant on concentric annuli (blue line is $z$-component), while at the ring-shaped domain walls between these annuli there is in-plane winding of the magnetization (orange line is the in-plane component).

For the sake of completeness, we remind that the general texture profile is $\{n_{x},n_{y},n_{z}\}\equiv\left\{\sin(f(r))\cos(q\theta),\sin(f(r))\sin(q\theta),\cos(f(r))\right\}$, with $r$ the radial coordinate and $\theta$ the polar angle, and the explicit profile for the bubble we choose as:

[TABLE]

where $g\left(x,w\right)=\left(1+\exp\left(-x/w\right)\right)^{-1}$ is an inverted Fermi function ($g(x,w)=n_{F}(-x,w)$), the $r_{0}=4w$ is a global offset, and $r_{i}=i\left(L-pw/2\right)/p$ is the position of the $i^{\rm th}$ domain wall, while $w$ is the radial width of the domain-wall.

We show in Supplementary Figure 10 the result of our radial tight-binding computations with the bubble profile of Supplementary Figure 9c. The momentum-resolved excitation spectrum, Supplementary Figure 10a, shows that the nearly-flat band and the effective gap survive; the wavefunction amplitude of the zero-energy mode (its partner is not shown) at angular momentum zero shows that the Majorana wavefunctions at the center and at the edge remain localized, although their detailed spatial profile has changed significantly.

In conclusion, if the bubble’s uniform magnetization domains are not too large in real space (the radial width of the domains in above calculation is 100 sites, corresponding to a lengthscale of nm to 10 nm), our results are robust. If the domains increase, one expects to reach a regime where the ferromagnetic nature of the domains dominates, and the effective gap closes. In magnetic bubble material, the size $w$ of a domain wall is estimated as $w=\sqrt{A/K}$ where $A$ is the micromagnetic exchange constant and $K$ the (effective) magnetocrystalline anisotropy Supplementary Reference BernandMantel:2018, . With the typical values $A\sim{\rm pJ}.{\rm m}^{-2}$ and $K\sim{\rm MJ}.{\rm m}^{-3}$ Supplementary Reference BernandMantel:2018, and Supplementary Reference Wiesendanger_target_exp, , we get $w\sim 1-10\,{\rm nm}$ which is consistent with our previous estimations for the radial winding number $p\lesssim 10$ and the skyrmion radius $R_{\rm sk}\sim 10-100\,{\rm nm}$.