Majorana bound states in magnetic skyrmions imposed onto a superconductor

TL;DR

This paper investigates how magnetic skyrmions in a superconductor can host Majorana bound states, proposing that skyrmion-vortex pairs with winding number one are promising for realizing and controlling MBS.

Contribution

It demonstrates that realistic skyrmions of winding number one, combined with vortices, can support Majorana bound states, advancing potential experimental control of MBS.

Findings

Majorana bound states are supported in skyrmion-vortex pairs with winding number one.

Higher-winding skyrmions are unlikely to form in chiral magnets.

Numerical results confirm the spatial profile of the MBS in the proposed system.

Abstract

We consider a superconducting film exchange-coupled to a close-by chiral magnetic layer and study how magnetic skyrmions can induce the formation of Majorana bound states (MBS) in the superconductor. Inspired by a proposal by Yang et al. [Phys. Rev. B 93, 224505 (2016)], which suggested MBS in skyrmions of even winding number, we explore whether such skyrmions could result from a merger of ordinary skyrmions. We conclude that the formation of higher-winding skyrmions is not realistic in chiral magnets. Subsequently, we present a possibility to obtain MBS from realistic skyrmions of winding number one, if a skyrmion-vortex pair is formed instead of a bare skyrmion. Specifically, we show that MBS are supported in a pair of a circular skyrmion and a vortex which both have a winding number of one. We back up our analytical prediction with results from numerical diagonalization and obtain the spatial profile of the MBS. In light of recent experimental progress on the manipulation of skyrmions, such systems are promising candidates to achieve direct spatial control of MBS.

Equations

$$n=\int \frac{{\text{d}}^{2}r}{4\pi }\mathbf{m}\cdot \left({\partial }_x\mathbf{m}\times {\partial }_y\mathbf{m}\right)$$

$${\mathbf{m}}_{\text{sk}}(r,\phi )=\left(\begin{array}{c}\cos f(\phi )\sin g(r)\\ \sin f(\phi )\sin g(r)\\ \cos g(r)\end{array}\right)$$

$$n={\left[\frac{f(\phi )}{2\pi }\right]}_0^{2\pi }{\left[\frac{-\cos g(r)}{2}\right]}_0^{\infty }$$

$$H$$

$$\mathcal{H}=-\left(\frac{1}{2m}{\nabla }^2+\mu \right){\tau }_z+\lambda {\mathbf{m}}_{\text{sk}}(\mathbf{r})\cdot \mathbf{\sigma }+\text{Re}(\Delta ){\tau }_x-\text{Im}(\Delta ){\tau }_y$$

$$L=-i{\partial }_{\phi }+\frac{n}{2}{\sigma }_z-\frac{b}{2}{\tau }_z$$

$${\Psi }_{\phi }^l(\phi )=e^{i\phi \left(l-\frac{n}{2}{\sigma }_z+\frac{b}{2}{\tau }_z\right)},$$

