Majorana bound states in magnetic skyrmions

Magnetic skyrmions are highly mobile nanoscale topological spin textures. We show, both analytically and numerically, that a magnetic skyrmion of an even azimuthal winding number placed in proximity to an $s$-wave superconductor hosts a zero-energy Majorana bound state in its core, when the exchange coupling between the itinerant electrons and the skyrmion is strong. This Majorana bound state is stabilized by the presence of a spin-orbit interaction. We propose the use of a superconducting trijunction to realize non-Abelian statistics of such Majorana bound states.

Figure 1

Radial wave function $\tilde{\mathrm{\Phi }}\left(r\right)={[{\tilde{u}}_{\uparrow }\left(r\right),{\tilde{u}}_{\downarrow }\left(r\right)]}^T$ of the MBS ${\chi }^{\mathrm{in}}$, with the solid (black) and dashed (red) lines denoting ${\tilde{u}}_{\uparrow }\left(r\right)$ and ${\tilde{u}}_{\downarrow }\left(r\right)$, respectively. We set $\mu =0$, $\alpha =1$ meV, ${\mathrm{\Delta }}_0=0.5$ meV, $R=100$ nm, $r_0=10$ nm, $n=2$, and $pR$ much larger than the shown range. The discontinuities in ${\tilde{\mathrm{\Phi }}}^{\prime }\left(r\right)$ at $r=r_0$ [less discernible for ${\tilde{u}}_{\downarrow }\left(r\right)]$ arise from singularities in $f^{\prime \prime }\left(r\right)$, Eq. (7).

Figure 2

Excitation spectrum (left panels) and the probability density $\rho \left(r\right)=r|{\mathrm{{\rm Y}}}^0{\left(r\right)|}^2$ of the lowest eigenstate in the $l=0$ sector (right panels, with the $u_{\uparrow }^0,u_{\downarrow }^0,v_{\uparrow }^0$, and $v_{\downarrow }^0$ components shown in yellow, blue, black, and red, respectively), for a system size $L=1\mu \mathrm{m}$ with Dirichlet boundary conditions and for realistic parameters $m=m_e$, $\mu =0$, $\alpha =1$ meV, ${\mathrm{\Delta }}_0=0.5$ meV, $\phi =0$, $R=100$ nm, $r_0=0$, and $n=2$. (a),(b) Skyrmion with $p=1$. (c),(d) Skyrmion with $p=10$. (e),(f) Skyrmion with $p=1$ in the presence of a SOI described by $H_{so}$, with $l_{so}=100$ nm. The energy of the lowest eigenstate in the $l=0$ sector returned by numerics is 3240 neV in (a), 20 neV in (c), and 201 neV in (e). Excitation energies for $l<0$ are not plotted, being identical to those for $l>0$. In right panels, wave functions are plotted for $r\le L/2$, to be separated from the degenerate state localized at the outer boundary. In (d) and (f), $|v_{\uparrow ,\downarrow }^0|=|u_{\uparrow ,\downarrow }^0|$ and only the hole components are plotted.

Figure 3

(a) Energy spectrum of $\mathcal{H}$ found numerically in the tight-binding model with hopping amplitude $t$ and lattice constant $a$, in the presence of a skyrmion with $n=2$ and $p=25$. Here, $m$ labels the eigenstates. The inset shows ten energies closest to the chemical potential. (b) The probability density profile of the lowest eigenstate at near-zero energy $E_M/t=1.6\times {10}^{-7}$. This electron state and its hole partner at $-E_M$ can be identified with two weakly overlapping MBSs (one at the center and one at the system edge) and are very well separated in energy from the remaining states. We have used $\alpha /t=1.2$, ${\mathrm{\Delta }}_0/t=0.4$, $\mu /t=0$, and $R/a=3.5$.

Figure 4

Braiding of two MBSs ${\chi }^A$ and ${\chi }^B$ (solid circles): (a) real-space manipulation, with the hollow circle indicating the temporary position of ${\chi }^A$, and (b) corresponding phase rotations in the MBS wave function, with the wavy line denoting a branch cut. The solid and dashed arrows show the trajectories of ${\chi }^A$ and ${\chi }^B$, respectively.

Figure 5

Comparison of analytical and numerical solutions for $u_{\uparrow }^0$ (upper panel) and $u_{\downarrow }^0$ (lower panel) components of the probability density $\rho \left(r\right)={r\left|\mathrm{\Phi }\left(r\right)\right|}^2$ of the MBS ${\chi }^{\mathrm{in}}$ in Fig. 2.

Figure 6

Localization length (defined as the inverse participation ratio over half of the space, ${\xi }^{-1}={\int }_0^{L/2}dr{\left[r{\mathrm{\Phi }}^{†}\left(r\right)\mathrm{\Phi }\left(r\right)\right]}^2$) of the lowest eigenstate in $l=0$ sector as a function of the exchange energy, for $m=m_e$, $\mu =0$, ${\mathrm{\Delta }}_0=0.5$ meV, $\phi =0$, $r_0=0$, $R=25$ nm, and $n=2$. The curves in different colors correspond to different system sizes, from $L=1\mu \mathrm{m}$ to $L=7\mu \mathrm{m}$ in step of 200 nm. The dotted arrow is a guide to the eye, denoting the convergence point fitted as $\alpha \approx 0.5218$ meV.

Figure 7

Probability densities of the lowest four eigenstates in $l=7$ sector in Fig. 2. (a) A subgap state localized at the outer boundary, $r=pR$. (b),(c) Subgap states localized near $r=0$. (d) An extended state.

Figure 8

Charge and spin expectation values, shown in color as compared with the color bars aside, of the subgap states in Fig. 2 (upper four-panel block) and Fig. 2 (lower four-panel block). The color bar represents a range of $\langle -0.1,0.1\rangle$ for both charge and spin, with full polarization corresponding to $\pm 1$. (a) Mean charge. (b) Mean spin along $\widehat{z}$. (c) Mean spin projection along the local skyrmion spin $\widehat{N}$. (d) Spectrum for the doubled system size $L=2\mu \mathrm{m}$, with other parameters unchanged. In all these figures, the upper limit on the energy is a numerical artifact, due to a limited number of eigenstates returned by the diagonalization routine. The true spectrum extends to much higher energies.