Universality of annihilation barriers of large magnetic skyrmions in chiral and frustrated magnets

TL;DR

This paper investigates the energy barriers that prevent large magnetic skyrmions from decaying in chiral and frustrated magnets, revealing a universal behavior for DMI-stabilized skyrmions and distinct core-based barriers for frustrated systems.

Contribution

It provides analytic expressions for saddle point textures and energy barriers of large skyrmions, highlighting the universal upper limit in DMI-stabilized cases and core-driven barriers in frustrated systems.

Findings

Universal upper limit for DMI-stabilized skyrmion barriers

Negative stiffness leads to core-based barriers in frustrated magnets

Analytic expressions for saddle point configurations

Abstract

Magnetic skyrmions are whirls in the magnetization whose topological winding number promises stability against thermal fluctuations and defects. They can only decay via singular spin configurations. We analyze the corresponding energy barriers of skyrmions in a magnetic monolayer for two distinct stabilization mechanisms, i.e., Dzyaloshinskii-Moriya interaction (DMI) and competing interactions. Based on our numerically calculated collapse paths on an atomic lattice, we derive analytic expressions for the saddle point textures and energy barriers of large skyrmions. The sign of the spin stiffness and the sign of 4th-order derivative terms in the classical field theory determines the nature of the saddle point and thus the height of the energy barrier. In the most common case for DMI-stabilized skyrmions, positive stiffness and negative 4th-order term, the saddle point energy approaches a universal upper limit described by an effective continuum theory. For skyrmions stabilized by frustrating interactions, the stiffness is negative and the energy barrier arises mainly from the core of a singular vortex configuration.

Tables

Table 1. Table 1: Coefficients of the higher order discretization scheme for a derivative ∂ α 𝐧 subscript 𝛼 𝐧 \partial_{\alpha}{{\mathbf{\boldsymbol{n}}}} as defined in Eq. ( 30 ).

$𝒪\left(a^{2N}\right)$	$d_1$	$w_{1,0}$	$w_{1,\pm 1}$	$w_{1,\pm 2}$	$w_{1,\pm 3}$	$w_{1,\pm 4}$
$N=1$	$2$	$0$	$\pm 1$
$N=2$	$12$	$0$	$\pm 8$	$\mp 1$
$N=3$	$60$	$0$	$\pm 45$	$\mp 9$	$\pm 1$
$N=4$	$840$	$0$	$\pm 672$	$\mp 168$	$\pm 32$	$\mp 3$

Table 2. Table 2: Coefficients of the higher order discretization scheme for a second derivative ∂ α 2 𝐧 superscript subscript 𝛼 2 𝐧 \partial_{\alpha}^{2}{{\mathbf{\boldsymbol{n}}}} as defined in Eq. ( 31 ).

$𝒪\left(a^{2N}\right)$	$d_2$	$w_{2,0}$	$w_{2,\pm 1}$	$w_{2,\pm 2}$	$w_{2,\pm 3}$	$w_{2,\pm 4}$
$N=1$	$1$	$-2$	$1$
$N=2$	$12$	$-30$	$16$	$-1$
$N=3$	$180$	$-490$	$270$	$-27$	$2$
$N=4$	$5040$	$-14350$	$8064$	$-1008$	$128$	$-9$

Table 3. Table 3: Coefficients of the higher order discretization scheme for a third derivative ∂ α 3 𝐧 superscript subscript 𝛼 3 𝐧 \partial_{\alpha}^{3}{{\mathbf{\boldsymbol{n}}}} as defined in Eq. ( 32 ) with N ′ = N − 1 superscript 𝑁 ′ 𝑁 1 N^{\prime}=N-1 .

$𝒪\left(a^{2N^{\prime }}\right)$	$d_2$	$w_{2,0}$	$w_{2,\pm 1}$	$w_{2,\pm 2}$	$w_{2,\pm 3}$	$w_{2,\pm 4}$	$w_{2,\pm 5}$
$N=2$	$2$	$0$	$\mp 2$	$\pm 1$
$N=3$	$8$	$0$	$\mp 13$	$\pm 8$	$\mp 1$
$N=4$	$240$	$0$	$\mp 488$	$\pm 338$	$\mp 72$	$\pm 7$
$N=5$	$30240$	$0$	$\mp 70098$	$\pm 52428$	$\mp 14607$	$\pm 2522$	$\mp 205$

Table 4. Table 4: Coefficients of the higher order discretization scheme for a 4th derivative ∂ α 4 𝐧 superscript subscript 𝛼 4 𝐧 \partial_{\alpha}^{4}{{\mathbf{\boldsymbol{n}}}} as defined in Eq. ( 32 ) with N ′ = N − 1 superscript 𝑁 ′ 𝑁 1 N^{\prime}=N-1 .

$𝒪\left(a^{2N^{\prime }}\right)$	$d_2$	$w_{2,0}$	$w_{2,\pm 1}$	$w_{2,\pm 2}$	$w_{2,\pm 3}$	$w_{2,\pm 4}$	$w_{2,\pm 5}$
$N=2$	$1$	$6$	$-4$	$1$
$N=3$	$6$	$56$	$-39$	$12$	$-1$
$N=4$	$240$	$2730$	$-1952$	$676$	$-96$	$7$
$N=5$	$15120$	$193654$	$-140196$	$52428$	$-9738$	$1261$	$-82$

Equations

$$\begin{array}{rl}E=& -J\sum _{⟨i,j⟩}{\mathbf{m}}_i\cdot {\mathbf{m}}_j-D\sum _{⟨i,j⟩}{\hat{d}}_{ij}\cdot \left({\mathbf{m}}_i\times {\mathbf{m}}_j\right)\\ & -{\mu }_{s}H\sum _i{m}_i^z-K\sum _i(m_i^z{)}^2.\end{array}$$

$$\mathbf{m}(\mathbf{r})=\frac{1}{\sqrt{V_q}}\int {\mathrm{d}}^{2}q\tilde{\mathbf{m}}(\mathbf{q})e^{i\mathbf{q}\cdot \mathbf{r}}$$

$$\begin{array}{rl}\delta E[\mathbf{m}]=\int {\mathrm{d}}^{2}r& \frac{\mathcal{J}}{2}(\nabla \mathbf{m}{)}^2-\mathcal{D}\mathbf{m}\cdot \left(\left(\hat{z}\times \nabla \right)\times \mathbf{m}\right)\\ & -{\mu }_0\mathcal{H}(m_z-1)-\mathcal{K}(m_z^2-1)\end{array}$$

$$\begin{array}{rl}\delta E/\mathcal{J}=\int {\mathrm{d}}^2\tilde{r}& \frac{1}{2}(\tilde{\nabla }\mathbf{n}{)}^2-\mathbf{n}\cdot \left(\left(\hat{z}\times \tilde{\nabla }\right)\times \mathbf{n}\right)\\ & -h(n_z-1)-\kappa (n_z^2-1)\end{array}$$

$$\delta E_{\text{□}}^{(4)}=-\int {\mathrm{d}}^{2}r\frac{J{a}^2}{24}{\left({\partial }_{\alpha }^2\mathbf{m}\right)}^2+\frac{Da}{6}\mathbf{m}\cdot \left(\left(\hat{z}\times {\hat{e}}_{\alpha }{\partial }_{\alpha }^3\right)\times \mathbf{m}\right)$$

$$\begin{array}{rl}\delta E_{\text{△}}^{(4)}=& -\int {\mathrm{d}}^{2}r\frac{\sqrt{3}J{a}^2}{32}{\left({\nabla }^2\mathbf{m}\right)}^2\\ & +\frac{\sqrt{3}Da}{8}\left({\nabla }^2\mathbf{m}\right)\cdot \left(\left(\hat{z}\times \nabla \right)\times \mathbf{m}\right)\end{array}$$

$$\begin{array}{rl}\frac{\delta E^{(4)}}{\mathcal{J}}=\int {\mathrm{d}}^{2}r& \frac{{\mathcal{K}}_4}{2}{\left({\tilde{\nabla }}^2\mathbf{n}\right)}^2+{\mathcal{K}}_3\left({\tilde{\nabla }}^2\mathbf{n}\right)\cdot \left(\left(\hat{z}\times \tilde{\nabla }\right)\times \mathbf{n}\right)\\ +& \frac{{\mathcal{K}}_4^{\prime }}{4}{\left({\left(\tilde{\nabla }\mathbf{n}\right)}^2\right)}^2+\frac{{\mathcal{K}}_4^{\prime \prime }}{4}{\left({\tilde{\partial }}_{\alpha }\mathbf{n}\cdot {\tilde{\partial }}_{{\alpha }^{\prime }}\mathbf{n}\right)}^2\end{array}$$