$$l\in \{\begin{array}{ll}\mathbb{Z}& \text{if }n+b\text{ even}\\ \mathbb{Z}+\frac{1}{2}& \text{if }n+b\text{ odd.}\end{array}$$

$${\mathcal{H}}^l(r)$$

$${\tilde{\mathcal{H}}}^l(r)=\frac{1}{\sqrt{r}}\left({\mathcal{H}}^l(r)+\frac{1}{2mr}{\partial }_r{\tau }_z-\frac{1}{8m{r}^2}{\tau }_z\right).$$

$${\mathcal{H}}_R$$

$${\mathcal{H}}_{\text{soc}}=-2i\alpha {\tau }_z{e}^{-i{n}^{\prime }\phi {\sigma }_z}\left(-{\sigma }_y{\partial }_r+\frac{1}{r}{\sigma }_x{\partial }_{\phi }\right)$$

$${\tilde{\mathcal{H}}}_{\text{soc}}^l(r)=\frac{\alpha }{\sqrt{r}}\left[2i{\sigma }_y{\tau }_z{\partial }_r+\frac{n^{\prime }-1}{r}i{\sigma }_y{\tau }_z+\frac{2l}{r}{\sigma }_x{\tau }_z+\frac{b}{r}{\sigma }_x\right]$$

$$m_{\beta }(\mathbf{r})=m_{\beta }({\mathbf{r}}_i)+(\mathbf{r}-{\mathbf{r}}_i)\cdot \nabla m_{\beta }({\mathbf{r}}_i),$$

$${\mathbf{D}}_{ij}=D{\hat{\mathbf{e}}}_z\times \frac{\left({\mathbf{r}}_j-{\mathbf{r}}_i\right)}{a}$$

$$E_{\text{DMI}}$$

$$E_{\text{DMI}}=\frac{\pi \eta D}{2a}\cos {\phi }_0\left[-R{\int }_{g(0)}^{g(\infty )}dg\frac{\sin (2g)}{2g}-{\int }_0^{\infty }drg(r)\right],$$

$$E_{\text{DMI}}=-\mathcal{K}\eta D\frac{R}{a}\cos {\phi }_0.$$

$$\mathbf{m}(\mathbf{r},d)=\frac{\overline{\mathbf{m}}(\mathbf{r},d)}{\left|\overline{\mathbf{m}}(\mathbf{r},d)\right|},$$

$$\overline{\mathbf{m}}(\mathbf{r},d)=$$

$$(1-w)\left[{\mathbf{m}}_{\text{SWS}}(x+d/2,y)+{\mathbf{m}}_{\text{SWS}}(x-d/2,y)\right]$$

$$+w{\mathbf{m}}_{\text{SWS}}(x,y)+{\mathbf{m}}_{\text{bg}},$$

$$w=\Theta (d_1-d){\left(\frac{d_1-d}{d_1}\right)}^2.$$

$$f_{\pm }(\phi )=\phi \pm \{\begin{array}{ll}0& \text{if }d>2d_0\\ \frac{\pi }{4}\left[1+\cos (\pi \frac{d-d_0}{d_0})\right]& \text{if }{d}_0<d\le 2d_0\\ \frac{\pi }{2}& \text{if }d\le d_0\end{array}.$$

$$\overline{\mathbf{m}}(\mathbf{r},d)=$$

$$(1-w)\left[{\mathbf{m}}_{\text{SWS}}^{+}(x+d/2,y)+{\mathbf{m}}_{\text{SWS}}^{-}(x-d/2,y)\right]$$

$$+w{\mathbf{m}}_{\text{DWS}}(x,y)+{\mathbf{m}}_{\text{bg}},$$

$$w=\frac{1}{2}\Theta (d_1-d)\left[1-\cos \left(\pi \frac{d_1-d}{d_1}\right)\right]$$

Full text

I Introduction

The search for Majorana bound states (MBS) in superconducting heterostructures has experienced massive theoretical and experimental efforts in recent years. The enormous interest in such MBS is owed to their intriguing properties, most notably the topological protection against local perturbations and the non-Abelian exchange statistics Alicea (2012); Beenakker (2013), which have fueled ambitions to use MBS as constituents of topological qubits in quantum computers Stern and Lindner (2013); Nayak et al. (2008).

By now, the most prominent platform for Majorana physics are semiconductor-superconductor hetero-nanowires Kitaev (2001); Lutchyn et al. (2010); Oreg et al. (2010), where a combination of strong spin-orbit coupling (SOC), superconductivity, and a Zeeman field yields an effective $p$-wave order parameter and a topological phase. Starting from the first measurements in 2012 Mourik et al. (2012); Rokhinson et al. (2012); Das et al. (2012); Deng et al. (2012), increasingly appealing experimental indications for the existence of MBS could be gained with improved nanowire architectures, as, e.g., in Refs. Deng et al., 2016; Gül et al., 2018. What remains difficult, though conceptionally possible, is to realize real-space braiding in nanowire-based systems.

Branching off from such nanowires, also linear arrangements of magnetic adatoms on superconductors (Shiba chains) have been studied Choy et al. (2011); Nadj-Perge et al. (2013); Pientka et al. (2013); Klinovaja et al. (2013); Braunecker and Simon (2013); Vazifeh and Franz (2013); Brydon et al. (2015); Peng et al. (2015); Christensen et al. (2016); Schecter et al. (2016), in which the internal magnetization replaces the external field. Furthermore, if the magnetic order in a chain or wire is helical rather than ferromagnetic, the winding of the magnetization induces an effective SOC, thus also replacing the need for intrinsic high-SOC materials. This interesting approach was soon expanded theoretically to two-dimensional Shiba lattices Nakosai et al. (2013); Pöyhönen et al. (2014); Röntynen and Ojanen (2015); Li et al. (2016), where non-collinear magnetic phases were shown to lead to a rich topological phase diagram. Signatures of a chiral topological phase have been measured recently in ferromagnetic islands on a superconductor Ménard et al. (2019).

In principle, similar ideas can be applied to bilayer systems where a superconducting substrate is covered by a thin magnetic film rather than decorated with isolated magnetic impurities. Non-collinear phases in magnetic films constitute by themselves a very active and fruitful field of research Hellman et al. (2017). In particular, chiral magnets with Dzyaloshinskii-Moriya interaction Dzyaloshinsky (1958); Moriya (1960) can exhibit phases with stable magnetic skyrmions Mühlbauer et al. (2009); Yu et al. (2010); Heinze et al. (2011). Skyrmions are particle-like textures characterized by a topological winding number $n$ Nagaosa and Tokura (2013), which is defined by

[TABLE]

(in terms of the unit vector $\mathbf{m}$). They give rise to emergent electrodynamics Schulz et al. (2012) and have been studied intensely, in particular, in context of potential applications in spintronics devices Wiesendanger (2016).

Given the technological maturity in the creation and control of skyrmions Finocchio et al. (2016); Fert et al. (2017); Everschor-Sitte et al. (2018), it seems natural to investigate in how far they can, if combined with superconductivity Vadimov et al. (2018), serve as carriers for MBS. However, relatively little work on magnetic skyrmions on superconductors has been done so far. For skyrmions on $s$-wave and $p$-wave superconductors, Yu-Shiba-Rusinov-like bound states have been considered Pershoguba et al. (2016); Pöyhönen et al. (2016). With regard to Majorana physics, effective $p$-wave phases with chiral edge modes have been predicted Chen and Schnyder (2015); Mascot et al. , whereas localized MBS may emerge in elongated skyrmions Güngördü et al. (2018) which resemble nanowires. For more complicated skyrmion textures, Yang et al. Yang et al. (2016) have derived the criterion that a pair of MBS can form for skyrmions of even winding number, with one state being localized at the skyrmion core and the partner state at the outer rim of the skyrmion. We will return to this criterion in more detail later.

The idea of having an MBS captured inside a skyrmion is very appealing, because the existing tools for skyrmion manipulation could potentially be employed to navigate single Majorana modes. It is particularly noteworthy that a recent experiment Hervé et al. (2018) has demonstrated the possibility to move individual skyrmions in the plane with a scanning tunneling microscope tip. In that way, one could envision explicit real-space braiding operations with MBS.

The goal of this paper is to point out a feasible way to generate MBS from magnetic skyrmions. Importantly, the required skyrmions must have a realistic spatial shape. We take the proposal Yang et al. (2016) of MBS in skyrmions of even winding number, in particular double-winding skyrmions (DWS), as a starting point. Skyrmions in chiral magnets are, in contrast, characterized by a winding number $n=1$ (perhaps $n=-1$), which we refer to as single-winding skyrmions (SWS). In the first part of our work, Sec. II, we investigate whether DWS could be formed by the fusion of two SWS in a chiral magnetic film, as the conservation of the total winding number in the system would insinuate. We show that, when a skyrmion merger is enforced, the destruction of one skyrmion is energetically favorable compared to the protection of the topological invariant. We note, though, that DWS may occur in systems where skyrmion formation is driven by different mechanisms, as dipolar interactions or frustration effects Yu et al. (2014); Leonov and Mostovoy (2015); Zhang et al. (2016); Ozawa et al. (2017); Hayami and Motome (2019), which we do not study in this paper.

In the second part, Secs. III and IV, we generalize the criterion for MBS presented in Ref. Yang et al., 2016 in order to cirumvent the requirement of a DWS, by combining an SWS with a vortex in the superconductor. Such skyrmion-vortex pairs have been shown to be stable under suitable conditions Hals et al. (2016). MBS being bound to vortex cores is a concept that has been known for many years Ivanov (2001) for two-dimensional $p$-wave superconductors. Topological superconductivity in two dimensions emerges, for instance, at topological insulator-superconductor interfaces Fu and Kane (2008); Linder et al. (2010), and signatures of MBS at vortices have indeed been reported Sun et al. (2016); Sun and Jia (2017); Liu et al. (2018). Bare vortices are, however, not easily individually controllable. In this regard, the skyrmion-vortex pairs that we address allow one to envision more flexible applications in quantum technology.

II Merger of single-winding skyrmions in chiral magnets

In this section, we investigate a merger of two SWS as a potential way to generate a DWS in a chiral magnet, having in mind the proposal in Ref. Yang et al., 2016. For simplicity, we do not include the proximitized superconductor at this point, because the relevant energy scale in the magnet, given by the exchange interaction, will be large compared to the induced superconductivity in the magnet.

The texture of a circular skyrmion can be expressed as

[TABLE]

in polar coordinates $(r,\varphi)$ relative to the skyrmion center, where $f(\varphi)$ and $g(r)$ denote the angular and radial profiles, respectively. In general, the winding number (1) is

[TABLE]

and one has to impose the boundary conditions $g(0)=\pi$ and $g(\infty)=0$. We assume an exponential radial profile, $g(r)=\pi e^{-r/R}$ with the skyrmion radius $R$, and the angular profile $f(\varphi)=n\varphi+\varphi_{0}$. If $n=1$, the offset angle $\varphi_{0}$ indicates Néel type ($\varphi_{0}=0,\pi$), Bloch type ($\varphi=\pm\pi/2$), or any intermediate configuration (all of which are topologically equivalent).

The energy of the magnetic configuration on a lattice is given by

[TABLE]

where the first term describes exchange coupling of nearest-neighbor sites, the second term is the Dzyaloshinskii-Moriya interaction Dzyaloshinsky (1958); Moriya (1960), the third term indicates easy-plane anisotropy if $K>0$, and the last term accounts for an external field. In Eq. (4), $\mathbf{m}_{i}$ denotes the unit vector parallel to the magnetic moment $\mathbf{M}_{i}=\mu\mathbf{m}_{i}$ at site $i$.

The Dzyaloshinskii-Moriya interaction favors a relative canting of neighboring magnetic moments and thereby enriches the phase diagram as compared to non-chiral magnets. While the ground state at $B=0$ typically exhibits spiral magnetic order, a skyrmion phase may occur at intermediate field strengths followed by collinear alignment at higher fields. Here, we take $\mathbf{D}_{ij}$ to be in-plane and orthogonal to the $\left<ij\right>$ bond (favoring Néel type winding). We choose a set of parameters in Eq. (4) which resembles a Co/Ru(0001) bilayer system Hervé et al. (2018) where stable SWS exist. This is illustrated in Fig. 1, where the energy of SWS and DWS is shown as a function of the skyrmion radius for multiple field strengths. An SWS has an energy minimum at an optimal radius for suitable field strengths. The DWS, in contrast, never develops such a minimum. This is understood as follows: Whenever $n\geq 2$, the rotation of the magnetization along the radial direction continuously varies with $\varphi$ between the preferred Néel type (where $\mathbf{D}_{ij}\cdot(\mathbf{m}_{i}\times\mathbf{m}_{j})>0$), the neutral Bloch type ($\mathbf{D}_{ij}\cdot(\mathbf{m}_{i}\times\mathbf{m}_{j})=0$), and the unfavorable (opposite) Néel type ($\mathbf{D}_{ij}\cdot(\mathbf{m}_{i}\times\mathbf{m}_{j})<0$). In total, these contributions cancel if $f$ is linear in $\varphi$ [see Appendix A for a calculation of the Dzyaloshinskii-Moriya energy of a skyrmion with the general shape given by Eq. (2)]. Although DWS do not form spontaneously, their energy profile is relatively flat at not too large $B$, such that a DWS might still have a sufficient lifetime if it can be created artificially, or even persist as a metastable state if the external field is switched off after the initial nucleation of skyrmions. Therefore we now turn to the DWS creation by a skyrmion merger.

When two SWS are merged, there are two conceivable scenarios: (i) a discontinuous process in which the overall winding number is allowed to change and (ii) the continuous merger resulting in a DWS due to the topological protection of $n$. Both cases are depicted in Fig. 2. In (i), a pointlike discontinuity appears, similar to the case studied in Ref. Milde et al., 2013, at the critical skyrmion distance $d_{0}=2R\text{ln}3$ [with our specific choice of $g(r)$]. The discontinuity can also be interpreted as an emergent antiskyrmion Takashima and Fujimoto (2014). To avoid this discontinuity in (ii), the textures of both SWS must be rotated such that their contributions cannot cancel in the overlap region.

It is clear that both processes involve an energy barrier. Therefore, any such process has to be enforced externally, e.g., by moving the SWS with a tunneling tip Hervé et al. (2018). In order to decide whether (i) or (ii) would be energetically favorable, we have implemented a specific instance of a continuous and a discontinuous merger (the ones shown in Fig. 2), starting from the same initial configuration. Details on the magnetic texture during the merger can be found in Appendix B. By means of Eq. (4), the energy can be tracked throughout each process. We present the result in Fig. 3. Note that the sections displayed in Fig. 2 do not represent the full size of the system on which the calculations were performed. We use the same parameters as in Fig. 1 and $R=73{.}6a$, corresponding to the energy minimum at $B=200\,\text{mT}$.

The main observation is that the barrier related to the discontinuity in (i) is small compared to the energy required in (ii). The reason is, qualitatively, that the rotation of the SWS textures in (ii) acts against the Dzyaloshinskii-Moriya term in the entire area of each SWS. In contrast, the “collision area” in (i) in the vicinity of the discontinuity is comparably small, containing only a few sites at which the exchange energy becomes large. To be more precise, the energy related to the Dzyaloshinskii-Moriya interaction is a numerical constant times $D(R/a)\cos\varphi_{0}$ (see Appendix A), where $\varphi_{0}$ changes from [math] to $\pm\frac{\pi}{2}$ in going from panel (a) to panel (d) in Fig. 2. For our parameters, this corresponds to $1600\,\text{meV}$ per skyrmion. This is consistent with the barrier height in Fig. 3, given that the overlap of the skyrmions is not taken into account in our estimate. In the discontinuous merger, on the other hand, the barrier height is determined by the exchange energy in an area $A$ around the discontinuity at $d_{0}$ in which $\mathbf{m}$ changes its direction. This area scales as $A\propto d_{0}^{2}$. The exchange energy density can be expressed as $\propto J(\nabla\mathbf{m})^{2}$, where the derivative of $\mathbf{m}$ scales with $1/d_{0}$. Hence, the exchange energy contribution from the vicinity of the discontinuity is a constant proportional to $J$ which does not scale with the skyrmion size.

In this section, we did not account for the proximitized superconducting layer. In principle, the superconductor can lead to additional energy contributions in the merger process. For instance, bound states similar to Yu-Shiba-Rusinov states exist in the SWS Pershoguba et al. (2016); Pöyhönen et al. (2016) and can induce antiferromagnetic interactions Yao et al. (2014). The energy scale of such effects is proportional to the superconducting gap. While these interactions can have important consequences on magnetic impurities on a superconductor (where they compete only with the Ruderman-Kittel-Kasuya-Yosida interaction), the direct magnetic exchange interaction of strength $J$ will be dominant in the case of a magnetic film covering the superconductor. Referring again to a Co/Ru(0001) bilayer, the superconducting transition temperature of ruthenium is $T_{c}=470\,\text{mK}$ Hulm and Goodman (1957), corresponding to $\Delta(T=0)\approx\frac{\pi}{e^{\gamma}}k_{B}T_{c}=0.07\,\text{meV}$, which is negligible compared to $J=13.1\,\text{meV}$ Hervé et al. (2018).

In total, we find that in an enforced SWS merger, the chiral magnetic film would favor the deletion of a skyrmion – notably despite its topological nature – over the formation of a DWS. Although our specific implementations of the two scenarios do not represent fully optimized paths, the result in Fig. 3 demonstrates the generic energy proportions, because the large discrepancy between small-area versus large-area energy cost cannot be overcome by small adjustments. The only possible exception would be extremely small skyrmions, which seem, however, unfavorable to localize a MBS. In essence, this rules out the possibility to generate MBS by DWS creation in the known Dzyaloshinskii-Moriya-driven skyrmion host materials. We turn therefore to cases where MBS are obtained from SWS in the following sections.

III Majorana criterion in a skyrmion-vortex pair

Here, we present a generalization of the argument given in Yang et al., 2016 to circumvent the necessity of DWS. Namely, we consider not a bare skyrmion on a superconductor, but the composite object of a skyrmion and a vortex in the superconductor Hals et al. (2016), thus incorporating the phase winding of the superconducting order parameter in addition to the skyrmion winding number. A similar argument solely based on the effect of phase winding was recently presented for superconducting hollow cylinders Lutchyn et al. ; Vaitiekenas et al. . Consider the Hamiltonian $H=\int d^{2}r\Psi^{\dagger}(r)\mathcal{H}\Psi(r)$ with

[TABLE]

in Nambu space, $\Psi=(c_{\uparrow},c_{\downarrow},c^{\dagger}_{\downarrow},-c^{\dagger}_{\uparrow})$, for electrons with effective mass $m$ at a chemical potential $\mu$ exchange-coupled to the magnetic skyrmion texture $\mathbf{m}_{\text{sk}}(\mathbf{r})$, in the presence of a superconducting $s$-wave gap $\Delta(\mathbf{r})$. The Pauli matrices $\tau_{x,y,z}$ and $\sigma_{x,y,z}$ act in particle-hole space and in spin space, respectively, and we have set $\hbar=1$.

We set $\mathbf{m}_{\text{sk}}(\mathbf{r})$ as in Eq. (2). The superconducting order parameter is spatially dependent to account for the vortex, $\Delta(\mathbf{r})=e^{ib\varphi}\Delta(r)$ with integer $b$ and $\Delta(r)=\Delta(1-e^{-r/R_{v}})$, where the vortex and the skyrmion cores are both located at the origin. This configuration can be energetically stable Hals et al. (2016); Baumard et al. (2019); Dahir et al. (2019). Then, the modified angular momentum operator

[TABLE]

commutes with the Hamiltonian, such that the eigenstates of the system can be separated into angular and radial parts, $\Psi^{l}(r,\varphi)=\Psi^{l}_{\varphi}(\varphi)\Psi^{l}_{r}(r)$, where $l$ denotes the eigenvalue of $L$. Following from Eq. (6), the angular part has the form

[TABLE]

such that $\Psi^{l}_{\varphi}(\varphi)=\Psi^{l}_{\varphi}(\varphi+2\pi)$ imposes the constraint on $l$ that

[TABLE]

Furthermore, the particle-hole operator $C=\sigma_{y}\tau_{y}K$ (with complex conjugation $K$) and $L$ have the commutation relation $[L,C]=2LC$, and consequently $LC\Psi^{l}=-lC\Psi^{l}$. Therefore, if $\Psi^{l}$ is an eigenstate of $C$ (a Majorana mode), then necessarily $l=0$. This restricts the emergence of MBS to skyrmion-vortex pairs where $n+b$ is even, as in the realistic case $n=b=1$. Thus, in such a pair, the existence of DWS is no longer required for the emergence of MBS.

The eigenproblem for $\mathcal{H}$ can be reduced to one dimension for the radial component, $H^{l}(r)\Psi^{l}_{r}(r)=\epsilon^{l}\Psi^{l}_{r}(r)$, with

[TABLE]

Because the radial probability density is $r|\Psi^{l}_{r}(r)|^{2}$, it is convenient to consider $\Phi^{l}(r)=\sqrt{r}\Psi^{l}_{r}(r)$. The Hamiltonian acting on $\Phi^{l}(r)$ reads

[TABLE]

The minimal Hamiltonian (5) could contain further terms that we have not taken into account so far. The presence of skyrmions implicitly requires SOC, which is missing in the Hamiltonian. Importantly, in the presence of SOC, the reasoning above remains valid only if the commutation with $L$ is preserved. This can be achieved if the usual Rashba Hamiltonian

[TABLE]

is replaced by a generalized SOC

[TABLE]

in which $n^{\prime}$ is chosen equal to $n$. Equation (12) describes a SOC in which the spin direction winds $n^{\prime}$ times as one encircles the origin. In result, it varies with the polar angle between the Rashba and Dresselhaus types whenever $n^{\prime}\neq 1$. In particular, $\mathcal{H}_{\text{soc}}$ is identical to the SOC considered in Ref. Yang et al., 2016 for $n^{\prime}=2$.  111Up to a momentum-independent term that also commutes with $L$. The SOC of Ref. Yang et al., 2016 can be written in polar form as $e^{-i\varphi\sigma_{z}}\mathcal{H}_{R}+\mathcal{H}_{R}e^{i\varphi\sigma_{z}}$. Unconventional SOC with a vortex-like structure was also suggested to lead to MBS in cold-atom systems in Ref. Sato et al., 2009. However, in that case only a complex phase winding $e^{im\varphi}$ with integer $m$ was introduced in the prefactor of the Rashba Hamiltonian, instead of a spin rotation, which is different from $\mathcal{H}_{\text{soc}}$. For arbitrary $n^{\prime}$, Eq. (12) represents a rather artificial type of SOC which is compatible with a skyrmion of the same winding number. Contrarily, in our proposal of a skyrmion-vortex pair with $n=b=1$, the SWS demands that also $n^{\prime}=1$, such that $\mathcal{H}_{\text{soc}}=\mathcal{H}_{R}$ recovers the ordinary Rashba Hamiltonian. After the reduction to the one-dimensional radial probability $\Phi^{l}(r)$, $\mathcal{H}_{\text{soc}}$ takes the form

[TABLE]

and can be added to Eq. (10). In total, electrons will be subject to the sum of the background SOC of Eq. (12) and the synthetic SOC generated by the spatially varying magnetization. We note that the criterion for the existence of MBS is independent of the presence or absence of (commutation-preserving) background SOC. The significance of $\mathcal{H}_{\text{soc}}$ is to reduce the overlap of the two Majorana modes at the skyrmion core and the rim of the system. To obtain the same effect based on synthetic SOC at $\alpha=0$, the skyrmion must have multiple spin flips in the radial direction Yang et al. (2016). Numerical results on the MBS will be shown in the next section.

Another remark is at hand concerning the role of SOC. One canonical example for the formation of localized MBS is a vortex in a two-dimensional spinless $p$-wave superconductor Alicea (2012); Beenakker (2013); Ivanov (2001), which may effectively result from $s$-wave pairing, Rashba SOC, and Zeeman splitting. In our language, such a vortex would correspond to $n=0$ and $b=1$, where obviously $n+b$ is odd. Yet, our criterion does not contradict the well-known results about MBS in that case. The reason is the requirement of Rashba SOC, which breaks the commutation with $L$ for $n=0$, such that our criterion is not valid. For a uniform out-of-plane Zeeman field 222Note that this is a stronger requirement than $n=0$, for which a non-winding in-plane component would still be allowed. without skyrmions, one can formulate a different criterion, namely, for a SOC of the general form Eq. (12) and a vorticity $b$, the operator $L^{\prime}=-i\partial_{\varphi}+\frac{1}{2}n^{\prime}\sigma_{z}-\frac{1}{2}b\tau_{z}$ commutes with the Hamiltonian. Though formally similar to $L$, $L^{\prime}$ is constructed based on the winding number of the SOC instead of a skyrmion. In analogy to the arguments presented above, the Majorana criterion $(n^{\prime}+b)\,\text{mod}\,2\overset{!}{=}0$ can then be derived. For Rashba SOC, where $n^{\prime}=1$, it follows that vortices with odd $b$ will harbor MBS (similar to the result in Ref. Lutchyn et al., ). The canonical case is recovered with $b=1$. Whenever a skyrmion and SOC coexist, they are only compatible in terms of commutation relations if $n=n^{\prime}$ (then, indeed $L=L^{\prime}$). Otherwise, the system does not have any rotational symmetry that could give rise to a localized MBS. We will present an example of such a case at the end of the next section.

As an interesting side remark, in Majorana systems based on either $L$ or $L^{\prime}$, equivalent formulations can often be found by suitable transformations. One can check that, once the Hamiltonian is reduced to the Majorana sector Sau et al. (2010); Yang et al. (2016), $n=b=0$ combined with unconventional SOC where $n^{\prime}=0$ is equivalent to $n=2$ and $b=0$ without SOC. The latter case was solved in Ref. Yang et al., 2016, see, e.g., Fig. 5 therein for the analytical Majorana wavefunctions. In Ref. Sato et al., 2009, it was shown that a vortex in the superconducting order parameter can be mapped to a complex-phase vortex in the SOC. Thus, for a uniform out-of-plane Zeeman field, the case $n^{\prime}=b=1$ is equivalent to $b=1$ and Rashba SOC with a vortex.

Besides SOC, we have not incorporated the electromagnetic vector potential in Eq. (5), thus orbital effects and screening supercurrents in response to the magnetic texture are not included. This is a valid approximation for the following reasons. First and foremost, any vector potential of the form $\mathbf{A}=A(r)\hat{\mathbf{e}}_{\varphi}$ (corresponding to the out-of-plane magnetic flux in a vortex) does not change the commutation relation $[H,L]=0$ and hence does not affect the Majorana criterion. Moreover, here we focus on a system where both the chiral magnetic layer and the superconducting film are thin. The supercurrent that emerges as a consequence of a domain wall in a superconductor-ferromagnet bilayer (we may view skyrmions as a limiting case of general domain wall structures in this regard by continuously shrinking a circular domain to a skyrmion core) will be bounded by $J_{\text{max}}=cM(d_{F}/\delta)$ Burmistrov and Chtchelkatchev (2005), where $c$ is the speed of light, $M$ the saturation magnetization, $d_{F}$ the thickness of the magnetic layer, and $\delta$ the width of the domain wall (here: related to the skyrmion radius). As we focus on $d_{F}\ll\delta$, where $d_{F}$ may correspond to monoatomic layer thickness Hervé et al. (2018) whereas $\delta$ will be on the order of 10 – 100 nanometers, we can expect that $J_{\text{max}}$ is insignificant in the thin-film limit. In this limit, vortices may not form spontaneously, but can still be induced by an external out-of-plane magnetic field Baumard et al. (2019). When skyrmions and vortices coexist, vortex pinning to the skyrmion center is energetically favorable Hals et al. (2016); Baumard et al. (2019); Dahir et al. (2019).

IV Results on Majorana bound states in skyrmion-vortex pairs

Now we present our results obtained for the eigenstates in skyrmion-vortex pairs by exact diagonalization of Eq. (10) [together with Eq. (13) whenever $\alpha\neq 0$]. The low-energy spectra are shown in Fig. 4. Our main result concerns the case with realistic winding numbers $n=1$ and $b=1$ and an exponential radial shape of both the skyrmion and the vortex in the presence of background SOC [Fig. 4(a)]. The skyrmion is chosen to be of Néel type. A pair of zero-energy states is found at $l=0$, where one of these states is localized at the core of the skyrmion, whereas the other one is located at the rim of the system. These are the MBS predicted by the criterion in the previous section. Thus, it is possible to nucleate a MBS in a skyrmion without a demand for higher winding numbers. The probability densities of the two zero modes are shown in Fig. 5(a).

Apart from the MBS, this system has further noteworthy features. First, a band of close-to-zero energy states appears at any eigenvalue of $L$. All of them are located at the rim of the system. A close-up of the spectrum (inset in Fig. 4(a)) reveals that this band exhibits a weak linear dispersion, thus none of the states at $l\neq 0$ are truly at zero energy. Rather, this band represents a chiral Majorana mode at the rim, where the tiny slope of the band with respect to angular momentum yields the group velocity when multiplied with the system radius. The presence of such a chiral mode is consistent with the results of Ref. Chen and Schnyder, 2015. A second feature is the emergence of further bound states at the core besides the MBS. In Fig. 5(b), two examples of such states are compared to the MBS. The components $(u_{\uparrow},u_{\downarrow},v_{\downarrow},-v_{\uparrow})$ of a MBS wavefunction generally have to satisfy the conditions $u_{\uparrow}=c^{*}v^{*}_{\uparrow}$, $u_{\downarrow}=c^{*}v^{*}_{\downarrow}$, where $c$ with $|c|=1$ is the eigenvalue of the particle-hole operator $C$. In contrast, all components are independent for the non-Majorana bound states. These additional bound states exist for two reasons. First, it is expected Pershoguba et al. (2016); Pöyhönen et al. (2016) that a skyrmion on a superconductor will generally induce bound states which relate to Yu-Shiba-Rusinov states in the limit of a very small skyrmion (a tiny skyrmion is conceptionally similar to a magnetic impurity). Second, it is known that also a vortex line in a superconductor causes bound states Caroli et al. (1964). Interestingly, the core states close to $l=0$ are reminiscent of a band, signaling an emergent chiral character. The slope has opposite sign compared with the band of rim states. The core and rim bands can be interpreted as the inner and outer chiral modes of a topological superconductor with the geometry of an annulus, where the hole stems from the vortex. In our calculations, bound states exist at various energies both above and below the effective superconducting gap. The localization of states outside the gap is not necessarily exponential. In Fig. 4, we have used the qualitative criterion that $70\,\%$ of the total probability of a state must be accumulated on $20\,\%$ of the radial extent of the system for this state to be classified as localized.

For comparison, Fig. 4(b) shows the spectrum obtained from a bare skyrmion in the absence of a vortex, i.e., $b=0$. No MBS are present. Note that $l$ takes half-integer values in this case, which is similar to the reversed situation $n=0$ and $b=1$. The outer chiral mode persists, but without a state at exactly zero energy. Also, localized states at the core are still present, though not as many as in the presence of a vortex, and without the formation of a band-like structure.

Next, we investigate the impact of the vortex radius, which will be of practical importance for experiments. To obtain MBS from skyrmion-vortex pairs, it is generally desirable that the vortex is small compared to the skyrmion: If electrons inside the skyrmion would essentially not feel a pairing potential because the vortex is too big, it is clear that the joint effect of the skyrmion texture and superconductivity would not be observable. In Fig. 4(a), $R_{v}$ was therefore set to $5$ compared to the skyrmion radius $R=300$. In panels (c) and (d), the radii are $R_{v}=200$, $R=300$ and $R_{v}=400$, $R=100$, respectively, with otherwise identical parameters. In both cases, MBS are still found, thus they persist over a large range of relative sizes of the skyrmion compared to the vortex. However, the localization of the inner MBS becomes weaker because of the reduced superconducting gap close to the core, see Figs. 5(c) and 5(d). In consequence, the overlap between the inner and outer MBS increases, thus the numerically obtained energy in the last case, $\epsilon_{\text{MBS}}\approx 8\times 10^{-8}$, is an order of magnitude larger than in the first panel (small $R_{v}$), where $\epsilon_{\text{MBS}}\approx 6\times 10^{-9}$.

For comparison with the earlier results of Ref. Yang et al., 2016, we have also considered skyrmion-vortex pairs in which the effective SOC stems from a winding of the magnetization along the radial direction, while $\alpha=0$. More specifically, $g(r)=\pi(1-r/R)$ if $r<25R$ and $g(r)=0$ otherwise, with $R=30$ (thus, the outer radius of the skyrmion texture is $25R$). In Fig. 4(e), the system is restricted to the skyrmion, whereas the skyrmion is embedded in a region of uniform magnetization in panel (f). In agreement with Ref. Yang et al., 2016, the radial winding of the magnetization is suitable to stabilize MBS. Panel (e) is reminiscent of the results in Ref. Yang et al., 2016, with a MBS at the core and one at the rim of the skyrmion. Notably, if the skyrmion has a surrounding, the outer mode is not localized at the rim of the system, but rather delocalized between the outer radius of the skyrmion and the rim of the system, cf. Figs. 5(e) and 5(f). This is ascribed to the absence of background SOC. Here, the effective SOC is only present inside the skyrmion, and no SOC emerges otherwise. In contrast, in the previously discussed cases with $\alpha\neq 0$, the uniform background SOC enforces the outer modes to reside at the system rim. This difference does not become apparent as long as the system size is restricted to the skyrmion itself. In addition, the surrounding of the skyrmion harbours numerous further bound states at energies below the effective gap, which form multiple weakly dispersive bands at different energies. The number of such states can be expected to increase further for larger system sizes. As in all other panels, multiple bound states near the core appear at various energies inside and outside the effective gap. We have also computed energy spectra for DWS without a vortex ($n=2$, $b=0$, not shown), which are qualitatively similar to panels (e) and (f).

Finally, we discuss the impact of SOC with respect to the rotational symmetry of the system, as expressed by the commutation with the operator $L$, cf. Eq. (6). If the symmetry is broken by a mismatch of $n$ and $n^{\prime}$, the reduction to a one-dimensional system fails. Therefore, we now diagonalize the Hamiltonian numerically in two dimensions. We consider two cases: (i) a skyrmion-vortex pair with $n=b=1$ and Rashba SOC as before, where $[\mathcal{H},L]=0$, and (ii) a DWS in the absence of a vortex ($n=2$, $b=0$) with Rashba SOC, such that $n\neq n^{\prime}$ and $[\mathcal{H}_{R},L]\neq 0$ breaks the commutation relation. In Fig. LABEL:Fig:Rashba, the probability density of the lowest-energy state is shown for both systems. In the first case, a symmetrically localized core state near zero energy is found, which is consistent with the existence of MBS in this case (only the inner state is shown in the figure) and qualitatively confirms our previous results. We note, though, that finite-size effects are much larger than in the one-dimensional calculations, owed to computational limitations in two dimensions. The energy of the MBS is therefore still $\approx 8\times 10^{-5}$ times the hopping amplitude in our calculations. In the second case, the lowest-energy state clearly reveals the broken symmetry. It is also not well-localized anymore and raises in energy. This state and its partner state (related through particle-hole symmetry $\mathcal{C}$) do not form a pair of spatially separated inner and outer states, as expected for MBS in skyrmions. These results demonstrate that the unconventional form of the background SOC $\mathcal{H}_{\text{soc}}$ in Eq. (12) is important when $n\neq 1$ and cannot be substituted by Rashba SOC. On the other hand, the compatibility of the skyrmion-vortex pair with $n=1$ with ordinary Rashba SOC represents a further advantage of our proposal compared to skyrmions of higher winding numbers. In Fig. LABEL:Fig:Rashba, we used a square lattice with $76\times 76$ sites and the numerical parameters, in units of the nearest-neighbor hopping, $\mu=0$, $\alpha=0.8$, $\Delta=0.5$, and a skyrmion radius of $R=12a$ (vortex radius $R_{v}=R/10$ in the first case).

V Summary

We have studied possibilities to obtain localized Majorana bound states (MBS) from skyrmions at a superconductor-chiral magnet interface. We have shown that the formation of skyrmions with winding number two, i.e. DWS, in in such systems is energetically unfavorable, even if one can enforce a merger of two circular skyrmions of winding number one. Instead, one of these skyrmions is deleted through a point-like discontinuity. Therefore, DWS as MBS carriers seem unrealistic.

Nevertheless, we have demonstrated that a composite object consisting of a skyrmion and a concentrical vortex can support a localized MBS at its core, if the winding numbers of the skyrmion and the vortex add up to an even number. This allows the creation of MBS with experimentally accessible skyrmions of winding number one. Importantly, this case is also compatible with ordinary Rashba spin-orbit coupling, which would break the rotational symmetry of the system otherwise. We have calculated the Majorana energies and their wavefunctions, see Fig. 4 and 5, showing the localization on the radial axis and the spatial separation of the inner and the outer mode. MBS at skyrmion cores indicate a route towards direct spatial control of non-Abelian quantum states in two dimensions.

We hope that our proposal will be realized in future experiments. Already now, magnetic bilayer systems hosting skyrmions are available. In addition, there is ongoing progress towards the observation of skyrmions below the superconducting transition temperature in such bilayers Mougel et al. (2019). In order to generate skyrmion-vortex pairs rather than bare skyrmions, the superconducting layer should preferably be a type-II superconductor. Alternatively, this layer could be fabricated sufficiently thin compared to the penetration depth of the material.

We expect that future work will address the stability of the MBS under distortions from the circular shape of the skyrmions. Shape distortions can be caused, e.g., by the geometry of the sample Hervé et al. (2018); Jin et al. (2017), or by dynamic oscillatory modes which might be activated when operations are performed with the skyrmions. Furthermore, real samples will harbour multiple skyrmions or skyrmion-vortex pairs. Hence, hybridization effects of the electronic states bound to different skyrmion-vortex pairs should be studied, as well as many-skyrmion arrangements in general.

Note added: Recently, we learned about Ref. Garnier et al., , which also addresses topological superconductivity induced by magnetic skyrmions.

Acknowledgements.

We thank Dmitry Aristov, Stefan Blügel, Igor Burmistrov, Marie Hervé, Markus Garst, Leonid Glazman, Jens Paaske, and Wulf Wulfhekel for useful and inspiring discussions. The work was supported by the Deutsche Forschungsgemeinschaft via the grants MI 658/9-1 (joint DFG-RSF project) and MI 658/7-2 (Priority Programme 1666 “Topological Insulators”). I.V.G. acknowledges support by the Russian Science Foundation through grant No. 17-12-01182.

Appendix A Dzyaloshinskii-Moriya energy of skyrmions

In this appendix, we provide a proof for the two statements from Sec. II of the main text that (i) the Dzyaloshinskii-Moriya contribution to the energy vanishes in a DWS, and (ii) the rotation of the magnetic texture in a SWS (by alteration of $\varphi_{0}$) is counteracting the Dzyaloshinskii-Moriya interaction, leading to an energy barrier proportional to $DR$.

For slowly varying magnetization (compared to the lattice constant $a$), $\mathbf{m}(\mathbf{r})$ can be expanded at site $i$,

[TABLE]

($\beta=x,y,z$) to get a continuum version of Eq. (4). Using the vector

[TABLE]

acting on a pair of nearest neighbors $\left<ij\right>$, the overall Dzyaloshinskii-Moriya energy for a magnetic texture as in Eq. (2) then takes the form

[TABLE]

Here, we denote by $\boldsymbol{\rho}$ all possible nearest-neighbor bond vectors in a given lattice type. The factor of $\frac{1}{2}$ corrects the double-counting of bonds in the integration. When $f(\varphi)=n\varphi+\varphi_{0}$ is inserted, the integration over $\varphi$ cancels all terms whenever $n\neq\pm 1$. In particular, a DWS has $E_{\text{DMI}}=0$, proving statement (i). For $n=1$,

[TABLE]

where we have used $\rho_{x}^{2}+\rho_{y}^{2}=a^{2}$ and introduced the coordination number $\eta$ of the lattice. A similar expression was found in Ref. Wang et al., 2018. The first integral always yields $\frac{1}{2}\text{Si}(2\pi)\approx 0{.}71$ for any monotonous $g(r)$ with correct boundary conditions. The second integral gives $\pi R$ for our specific choice of $g(r)$. In total,

[TABLE]

with the numerical constant $\mathcal{K}\approx 6.05$. With respect to a rotation of the skyrmion texture, i.e. changing $\varphi_{0}$, the Dzyaloshinskii-Moriya energy is proportional to $\cos\varphi_{0}$, such that $\varphi_{0}=0$ (Néel), $\varphi_{0}=\frac{\pi}{2}$ (Bloch) and $\varphi_{0}=\pi$ (anti-Néel) correspond to minimal, zero, and maximal energy, respectively, if $n=1$, and vice versa if $n=-1$. This proves statement (ii). With the parameters used in Sec. II on the triangular lattice ($\eta=6$), one obtains that rotating the texture of one skyrmion by $\pm\frac{\pi}{2}$ changes the Dzyaloshinskii-Moriya energy in the system approximately by $1600\,\text{meV}$.

Appendix B Details on the merger implementation

Here, we describe how the continuous and discontinuous mergers of two SWS, cf. Fig. 2, have been implemented to obtain the result shown in Fig. 3. First, we introduce the notion of a skyrmion field $\mathbf{m}_{\mathrm{sk}}(x,y)$, which is obtained from Eq. (2) by substracting the uniform background $\mathbf{m}_{\mathrm{bg}}=\hat{\mathbf{e}}_{z}$. The boundary conditions for $g(r)$ imply that the skyrmion field vanishes at infinity. We denote the SWS and DWS fields $\mathbf{m}_{\text{SWS}}$ and $\mathbf{m}_{\text{DWS}}$, respectively. A configuration with more than one skyrmion is now obtained by summing $\mathbf{m}_{\textrm{bg}}$ and multiple skyrmion fields and subsequently normalizing to $\left|\mathbf{m}(\mathbf{r})\right|=1$. For the merger process, we initially place two identical Néel SWS (with optimal radius, cf. Fig. 1) on the $x$ axis at a distance $d$, which is gradually decreased.

For the discontinuous merger, the overall magnetic texture is

[TABLE]

with the not normalized vector field

[TABLE]

and the weight

[TABLE]

By means of $w$, we smoothly interpolate between the superposition of two skyrmions and a single skyrmion for small $d$ ($d_{1}$ is specified below) in order to ensure that the final skyrmion is identical to each of the initial ones. Otherwise, the final state would also be a SWS, but with a different radial profile.

This process starts with two SWS, thus $n_{\text{tot}}=2$, but terminates with a single SWS, $n_{\text{tot}}=1$. The change of the total winding number happens at a critical distance $d_{0}$, where Eq. (19) has a discontinuity at $(x,y)=(0,0)$. Namely, $\overline{\mathbf{m}}(0,d_{0})=0$, such that the normalization fails. We find that $d_{0}=2R\text{ln}3$ for the exponential radial profile. It is easy to show that $\overline{\mathbf{m}}\neq 0$ in any other case. In Eq. (21), $d_{1}$ was set to $0{.}8d_{c}$.

In the continuous merger, where $n_{\text{tot}}$ is conserved, a point where $\overline{\mathbf{m}}=0$ has to be avoided by rotating the two SWS at $d>d_{0}$. More precisely, the offset angle $\varphi_{0}$ in $f(\varphi)$ is smoothly changed from [math] to $\pm\frac{1}{2}\pi$, with $+$ ($-$) for the left (right) SWS:

[TABLE]

The angles at $d<d_{0}$ are chosen to be precisely $\pm\frac{1}{2}\pi$ in order for the resulting skyrmion to have the rotational symmetry described by $L$, Eq. (6). We show the offset angle of $f_{+}(\varphi)$ in Fig. 7 for clarity. The respective fields are labeled $\mathbf{m}_{\text{SWS}}^{\pm}$.

The vector field prior to normalization now reads

[TABLE]

and here we use the weight

[TABLE]

to smoothly transition into the final state, with $d_{1}=0.8d_{0}$ as before. The final-state skyrmion field $\mathbf{m}_{\text{DWS}}$ is slightly shrinked in comparison to the initial SWS, $R_{\text{DWS}}\overset{\sim}{=}\frac{2}{3}R_{\text{SWS}}$ to reduce the energy of the state.

The details of this specific implementation of the skyrmion mergers are to some degree arbitrary, and some optimization with respect to the energy throughout the process could be done. This effort seems futile, though, as the modification of $\varphi_{0}$ in the continuous case must persist. This feature is mainly responsible for the energy difference between the two scenarios, as discussed in the main text, such that most modifications of the implementation would not change the physical result, namely, the deletion of one SWS in an enforced skyrmion merger. In order to reverse the energy proportions, the skyrmion would have to be shrinked in size: the barrier in the discontinuous merger is approximately constant, whereas it is proportional to $R$ in the continuous merger (see previous appendix). Figure 3 would suggest to consider skyrmions which are roughly five times smaller in order to equalize the energy cost in both processes. In practice, however, $R$ is not an independent parameter but has a nontrivial dependence on the coefficients in Eq. (4). Therefore, such an estimate cannot be extracted from our results.