Tunable skyrmion-skyrmion binding on the surface of a topological insulator

TL;DR

This paper demonstrates that skyrmions on the surface of a magnetic topological insulator can attract and form bound states, with the interaction tunable via chemical potential, magnetic field, and anisotropy, unlike in conventional ferromagnets.

Contribution

It introduces a mechanism for skyrmion binding based on electronic orbital hybridization and shows how this can be controlled in topological insulator systems.

Findings

Skyrmions can form bound states due to electronic hybridization.

Binding can be tuned by chemical potential, magnetic field, and anisotropy.

Phase diagram includes stable, unstable, and metastable skyrmion pairs.

Abstract

We show that skyrmions on the surface of a magnetic topological insulator may experience an attractive interaction that leads to the formation of a skyrmion-skyrmion bound state. This is in contrast to the case of skyrmions in a conventional chiral ferromagnet, for which the intrinsic interaction is repulsive. The origin of skyrmion binding in our model is the molecular hybridization of topologically protected electronic orbitals associated with each skyrmion. Attraction between the skyrmions can therefore be controlled by tuning a chemical potential that populates/depopulates the lowest-energy molecular orbital. We find that the skyrmion-skyrmion bound state can be made stable, unstable, or metastable depending on the chemical potential, magnetic field, and easy-axis anisotropy of the underlying ferromagnet, resulting in a rich phase diagram. Finally, we discuss the possibility to realize this effect in a recently synthesized Cr doped ${\left(\mathrm{Bi}_{2-y}\mathrm{Sb}_{y}\right)}_{2}\mathrm{Te}_3$ heterostructure.

Figures

Tables

Table 1. Table 1: Magnetic-free-energy parameters and single-skyrmion length scales for the numerical results used throughout this work. The single-skyrmion length scales, ξ 𝜉 \xi and R 𝑅 R , are determined from B 𝐵 B and K 𝐾 K through Eqs. ( 9 ) and ( 10 ). We choose B 𝐵 B and K 𝐾 K within the stable region of the skyrmion phase diagram [see Fig.  2 (a)] to tune the ratio ξ / R 𝜉 𝑅 \xi/R . We set λ 0 = R subscript 𝜆 0 𝑅 \lambda_{0}=R when calculating the electronic states to ensure that the skyrmions host discrete sets of orbitals.

Color	$\kappa \xi$	$\kappa R$	$B/B_s$	$K/B_s$
$\bullet$	0.73	2.85	0.1	0.9
$\bullet$	0.69	2.35	0.1	1.0
$\bullet$	0.66	1.94	0.1	1.1
$\bullet$	0.71	1.84	0.2	0.9
$\bullet$	0.81	2.14	0.25	0.64
$\bullet$	0.89	2.14	0.35	0.45

Equations

$$f$$

$$F_M$$

$$\kappa$$

$$B_s$$

$$\text{functionalderivative}{F}_M\mathbf{m}=$$

$$\begin{gathered} \mathbf{m}\text{quantity}(\mathbf{r})=\text{matrixquantity}(\sin \text{quantity}[\Theta \text{quantity}(r)]\cos \text{quantity}(W\varphi +{\varphi }_0) \\ \sin \text{quantity}[\Theta \text{quantity}(r)]\sin \text{quantity}(W\varphi +{\varphi }_0) \\ \cos \text{quantity}[\Theta \text{quantity}(r)]), \end{gathered}$$

$$\text{functionalderivative}{F}_M\Theta =0,$$

$$\Theta \text{quantity}(r)\propto {\mathcal{K}}_1\text{quantity}(\frac{r}{\xi }),$$

$$\xi =\frac{1}{\kappa }\sqrt{\frac{B_s}{B+2K}},$$

$$\Theta \text{quantity}(R)=\frac{\pi }{2}$$

$$\begin{gathered} H_{\mathrm{MTI}}=\hslash v\sum _{\mathbf{k}}{c}_{\mathbf{k}}^{†}{\text{quantity}(\mathbf{k}\times \mathbf{\sigma })}_z{c}_{\mathbf{k}} \\ -\Delta \int d^2\mathbf{r}{\psi }_{\mathbf{r}}^{†}\mathbf{m}\text{quantity}(\mathbf{r})\cdot \mathbf{\sigma }{\psi }_{\mathbf{r}}, \end{gathered}$$

$${\lambda }_0=\frac{\hslash v}{\Delta },$$

$$\displaystyle\matrixquantity(-{m}_{z}\quantity(r)&-\lambda_{0}\derivative{r}-\lambda_{0}\frac{j+\frac{1}{2}}{r}\\ \lambda_{0}\derivative{r}-\lambda_{0}\frac{j-\frac{1}{2}}{r}&{m}_{z}\quantity(r))\boldsymbol{\chi}_{j}\quantity(r)=\frac{E_{j}}{\Delta}\boldsymbol{\chi}_{j}\quantity(r).$$

$$\begin{gathered} {\mathbf{\chi }}_j\text{quantity}(r)\sim \text{matrixquantity}(\sqrt{1-\frac{E_j}{\Delta }}{\mathcal{K}}_{j-\frac{1}{2}}\text{quantity}(\frac{r}{{\lambda }_j}) \\ \sqrt{1+\frac{E_j}{\Delta }}{\mathcal{K}}_{j+\frac{1}{2}}\text{quantity}(\frac{r}{{\lambda }_j})), \end{gathered}$$

$${\lambda }_j=\frac{{\lambda }_0}{\sqrt{1-{\text{quantity}(\frac{E_j}{\Delta })}^2}}.$$

$$H^{\text{Mol}}=\sum _{i=1,2}{E}_0{c}_i^{†}{c}_i+t\text{quantity}(c_1^{†}{c}_2+c_2^{†}{c}_1),$$

$$t=\text{matrixelement}{\psi }^1{M}_1{\psi }^2=\text{matrixelement}{\psi }^2{M}_1{\psi }^1,$$

$$M_i=-\Delta \text{quantity}[{\mathbf{m}}^{\mathrm{sk}}\text{quantity}(\mathbf{r}-{\mathbf{r}}_i)-\hat{\mathbf{z}}]\cdot \mathbf{\sigma }$$

$$E_{\pm }=E_0\pm \text{absolutevalue}t\text{quantity}(x).$$

$$\Phi \text{quantity}(x)=$$

$$\mathcal{Z}=$$

$$E_{-}\text{quantity}(x_{\mu })\equiv \mu .$$