$$\begin{array}{rl}& \frac{\delta E}{\mathcal{J}}\approx \int {\mathrm{d}}^{2}x\frac{1}{2}{\left(\nabla \mathbf{n}\right)}^2-\frac{R_s}{\xi }\mathbf{n}\cdot \left(\hat{z}\times \nabla \right)\times \mathbf{n}\\ & +{\left(\frac{\xi }{R_s}\right)}^2\left(\frac{{\mathcal{K}}_4}{2}{\left({\nabla }^2\mathbf{n}\right)}^2+\frac{{\mathcal{K}}_4^{\prime }}{4}((\nabla \mathbf{n}{)}^2{)}^2+\frac{{\mathcal{K}}_4^{\prime \prime }}{4}({\partial }_{\alpha }\mathbf{n}\cdot {\partial }_{{\alpha }^{\prime }}\mathbf{n}{)}^2\right)\\ & +\mathcal{O}\left(\left(\frac{R_s}{\xi }\right){}^2,{\mathcal{K}}_3\frac{\xi }{R_s}\right).\end{array}$$

$$\frac{\delta E}{\mathcal{J}}\approx \int {\mathrm{d}}^2\tilde{\mathbf{r}}\frac{1}{2}(\nabla \mathbf{n}{)}^2,$$

$${\mathbf{n}}_{\text{BP}}(\tilde{\mathbf{r}})=\hat{r}\sin {\theta }_0(\tilde{r})+\hat{z}\cos {\theta }_0(\tilde{r}),{\theta }_0(\tilde{r})=2\arctan \left(\frac{{\tilde{r}}_{\text{BP}}}{\tilde{r}}\right)$$

$$\delta E_0[{\mathbf{n}}_{\text{BP}}]/\mathcal{J}=4\pi$$

$$\frac{\delta E}{\mathcal{J}}$$

$${\bar{\mathcal{K}}}_4$$

$${\xi }_{\text{SP}}$$

$$\frac{\delta E_{\text{SP}}}{\mathcal{J}}$$

$$\begin{array}{rl}E=& -J_1\sum _{⟨i,j⟩}{\mathbf{m}}_i\cdot {\mathbf{m}}_j+J_2\sum _{⟨⟨i,j⟩⟩}{\mathbf{m}}_i\cdot {\mathbf{m}}_j\\ & -{\mu }_{s}H\sum _i{m}_i^z-K\sum _i(m_i^z{)}^2,\end{array}$$

$$\begin{array}{rl}\delta E[\mathbf{m}]/J_1=\int {\mathrm{d}}^{2}r& -\frac{{\mathcal{I}}_1}{2}(\nabla \mathbf{m}{)}^2+\frac{{\mathcal{I}}_2}{2}({\nabla }^2\mathbf{m}{)}^2\\ & -{\mu }_0\mathcal{H}(m_z-1)-\mathcal{K}(m_z^2-1).\end{array}$$

$${\mathcal{I}}_1=\sqrt{3}\left(3\frac{J_2}{J_1}-1\right),{\mathcal{I}}_2=\frac{a^2\sqrt{3}}{16}\left(9\frac{J_2}{J_1}-1\right),$$

$$\xi =\sqrt{\frac{{\mathcal{I}}_2}{{\mathcal{I}}_1}}=\frac{a}{4}\sqrt{\frac{J_1-9J_2}{J_1-3J_2}}$$

$$\begin{array}{rl}\delta E/J_1={\mathcal{I}}_1\int {\mathrm{d}}^2\tilde{r}& -\frac{1}{2}(\tilde{\nabla }\mathbf{n}{)}^2+\frac{1}{2}({\tilde{\nabla }}^2\mathbf{n}{)}^2\\ & -h(n_z-1)-\kappa (n_z^2-1),\end{array}$$

$$\delta E/J_1={\mathcal{I}}_{1}C(h,\kappa ).$$

$$\frac{\delta E}{J_1}$$

$$+{\mathcal{I}}_1{C}^{(ii)}(h,\kappa )+e^{(iii)}(J_1,J_2)$$

$$=e_0-\frac{\pi }{2}{\mathcal{I}}_1\ln \left(\frac{1}{{\mathcal{I}}_1}\right)+{\mathcal{I}}_1{e}_1(h,\kappa )$$

$$e_0\approx 2.8156$$

$$0\le \delta E{|}_{M_z^S}-\delta E^{\text{SP}}\le -({\mu }_{s}H+2K)\Delta M_z\sim J_1(h+\kappa ){\mathcal{I}}_1,$$

$$D/J=\frac{a}{\xi }.$$

$${\mu }_{0}H/J=\frac{a^2}{{\xi }^2}h,K/J=\frac{a^2}{{\xi }^2}\kappa$$

$${\mu }_{0}H/J=\frac{3a^2}{2{\xi }^2}h,K/J=\frac{3a^2}{2{\xi }^2}\kappa$$

$$J_2/J_1=\frac{1-16(\xi /a{)}^2}{3(3-16(\xi /a{)}^2)}.$$

$$\begin{array}{rl}{\mu }_{0}H/J_1& =\frac{3h}{(\xi /a{)}^2(-3+16(\xi /a{)}^2)},\\ K/J_1& =\frac{3\kappa }{(\xi /a{)}^2(-3+16(\xi /a{)}^2)}.\end{array}$$

$${\mathcal{I}}_1=\frac{2\sqrt{3}}{-3+16(\xi /a{)}^2}\sim (a/\xi {)}^2,$$

Full text

I Introduction

A magnetic skyrmion is a magnetic texture where the spin-orientation varies smoothly over many lattice constantsBogdanov and Yablonskii (1989); Mühlbauer et al. (2009); Nagaosa and Tokura (2013); Wiesendanger (2016); Fert et al. (2017); Everschor-Sitte et al. (2018). Its distinctive feature is the topological winding number: the magnetization of each skyrmion winds once in all possible directions. Arguably, the most fundamental property of a skyrmion is its stability: it is protected by an energy barrier from its destruction by thermal fluctuationsBogdanov and Hubert (1999); Büttner et al. (2018); Wang et al. (2018); Bernand-Mantel et al. (2018). The size of this energy barrier is therefore an important quantity as it determines in what sense a skyrmion is ”topologically protected”. From the application perspective, the energy barrier determines not only the thermal stability of skyrmions but, for example, also how much energy is needed to create skyrmions.

Smooth deformations of the magnetic texture cannot change the topology, therefore the creation or destruction of a skyrmion can only occur via singular spin configurations, where the spin varies rapidly on the atomic scaleNagaosa and Tokura (2013); Milde et al. (2013). A central question both of fundamental and practical importance is whether the energy barrier is ultimately determined by this singular spin configuration and, hence, the non-universal local physics on the length scale of a lattice constant, or whether it is governed by the universal physics encoded in (classical) continuum field theories. If the first scenario holds, it is often very difficult to calculate the relevant energy barriers: most experimental systems are metallicNagaosa and Tokura (2013); Everschor-Sitte et al. (2018) and a full quantum mechanical calculation of an unstable high-energy state, necessarily including the electronic degree of freedom, would be very challenging. We will show that for a wide class of skyrmion systems the second scenario holds. Despite the singular nature of the skyrmion destruction path, continuum field theories can be used to obtain quantitative (and even analytical) values for the energy barriers. Here we argue that forth-order derivative terms, which are usually not considered when using, e.g., a micromagnetic model, determine that nature of the saddle point and thus the height of the energy barrier in an essential way. A continuum field theory or, equivalently, a variant of micromagnetic simulations, which includes effects of 4th-order derivatives, can be used to calculate the energy barriers provided that the spin-stiffness is positive and that 4th-order terms have a negative sign. We will investigate two rather different types of skyrmions which are destroyed via two completely different singular spin configurations: skyrmions stabilized by Dzyaloshinskii-Moriya interactions (DMI) and skyrmions in frustrated magnets.

Due to the importance of skyrmion stability, many studies have investigated the creation and annihiliation of skyrmions and the height of the energy barrier Bessarab et al. (2015); Hagemeister et al. (2015); Rybakov et al. (2015); Rohart et al. (2016); Lobanov et al. (2016); Siemens et al. (2016); Yin et al. (2016); Stosic et al. (2017); Cortés-Ortuño et al. (2017); Bessarab et al. (2018); von Malottki et al. (2017); Uzdin et al. (2018a, b); Leonov and Inoue (2018); Varentsova et al. (2018); Desplat et al. (2018); Müller et al. (2018); Haldar et al. (2018); Derras-Chouk et al. (2018); Garanin et al. (2018); Büttner et al. (2018); Wang et al. (2018); Bernand-Mantel et al. (2018); Desplat et al. (2019); Derras-Chouk et al. (2019); von Malottki et al. (2019); Meyer et al. (2019); Potkina et al. (2019); Vlasov et al. (2019). From a methods perspective, perhaps the most important numerical approach to calculate energy barriers in classical spin systems is the geodesic nudged elastic band (GNEB) method pioneered by Bessarab et al. Bessarab et al. (2015), where one energetically optimizes spin configurations which smoothly connect an initial state (e.g. a magnetic skyrmion) and a final state (e.g. the ferromagnet). The highest-energy intermediate state, the saddle point, determines the energy barrier. Most results are availiable for skyrmions stabilized by DM interaction. Here it was found that, typically, skyrmions decay via shrinking to the size of a lattice constantRohart et al. (2016) if they cannot escape via the edge of the systemCortés-Ortuño et al. (2017); Bessarab et al. (2018) or decay via quantum mechanical tunnellingVlasov et al. (2019). Furthermore, two distinct concepts were identified that lead to higher energy barriers and, hence, more stable skyrmions: larger energy barriers were obtained for larger skyrmionsVarentsova et al. (2018); Potkina et al. (2019) while the energy barrier of smaller skyrmions was found to be increased by frustrating next-nearest neighbor exchange interactionsvon Malottki et al. (2017). For a particular case of the latter scenario, in contrast to the symmetric shrinking path, an asymmetric decay path was reportedMeyer et al. (2019) where the skyrmion annihilates via a singular defect at an off-center position. Similar paths with off-centered and on-centered defects were reported by Desplat et al.Desplat et al. (2019) for a model system where skyrmions were stabilized by frustrated interactions on a square lattice without DMI. Moreover, other methods have been employed to calculate the energy barrier for the decay of a skyrmion, including stochastic simulationsSchütte et al. (2014); Hagemeister et al. (2015); Siemens et al. (2016); Rohart et al. (2016); Yin et al. (2016); Derras-Chouk et al. (2019) or a forced decay by increasing the external fieldDerras-Chouk et al. (2018, 2019). Also attempts have been made to predict the energy barrier by approximate analytical solutions of the micromagnetic model which can be used to gather information about the isolated skyrmionBüttner et al. (2018); Bernand-Mantel et al. (2018). However, as we show in this paper, the fourth order gradient corrections beyond the standard micromagnetic model are essential to not only describe the energy barrier but also determine whether the physics of the barrier are universal or depend on the complicated microscopic details. A recent study by Derras-Chouk, Chudnovsky, and Garanin Derras-Chouk et al. (2019), which was performed in parallel to our investigations, recognizes the importance of 4th-order terms for the calculation of thermal activation barriers. We will discuss similarities and differences to the analysis of Ref. 33 below.

Experimentally, the most direct way to measure the energy barrier $\Delta E$ is to investigate the temperature dependence of the lifetime of skyrmionsWild et al. (2017); Wilson et al. (2019), which is expected to be proportional to $e^{-\Delta E/T}$. For skyrmions in rather thick films of Cu2OSe3O it was shownWild et al. (2017) that the activation barriers are determined by the emergence of a Bloch-point, i.e., a three-dimensional topological defect which caps the end of a skyrmion string. Furthermore, the pre-exponential factor of the Néel-Arrhenius law was shown to induce major corrections of $40$ orders of magnitude to the lifetime. For thin layers, however, this huge effect is not confirmed experimentally.

In the following, we will first consider a microscopic model for skyrmions stabilized by Dzyaloshinskii-Moriya interactions and investigate under which condition and in which sense its collapse can be described by a continuum field theory. A second chapter will focus on skyrmions in frustrated magnets where skyrmion annihilation is governed by a completely different process.

II Skyrmions stabilized by Dzyaloshinskii-Moriya interactions

II.1 A minimal atomistic model

In the following, we will consider a minimal model with only nearest neighbor interactions for a system which stabilizes magnetic skyrmions by (interfacial) Dzyaloshinskii-Moriya interactions in a system without inversion symmetry. Ferromagnetic nearest-neighbor spin-orbit interactions of strength $J$ favour a parallel alignment of spins, while the interfacial interfacial Dzyaloshinskii-Moriya interaction (DMI) parametrized by $D$ induces inhomogeneous magnetic textures. Furthermore, we consider the effects of an external magnetic field $H$, and a uni-axial anisotropy $K$. The energy for magnetic moments ${{\mathbf{\boldsymbol{m}}}}_{i}\!=\!{{\mathbf{\boldsymbol{M}}}}_{i}/M$ located on lattice sites ${{\mathbf{\boldsymbol{r}}}}_{i}$ then reads

[TABLE]

Here, $\langle i,j\rangle$ are nearest neighbors and every bond is counted once. The DMI-vector $\hat{d}_{i\!j}\!=\!\hat{z}\!\times\!\frac{{{\mathbf{\boldsymbol{r}}}}_{i}\!-\!{{\mathbf{\boldsymbol{r}}}}_{j}}{|{{\mathbf{\boldsymbol{r}}}}_{i}\!-\!{{\mathbf{\boldsymbol{r}}}}_{j}|}$ fixes the sense of rotation for spin spirals, resulting in Néel-type skyrmions, see Fig. 1a. We will consider both square and triangular lattices in the following. In both cases, the lattice constant is given by $a$.

For calculations of the minimal energy path, we use the geodesic nudged elastic band methodBessarab et al. (2015) (GNEB), see Sec. B.1. This method allows to construct a minimal energy path (MEP) linking the skyrmion state to a ferromagnetic state. The maximum of the MEP defines the activation barrier for skyrmion annihiliation or creation. As an example, a MEP for a given parameter set on a triangular lattice is shown in Fig. 1. The snapshots of the spin-configuration obtained along the path, Fig. 1a-f, show that the decay can be described as the shrinking of the skyrmion along a so-called collapse path Bessarab et al. (2015); Rohart et al. (2016); Lobanov et al. (2016); Cortés-Ortuño et al. (2017); Bessarab et al. (2018); Varentsova et al. (2018); Desplat et al. (2018). Our goal will be to obtain an analytic understanding of the value of the energy at the saddle point configuration, shown in Fig. 1c and as a close-up in Fig. 1f, respectively.

Our numerical studies confirm the collapse path mechanism also for other parameter sets on both the square and triangular lattice geometries. The resulting energy barriers for the annihilation of a skyrmion are summarized in Fig. 2. In agreement with the previous study by Varentsova et al.Varentsova et al. (2018), we observe that the annihilation barrier becomes larger when the Dzyaloshinskii-Moriya interaction $D/J$ is smaller. As $J/D$ determines the radius of the skyrmion, this also implies that larger skyrmions are more stable.

In the following, we will derive an expression for the asymptotic behavior of the energy barrier.

II.2 An effective continuum theory for large skyrmions

In order to study the asymptotic limit for very large skyrmions, we derive an effective continuum theory for the atomistic energy functional, Eq. (1). Therefore, we write the magnetization on the discrete lattice with lattice constant $a$ in a Fourier representation

[TABLE]

and subsequently apply an approximation for smooth variations of the magnetization, $|{{\mathbf{\boldsymbol{q}}}}|\ll a^{-1}$, where we, in a first step, consider all terms up to quadratic order $\mathcal{O}(q^{2})$. After applying the continuous back-transformation into real-space again, the leading-order gradient expansion of the energy functional reads

[TABLE]

with the effective continuum interaction constants $\mathcal{J}\!=\!J$, $\mathcal{D}\!=\!D/a$, $\mathcal{H}\!=\!H/a^{2}$, and $\mathcal{K}\!=\!K/a^{2}$ for the square lattice, or $\mathcal{J}\!=\!\sqrt{3}J$, $\mathcal{D}\!=\!\sqrt{3}D/a$, $\mathcal{H}\!=\!2H/(\sqrt{3}a^{2})$, and $\mathcal{K}\!=\!2K/(\sqrt{3}a^{2})$ for the triangular lattice, respectively. Note that this second order effective theory is independent of the geometry of the lattice and obeys a continuous symmetry for rotations around the $\hat{z}$-axis. In particular, the triangular lattice with $\mu_{s}H=1.50D^{2}/J$ and the square lattice with $\mu_{s}H=1.00D^{2}/J$ assume the same solution in this limit and all skyrmion solution are circular symmetric.

From a dimensional analysis, we find that the continuum model has the intrinsic length scale $\xi=\mathcal{J}/\mathcal{D}=aJ/D$. By rewriting the spatial distances in units of this length scale, ${{\mathbf{\boldsymbol{r}}}}=\xi\tilde{{{\mathbf{\boldsymbol{r}}}}}$, one obtains the dimensionless energy functional

[TABLE]

with $h=\mu_{0}\mathcal{J}\mathcal{H}/\mathcal{D}^{2}$ and $\kappa=\mathcal{J}\mathcal{K}/\mathcal{D}^{2}$ where for convenience we use ${{\mathbf{\boldsymbol{n}}}}(\tilde{{{\mathbf{\boldsymbol{r}}}}})={{\mathbf{\boldsymbol{m}}}}({{\mathbf{\boldsymbol{r}}}})$.