$$\Phi =\{\begin{array}{ll}{E}_{-}\text{quantity}(x)-\mu & x<x_{\mu }\\ 0& x>x_{\mu }\end{array}.$$

$$F(x)=F_M(x)+\Phi (x).$$

$$F_M^{\prime }(x_B)+E_{-}^{\prime }(x_B)=0.$$

$$\mu >E_0.$$

$$F(x_B)<F(x\to \infty ),$$

$$F(x_B)>F(x\to \infty ),$$

$$x_{\mu }<x_B,$$

$$F_M(x)\sim F_0{\mathcal{K}}_1\text{quantity}(\frac{x}{\xi }),$$

$$t(x)\sim t_0{\mathcal{K}}_0\text{quantity}(\frac{x}{\lambda }),$$

$$F_M(x)$$

Full text

I Introduction

The surface of a strong three-dimensional topological insulator hosts a set of two-dimensional surface states protected by topology, provided that the surface does not break the protecting symmetries.Hasan and Kane (2010); Bansil et al. (2016); Chiu et al. (2016) Within the bulk gap, these states are characterized by a chiral Dirac cone dispersion. When time-reversal symmetry is broken, these surface states are no longer protected and may be gapped. In a magnetic topological insulator (MTI), magnetic moments result in broken time-reversal symmetry. The magnetic moments may couple to each other through direct or indirect [Ruderman-Kittel-Kasuya-Yosida (RKKY)] exchange and form an ordered state.Liu et al. (2009); Zhu et al. (2011) These moments may also couple to the electronic subsystem through a Zeeman-like term proportional to the local magnetization. In the case of uniform out-of-plane ferromagnetic order, the magnetization gives the surface Dirac electrons a finite mass, resulting in a gap in the surface-state spectrum.

The gapped surface states can be effectively described by a massive Dirac model. Therefore, any sign change of the Dirac mass leads to localized Jackiw-Rebbi modes Tokura et al. (2019); Jackiw and Rebbi (1976); Jackiw (1984). Consequently, when the massive Dirac electrons are coupled to ferromagnetic moments there will be a set of protected one-dimensional edge states associated with the boundary between magnetic domains with opposite magnetization. Such edge states are responsible for the anomalous quantum Hall effect,Chang et al. (2013) and are thought to give rise to butterfly hysteresis in the magnetotransport properties of magnetic topological insulators.Nakajima et al. (2015); Tiwari et al. (2017) Magnetic skyrmions, topological defects in the magnetization, give rise to similar topologically protected electronic states at the skyrmion perimeter.Hurst et al. (2015)

Magnetic skyrmions are localized topologically stable configurations of the magnetization for chiral ferromagnets.Nagaosa and Tokura (2013); Rößler et al. (2006) We consider skyrmions in a planar magnetic system having an easy-axis anisotropy and in the presence of a finite applied magnetic field perpendicular to the surface (along $\mathbf{\hat{z}}$), both of which tend to stabilize the skyrmion phase.Rößler et al. (2006) Far from a skyrmion, the magnetization, $\mathbf{{m}}$, is uniformly aligned with the applied magnetic field, $\mathbf{{m}}\parallel\mathbf{\hat{z}}$. At the center of a skyrmion, the magnetization is anti-aligned with the applied magnetic field, $\mathbf{{m}}\parallel-\mathbf{\hat{z}}$. Across the perimeter of a skyrmion, moving radially outward, $m_{z}$ smoothly interpolates between these boundary conditions as the magnetization vector tilts in-plane, with the in-plane component of magnetization perpendicular to the to the radial direction $\hat{r}$ (for Bloch-type skyrmions) [Fig. 1​(a)]. The magnetization vector $\mathbf{{m}}\quantity(\mathbf{r})$ at position $\mathbf{r}$ wraps the unit sphere as $\mathbf{r}$ traverses the entire plane. This wrapping is a manifestation of the skyrmion’s topological nature; the magnetization associated with a skyrmion cannot be smoothly deformed to a uniform ferromagnetic state. This fact, coupled with a ferromagnetic exchange interaction that prohibits sharp changes in the magnetization, makes skyrmions exceptionally stable. This stability is present even if the skyrmion is higher in energy than the uniform ferromagnetic phase. In Ref. [Hurst et al., 2015], Hurst et al. consider a topological insulator surface coupled to a planar magnet that hosts a skyrmion. They find a discrete set of protected orbitals bound to the skyrmion [see, e.g., Fig. 1​(b)]. Similar electronic states have also been investigated in related systems.Ferreira and Loss (2013); Uchoa et al. (2015)

In this paper, we consider the interaction between a pair of skyrmions on the surface of a MTI. In the absence of the Dirac surface states, a pair of skyrmions in a chiral ferromagnet experiences a mutually repulsive interaction at all distances.Lin et al. (2013) This repulsive interaction decays exponentially with the inter-skyrmion separation over the magnetic healing length, $\xi$. The healing length is controlled by the applied magnetic field and easy-axis anisotropy. Once the Dirac electronic system is coupled to the magnetic system, skyrmion-bound electronic states form. Their wave functions overlap and hybridize when two skyrmions are close. If an electron occupies the lowest-lying hybridized electronic state, there will be an attractive contribution to the skyrmion-skyrmion interaction from the associated molecular binding energy. This attractive interaction also decays exponentially with separation and is governed independently by the skyrmion-bound orbital decay length. Occupation of the lowest-energy electronic orbital, and hence the attractive interaction, can be controlled by coupling the electronic system to a reservoir and tuning its chemical potential. Frustrated magnetic systems, with competing antiferromagnetic and ferromagnetic couplings, may show an oscillating repulsive/attractive interaction between skyrmions as a function of distance, even in the absence of an electronic system.Leonov and Mostovoy (2015); Kharkov et al. (2017) In this work, in contrast, our focus is on chiral ferromagnetic systems, where the magnetic interaction is purely repulsive.

A central result of this work is a zero-temperature skyrmion-skyrmion binding phase diagram, displaying the stability of bound skyrmion pairs as a function of the electrons’ chemical potential and magnetic healing length. The healing length may be controlled through the applied magnetic field and the easy-axis anisotropy. We also extend these results to low, but finite, temperature and establish concrete conditions to realize this effect experimentally. Understanding and controlling magnetic skyrmion binding through the electronic system may be important for both classical skyrmionic devicesFert et al. (2017) and for possible qubit implementations.Ferreira and Loss (2013)

The remainder of this paper is structured as follows. In Section II, we review the theory of individual magnetic skyrmions. Evaluating the magnetic free energy of a system containing two skyrmions as a function of skyrmion-skyrmion separation reveals a short-range repulsive interaction. In Section III, we consider the electronic subsystem. Topologically protected electron orbitals bound to single skyrmions hybridize to form molecular orbitals in a two-skyrmion system. We determine the grand potential of this system in contact with an electronic reservoir as a function of skyrmion-skyrmion separation and conclude that the electronic system gives rise to an attractive interaction. In Section IV, we analyze the total skyrmion-skyrmion interaction as a function of tunable parameters, finding a phase diagram for skyrmion-skyrmion binding at zero temperature, and then addressing the case of low, but finite, temperature. Finally, in Section V, we discuss the possibility to realize this effect in $\textrm{Cr}_{x}{(\textrm{Bi}_{1-y}\textrm{Sb}_{y})}_{2-x}\textrm{Te}_{3}/{(\textrm{Bi}_{1-y}\textrm{Sb}_{y})}_{2}\textrm{Te}_{3}$ and justify the approximations we have made in the context of this material.

II Short-range magnetic repulsion

In this section, we consider the magnetic subsystem in isolation to determine the range and character of the intrinsic repulsive skyrmion-skyrmion interaction. Starting from a standard free energy for a planar chiral magnet, with the addition of an easy-axis anisotropy, we determine the single-skyrmion magnetization. The skyrmion is then characterized by two length scales: the radius, $R$, and the healing length, $\xi$ [see Fig. 1​(a)]. Next, we numerically determine the magnetic free energy of a system of two skyrmions as a function of inter-skyrmion separation. This analysis reveals that the magnetic repulsive interaction decays exponentially with inter-skyrmion separation over the healing length, $\xi$. If $\xi$ is sufficiently short, the repulsion will be overcome by a longer-range attractive interaction due to the electronic subsystem.

The standard free energy for a two-dimensional chiral magnet stabilizing skyrmions includes the ferromagnetic exchange interaction, the Dzyaloshinskii-Moriya interaction (DMI), and a perpendicular applied magnetic field.Nagaosa and Tokura (2013) We consider this magnetic free energy with an additional term accounting for an easy-axis anisotropy, which may stabilize skyrmions in the absence of an applied magnetic field. The free energy density, $f$, and total magnetic free energy, $F_{M}$, are

[TABLE]

Here, $J>0$ is proportional to the exchange constant, $D$ arises from the DMI,Dzyaloshinsky (1958); Moriya (1960) $B$ is proportional to the out-of-plane applied magnetic field, and $K$ is the easy-axis anisotropy. The magnetic system has a natural inverse length scale $\kappa$, and a natural energy scale $B_{s}$:

[TABLE]

We neglect fluctuations in $\absolutevalue{\mathbf{{m}}}$, which may be penalized by terms quadratic and quartic in $\absolutevalue{\mathbf{{m}}}$ (not explicitly included here). Without loss of generality, we set $\absolutevalue{\mathbf{{m}}}=1$. For $D=0$ and $B>0$, the ground state of the magnetic system is uniform, with $\mathbf{{m}}=\mathbf{\hat{z}}$ throughout the plane. For finite $D$ and sufficiently large $B$ and/or $K$, skyrmions may exist as locally stable features in the magnetization.Rößler et al. (2006) The presence of a skyrmion in the ferromagnetic background may either increase or decrease the magnetic free energy depending on $B$ and $K$. A discussion of the phase diagram of chiral magnets governed by Eq. (1) is given in Refs. [Güngördü et al., 2016; Banerjee et al., 2014].

For sufficiently large $B$ and/or $K$, a magnetic skyrmion minimizes the magnetic free energy, i.e., it solves

[TABLE]

with the boundary conditions $\mathbf{{m}}\quantity(r=0)=-\mathbf{\hat{z}}$ and $\mathbf{{m}}{\quantity(r\rightarrow\infty)}=\mathbf{\hat{z}}$ for positive $B$. These boundary conditions ensure that the skyrmion is centered at $r=0$, and that the magnetization far from the skyrmion approaches the value it would have in the uniform ferromagnetic phase. The magnetization profile of the skyrmion may be parameterized as

[TABLE]

where $W$ is the winding number of the skyrmion, and $\phi_{0}$ determines its chirality. For $D>0$, the specific form of the DMI considered in Eq. (1) stabilizes spiral-like (Bloch) skyrmions, with a definite chirality, $\phi_{0}=+\ifrac{\pi}{2}$,111More generally, the DMI for a two-dimensional thin film can be written as $f_{DM}=D_{1}\mathbf{m}\cdot\left(\nabla\times\mathbf{m}\right)+D_{2}\left[\left(\mathbf{m}\cdot\nabla\right)m_{z}-m_{z}\left(\nabla\cdot\mathbf{m}\right)\right]$. For concrete calculations in this manuscript, we have chosen $D_{1}=D>0$, $D_{2}=0$, leading to stable Bloch-type skyrmions with $\phi_{0}=+\pi/2$. In the more general case ($D_{1}\neq 0,\,D_{2}\neq 0$), some other (fixed) value of $\phi_{0}$ will be stabilized. For any value of $\phi_{0}$, the qualitative arguments leading to magnetic repulsion hold [see the discussion following Eq. (10), below]. Because the electronic states are topological in origin, we expect their presence to be robust to changes in the in-plane magnetization, so the electronic interaction is also robust to changes in the type of skyrmion (value of $\phi_{0}$), although the wavefunctions may be subject to a more complicated wave equation, Eq. (13). We therefore expect our qualitative analysis of skyrmion binding to hold for skyrmions of Néel-type ($\phi_{0}=0$), Bloch-type ($\phi_{0}=\pi/2$), and intermediate type [$\phi_{0}\in(0,\pi/2)$]. and a definite winding number, $W=+1$ (see, e.g., Ref. Han, 2017).222This is in contrast with the case of a centrosymmetric frustrated magnetic system, not considered here, where skyrmions of either winding number, $W=\pm 1$, may be stable.Lin and Hayami (2016) An individual skyrmion satisfies the equation

[TABLE]

with the boundary condition $\Theta\,\quantity(r=0)=\pi$, $\Theta\quantity(r\rightarrow\infty)=0$. Away from the skyrmion core, for $\Theta\ll 1$ and $\kappa r\gg 1$, Eq. (7) is solved asymptotically by

[TABLE]

where $\mathcal{K}_{n}$ is the modified Bessel function of the second kind. The healing length,

[TABLE]

sets the scale for decay of the in-plane magnetization far from the skyrmion. The asymptotic solution given in Eq. (8) can be further simplified to $\Theta\quantity(r)\propto\sqrt{\ifrac{\xi}{r}}\exp(-\ifrac{r}{\xi})$ in the limit $r\gg\xi$. The skyrmion radius, $R$, may be found numerically by inverting

[TABLE]

for $\Theta\quantity(r)$ satisfying Eq. (7). At the skyrmion radius, $r=R$, the magnetization is entirely in-plane. The magnetic texture of an isolated skyrmion is plotted in Fig. 1​(a).

Although Eq. (5) has a skyrmion solution for any finite $B$ or $K$, this solution does not describe a minimum in the free energy for sufficiently small $B$ and $K$. Below a critical line in the $B,K$ plane, skyrmions are unstable to a transition towards a spin-spiral phase.Han (2017) We numerically determine this phase boundary, plotted in Fig. 2(a) (see also Table 1) through the procedure described in Appendices A and  B. Above this phase boundary, localized skyrmions exist as stable minima of the magnetic free energy. A more extensive review of the theory of magnetic skyrmions is given in Refs. [Nagaosa and Tokura, 2013; Han, 2017].

We now consider a system of two skyrmions separated by a distance $x$. Our goal is to establish that the intrinsic inter-skyrmion interaction is repulsive and short-range, with length scale $\xi$. This form of the intrinsic interaction is expected based on the following reasoning. Far from the skyrmion core, the magnetization is mostly aligned with the applied field save for a small in-plane component. Two skyrmions repel each other because they favor opposing in-plane magnetization in the region between them.Lin et al. (2013) The magnitude of this in-plane component is determined by the length scale $\xi$. We therefore expect a short-range repulsive interaction between skyrmions with length scale $\xi$.

To verify this picture, we treat this problem numerically on a lattice in a similar manner to Ref. [Lin et al., 2013]. We prepare an initial magnetic configuration with two skyrmions separated by a specified distance. We then freeze a magnetic moment within each skyrmion core and minimize the magnetic free energy with respect to variations in the remaining moments. Given this relaxed magnetization, we calculate the magnetic free energy and establish the final inter-skyrmion separation, $x$ for the relaxed configuration. By repeating this calculation with different initial inter-skyrmion separations, we determine the magnetic free energy as a function of $x$. Further details are given in Appendix A. The results of this calculation are given in Fig. 3. We find that the intrinsic magnetic inter-skyrmion interaction is well fit by a modified Bessel function (of first order and second kind) with the separation, $x$, scaled by $\xi$, consistent with Eq. (8). This agreement holds over several orders of magnitude of interaction strength. For sufficiently large inter-skyrmion separations, errors in the numerical calculation of $F_{M}(x)$ become comparable to the interaction energy [$F_{M}(x)-F_{M}(\infty)$], and we can no longer make numerical predictions. These results are entirely consistent with the analogous calculation presented in Ref. [Lin et al., 2013], where it was shown that the magnetic healing length determines the range of the repulsive skyrmion-skyrmion interaction for $K=0$.

In the analysis given above, we have neglected the long-range magnetic dipolar interaction. This can be justified for skyrmions with a sufficiently small radius $R$ (reducing the total moment associated with each skyrmion) and for a moderate separation $x$. In Sec. V, we provide a more detailed justification of this approximation in the context of a candidate material, $\textrm{Cr}_{x}{(\textrm{Bi}_{1-y}\textrm{Sb}_{y})}_{2-x}\textrm{Te}_{3}/{(\textrm{Bi}_{1-y}\textrm{Sb}_{y})}_{2}\textrm{Te}_{3}$.

III Attractive interaction

This section addresses the attractive interaction effected by the electronic system. First, we consider the electronic structure of the MTI surface resulting from the hybridization of single-skyrmion orbitals. Next, we integrate out the electronic subsystem, accounting for an electronic reservoir at fixed chemical potential, to determine the grand potential as a function of skyrmion-skyrmion separation. This analysis leads us to conclude that the MTI surface states can give rise to an attractive skyrmion-skyrmion interaction.

In the absence of a magnetic subsystem, a three-dimensional topological insulator surface is characterized by an odd number of spin-orbit-coupled Dirac cones. We consider a single Dirac cone centered at $\mathbf{k}=0$, which is coupled to the magnetic system through a Zeeman-like term proportional to the local magnetization:

[TABLE]

where $c_{\bf k}$ is a spinor of electronic annihilation operators in momentum and $\psi_{\bf r}={1\over\sqrt{A}}\sum_{\bf k}\exp(i{\bf k\cdot r})c_{\bf k}$ is its real space counterpart defined with the system’s area $A$. Furthermore, $v$ is the Fermi velocity, and $\Delta$ is proportional to the exchange interaction between the surface electrons and the magnetic system. These parameters imply a natural length scale:

[TABLE]

and a natural energy scale, $\Delta$.

The MTI surface Hamiltonian, Eq. (11), neglects the magnetic vector potential. This is justified if the $U(1)$ phase acquired by electrons in the skyrmion-bound orbitals is small, i.e., the flux through the skyrmion area is smaller than the flux quantum. This condition can be written as $R^{2}\tilde{B}\ll\Phi_{0}$ where $\tilde{B}\propto B$ is the dimensionful applied magnetic field. For a skyrmion of radius $R=50$ nm, this condition is satisfied for $\tilde{B}\ll 1$ T. Equation (11) also neglects the direct Zeeman coupling of the surface electrons to the applied magnetic field. This is justified if the Zeeman-like exchange coupling to the magnetization is large compared to the Zeeman energy. In Sec. V, we show that this limit is experimentally relevant.

For a uniform magnetization and finite $\Delta$, ${m}_{z}$ gives the surface Dirac electrons a finite mass, while the in-plane components simply shift the Dirac point from $\mathbf{k}=0$. The resulting massive Dirac cone is characterized by a finite Berry curvature, with sign determined by the sign of $m_{z}$. Thus, gapless states will be found where $m_{z}$ changes sign.Jackiw and Rebbi (1976); Jackiw (1984) These states are given below as the eigenstates of Eq. (11) with a magnetization given by the solution of Eq. (7). Since the skyrmion magnetization has rotational symmetry, the energy eigenstates are a product of angular and radial functions:

[TABLE]

where $j$ is the half-integer total-angular-momentum quantum number, and $\boldsymbol{\chi}_{j}\quantity(r)=\quantity(\chi_{j}^{\uparrow}\quantity(r),\chi_{j}^{\downarrow}\quantity(r))^{T}$ is the radial wave function. If the in-plane component of the magnetization is divergenceless, as in the case of Bloch skyrmions, we may work in a gauge where it does not enter the radial wave equation:Hurst et al. (2015)

[TABLE]

For ${r-R}{}\gg\xi$, the magnetization approaches $\mathbf{\hat{z}}$ and the asymptotic behavior of the bound-state wavefunction is

[TABLE]

where

[TABLE]

For $R\gtrsim\lambda_{0}$ there is a discrete set of states within the Dirac mass gap with $\bra{\mathbf{r}}\ket{\psi_{j}}$ localized to the skyrmion. Away from the skyrmion, these wavefunctions decay with length scale $\lambda_{j}$. Hurst et al. (Ref. Hurst et al., 2015) provide an explicit solution for the wavefunctions in the limit $\xi=0$ as well as a transcendental equation for the electronic spectrum as a function of $R$ (similar solutions are also given for related models in Refs. Ferreira and Loss, 2013; Uchoa et al., 2015). We determine the bound states exactly on a lattice for the full single-skyrmion magnetization, $\mathbf{{m}}^{\mathrm{sk}}\quantity(\mathbf{r})$. The magnetization profile, $\mathbf{{m}}^{\mathrm{sk}}\quantity(\mathbf{r})$, is determined from the solution to Eq. (7) using the approach discussed in Appendix C. In Fig. 4, we plot the spectrum from Ref. [Hurst et al., 2015] as a function of $R$ (black lines) together with the numerically determined spectrum for a skyrmion with $\xi=0.256\,R$ (blue dots).

In a system hosting two skyrmions, the single-skyrmion orbitals will overlap and hybridize to form molecular orbitals. We focus on the lowest-energy electronic state with $j=1/2$ in each skyrmion and construct a two-level Hamiltonian to describe the orbital hybridization:

[TABLE]

where $E_{0}$ is the energy of the lowest-energy single-skyrmion orbital and $c_{i}$ annihilates an electron in the lowest-energy single-skyrmion orbital of skyrmion $i$. The tunneling amplitude is given by

[TABLE]

where $\ket{\psi^{i}}$ are the lowest-energy ($j=\ifrac{1}{2}$ for $R=\lambda_{0}$ and $v>0$) single-skyrmion orbitals centered at $\mathbf{r}_{i}$, defined through Eq. (III), and

[TABLE]

is the ‘atomic potential’ of skyrmion $i$ relative to the ferromagnetic background. The tunneling amplitude, $t\quantity(x)$, decays exponentially for ${x-2R}{}\gg\lambda$ where $\lambda\equiv\lambda_{j}$ for $j$ corresponding to the lowest-energy single-skyrmion orbital. This parametric dependence comes from the asymptotic behavior of the single-skyrmion orbitals defined in Eq. (14).

Equation (16) is diagonalized by molecular orbitals with energies

[TABLE]

We neglect the doubly-occupied state for a skyrmion pair at a sufficiently small separation $x$ and for a weakly-screened long-range Coulomb interaction. In Sec. V, we justify this approximation in the context of recently synthesized Cr doped ${\quantity(\mathrm{Bi}_{2-y}\mathrm{Sb}_{y})}_{2}\mathrm{Te}_{3}$ heterostructures.

Figure 5 presents the binding energy as a function of $x$ for skyrmions with radius $R=\lambda_{0}$. The binding energy is the difference between the lowest-lying molecular orbital, $E_{-}\quantity(x)$, and the single-skyrmion orbital energy, $E_{0}=E_{-}\quantity(\infty)$. The binding energy indeed decays with length scale $\lambda$, and is well fit by $t_{0}\mathcal{K}_{0}\quantity(\frac{x}{\lambda})$, as expected from Eq. (17), and the asymptotic single-skyrmion wavefunctions, given in Eq. (14).

The effective interaction between a skyrmion pair is determined by integrating out the electronic degrees of freedom. The electronic state is determined through the grand potential for a fixed separation $x$:

[TABLE]

Where we set Boltzman’s constant to 1 and assume that the electronic subsystem is in contact with a reservoir at temperature $T$ and chemical potential $\mu$. The molecule is restricted to support $N=0,1$ electrons, consistent with the discussion following Eq. (19).

The $T=0$ behavior of $\Phi$ is sketched in Fig. 6. For $\mu>E_{0}$, the lowest-energy molecular orbital is occupied at all $x$ and the grand potential is $\Phi=E_{-}\quantity(x)-\mu$. For $\mu<E_{0}$, the lowest-energy molecular orbital crosses the chemical potential at a separation $x_{\mu}$, above which it becomes unoccupied. $x_{\mu}$ is defined by

[TABLE]

For a monotonically increasing energy $E_{-}\quantity(x)$ at $T=0$, the grand potential is therefore:

[TABLE]

Since $-\Phi^{\prime}(x)<0$, at $T=0$, the electronic system effects an attractive force for $x<x_{\mu}$ or when $\mu>E_{0}$. If this attractive interaction overcomes the short-range repulsive magnetic interaction discussed in Sec. II, a stable bound skyrmion molecule will form.

IV Pair-binding phase diagram

The total free energy is given by the sum of the contributions from the electronic system and the magnetic system,

[TABLE]

Based on the possible functional forms of $F(x)$ we identify four regimes whose nature and boundaries will be discussed below. The phase diagram is depicted in Fig. 7(a) and a representative free-energy plot for each phase is shown in Figs. 7(b)-(e). We denote the regimes of our model as:

(i) the ‘stable’ regime, where the free energy has a unique global minimum at finite skyrmion-skyrmion separation. We denote this separation $x_{B}$ and refer to it as the bond length. The stable regime is colored green in Fig. 7.

(ii) the ‘stable-hysteretic’ regime, where, in addition to the global minimum at $x_{B}$, there is a local minimum as $x\to\infty$. This regime is colored in beige in the Fig. 7.

(iii) the ‘metastable-hysteretic’ regime, where there is a local minimum at $x_{B}$ and a global minimum at infinite separation. This phase is blue in Fig. 7.

(iv) the ‘unstable’ regime, where the free energy is a monotonically decreasing function of separation, and there is no bound state. This phase is colored white in Fig. 7.

IV.1 Zero temperature

At zero temperature, the electronic free energy is defined piecewise through $\Phi=\min\left[E_{-}(x)-\mu,0\right]$, as shown in Eq. (23) and Fig. 6(b). For any pair of $\mu$, $\xi$, the phase at zero temperature is determined from the balance of the repulsive (magnetic) and attractive (electronic) forces at a finite separation. We recall that the magnetic contribution to the free energy, $F_{M}$, is always repulsive and decays quickly over a length scale $\xi$, while the electronic contribution is attractive whenever the lower orbital state is populated and is governed by the length scale $\lambda$. Moreover, we note that, for a bound state to form, the repulsion should be shorter-range than the attraction, $\xi<\lambda$. We therefore scan $\xi<\lambda$ and find the phase boundaries for each $\xi$ as a function of the chemical potential, $\mu$.

(i) The stable phase—if $\mu$ is above the lowest molecular orbital energy, $E_{-}$, for all separations $x$, i.e., $\mu>E_{0}$, then the orbital is always occupied and $\Phi(x)=E_{-}(x)-\mu$ for all $x$. We therefore have both repulsion and attraction everywhere and a balance of forces occurs at $x_{B}$ satisfying

[TABLE]

This stable phase is characterized by a free energy with a unique minimum at $x_{B}$. In Fig. 7, the regime is bounded from below by

[TABLE]

(ii) The stable-hysteretic regime—in addition to the stable regime, there are two other regimes in which a bound state occurs. As discussed above, when the chemical potential crosses $E_{-}(x)$ at some separation $x_{\mu}$, there is an attractive force for $x<x_{\mu}$ but no attraction for $x>x_{\mu}$. If the forces balance at a separation where the lower-energy orbital is populated, i.e., $x_{B}<x_{\mu}$, there will be a skyrmion bound state. The stable-hysteretic regime and the stable regime are distinguished by the behavior of the free energy at large separations. Because there is no attractive interaction for $x>x_{\mu}$ in the stable-hysteretic regime, there is a free-energy minimum as $x\to\infty$. Between the to minima there is a cusp-like energy barrier at $x_{\mu}$ as depicted in Fig. 7(c). In the stable-hysteretic phase the minimum at $x_{B}$ arises at a lower free energy than the minimum at $x\to\infty$,

[TABLE]

so the bound state is stable. However, two skyrmions prepared at large separation experience a repulsive interaction and so remain unbound. This is the hysteretic behavior alluded to in the name of the phase.

(iii) The metastable-hysteretic regime—the metastable-hysteretic regime is distinguished from the stable-hysteretic regime by the relative energies of the bound and unbound configurations. In this regime the unbound configuration globally minimizes the free energy,

[TABLE]

but the molecular orbital is occupied at the bond length, $x_{B}<x_{\mu}$. The free-energy minimum at $x_{B}$ is local and the bound state is therefore metastable. The lower bound of this phase is found when the molecular orbital is depopulated at $x_{B}$, $x_{B}=x_{\mu}$, and the attraction can no longer overcome the repulsion at any separation.

(iv) The unstable regime—in this regime, the repulsive interaction dominates the free energy for all separations. In this case the only free-energy minimum occurs as $x\to\infty$. This occurs when the molecular orbital is unoccupied at $x_{B}$, as defined in Eq. (25), so $x_{B}$ is not a minimum of the free energy. The unstable phase is defined by the inequality:

[TABLE]

and the phase boundary is found by comparing the two lengths.

In Fig. 7(a), we estimate the phase boundaries as described above using the asymptotic behaviour of $F_{M}(x)$ and $\Phi(x)$. In Fig. 3, we show that for $x\gg\xi,2R$, the magnetic free energy, $F_{M}(x)$ has the asymptotic form

[TABLE]

where $F_{0}$ characterizes the strength of the repulsive interaction. Similarly, in Fig. 5, we show that the tunneling amplitude, $t(x)$ at $x\gg\lambda$, has the asymptotic form

[TABLE]

where $t_{0}$ characterizes the strength of the tunnel splitting.

We use the above asymptotic forms to find the bond length $x_{B}$, the scale $x_{\mu}$, and the free energy at the minimum $F(x_{B})$ and use these quantities to draw the phase diagram in Fig. 7. Moreover, it is possible to approximate the functional form of the free energy further by taking the asymptotic limit of the above Bessel functions:

[TABLE]

Balancing the derivatives of the two expressions above and neglecting subleading terms in $\lambda/x$ and $\xi/x$, we obtain:

[TABLE]

This limit is consistent with our assumption of $x\gg\lambda>\xi$ when

[TABLE]

Equation (34) shows that the bond length becomes shorter as the ratio $\xi/\lambda$ is decreased. The minimum of the free energy can be found by substituting Eq. (34) into the approximated form of $F(x)$ using Eqs. (30), (31) or Eqs. (32), (33). For $x\gg\xi,\lambda$, $E_{-}(x)$ and $F_{M}(x)$ are approximately exponential giving

[TABLE]

This helps simplify the expression for $F(x_{B})$ as the forces [$-F_{M}^{\prime}(x)$, $-E_{-}^{\prime}(x)$] are equal in magnitude and opposite in sign at the minimum $x_{B}$. We may therefore write

[TABLE]

which demonstrates that, as expected, the bound-configuration minimum is deeper for smaller ratios of $\xi/\lambda$.

One should note that these estimates are based on the asymptotic behavior of $t(x)$ and $F_{M}(x)$ and hence are accurate only in the limits mentioned above. However, we stress that a bound state may exist well beyond these limits. For example, in the limit where $\xi\ll\lambda$ we may view the magnetic repulsion as a hard-shell repulsion. Therefore, the bound state would form at $x=2R$ where the repulsion drops to zero while the attractive force is given by $-E_{-}^{\prime}(2R)\neq 0$.

IV.2 Finite temperature

At finite temperature the free energy changes due to thermal occupation of the molecular orbital states as well as thermal fluctuations in the skyrmion separation $x$ and possible changes in the magnetic free energy. In this section we include the temperature only through its effect on the population of the molecular orbitals. Beginning in the stable regime at zero temperature, we determine the conditions for which there may still be a bound state when the temperature is raised. This may be addressed qualitatively by plotting $F(x)$ as defined in Eq. (24), with the finite-temperature $\Phi$ given by Eq. (20), see, e.g., Fig. 8.

We assess the effect of temperature on the bound state in the following way. We assume a large chemical potential such that one of the molecular orbitals is always occupied, $\mu-E_{\pm}\gg T$ (but we still neglect double occupancy). At finite temperature, the average occupations of the states with energies $E_{-}$ and $E_{+}$ shift away from the zero-temperature limits of 1 and 0. The electronic free energy, $\Phi(x)=\bar{E}-TS$, has contributions from both the entropy $S$ and the average energy $\bar{E}$. When $T\gg E_{+}-E_{-}$, the populations of the two states are comparable, giving $S\to\ln 2$, and the average energy also becomes independent of $x$: $\bar{E}\to\left(E_{+}(x)+E_{-}(x)\right)/2=E_{0}$. Since $\Phi$ becomes $x$-independent, the attractive force is suppressed as the temperature is raised. We therefore define the temperature scale $T^{*}$, controlled by the typical molecular energy scale at the free-energy minimum:

[TABLE]

For $T\ll T^{*}$, the population of the $E_{-}$ state is significantly larger than that of $E_{+}$, the zero-temperature analysis applies, and the bound state persists. Above this temperature, the electronic contribution to the free energy is suppressed. For $T>T^{*}$ the suppression of the attractive force is linear in $t(x)/T$ and therefore one may find a free-energy minimum even above $T^{*}$. Since one of the molecular states is always occupied under the conditions described above, we can write the partition function

[TABLE]

and consequently find the attractive force,

[TABLE]

At high temperature, the attractive force $f(x)$ is suppressed, and decays exponentially at half the length scale ($\lambda/2$) relative to the decay length of the low-temperature attractive force ($\lambda$):

[TABLE]

Although it is suppressed by temperature, the above force may still overcome the magnetic repulsion at large separations provided $\xi<\lambda/2$. For $\xi\ll\lambda/2$, we expect to see a free-energy minimum at $T>T^{*}$ which becomes shallow as the temperature is increased. For $\lambda>\xi\gtrsim\lambda/2$, we expect the minimum to vanish above $T^{*}$. The two types of behavior can be seen in Fig. 8.

V Potential realization

A promising candidate material to host skyrmions that may realize this effect is the topological insulator heterostructure $\textrm{Cr}_{x}{(\textrm{Bi}_{1-y}\textrm{Sb}_{y})}_{2-x}\textrm{Te}_{3}/{(\textrm{Bi}_{1-y}\textrm{Sb}_{y})}_{2}\textrm{Te}_{3}$. In Ref. [Yasuda et al., 2016], Yasuda, et al. argue that this heterostructure hosts a skyrmion ground state for $\tilde{B}\lesssim 0.1$ T at $T\lesssim 10$ K. Band structure calculations suggest that this material has an electronic gap of $\Delta\approx 10$ meV and a Fermi velocity of $\hbar v\approx 500$ meV$\cdot$nm.Yasuda et al. (2016) This suggests that skyrmions of radius $R\gtrsim 50\,\mathrm{nm}$ would host skyrmion bound states in this material. Accurate characterization of the actual skyrmion size, e.g., by magnetic atomic force microscopy, would be necessary to ensure that this condition is fulfilled. For a tunnel coupling $t\quantity(2R)=1\,\mathrm{meV}$, and skyrmion-skyrmion binding energy $\simeq 0.1\,\mathrm{meV}$, skyrmion-skyrmion bound states should form when $T\lesssim 1$ K.

Throughout this paper, we have neglected the repulsive dipolar interaction between skyrmions as well as the magnetic vector potential. We have also assumed that double-occupancy of the molecular orbitals is suppressed by a weakly-screened Coulomb interaction. Here we argue that these approximations are justified for the heterostructure considered in Ref. [Yasuda et al., 2016] for skyrmions of radius $R=\lambda_{0}=50$ nm and with $t\quantity(2R)=1$ meV. Such skyrmions will have $E_{j}/\Delta\approx-0.6$ for their lowest-energy single-skyrmion orbitals (see Fig. 4), and will therefore have an orbital decay length, defined through Eq. (15), $\lambda\approx 60$ nm.

The repulsive dipolar interaction between skyrmions will dominate both the short-range repulsive interaction considered in Sec. II and the short-range attractive interaction examined in Sec. III at large separations. We approximate the dipolar interaction energy as $E_{D}\approx\frac{\mu_{0}\absolutevalue{\mathbf{d}}^{2}}{4\pi x^{3}}$ where $\mathbf{d}\approx\frac{\pi\alpha gsR^{2}l}{V}\mu_{B}\mathbf{\hat{z}}$ is the dipole moment due to the core of a skyrmion. We estimate $\alpha\approx 0.6$ as the number of dopants per unit cell,Yasuda et al. (2016) $V\approx 0.59$ nm3 as the unit cell volume,Jenkins et al. (1972) $l\approx 2$ nm as the depth of the magnetic layer,Yasuda et al. (2016) and $gs\approx 1$ as the magnitude of the moment of the dopant atoms. The force due to the dipolar interaction balances the force due to the attractive interaction at $x_{D}$ satisfying $\Phi^{\prime}\quantity(x_{D})=E_{D}^{\prime}\quantity(x_{D})$, with $\Phi\quantity(x)$ defined in Eq. (20). Under the above assumptions, this equality is satisfied for $x_{D}\approx 450$ nm at $T=20$ mK. Below this separation, the repulsive dipolar force may be safely neglected.

In Eq. (11), we have neglected the contribution of the magnetic vector potential by arguing that the $U(1)$ phase acquired by an electron in a skyrmion-bound orbital of radius $R=50$ nm is small for $\tilde{B}\ll 1$ T. This is easily satisfied for $\tilde{B}\lesssim 0.1$ T, the field strength at which the skyrmion crystal reported in Ref. [Yasuda et al., 2016] is stable.

Finally, we assumed that the molecular orbitals support only $N=0$ or $1$ electrons because of a large charging energy. The doubly-occupied state will be irrelevant if the chemical potential is well below the sum of the excited molecular orbital energy and the charging energy, i.e., if $E_{+}\quantity(x)+U\quantity(x)-\mu\gg T$, where $U\quantity(x)$ is the charging energy. This is naturally satisfied at all separations if $E_{0}-\mu\gg T$. Thus, the majority of the phase diagram, Fig. 7(a), is valid even for a vanishing Coulomb interaction, provided that the temperature is small compared to the binding energy. For a weakly-screened Coulomb interaction, $U\quantity(x)=\frac{e^{2}}{4\pi\epsilon}\frac{\exp(-\ifrac{x}{\lambda_{D}})}{x}$, where $\lambda_{D}$ is the Debye length and $\epsilon$ is the permittivity, the doubly-occupied state may be neglected for separations smaller than $x_{C}$ satisfying $E_{+}\quantity(x_{C})+U\quantity(x_{C})-\mu=T$. Assuming $\mu=E_{0}$, $\lambda_{D}=500$ nm, $\epsilon=\epsilon_{0}$ (i.e., $\frac{e^{2}}{4\pi\epsilon}=1\,e$V-nm), and $T=20$ mK, we find $x_{C}\approx 2.5$ $\mu$m.

Based on the length scales estimated above, our analysis will be applicable provided that the skyrmions are confined to a sample of size $L<\mathrm{min}\left(x_{C},x_{D}\right)\approx 450$ nm. Larger samples may still host stable bound states, but the neglected effects will dominate the skyrmion-skyrmion interaction at large separations.

In addition to the promising MTI candidate material studied in Ref. Yasuda et al., 2016, there may also be other systems where it is possible to realize this effect, e.g. a TI device functionalized with a magnetic top layer,Gong et al. (2017); Huang et al. (2017) or even in a 2-dimensional topological insulator (2DTI) setting. This can be realized in a 2DTI where there are two massive Dirac cones of unequal gap sizes (such as the Haldane modelHaldane (1988) with sublattice asymmetrySemenoff (1984)). If the coupling to the magnetic system changes the mass of the two valleys in the same way, there would be a regime in which a sign change in the magnetization amounts to a sign change in the Dirac mass of only one valley. In the presence of skyrmions this would lead to electronic bound states of the kind discussed here.

The DMI assumed here is essential to stabilize skyrmions. This interaction is allowed in the bulk only for non-centrosymmetric materials. However, an interface may also break inversion symmetry leading to a robust DMI (see, e.g., Ref. Soumyanarayanan et al., 2016 for a review). A magnetic thin film on top of a topological insulator may then naturally lead to the required DMI at the interface,Zarzuela and Tserkovnyak (2017); Zarzuela et al. (2018) without the need to specialize to non-centrosymmetric materials.

VI Summary and conclusion

In summary, we have shown that a pair of skyrmions on the surface of a magnetic topological insulator may experience a mutually attractive interaction. In the absence of the topological insulator’s Dirac surface electrons, the skyrmions will experience a repulsive interaction. The contribution of the electronic system may be switched on and off by tuning the chemical potential. For large chemical potential, the skyrmions form a bound state, while for low chemical potential they remain unbound. For intermediate chemical potential, both the bound and unbound states are locally stable. A chemical-potential sweep from low to high and back should therefore lead to hysteretic binding and unbinding of the skyrmion pair. This hysteretic binding may find use in skyrmionic devices,Fert et al. (2017) either as means of information storage, or as means to couple electronic and skyrmionic degrees of freedom. It may also be possible to generate GHz spin waves without microwave-frequency magnetic fields via an AC gate voltage. Such a magnon source could allow direct coupling of conventional and spin-wave circuits. The conditions laid out here for skyrmion binding could furthermore be used to realize stable qubits from the resulting molecular states.Ferreira and Loss (2013) Beyond the predictions made here, the hybridization of skyrmion-bound orbitals may lead to novel magnetoelectric effects within the skyrmion crystal phase, and the contribution of the electronic free energy may modify the skyrmion crystal phase boundary.

There is significant experimental effort in growing and characterizing magnetic topological insulators capable of supporting quantized anomalous Hall conductance and related domain-wall-bound-state phenomena. Current state-of-the-art Cr doped ${\quantity(\mathrm{Bi}_{2-y}\mathrm{Sb}_{y})}_{2}\mathrm{Te}_{3}$ heterostructures show signs of a skyrmion crystal phase,Yasuda et al. (2016); Liu et al. (2017); Jiang et al. (2019) while the recent discovery of easy-axis ferromagnetic van der Waals materials may lead to a new class of magnetic topological insulators.Gong et al. (2017); Huang et al. (2017) Magnetic topological insulators proximity-coupled to a ferromagnetic thin film have already been shown to host skyrmions.Zhang et al. (2018) Observing the effect considered in this work may be a natural next step toward demonstrating that skyrmion physics on the surface of a topological insulator leads to new and exciting effects.

Acknowledgements.

This work was enabled in part by support provided by NSERC, FRQNT, INTRIQ, CIFAR, Nordea Fonden, the FRQNT doctoral scholarship, and Compute Canada.

Appendix A Numerical solution of Eq. (5)

In the main text, we introduce a phenomenological form of $F_{M}\quantity(x)$, the magnetic free energy for two skyrmions separated by distance $x$. Since we have no exact analytical solution for the free energy, $F_{M}(x)$, or generally Eq. (5) with arbitrary boundary conditions, we must numerically determine $\mathbf{{m}}\quantity(\mathbf{r})$ and $F_{M}\quantity(x)$ to justify our phenomenological model. We outline our numerical approach to minimize the magnetic free energy here.

To minimize $F_{M}$, we evolve $\mathbf{{m}}$ with the partial differential equation

[TABLE]

where $\tau$ is the simulation time and

[TABLE]

Eq. (43) is simply the component of the Landau-Lifschitz-Gilbert equationGilbert (2004) that leads to relaxation. In the long-time limit, when $\partialderivative{\mathbf{{m}}}{\tau}=0$, the magnetization has been evolved to a configuration that minimizes $F_{M}\quantity[\mathbf{{m}}]$ under the constraint $\absolutevalue{\mathbf{{m}}}=1$. We implement boundary conditions by setting $\mathbf{H}\quantity(\mathbf{r})=0$ for $\mathbf{r}$ at boundaries and choosing the correct initial condition.

To find $F_{M}\quantity(x)$, we initialize the magnetization in an approximate configuration in the two-skyrmion topological sector that should be close to the ground state. We separate the magnetization into a left and a right region and use the numerical solution of Eq. (7) to prepare a single skyrmion centered at $-\ifrac{x^{\prime}}{2}$ in the left region. Similarly, we prepare a skyrmion at $\ifrac{x^{\prime}}{2}$ in the right region. Then we fix $\mathbf{H}\quantity(\pm\ifrac{x^{\prime}\hat{\mathbf{x}}}{2})=0$ to pin the skyrmion cores and evolve the magnetization under Eq. (43) until dynamics cease. From the resulting magnetization we can numerically evaluate Eq. (2) to find $F_{M}$. However, the pinned skyrmion-skyrmion separation, $x^{\prime}$, is not the distance between the skyrmion centers. Since the skyrmions repel, they will move away from each other under Eq. (43) until the pinned magnetic moments reach the skyrmion perimeters. We determine the location of the skyrmion centers, $\mathbf{r}_{\pm}$, by fitting the final $m_{z}$ with the ansatz

[TABLE]

where $m_{z}=\cos\Theta$. This fitting procedure is justified in the limit $\ifrac{\quantity(x-2R)}{\xi}\gg 1$.

Appendix B Numerical determination of single-skyrmion instability phase boundary

In the Sec. II, we argue that the magnetic free energy is stationary for a skyrmion texture, i.e. Eq. (5) is solved by Eq. (18), with $\Theta$ given by the solution of Eq. (7). While this is true, for sufficiently low $\quantity(B,K)$, the single-skyrmion texture does not minimize $F_{M}$. It is instead a saddle-point solution. In this region of parameter space, skyrmions are unstable to transitions towards a spin-spiral phase. This phase boundary is well known in the literature,Banerjee et al. (2014) and is identified by our numerical simulations. For $\quantity(B,K)$ below the instability phase boundary, we find that Eq. (43), which minimizes $F_{M}$ under the constraint $\absolutevalue{\mathbf{{m}}}=1$, takes an initial single-skyrmion configurations to a final spiral configurations. To determine the phase boundary, we determine the lifetime of the skyrmion configuration, prepared using the numerical solution of Eq. (7), under evolution with Eq. (43). We define the lifetime to be the time elapsed before the magnetization deviates from cylindrical symmetry. Specifically, we calculate the time $\tau^{*}$ at which the $z$ component of the magnetization deviates from its azimuthal average by a small threshold:

[TABLE]

where

[TABLE]

The skyrmion lifetime diverges at the instability phase boundary. We fit $\tau^{*}\quantity(B,K)$ for $B,K$ within the unstable phase using

[TABLE]

to determine $K^{*}\quantity(B)$, the critical anisotropy for a given applied field. The phase boundary plotted in Fig. 7 is determined by these $K^{*}\quantity(B)$ for varying $B$. This phase boundary is consistent with the phase diagram presented in Ref. [Banerjee et al., 2014].

Appendix C Numerical calculation of in-gap electronic states

To justify our tight-binding model and verify the analytical single-skyrmion orbitals, we determine the electronic surface states numerically on a lattice. We can construct a lattice model whose low-energy excitations are governed by Eq. (11). Following Ref. [Marchand and Franz, 2012], we consider an electronic system governed by

[TABLE]

where

[TABLE]

The term proportional to $\eta$ in Eq. (51) gaps the lattice model’s Dirac cones at the edge of the Brillouin zone, leaving low-energy excitations only around $\mathbf{k}=0$. Thus for $\eta>\Delta$, the low-energy, long-wavelength, physics of this model should be given by Eq. (11).

We used the python package KwantGroth et al. (2014) to determine the electronic structure associated with the magnetization generated from the simulations described in Appendix A.