The scaling analysis reveals two important features of skyrmions in the continuum limit: (i) for the same set of effective parameters ($h$,$\kappa$) they are described by the same solution ${{\mathbf{\boldsymbol{n}}}}(\tilde{{{\mathbf{\boldsymbol{r}}}}})$ and (ii) the size of the skyrmion is linear in $\xi=\mathcal{J}/\mathcal{D}$ for fixed $h$ and $\kappa$. In the following, we will investigate corrections to the continuum limit and derive an expression for the saddle point of the collapse path for large skyrmions.

II.2.1 Corrections to the continuum theory

The above continuum theory, Eq. (4), was derived only up to quadaratic order. Corrections will, however, be important, especially when the skyrmion shrinks along the collapse path. We can derive the leading order correction terms $\delta E^{(4)}$ to the quadratic continuum theory by repeating the derivation sketched above where we take into account all terms up to fourth order $\mathcal{O}(q^{4})$ instead.

In contrast to the $\mathcal{O}(q^{2})$-theory, Eq. (3), the fourth order terms depend on the lattice geometry. For the square lattice, the expression takes the form

[TABLE]

with implicit summation over the spatial index $\alpha=x,y$ and we use the unit vectors $\hat{e}_{x}=\hat{x}$ and $\hat{e}_{y}=\hat{y}$. This expression inherits the discrete $90^{\circ}$ rotation symmetry of the square lattice. The correction term for the triangular lattice geometry, in turn, reads

[TABLE]

with the Laplace operator $\nabla^{2}=\partial^{2}_{x}+\partial^{2}_{y}$. Note that here the rotation symmetry is still continuous.

Using the rescaled coordinates of Eq. (4), the leading corrections to the continuum theory of a circluarly symmetric skyrmion can in general be written as

[TABLE]

where the (negative) values of $\mathcal{K}_{4}=-\frac{\mathcal{D}^{2}a^{2}}{16\mathcal{J}^{2}}$ and $\mathcal{K}_{3}=-\frac{\mathcal{D}^{2}a^{2}}{8\mathcal{J}^{2}}$ are accidentially the same for the square and the triangular lattice when expressed in terms of the coupling constants of the continuum theory. The prefactor of the quartic 4-spin terms vanish for our microscopic models, $\mathcal{K}_{4}^{\prime}=\mathcal{K}_{4}^{\prime\prime}=0$. We have included them neverthess here as they may exist in real materials, especially in metallic compounds Heinze et al. (2011). As they have the same scaling properties as the $\mathcal{K}_{4}$ term, they turn out to be equally important for the calculation of the activation barriers. For the square lattice, the rotationally invariant term $(\!\tilde{\nabla}^{2}{{\mathbf{\boldsymbol{n}}}}\!)^{2}$ has been obtained from $(\partial_{\alpha}^{2}{{\mathbf{\boldsymbol{m}}}})^{2}$ by averaging over spatial directions. A corresponding correction term proportional to $(\partial_{\alpha}^{2}{{\mathbf{\boldsymbol{n}}}})^{2}-\frac{3}{2}(\!\tilde{\nabla}^{2}{{\mathbf{\boldsymbol{n}}}})^{2}$ has been omitted in Eq. (7) as its expecation value vanishes for rotationally invariant textures. Therefore it does not contribute to linear order perturbation theory and gives only subleading corrections.

Using that $|\mathcal{K}_{3/4}|\ll 1$ for skyrmions, which are much larger than the lattice constant, $\mathcal{J}/\mathcal{D}\gg a$, we can calculate skyrmion energy up to linear order in $\mathcal{K}_{3/4}$, i.e., up to order $(a/\xi)^{2}$ by first computing the rotationally-symmetric minimum of the leading-order energy functional Eq. (4), ${{\mathbf{\boldsymbol{n}}}}_{s}(\tilde{{{\mathbf{\boldsymbol{r}}}}})=\hat{r}\sin\theta_{s}(\tilde{r})+\hat{z}\cos\theta_{s}(\tilde{r})$. The leading-order correction due the resulting energy is then computed by evaluating $\delta E^{(4)}[{{\mathbf{\boldsymbol{n}}}}_{s}(\tilde{{{\mathbf{\boldsymbol{r}}}}})]$, Eq. (7), for this skyrmion texture. A comparison of this perturbative result and the result of the lattice model, Eq. (1), is shown in Fig. 3. As expected, the two distinct approaches are in very good agreement and the corrections to the leading order continuum description are suppressed by the tiny factor $(a/\xi)^{2}\sim(\mathcal{J}/\mathcal{D})^{2}$ for skyrmions which are much larger than the lattice spacing. As in most experimental systems the relativistic spin-orbit interaction is much weaker than the magnetic exchange interaction, hence, this is also the experimentally relevant limit.

We can, however, expect that the corrections are much more important for the computation of the energy barrier along the collapse path of the skyrmion.

II.2.2 The saddle point of the collapse path

In Fig. 1 we have shown that the skyrmion shrinks before it collapses as has been discussed in numerous previous studies Bessarab et al. (2015); Rohart et al. (2016); Lobanov et al. (2016); Cortés-Ortuño et al. (2017); Bessarab et al. (2018); Varentsova et al. (2018); Desplat et al. (2018). The main reason for the collapse is that the destruction of a skyrmion necessarily involves a singular, vortex-like spin configuration. As discussed in more detail in Sec. III, a vortex-like solution of radius $R_{v}\gg a$ is associated with an energy cost approximately given by $\pi\mathcal{J}\ln[R_{v}/a]$. This energy cost can be minimized by making $R_{v}$ smaller and smaller. Along the collapse path, the skyrmion therefore shrinks in size and the vortex is formed for $R_{v}\sim a$. We will show below that the saddle point, i.e. the energetic bottleneck along the collapse path, is for a wide class of systems obtained from a very small skyrmion and not from the vortex solution. Moreover, for sufficiently large initial skyrmions this small skyrmion is still much larger than the underlying atomic lattice and, thus, can be described by a continuum theory.

We first analyze skyrmions with radius $R_{s}$ for $a\ll R_{s}\ll\xi$. This is done most conveniently by introducing dimensionless coordinates ${{\mathbf{\boldsymbol{x}}}}$ with ${{\mathbf{\boldsymbol{r}}}}=R_{s}{{\mathbf{\boldsymbol{x}}}}$ or, equivalently, $\tilde{{{\mathbf{\boldsymbol{r}}}}}=\frac{R_{s}}{\xi}{{\mathbf{\boldsymbol{x}}}}$. In these units we find

[TABLE]

The DM interaction is suppressed by the small prefactor $R_{s}/\xi$ and also the higher-order gradient terms are small as $\mathcal{K}_{4}(\xi/R_{s})^{2}\sim(a/R_{s})^{2}\ll 1$ for $R_{s}\gg a$. Therefore, the main contribution to the energy arises for $a\ll R_{s}\ll\xi$ from the scale-invariant term

[TABLE]

where we now use again the dimensionless variables $\tilde{r}$ of Eq. (4) and (7).

A solution which minimizes the energy (9) for fixed skyrmion winding number was found by Belavin and Polyakov Belavin and Polyakov (1975). This Belavin-Polyakov skyrmion ${{\mathbf{\boldsymbol{n}}}}_{\text{\tiny{BP}}}(\tilde{{{\mathbf{\boldsymbol{r}}}}})$ can be expressed as

[TABLE]

where $\tilde{r}_{\text{\tiny{BP}}}=r_{\text{\tiny{BP}}}/\xi$ is the radius of this skyrmion, i.e., where the magnetization is in-plane. The resulting energy of the Belavin-Polyakov skyrmion with respect to the polarized state evaluates to

[TABLE]

and is independent of the skyrmion radius $r_{\text{\tiny{BP}}}$. We therefore have to compute the leading order correction by evaluating all terms in Eq. (8) using the Belavin-Polyakov solution.

The evaluation of the correction from the DM interaction turns out to be non-trivial. A straightforward integration gives $\int\!d^{2}\tilde{r}\,\,{{\mathbf{\boldsymbol{n}}}}_{\text{\tiny{BP}}}\!\cdot\!(\hat{z}\!\times\!\tilde{\nabla})\!\times\!{{\mathbf{\boldsymbol{n}}}}_{\text{\tiny{BP}}}=4\pi r_{\text{\tiny{BP}}}$. While this integral is well converging, one has to take into account a mathematical subtlety (apparently overlooked in Ref. [Derras-Chouk et al., 2019]) which arises because the Belavin-Polyakov decays very slowly (proportional to $1/\tilde{r}$) for large distances towards the ferromagnetic state. In the presence of anisotropy terms or external magnetic fields, the real skyrmion decays, however, faster in this limit. This effect can be described by multiplying the $\theta_{0}(\tilde{r})$ in Eq. (10) with a smooth cutoff function $f_{c}(\tilde{r})$ with $f_{c}(\tilde{r}\ll\tilde{r}_{c})=1$ and $f_{c}(\tilde{r}\gg\tilde{r}_{c})=0$, where $\tilde{r}_{c}\gg r_{\text{\tiny{BP}}}$ is a cutoff length much larger than the radius of the Belavin-Polyakov skyrmion. With this cutoff, the integral obtains an extra contribution $-4\pi\tilde{r}_{\text{\tiny{BP}}}\int_{0}^{\infty}dx\,f_{c}^{\prime}(x)=4\pi r_{s}(f_{c}(0)-f_{c}(\infty))=4\pi r_{\text{\tiny{BP}}}$ in this limit. While this contribution vanishes for the idealized Belavin-Polyakov skyrmion ($f_{c}=1$), it changes the value of the integral from $4\pi r_{\text{\tiny{BP}}}$ to $8\pi r_{\text{\tiny{BP}}}$ for any finite cutoff function.

Taking this effect into account, we obtain for $a\ll r_{\text{\tiny{BP}}}\ll\xi$

[TABLE]

A similar formula (up to prefactors as discussed above) was obtained very recently in Ref. [Derras-Chouk et al., 2019]. The property of this function depends on the sign of $\bar{\mathcal{K}}_{4}$ as is shown in Fig. 4. We will first discuss the case of negative $\bar{\mathcal{K}}_{4}$ realized in the microscopic model of Eq. (1). For negative $\bar{\mathcal{K}}_{4}$, the energy $\delta E$ has a maximum as a function of $\tilde{r}_{\text{\tiny{BP}}}$ which defines the saddle point of the collapse path, see Fig. 4. We denote the size of the Belavin-Polyakov skyrmion at the saddle point by $\xi_{\text{\tiny{SP}}}$ and its energy relative to the ferromagnetic state by $\delta E_{\text{\tiny{SP}}}$ and obtain for $\bar{\mathcal{K}}_{4}<0$

[TABLE]

where we inserted the value of $\bar{\mathcal{K}}_{4}$ for the microscopic model, Eq. (1), and used $\xi=\mathcal{J}/\mathcal{D}=aJ/D$ as above. Note that for large skyrmions with $D\ll J/6$, we find that $\xi_{\text{\tiny{SP}}}\gg a$ which formally justifies the use of the continuum model for the calculation of the saddle point. For larger skyrmions, the height of the barrier approaches the universal limit $4\pi\mathcal{J}$. The importance of this energy scale (without a quantitative discussion of the range of validity and of corrections) has also been realized in a number of previous studiesKiselev et al. (2011); Leonov et al. (2016); Büttner et al. (2018); Bernand-Mantel et al. (2018).

In Fig. 5 we compare the analytical formula, Eq. (13), (grey line) to numerical results obtained from the GNEB method in the original atomistic model, Eq. (1), for various values of anisotropy and magnetic field and for both the square and the triangular lattice. For $D\ll J/6$, a very good agreement is obtained, fully confirming the validity of the analytical formula for $\bar{\mathcal{K}}_{4}<0$. The comparison to the numerics clearly show that the energy bottleneck arises in a regime where the continuum theory is valid.

Note that the energy of the saddle point depends for large skyrmions only very weakly on the magnetic field and the anisotropy as they provide only subleading corrections to the saddle point energy. In contrast, the energy of the skyrmion does depend strongly on those parameters, see Fig. 3, but is almost independent of $\bar{\mathcal{K}}_{4}$. As the annihilation barrier is given by the difference of saddle-point and skyrmion energy, the field dependence of the barrier height is dominated by the field dependence of the skyrmion energy.

For a positive sign of $\bar{\mathcal{K}}_{4}$, the energy of the skyrmion as function of its radius does not have a maximum for $r_{\text{\tiny{BP}}}>a$, see Fig. 4. In this case, the energy of the saddle point can become larger than $4\pi\mathcal{J}$ and determined by the microscopic physics at the length scale of the lattice constant $a$. In a metal, for example, one can expect that details of the electronic states scattering from a highly non-colinear local spin configuration will determine the energy barrier. Within our microscopic model, positive values of $\bar{\mathcal{K}}_{4}$ can be obtained by adding a frustrated antiferromagnetic next-nearest neighbor interaction to the microscopic model (details are discussed in the next section). In Fig. 6 we show the resulting energy barriers calculated with the GNEB method for frustrating next-nearest neighbor interactions which drive $\tilde{\mathcal{K}}_{4}$ positive. For sufficient strong frustration, the energy barrier becomes larger than $4\pi\mathcal{J}$, see Fig. 6.

III Skyrmions in frustrated Heisenberg ferromagnets

Recent theoretical studies proposed competing interactions as an alternative mechanism to stabilize magnetic skyrmions in symmetric magnetic layersIvanov et al. (1990); Okubo et al. (2012); Leonov and Mostovoy (2015); Lin and Hayami (2016); Hayami et al. (2016). Under certain conditions, skyrmions can exist in these systems as stable quasi-particles without the need for broken lattice symmetry and spin-orbit coupling. The lack of broken inversion symmetry results in equal energies for Néel and Bloch skyrmions and all intermediate helicities and, furthermore, also for anti-skyrmions.

In the following, we will determine the minimal energy transition path from a skyrmion to the polarized state or an anti-skyrmion and develop an approximate expression for large skyrmions.

III.1 A minimal atomistic model with inversion symmetry

For the following analysis of the stability of skyrmions in symmetric magnets, we consider a minimal model as proposed in Ref. [Leonov and Mostovoy, 2015]: Classical Heisenberg spins ${{\mathbf{\boldsymbol{m}}}}_{i}={{\mathbf{\boldsymbol{M}}}}_{i}/M$ on a triangular lattice interact ferromagnetically with their nearest neighbors and antiferromagnetically with their next nearest neighbors. Hence the system is frustrated as both interactions can never be optimized simultaneously and, as a compromise, co-planar spirals can form. In addition, the external magnetic field ${{\mathbf{\boldsymbol{H}}}}=H\hat{z}$ and a uni-axial anisotropy $K$ introduce a preferred axis for the magnetization and stabilize individual skyrmions. The energy of the magnetic momenta ${{\mathbf{\boldsymbol{m}}}}_{i}$ located on lattice sites ${{\mathbf{\boldsymbol{r}}}}_{i}$ then reads

[TABLE]

where we have adjusted the signs such that $J_{1}>0$ describes a ferromagnetic nearest-neighbor, while $J_{2}>0$ encodes the antiferromagnetic next-nearest-neighbor coupling. The existence of metastable skyrmions requiresLeonov and Mostovoy (2015) that $J_{2}>\frac{1}{3}J_{1}$.

An example of a minimal energy path (MEP) obtained with the geodesic nudged elastic band method (GNEB) for the decay of a skyrmion into the polarized state is shown in Fig. 7 for $J_{2}/J_{1}=31/90$, $\mu_{s}H/J_{1}=1/450$, and $K/J_{1}=1/900$. In contrast to the above discussed model with DMI, Eq. (1), the energy functional of the centro-symmetric magnetic is symmetric under rotations of the magnetic momenta around the $\hat{z}$-axis, hence the helicity of a skyrmion is a zero mode and leads to interesting dynamicsZhang et al. (2017); Leonov and Mostovoy (2017); Lohani et al. (2019). However, this additional degree of freedom does not influence the minimal energy path. Thus, for convenience, the skyrmion in Fig. 5 is displayed with the same helicity as in Sec. II but the result is the same for Bloch-type or any other helicity.

As can be clearly seen from the real-space images along the MEP, Figs 7a-e, the GNEB method now finds a completely different decay mechanism compared to the DMI-stabilized skyrmion, Figs 1a-e. Instead of shrinking the entire skyrmion to the size of a lattice constant, the spins in the inner region with $m^{z}\leq 0$ turn in-plane, forming a large vortex, while the outer region remains almost constant. The saddle point, Figs 7c or f, is assumed when the core of the skyrmion is approximately in-plane with a $120^{\circ}$ vortex-like singularity centered on a triangular plaquette. The transition from the saddle point to the polarized state then occurs via first rotating the spins around the core of the vortex out of plane and then smoothly turning the remaining texture in the $\hat{z}$-direction. The energy along the MEP is shown in the lower panel of Fig 7. The energy of the saddle point is an order of magnitude larger than the energy of the skyrmion, in contrast to the DMI-skyrmion, Fig. 1, whose energy is of the same order as the energy of the saddle point. However, a comparison in units of the nearest neighbor exchange interaction $J$ or $J_{1}$ reveals that the energy of the saddle point is, in fact, smaller than for the DMI-skyrmions.

In the following we will analyze the energy barrier set by the saddle point and find that the skyrmion and the saddle point are governed by different physics which leads to distinct energy scales for both. This is in contrast to the previous section, where both the skyrmion and the saddle point are well described within a micromagnetic approximation.

III.2 An effective continuum theory for large skyrmions

In order to study the energy of the saddle point of the MEP for a large skyrmion, we can derive an effective energy functional with the same small-$q$ expansion as in Sec. II.2. In the absence of DMI, skyrmions are stabilized by a $(\nabla{{\mathbf{\boldsymbol{m}}}})^{2}$ contribution to the energy with negative sign. To obtain a well-defined field theory bounded from below, one has to consider also the next order $\mathcal{O}(q^{4})$ as already computed in Eq. (6). The resulting $\mathcal{O}(q^{4})$ continuum approximation of the energy functional reads

[TABLE]

Note that we use different conventions compared to the previous section. Energies are measured in units of $J_{1}$ and the coupling constants are obtained as

[TABLE]

while $\mathcal{H}=2H/(\sqrt{3}a^{2}J_{1})$ and $\mathcal{K}=2K/(\sqrt{3}a^{2}J_{1})$. Importantly, skyrmions are only stabilized Leonov and Mostovoy (2015) for negative spin stiffness $-J_{1}\mathcal{I}_{1}$, i.e., for $\mathcal{I}_{1}>0$ or $J_{2}>\frac{1}{3}J_{1}$. In this regime, we obtain $\mathcal{I}_{2}>0$. Both $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ depend only on $J_{2}/J_{1}$ and are related by $\mathcal{I}_{2}=a^{2}(\sqrt{3}/8+3\,\mathcal{I}_{1}/16)$.

From a dimensional analysis we can extract the length scale $\xi$

[TABLE]

which is tuned to vary the size of the skyrmion. Large skyrmions with $\xi\gg a$ are stabilized for small frustration $\mathcal{I}_{1}\to 0$, i.e., by fine-tuning the nearest and next-nearest neighbor couplings $J_{2}\to\frac{1}{3}J_{1}$. Experimentally, this is more complicated to achieve than large DMI-stabilized skyrmions, as the latter only require a small spin-orbit coupling, realized in most materials.

A rescaling of the energy of the skyrmion using dimensionless units with ${{\mathbf{\boldsymbol{r}}}}=\xi\,\tilde{{{\mathbf{\boldsymbol{r}}}}}$ gives

[TABLE]

with $h=\mu_{0}\mathcal{H}\,\mathcal{I}_{2}/\mathcal{I}_{1}^{2}$ and $\kappa=\mathcal{K}\,\mathcal{I}_{2}/\mathcal{I}_{1}^{2}$.

For a fixed value of the rescaled parameters ($h$,$\kappa$) the energy of a skyrmion is therefore linear in $\mathcal{I}_{1}\sim(a/\xi)^{2}$ and vanishes for larger and larger skyrmions,

[TABLE]

where $C(h,\kappa)$ is a numerical constant depending on the parameters $h$ and $\kappa$. This linear dependence on $\mathcal{I}_{1}$ is confirmed in Fig. 8 where we plot the energy of the skyrmion in the atomistic model, Eq. (14), as function of the frustration for various fixed values of the rescaled external field $h$ and rescaled anisotropy $\kappa$.

III.2.1 The saddle point of the minimal energy path for large skyrmions

In Sec. II.2.2, we have argued that the DMI-stabilized skyrmion has to shrink before it is destroyed along its minimal energy path because a vortex of radius $R$ is associated with an energy cost of $\mathcal{J}\pi\ln(R/a)$ which can be minimized by reducing the size of $R$. For a skyrmion stabilized by frustration, the positiv constant $\mathcal{J}$ is replaced by the negative constant $-J_{1}\mathcal{I}_{1}$. Therefore there is no reason to expect a shrinking of the skyrmion. Instead, the decay occurs via the formation of a vortex with a size comparable to the skyrmion size, fully consistent with our numerical results, see Fig. 7. Fig. 7c and Fig. 7f show that the saddle point configuration is indeed a vortex state.

To estimate the energy cost of a vortex, we decompose the skyrmion into three parts: (i) an outer region $\tilde{r}>\tilde{R}$ where the magnetization smoothly turns into the polarized state, (ii) the region close to the center $\tilde{r}_{0}<\tilde{r}<\tilde{R}$ where the magnetization is approximately in-plane and pointing outwards ${{\mathbf{\boldsymbol{n}}}}(\tilde{{{\mathbf{\boldsymbol{r}}}}})=\hat{r}$ and (iii) the core of the vortex $\tilde{r}<\tilde{r}_{0}=r_{0}/\xi$, where $r_{0}$ is a length scale of the order of the lattice constant, $r_{0}\sim a$. In region (iii), one has to take into account the discrete lattice and the microscopic interactions as the singular nature of the vortex impedes a continuum description. Its energy $e^{(iii)}(J_{1},J_{2})$ is mainly determined by the large interactions $J_{1}$ and $J_{2}$ and only weakly affected by the small anisotropy or external field. In turn, the magnetization in the regimes (i) and (ii) varies smoothly and hence we can use the continuum model, Eq. (18). From the outer region (i) we obtain a term $\mathcal{I}_{1}C^{(i)}(h,\kappa)$ which is linear in the small prefactor $\mathcal{I}_{1}$ and depends on $h$ and $\kappa$. More interesting is region (ii) as the gradient terms give contributions which are singular as a function of $r_{0}$. Collecting all terms we find

[TABLE]

where $e_{0}=\frac{\sqrt{3}\,\pi}{16\,\tilde{r}_{0}^{2}}+e^{(iii)}(J_{1},J_{2})$ is a constant independent of $\mathcal{I}_{1}$ in the limit of small $\mathcal{I}_{1}$. $e_{0}J_{1}$ is the energy of the saddle point in the limit $\mathcal{I}_{1}=0$. It cannot be calculated in any continuum theory as it captures the singular center of the vortex but, instead, it can be computed from from the atomistic energy functional Eq. (14) by considering a vortex centered on a triangular plaquette for $\mathcal{I}_{1}=0$, i.e., for $J_{2}=J_{1}/3$ and $H=K=0$. Note that for these parameters the energy of an infinite-size vortex is finite (as $\mathcal{I}_{1}$ vanishes), while its energy diverges for $J_{2}\neq J_{1}/3$. A very good quantitative estimate for $e_{0}$ with a precision of about $1\%$ can simply be obtained by calculating the energy of the in-plane vortex configuration ${{\mathbf{\boldsymbol{m}}}}_{i}=\frac{{{\mathbf{\boldsymbol{r}}}}_{i}-{{\mathbf{\boldsymbol{r}}}}_{c}}{|{{\mathbf{\boldsymbol{r}}}}_{i}-{{\mathbf{\boldsymbol{r}}}}_{c}|}$, where $r_{c}$ is the center of a triangular plaquette. From this calculalation we obtain $e_{0}=2.8482$. A simple optimization of this solution obtained by minimizing the energy while keeping the spins in-plane yields with very high precision the value

[TABLE]

for classical spins on a triangular lattice.

In Fig. 9 we compare the analytical formula (20) to GNEB calculations of the energy barrier for various sets of parameters. In the large-skyrmion limit, $\mathcal{I}_{1}\to 0$, we recover the result $e_{0}J_{1}$ from Eq. (21) with high precision. To describe the deviations from this asymptotic value, we can use Eq. (20), where we treat $e_{1}(h,\kappa)$ for fixed $h$ and $\kappa$ as the only fitting parameter. Here the logarithmic corrections characteristic of the vortex solution are crucial for the energy of the saddle point for finite but small $\mathcal{I}_{1}$, see Eq. (20). We thus obtain a highly accurate quantitative description of the energy barrier for $\mathcal{I}_{1}\lesssim 0.1$, fully corrobating the physical picture of a large-vortex collapse path for (large) skyrmions stabilized by frustrating interactions.

Our analytical treatment was originally based on the picture of a decay path of the skyrmion which is rotationally symmetric. As discussed in more detail in Appendix C, for large anisotropies the decay path of the skyrmion obtained from GNEB calculations is highly asymmetric Meyer et al. (2019) as it is more favourable to locate the vortex in a region where skyrmion-spins are in-plane. As the transition still occurs via a large vortex, this does, however, not invalidate our analysis and Eq. (20) is still valid for small $\mathcal{I}_{1}$. Only the numerical value of the constant $e_{1}(h,\kappa)$ is affected by the asymmetric decay path.

III.2.2 Decay into an anti-skyrmion

The decay of a skyrmion into a polarized state is not the only decay mechanism which can be of importance. Other mechanisms include the bimeron-instabilityEzawa (2011) where the skyrmion elongates arbitrarily which can trigger the duplication of a skyrmionMüller et al. (2018). For large uni-axial anisotropy when the decay is asymmetric, see Sec. C, we also observe an intermediate skyrmion-like state with zero winding number similar to Ref. 35 which, however, decays immediately. In the frustrated system without DMI discussed here, the lack of a handedness leads to the degeneracy of skyrmions and anti-skyrmions. Furthermore, they not only have the same energy but also the same magnetization which enables quantum mechanical tunnellingLohani et al. (2019).

For the calculation of the minimal energy path from a skyrmion to an anti-skyrmion we prepare an initial path which exhibits a maximum with two simultaneous vortices, see Fig. 10f. This artificial maximum has approximately twice the energy of a single singularity. During the application of the GNEB method, see Sec. B.1, the initial maximum decays into two separate maxima with one vortex each. The intermediate minimum turns out to be the fully polarized state. Therefore, the total minimal energy path is the decay of the skyrmion into a ferromagnetic state followed by the creation of the anti-skyrmion in a reverse process. The resulting path is shown in Fig. 10 (red) for $J_{2}/J_{1}=9/25$, $\mu_{s}H/J_{1}=3/250$, $K/J_{1}=3/500$ and can be constructed from the individual decay/creation paths discussed in Sec. III.2.

Since the minimal energy path is via an intermediate polarized state, we conclude that a direct transition from a skyrmion to a antiskymrion by thermal fluctuations is very unlikely and of no practical importance. Assuming that skyrmion and antiskyrmion have an energy which is smaller than the energy of the ferromagnet, it is more likely that skyrmions (or antiskyrmions) proliferate spontaneously. This conclusion holds at least in the absence of defects and boundaries.

III.2.3 Conservation of magnetization

The Hamiltonian of Eq. (14) describes formally a system where the total magnetization in z-direction $M_{z}=\sum_{{\mathbf{\boldsymbol{r}}}}m_{z}({{\mathbf{\boldsymbol{r}}}})$ is a conserved quantity. In a real material, the magnetization is never exactly conserved due to (weak) spin-orbit coupling or due to spin-lattice relaxation processes involving spin-flip terms and phonons. In numerical simulations of the Landau-Lifshitz-Gilbert equation, this is usually taken into account in a phenomenological way by introducing the Gilbert damping $\alpha$. The presence of spin-flip processes in the thermal fluctuations justifies that we have ignored spin conservation when calculating the the relevant energy barriers.

There are, however, situations in which one would like to know the effective energy barrier in the presence of the magnetization constraint. This barrier will, for example, be important for the quantum-mechanical tunneling process from a skyrmion to an antiskyrmion in a frustrated magnet with little or no spin-orbit interactionsLohani et al. (2019). Even in the classical limit, there can be the situation where a process, which does not requiring spin-relaxation and has a larger energy barrier, competes with a process with a lower barrier requiring spin relaxation.

In general, any constraint to a system leads to an increase of its ground-state energy. Similarly, the energy of the saddle point, $\delta E|_{M_{z}^{S}}$, where the magnetization is fixed to the magnetization $M_{z}^{S}$ of the skyrmion, will be higher than the saddle point energy, $\delta E^{\text{\tiny{SP}}}$, without that constraint, $\delta E^{\text{\tiny{SP}}}\leq\delta E|_{M_{z}^{S}}$. The vortex-style saddle point has a higher magnetization than the initial skyrmion state, $\Delta M_{z}=M_{z}^{S}-M_{z}^{\text{\tiny{SP}}}<0$. Asumming an infinitely large system, we obtain an upper bound for the energy difference of the saddle points, $\Delta E|_{M_{z}^{0}}-\delta E^{\text{\tiny{SP}}}$ by using a spin configuration which is identical to the saddle point configuration close to the origin, but it absorbs the excess magnetization $\Delta M_{z}$ by slightly tilting the ferromagnetic magnetization far away from the origin. This costs the energy $-(\mu_{s}H+2K)\Delta M_{z}$ and we therefore find

[TABLE]

where we used that $|\Delta M_{z}|\sim\xi^{2}\sim a^{2}/\mathcal{I}_{1}$. At least in the limit of small $\mathcal{I}_{1}$, the energy barrier $\delta E|_{M_{z}^{S}}$ for processes which conserve magnetization is therefore only slightly higher than the barriers $\delta E^{\text{\tiny{SP}}}$ calculated in Fig. 9.

IV Discussion

We have calculated numerically the barriers for annihilation and creation of successively larger skyrmions in monolayers of chiral magnets and frustrated magnets. Three different scenarios for the destruction of a skyrmion state emerged. DMI-stabilized skyrmions in chiral magnets decay by shrinking. Depending on the sign of the 4th-order derivative terms, $\bar{\mathcal{K}}_{4}$, the saddle-point spin configuration which determines the energy barrier is either smooth (i) or singular (ii) on the length scale set by the lattice constant.

Most common seems to be case (i), $\bar{\mathcal{K}}_{4}<0$, where the saddle-point energy of a large DMI skyrmion is simply given by $4\pi\mathcal{J}$, the energy of the famous Belavin-Polyakov skyrmion in a system with spin-stiffness $\mathcal{J}$. Corrections proportional to $(D/J)^{2/3}$ always reduce the saddle point energy. In this case, one can use micromagnetic simulations to calculate the energy barrier but to obtain the correction to the universal value $4\pi\mathcal{J}$, one has to include 4th-order derivative terms in the model, see Appendix B.2. While the spin stiffness $\mathcal{J}$ can be measured directly using Neutron scattering and is routinely calculated using ab-initio methods, a measurement (or ab-initio calculation) of $\bar{\mathcal{K}}_{4}$ is needed as an input to the micromagnetic simulations, which is often much more challenging to obtain.

In case (ii), $\bar{\mathcal{K}}_{4}>0$, the saddle point energy of a DMI skyrmion is determined by microscopic physics on the length scale of the lattice spacing. The saddle point energy is typically larger than $4\pi\mathcal{J}$ in this case and determined by a spin configuration where spin orientation varies rapidly on the length scale of a lattice constant. This case is, for example, realized when frustrating interactions reduce the size of the spin-stiffness $\mathcal{J}$ while keeping the energy of a local vortex configuration high. While one can treat such a situation easily by brute-force numerics in idealized classical-spin models, it is very difficult to model this situation for an insulating spin-$1/2$ quantum spin system or for metallic compounds. In the latter case it would be interesting to combine the GNEB approach with an electronic ab-initio calculation.

A third case is realized in Heisenberg ferromagnets where frustration induces a negative spin stiffness, which provides an alternative stabilization mechanism for skyrmions. In this case, the energy barrier is determined by a vortex solution which has a size of the order of the skyrmion radius. Here the main contribution to the energy arises from the singular core of the vortex state, which again depends on microscopic details on the length scale of the lattice spacing. For classical spin models and large skyrmions, it was possible to get an accurate value for the energy barrier by simply considering an infinite vortex configuration in a model fine-tuned to have a vanishing spin-stiffness.

In our study, we have only considered a magnetic monolayer. It would be very interesting to extend the analysis to multi-layer compounds or even three dimensional systems. For very thin systems and a few Heisenberg-coupled monolayers, we can expect that one can simply multiply all energies by the number of layers. For thick system, in contrast, the nature of the saddle point configuration changes. The relevant singular spin configuration is in this case a Bloch point (or a pair of Bloch points)Milde et al. (2013); Schütte and Rosch (2014); Rybakov et al. (2015); Wild et al. (2017). In this case, the core-contribution of the energy of order of the exchange coupling $J$ is expected to be much smaller (see, e.g., supplement of Ref. [39]) than contributions from the surrounding of order $J^{2}/D$.

In conclusion, we have studied the topological protection of two-dimensional skyrmions. For a large class of 2d systems, we could identify the precise conditions required to obtain a universal energy barrier approximately given by $4\pi\mathcal{J}$. For other cases, however, microscopic details become more important which make the computation of the topological energy barrier a challenging task for real materials.

V Acknowledgments

J.M. thanks G. P. Müller, P. Bessarab, N. Kiselev, M. Hoffmann, S. von Malottki, S. Meyer, and S. Blügel for the fruitful discussions. This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) CRC 1238 project C04, project number 277146847. We furthermore thank the Regional Computing Center of the University of Cologne (RRZK) for providing computing time on the DFG-funded High Performance Computing (HPC) system CHEOPS as well as support.

Appendix A Parameters

We chose the interaction parameters in our calculations by first fixing the parameters of the continuum theory, which allows to compare data sets described by the same (lowest order) continuum theory.

A.1 Non-symmetric magnets

For non-symmetric magnets with DMI, Sec. II, we can map the atomistic model, Eq. (1), onto the continuum model in dimensionless units, Eq. (4). The length scale $\xi$ of the system fixes the DMI constant $D$:

[TABLE]

The dimensionless magnetic field $h=\mu_{0}\mathcal{J}\mathcal{H}/\mathcal{D}^{2}$ and uni-axial anisotropy $\kappa=\mathcal{J}\mathcal{K}/\mathcal{D}^{2}$ are fixed for a series of parameters. Hence, the atomistic interaction parameters used in the simulations are

[TABLE]

on the square lattice and

[TABLE]

on the triangular lattice. The data shown in Sec. II is computed for the dimensionless parameters $(h,\kappa)=(1.1,0)$ (dark blue), $(1,0)$ (blue and red), $(0.75,0)$ (light blue), and $(0,1.3)$ (yellow). The scales $\xi$ were chosen as $\xi/a=s/0.8$ with $s=1,2,4,8,16,32,64$, except for the anisotropy-stabilized skyrmion where the data for $s=64$ did not converge in a reasonable time.

A.2 Centro-symmetric magnets

For symmetric magnets, Sec. III, the atomistic model, Eq. (14), maps onto the continuum model in dimensionless units in Eq. (18). Here, the length scale $\xi$ is given in Eq. (17) which yields an expression for the next-nearest neighbor interaction $J_{2}$:

[TABLE]

For fixed dimensionless field $h=\mu_{0}\mathcal{H}\,\mathcal{I}_{2}/\mathcal{I}_{1}^{2}$ and anisotropy $\kappa=\mathcal{K}\,\mathcal{I}_{2}/\mathcal{I}_{1}^{2}$, the corresponding atomistic parameters then read

[TABLE]

The data shown in Sec. III is computed for the dimensionless parameters $(h,\kappa)=(\frac{21}{80},\frac{7}{80})$, $(\frac{7}{40},\frac{7}{80})$, $(\frac{7}{40},\frac{7}{160})$, and $(\frac{7}{160},\frac{35}{160})$ as indicated in the plots. The scales $\xi$ were chosen as $\xi/a=\frac{\sqrt{7}}{4}s$ with $s=1\frac{3}{4},2,2\frac{1}{2},3,4,5,6,8,12$ with the exception of $(h,\kappa)=(\frac{7}{160},\frac{35}{160})$, where $s=12$ did not converge in a reasonable time. Note, that the frustration parameter scales as

[TABLE]

therefore, the data in Sec. III is not plotted as a function of the inverse size of the skyrmion, in constrast to the data presented in Sec. II.

Appendix B Methods

B.1 Methods: Minimal energy path calculations

For the calculation of minimal energy transition paths we use the geodesic nudged elastic band (GNEB) method with additional climbing image (CI) as described in Ref. 12. In some figures, however, we use an advanced method for the calculation of the spring forces between adjacent images as explained below.

We first prepare a single skyrmion in a polarized background with periodic boundary conditions and relax the texture into a local energy minimum. The sample size is chosen large enough so that the self-interaction effects over the periodic boundaries are negligible. An interpolation between this skyrmion and the polarized state (or the anti-skyrmion) over a total of $40$ layers serves as the starting point for the GNEB method. The total geodesic distances $\lambda_{i}$ between two adjacent layers $i-1$ and $i$ is defined as

[TABLE]

where $\phi_{i}({{\mathbf{\boldsymbol{r}}}})$ is the rotation angle (or geodesic distance) of the spin ${{\mathbf{\boldsymbol{n}}}}({{\mathbf{\boldsymbol{r}}}})$ at site ${{\mathbf{\boldsymbol{r}}}}$ from layer $i-1$ to $i$. The reaction coordinate $\Lambda(i)$ is the sum of all distances $\lambda_{n}$, $n<i$, and serves as the measure for the evolution of the texture in Figs. 1,7,9. Following Ref. 12, the spring forces between adjacent images are proportional to the distance $\lambda_{i}$. However, for the calculation of some data sets we decided to modify this part of the method by considering the actual length of the nudged band instead. A cubic interpolation of the band in the combined space of energy and reaction coordinate is already implemented as a part of the method and it is used for plotting the continuous interpolations in Figs. 1,7,9. By using a numerical approximation of the length of this curve we achieve an improved distribution of images along the nudged elastic band and therefore observe an improved convergence behavior. A calculation of a minimal energy path is converged when the slope of the elastic band in the climbing image vanishes. In our simulations this is assumed when the total projected force on the climbing image is below $10^{-10}J$ or $10^{-10}J_{1}$, respectively.

B.2 Methods: High-accuracy numerical calculations of the continuum limit

We have shown that 4th order gradient terms play an important role for the topological energy barrier of 2d magnetic skyrmions. For the discretization of a continuum model, it is therefore essential to properly model higher-order terms. A second motivation for an improved discretization of continuum models is an improvement of numerical precision. For example, the calculation of the leading $\mathcal{O}(q^{4})$ perturbative corrections to the $\mathcal{O}(q^{2})$ continuum model in Sec. II.2.1 builds on the exact skyrmion solution. However, an analytic expression of this solution is not known. We therefore calculate numerically with very high precision the skyrmion texture which minimizes the energy of the dimensionless continuum model, Eq. (4).

Usually, this is achieved by discretizing the continuum model on a square or triangular lattice with lattice constant $a$ and replacing derivatives with nearest neighbors finite difference approximations. This standard approach, however, yields exactly the same formalism as the atomistic model and therefore the same $\mathcal{O}(a^{2})$ scaling of the numerical error. Therefore, when we try to calculate the texture with high numerical precision, the runtime of the simulation soon becomes critical: only a small factor $4$ in accuracy requires twice as many simulated spins per spatial direction, resulting in a factor $4$ longer computational time per iteration step. In addition, if we assume a constant step-width $\delta t$ for the relaxation algorithm, it is limited by $\delta t/a^{2}\leq c$. We here use $\delta t/a^{2}=0.1$ for most calculations. The total runtime for the same level of convergence is therefore enhanced by a factor $16$.

In order to achieve faster convergence, we discretize the continuous model on a square lattice with lattice spacing $a$. We replace simple derivatives with higher order stencils which take into account $N$ neighbors on every side:

[TABLE]

where $\alpha=x,y$ is the spatial direction. The second order derivatives are substituted by

[TABLE]

For the discretization of the 4th order correction terms, see Sec. II.2.1, we only need the square of second derivatives. If, however, one wishes to calculate the dynamics of this system, e.g., in order to relax the magnetic texture, then also higher order derivatives and their discretizations become inevitable. In order to achieve the same scaling of the numerical error, these higher order contributions require more neighbors to be included:

[TABLE]

The weights $w_{i,j}$ and $w_{i,j}$ and the denominators $d_{i}$ are given in Tabs. 1, 2,3, and4 respectively. Furthermore, a discussion for open boundary conditions and the coefficients for the micromagnetic model including DMI can be found in Ref. 55.

This approximation scheme yields an improved $\mathcal{O}(a^{2N})$ convergence while the total runtime for a fixed number of spins only increases by a factor $N$. Here, we use a discretization of $a=0.05$ on a square lattice which, following Ref. 55, leads to a total error in the energy of the skyrmion of the order of $10^{-10}J$ while the standard $\mathcal{O}(a^{2})$ scheme yields an error or the order $10^{-2}J$ for the same discretization. Note that the numerical error of the $\mathcal{O}(a^{2})$ scheme is, in fact, the size-dependent energy of the atomistic model, which is also derived as a result of the perturbative analysis in Sec. II.2.1.

Appendix C Asymmetric decay of skyrmions

For skyrmions in a symmetric magnet without Dzyaloshinskii-Moriya interaction we observe that the initially symmetric collapse path becomes asymmetric for certain parameter regions $(h,\kappa)$, see Fig. 11, similar to the decay reported in Ref. 35. Both the symmetric, Fig. 7, and the asymmetric collapse path have in common that they let not shrink the skyrmion radius to zero and therefore require a vortex singularity. In large skyrmions where the energy is dominated by the exchange interactions, the vortex core consists of three neighboring spins lying in a common plane with a $120^{\circ}$ angle relative to each other with subleading corrections due to field and anisotropy. In the symmetric scenario this vortex core is depicted in Fig. 7f where the common plane is the xy-plane. In the asymmetric scenario one of the spins in the vortex core points in the $\hat{z}$-direction.

The total magnetization of a vortex core is independent of the orientation of the common plane of spins, therefore the external magnetic field $h$ does not favor any particular direction. The uni-axial anisotropy $\kappa$, however, is proportional to $m_{z}^{2}$ and therefore favors vortex cores where one of the spins is aligned out of plane. Therefore we only observe the asymmetric decay for the systems with a relatively large $\kappa=\frac{35}{160}$. Note that, furthermore, the discussion in the limit of large skyrmions, Sec. III.2, is also valid for the asymmetric decay as the energy contributions for the saddle point only arise from the vortex which is in the limit $\xi\to\infty$ not only independent of the external field and anisotropy but also of the orientation of the vortex in the plane. Furthermore this implies that in this extreme limit the saddle point is not uniquely symmetric or asymmetric but can be any which has a negative impact on the convergence of large systems.