Paths to annihilation of first and second-order (anti)skyrmions via (anti)meron nucleation on the frustrated square lattice

TL;DR

This paper investigates how small skyrmions and antiskyrmions on a frustrated square lattice annihilate through meron or antimeron nucleation, revealing complex energy landscapes and validating theoretical lifetime predictions with simulations.

Contribution

It introduces detailed mechanisms of skyrmion and antiskyrmion annihilation on a frustrated lattice, including energy barriers and nucleation pathways, with validation against Langevin simulations.

Findings

Annihilation occurs via meron or antimeron nucleation.

Complex energy landscape with multiple stable and unstable states.

Excellent agreement between transition state theory and simulations.

Abstract

We study annihilation mechanisms of small first- and second-order skyrmions and antiskyrmions on the frustrated $J_1-J_2-J_3$ square lattice with broken inversion symmetry (DMI). We find that annihilation happens via the injection of the opposite topological charge in the form of meron or antimeron nucleation. Overall, the exchange frustration generates a complex energy landscape with not only many (meta)stable and unstable local energy solutions, but also many possible paths connecting them. Whenever possible, we compute the activation energy and attempt frequency for the annihilation of isolated topological defects. In particular, we compare the average lifetime of the antiskyrmion calculated with transition state theory with direct Langevin simulations, where an excellent agreement is obtained.

Tables

Table 1. Table 1 : Internal energy barrier Δ ​ E Δ 𝐸 \Delta E and rate prefactor f 0 subscript 𝑓 0 f_{0} for the annihilation of a skyrmion, an antiskyrmion, and a second-order skyrmion with different values of the reduced applied field b 𝑏 b and reduced DMI d 𝑑 d .

mechanism	$b$ ($J_1$)	$d$ ($J_1$)	$\Delta E$ ($J_1$)	$f_0$ (GHz)
		skyrmion
antimeron nucl.	0.2	0.03	1.27	109.9 $\times$ 4
iso. collapse	0.3	0.07	1.62	3639.5
iso. collapse	0.7	0.2	1.80	1247.5
		antiskyrmion
meron nucl. SP1	0.2	0.03	0.022	1.4 $\times$ 2
meron nucl. SP2	0.2	0.03	0.43	10.6 $\times$ 4
		2nd. order sk.
antimeron nucl.	0.14	0.005	0.0011	-
sk. pair division	0.14	0.005	0.062	-

Equations

$$E=-\sum_{\langle ij\rangle}J_{ij}\mathbf{m}_{i}\cdot\mathbf{m}_{j}-\sum_{\langle ij\rangle}\mathbf{D}_{ij}\cdot\big{(}\mathbf{m}_{i}\times\mathbf{m}_{j}\big{)}-K\sum_{i}m_{z,i}^{2}-B_{z}\sum_{i}m_{z,i},$$

$$N_s=\int \mathrm{d}{\mathbf{r}}^2{\rho }_s=-n,$$

$${\rho }_s=\frac{1}{4\pi }\mathbf{m}\cdot ({\partial }_x\mathbf{m}\times {\partial }_y\mathbf{m}).$$

$$\frac{\mathop{}\!\mathrm{d}\mathbf{m}_{i}}{\mathop{}\!\mathrm{d}\bar{t}}=\frac{1}{\alpha}\mathbf{m}_{i}\times\big{[}\mathbf{b}_{\text{eff}}+\mathbf{b}_{\text{fl}}(\bar{t})\big{]}-\mathbf{m}_{i}\times\big{(}\mathbf{m}_{i}\times\big{[}\mathbf{b}_{\text{eff}}+\mathbf{b}_{\text{fl}}(\bar{t})\big{]}\big{)}.$$

$$\begin{split}&\langle b_{\text{fl},j}(\bar{t})\rangle=0,\\ &\langle b_{\text{fl},j}(\bar{t})b_{\text{fl},k}(\bar{t}^{\prime})\rangle=2\Big{(}\frac{\alpha^{2}}{1+\alpha^{2}}\frac{k_{B}T}{J_{1}}\Big{)}\delta_{jk}\delta(\bar{t}-\bar{t}^{\prime}),\\ \end{split}$$

Full text

I Introduction

Magnetic skyrmions are solitonic, particlelike, topologically nontrivial magnetic textures. The existence of (meta)stable skyrmionic solutions in a system requires the introduction of a characteristic length-scale via competing interactions. This is typically achieved in noncentrosymmetric magnets with the Dzyaloshinskii-Moriya interaction (DMI) Dzyaloshinskii (1958); Moriya (1960). That particular type of solutions, commonly referred to as chiral skyrmions, were theoretically investigated in the 1990s Bogdanov and Yablonskii (1989); Bogdanov and Hubert (1994). In magnets with inversion symmetry, skyrmions can be stabilized via dipolar interactions (skyrmion bubbles) Chikazumi and Graham (2009); Malozemoff and Slonczewski (1979); Yu et al. (2012); Finazzi et al. (2013); Nagao et al. (2013); Koshibae and Nagaosa (2016), as well as frustrated exchange couplings Okubo et al. (2012); Nagaosa and Tokura (2013); Lin and Hayami (2016); Rózsa et al. (2017); Zhang et al. (2017). Systems with frustrated exchange are particularly interesting, as they allow many different kinds of topological defects to coexist. The frustrated $J_{1}-J_{2}$ system on the hexagonal lattice was extensively investigated in Ref. Leonov and Mostovoy (2015) and, in general, the frustrated hexagonal lattice is the most commonly studied Okubo et al. (2012); Dupé et al. (2014); von Malottki et al. (2017); Rózsa et al. (2017); Ritzmann et al. (2018). Skyrmions decouple from the lattice in the long wavelength limit and therefore, in that case, the choice of the lattice does not matter Lin and Hayami (2016). Collapse mechanisms of skyrmions and antiskyrmions in chiral thin films with frustrated exchange were previously investigated Bessarab et al. (2018); von Malottki et al. (2017), and appear similar to the isotropic-type of collapse calculated in non-frustrated chiral systems Bessarab et al. (2015); Lobanov et al. (2016); Uzdin et al. (2017); Cortés-Ortuño et al. (2017); Stosic et al. (2017); Desplat et al. (2018). A third mechanism involving the injection of singularity with opposite topological charge was reported in Cortés-Ortuño et al. (2017), but the path appears to go over a higher order saddle point. Additionally, this mechanism seems to emerge for larger skyrmion sizes, yet it is perhaps the closest to the type of mechanisms we report in this work.

In this article, we look at paths to annihilation of first- and second-order skyrmions and antiskyrmions on the frustrated $J_{1}-J_{2}-J_{3}$ square lattice previously studied by Lin et al. Lin and Hayami (2016). We relax skyrmionic solutions spanning over only a few sites in diameter that do not exhibit translational invariance on the lattice, and instead experience pinning at particular lattice points. In the present system, the solutions are too small to decouple from the underlying lattice. In this context, we obtain new types of collapse mechanisms which have, to date, not been reported. Interfacial DMI is added to break the invariance with respect to the rotation of helicity (i.e. the “spin” degree of freedom Lin and Hayami (2016)), as well as the degeneracy between skyrmions and antiskyrmions.

The paper is organized as follows. In Sec. II, we present the model Hamiltonian and some useful definitions, as well as a brief zoology of some of the different topological defects that can be stabilized in the system. In Sec. III, we firstly give regions of existence of metastable first- and second-order skyrmions and antiskyrmions in parameter space. We then present different mechanisms by which they annihilate, and we confirm them by means of direct Langevin simulations. We subsequently compute attempt frequencies associated with the different mechanisms via transition state theory calculations. Finally, the results are discussed in Sec. IV.

II Model

We simulate a bidimensional square lattice of $N$ magnetic spins $\{\mathbf{m}_{i}\}$ with a constant magnitude that we set to unity. We use the classical Heisenberg model Hamiltonian,

[TABLE]

where $J_{ij}$ is the strength of the isotropic exchange coupling extended to third nearest neighbors ($J_{1}-J_{2}-J_{3}$), $\mathbf{D}_{ij}$ is the Dzyaloshinskii vector restricted to first nearest neighbors, $K$ is the perpendicular, uniaxial anisotropy constant, and $B_{z}$ is the perpendicular applied magnetic field. The DMI is interfacial and favors Néel skyrmions, in which case the Dzyaloshinskii vector is defined as $\mathbf{D_{ij}}=D\mathbf{r}_{ij}\times\mathbf{e}_{z}$, where $\mathbf{r}_{ij}$ is the in-plane direction between sites $i$ and $j$ Thiaville et al. (2012). We introduce the reduced parameters: $j_{\alpha}=J_{\alpha}/J_{1},(\alpha=1,2,3)$ ; $d=\lvert\mathbf{D}_{ij}\lvert/J_{1}$; $k=K/J_{1}$ ; $b=B_{z}/J_{1}$. We set $k=0.1$ for the rest of this study. The isotropic exchange coupling parameters allow a spin-spiral ground state to be realized in the absence of other interactions, and are the following Lin and Hayami (2016): $j_{2}=-0.35$ ; $j_{3}=-0.15$. Additional parameter values are given in Appendix A. We define the topological charge Belavin (1975),

[TABLE]

in which $\rho_{s}$ is the toplogical charge density defined as,

[TABLE]

The quantity $n$ in Eq. (2) is the winding number, or vorticity, and counts the number of times the spin configuration wraps the unit sphere $S^{2}$. Solutions with a positive winding number are called skyrmions, and solutions with a negative winding number are called antiskyrmions.Contrary to chiral skyrmions stabilized by DMI where the type of the DMI favors one particular class of solutions, exotic spin textures with various winding numbers and helicities can coexist in frustrated magnets Leonov and Mostovoy (2015); Lin and Hayami (2016); Rózsa et al. (2017); Zhang et al. (2017) [Fig. 1]. Nonetheless, the possibility to stabilize structures with an arbitrary integer topological charge in chiral magnets (“skyrmionic sacks”) has recently been explored Rybakov and Kiselev (2019). In the absence of inversion symmetry breaking, skyrmions [Figs. 1(b), 1(c)] and antiskyrmions [Fig. 1(a)] with opposite winding numbers are degenerate in energy. Additionally, magnets with inversion symmetry possess a spin degree of freedom (where “spin” refers to the spin of the skyrmion as a particle), associated with a rotation of helicity. Bloch-type skyrmions [Fig. 1(b)], Néel-type skyrmions [Fig. 1(c)], and all intermeditate helicity states are therefore degenerate in energy Lin and Hayami (2016).

III Results

III.1 Coexistence of skyrmion and antiskyrmion solutions

We introduce a weak DMI to the system in order to break invariance with respect to the rotation of helicity, as well as the degeneracy between skyrmions and antiskyrmions. We look for a set of parameters that allow metastable skyrmion and antiskyrmion solutions to co-exist. To do so, we initialize the system close to a(n) (anti)skyrmion state and run an iterative energy minimization procedure similar to the one described in Aharoni (2000). If a(n) (anti)skyrmion solution exists, it is relaxed. We vary the perpendicular applied field $b$ in $[0,0.4]$, and the DMI coupling constant $d$ in $[0,0.07]$. We obtain the diagrams of Fig. 2 for skyrmions [Fig. 2(a)] and antiskyrmions [Fig. 2(b)].

At low field or high DMI, a spiral state is relaxed. Above a critical applied field, the spin-spiral state becomes the fully polarized FM state via a second-order phase transition Lin and Hayami (2016). From that state, metastable skyrmion solutions can be stabilized. In the presence of a perpendicular easy-axis anisotropy, the skyrmion lattice (SkX) becomes thermodynamically stable at intermediate applied field Lin and Hayami (2016). In this region, isolated skyrmions are more favorable than the FM state. They are metastable in the sense that they are local energy solutions, but not the ground state. At finite temperature, the isolated skyrmion state is rapidly destroyed by the nucleation of many topological defects in its vicinity. By contrast, the metastable solutions we are interested in are isolated skyrmionic defects as excitations of the FM ground state, which we mark as orange dots on Fig. 2. Since the symmetry of the DMI favors Néel skyrmions, degeneracies are lifted, and antiskyrmions become progressively less stable than their skyrmion counterparts as the DMI increases. They also relax diagonally on the square lattice, such that their Bloch-type axes are along the first-neighbor axes, which are the only ones that are DMI-coupled. This is due to the fact that one of the Néel-type axes of the antiskyrmion exhibits the opposite chirality to the one favored by the DMI. This introduces frustration in the orientation of the antiskyrmion solutions, which is not present in first-order, radially-symmetric skyrmion solutions.

Finally, Figs. 2(c) and 2(d) show the range of existence of second-order skyrmion and antiskyrmion solutions, respectively. They exist as stable solutions at low $d$ in the skyrmion lattice phase. By increasing the DMI at constant applied field, a bound skyrmion pair is relaxed instead of the second-order skyrmion, but the second-order antiskyrmion solution persists. Metastable solutions are found within a small window of the FM phase at sufficiently low field and DMI. We find second-order skyrmion and antiskyrmions to be quasi-degenerate in energy whenever they both exist, with small differences most likely caused by the underlying lattice. However, the range of existence of the antiskyrmion appears to be wider, as the first-order antiskyrmions are unfavored by the DMI.

III.2 Paths to annihilation

Transition mechanisms correspond to minimum energy paths (MEPs) on the $2N$-dimensional energy surface. We compute them via geodesics nudged elastic bands (GNEB) Bessarab et al. (2015) calculations with a climbing image (CI-GNEB) Henkelman et al. (2000) scheme to precisely identify the first-order saddle point (SP) on the path. Successive states of the system along the reaction coordinate are referred to as images. For each mechanism, we plot the spin maps at selected images and the corresponding topological charge density as defined in Eq. (3) [Figs. 3, 5, 7]. The total energy along the path is plotted in Figs. 4, 6, and 8 in units of the isotropic exchange coupling constant between first neighbors, $J_{1}$, as well as the total topological charge at each image, as defined in Eq. (2).

Skyrmion.

The first set of parameters we examine is $(b,d)=(0.2,0.03)$. It is a region of parameter space where both metastable skyrmions and antiskyrmions solutions exist. In this context, we observe an annihilation mechanism which is different from the usually reported isotropic type of collapse [Fig. 3(a)]. One half of the skyrmion unfolds, and a half-antiskyrmion, or antimeron, is nucleated in its place. The injection of the opposite topological charge in the system constitutes the first-order saddle point on the transition path [Fig. 4]. This state possesses four different realizations, equivalent to rigid $\pi/2$ rotations of the whole sample. The remaining meron and antimeron then naturally annihilate. We found that the isotropic type of collapse also exists in this system, but it involves a higher order SP and therefore constitutes a less probable route. The second set of parameters is $(b,d)=(0.3,0.07)$. In the space of parameters, it is situated just above the limit at which antiskyrmion solutions no longer exist. In this case, the skyrmion undergoes a regular collapse [Fig. 3(b)] involving a first-order SP (image 6), but the core displays a weak asymmetry in a way that is reminiscent of Fig. 3(a). Since we are close to the region of antiskyrmions’ metastability, this mechanism can be considered as intermediate between antimeron nucleation and isotropic collapse. Finally, we select a region of parameter space with high DMI, far from the metastability region of antiskrymions: $(b,d)=(0.7,0.2)$. As expected, we obtain a perfectly isotropic collapse [Fig. 3(c)] without the injection of an opposite topological charge. In all cases, the SP configuration corresponds to an almost zero topological charge, and the variation of the topological charge along the path is more abrupt with increasing DMI.

Antiskyrmion.

We revert back to the previous set of parameters, $(b,d)=(0.2,0.03)$, which allows metastable antiskyrmions to coexist with skyrmions. We look at the path to annihilation of an antiskyrmion and we relax two different mechanisms passing through a first-order SP [Figs. 5(b) and 5(a)]. In both cases, the annihilation process begins with the rotation of the antiskyrmion, such that its Néel-type axes are aligned along the first-neighbor axes (images 1-7). This configuration constitutes an energy maximum as it is unfavored by the DMI, and the antiskyrmion can reach this state either as a full antiskyrmion with full integer topological charge (SP2 on Figs. 5(a) and 6(a)), or the spins may begin to unwind, which leads to a decrease in the topological charge (SP1 on Figs. 5(b) and 6(b)). Past the SP, similarly to the annihilation mechanism of its skyrmion counterpart, the topological charge drops as a meron is nucleated and annihilates with the remaining antimeron. SP1 exists in four possible realizations ($\pi/2$-rotations), while SP2 exists in two ($\pi$-rotation). It is interesting to note that the fact that the antiskyrmion is frustrated in its orientation on the lattice suffices to distinguish between a metastable solution and a SP. Here, the coupling to the spin lattice entails that the antiskyrmion, as well as the second-order skyrmion and all non-radially symmetric solutions behave as strongly correlated spin-ice Yavors’kii et al. (2008).

Second-order skyrmion.

We select paramaters that allow both the first- and second-order skyrmions to be metastable: $(b,d)=(0.14,0.005)$. We initially set the final state of the GNEB calculation to a Néel skyrmion, and we obtain the mechanism in Fig. 7(a), with the corresponding energy profile shown on Fig. 8(a). The second-order skyrmion first rotates on the lattice and reaches the SP (image 3). The nucleation of an antimeron along the Néel axis follows (image 4) and the energy drops, accompanied by an abrupt change in the topological charge. The antimeron subsequently annihilates with part of the second-order skyrmion, leaving a first-order Néel skyrmion in its place. Another route to annihilation for the second-order skyrmion corresponds to its division into a bound skyrmion pair, as shown on Fig. 7(b). The bound skyrmion pair is a metastable solution, since the interaction potential of skyrmions in frustrated magnets is non-monotonic as a function of distance, and is found in turn positive (repulsive) and negative (binding) Leonov and Mostovoy (2015); Lin and Hayami (2016). This mechanism does not involve any significant change in the topological charge of the total system [Fig. 8(b)], but instead a redistribution of the topological charge density. Both of the above mechanisms possess four realizations, corresponding to $\pi/4$ rigid rotations of the sample. As for the fate of the skyrmion pair, the binding potential is quite shallow compared to the activation barrier for annihilation (from the results of Lin and Hayami (2016), we can estimate the unbinding barrier to be $\sim 10^{-4}$ $J_{1}$). Therefore, the skyrmions will in most cases separate before annihilating individually.

Similar processes where reported in Zhang et al. (2017) from dynamics simulations, in which the second-order skyrmion division was triggered by both a driving current at zero-temperature, and thermal fluctutations at finite temperature. The finite temperature simulations also showed the decay of the second-order skyrmion into a first-order skyrmion.

III.3 Langevin simulations

In order to confirm the previously calculated MEPs, we perform direct atomistic Langevin simulations at low temperature. Details on the method can be found in Appendix B. We observe annihilation mechanisms for which we show snapshots of the spin configuration in Fig. 9. For the antiskyrmion at $T=$ 80 K and $(b,d)=(0.2,0.03)$ [Fig. 9a], we observe the meron-nucleation type of annihilation that we reported on Fig. 5. Since we have short lifetimes within the nanosecond (ns) timescale, we can estimate the antiskyrmion’s stability directly from Langevin simulations. We use the topological charge of the system to track the collapse of the antiskyrmion. We record 390 collapses and obtain an average lifetime $\tau_{\text{DL}}(80$K$)=0.49$ ns, with a standard deviation of the order of the lifetime, $\sigma=0.46$ ns.

For the second-order skyrmion and antiskyrmion, as they overall exhibit lower activation energies, we perform simulations at $T=50$ K, and $(b,d)=(0.14,0.005)$. The division of the second-order skyrmion into a bound skyrmion pair seems to occur along at least two different paths shown on Fig. 9b and 9c, in which the latter corresponds to the MEP on Fig. 7(b). In Fig. 9b, the division occurs through excitation of one half of the second-order skyrmion, which eventually leads to the separation of the two halves. On the other hand, in Fig. 9c, we witness a more symmetric kind of division. For the same set of parameters, we also observe the decay into a single skyrmion [Fig. 9d, which corresponds to the path on Fig. 7(a)]. Finally, on Fig. 9e, we show the division of a second-order antiskyrmion into a pair of first-order antiskyrmions.

III.4 Annihilation rates

The rate of thermally activated processes with an internal energy barrier $\Delta E$ can be described by the Arrhenius law, $f(T)=f_{0}e^{-\Delta E/k_{B}T}.$ We use a form of Langer’s theory for the decay of metastable states Langer (1969) adapted to magnetic spin systems Coffey et al. (2001); Fiedler et al. (2012); Desplat et al. (2018) to compute individual rate prefactors $f_{0}$ for each of the mechanisms described previously, while the CI-GNEB scheme directly gives us a precise value for the activation energy. More details on the methods in this paragraph are given in Ref. Desplat et al., 2018. The results are gathered in Table 1. To account for all equivalent realizations of a given mechanism (rigid rotations of the spin configurations at the SP), the prefactor $f_{0}$ is multiplied by a factor two or four when needed, in accordance with Sec. III.2.

Based on these results, we can estimate the lifetime of the antiskyrmion calculated with Langer’s theory (with the subscript TST, for transition state theory), and compare it with the average lifetime previously obtained from direct Langevin simulations. Assuming independent processes, we have $f_{\text{TST}}(80$K$)=f_{01}e^{-\beta\Delta E_{1}}+f_{02}e^{-\beta\Delta E_{2}}=\tau_{\text{TST}}^{-1}(80$K$)=2.08$ GHz and $\tau_{\text{TST}}(80$K$)$ = 0.48 ns. Although we are not completely within the scope of Langer’s theory ($\beta\Delta E_{1}=0.3$, as opposed to the recommended $\beta\Delta E\geq 5$ Coffey et al. (2001)), this shows an excellent agreement with the result from direct Langevin, $\tau_{\text{DL}}(80$K$)$ = 0.49 ns.

In the case of the second-order skyrmion’s decay into a first-order skyrmion, we find a Goldstone mode at both the metastable and transition states, which do not seem to clearly correspond to a translational or rotational invariance, hence we cannot give a definite value for the attempt frequency in this case.

IV Discussion and Conclusion

IV.1 Topological transitions: taking a ball out of a net

Drawing inspiration from Ref. Rohart et al., 2016, we project images along the MEPs onto the space of configurations – the unit sphere [Fig. 10]. Vertices correspond to the tip of the magnetic vectors with their origin in the center of the sphere, and the edges represent the exchange coupling between first neighbors. The ferromagnetic background (spins along $+z$) is found at the north pole of the sphere, while the core of a skyrmion typically points along the south pole (along $-z$). The view is set just below the south pole, looking axially towards the $+z$ direction. Annihilating a topological defect to recover the uniformly magnetized state means bringing the mesh back to the north pole by moving the vertices on the sphere. This can be thought of as taking a ball (solid blue sphere) out of a net (dark blue mesh) by deforming the net Rohart et al. (2016). Annhilations that occur via the injection of the opposite topological charge [Figs. 10(a) and 10(c)] involve a rearranging of the spin-mesh in a way that allows the entire net to be removed along a single direction. By contrast, an isotropic collapse consists in progressively spreading the mesh open by an equal amount along all directions [Fig. 10(b)]. Finally, a second-order skyrmion corresponds to a kind of double net [Fig. 10(d)]. Its decay into a first-order skyrmion via antimeron injection consists in the removal of one of the nets, again along a given direction, so that only a single net remains.

IV.2 Conclusion

In regions of parameter space for which both metastable skyrmion and antiskyrmion solutions can be realized, the most probable paths to annihilation – that is, the ones that go over a first-order saddle point – seem to the paths involving the injection of the opposite topological charge into the system, in the form of the nucleation of merons and antimerons. The injection of the opposite charge usually happens just after the saddle point configuration, or at the saddle point in the case of the skyrmion. The reason these peculiar paths exist in the present system seem to be the interplay of the frustrated exchange, and the small sizes of the topological defects (small number of magnetic sites involved in the spin textures) that can be stabilized, which is directly linked with the period of the spin spiral ground state determined by the exchange frustration. Alternatively, the division of the second-order skyrmion into a bound skyrmion pair involves no change in the total toplogical charge of the system, only a redistribution of the charge density, and is also a valid path.

Overall, the skyrmionic solutions in this system are found to be rather unstable, as they exhibit low internal energy barriers, combined with many possible paths to collapse with attempt frequencies in the range of several GigaHertz (GHz). The existence of many possible paths is a direct consequence of the exchange frustration, which drastically complexifies the energy landscape and is responsible for the emergence of many (meta)stable and saddle point states, as well as many possible paths connecting them. The low stability allowed for a direct comparison of the average lifetime of the antiskyrmion at $80$ K computed from direct Langevin simulations, with transition state theory calculations, and we obtained a very good agreement. Nevertheless, in this context, the direct use of transition state theory to compile an overall reliable lifetime for any given structure seems ill-advised, as it is difficult to account for all possible mechanisms. As for the low activation energies and not particularly low attempt frequencies, we can once more relate them to the skyrmions being very small, and possessing only a few internal modes Desplat et al. (2018); von Malottki et al. (2019).

Acknowledgements.

This work was supported by the Horizon 2020 Framework Programme of the European Commission, under Grant agreement No. 665095 (MAGicSky), and the Agence Nationale de la Recherche under Contract No. ANR-17-CE24-0025 (TOPSKY). The authors thank B. Dupé for enlightening discussions, and D. Suess for previous help with code testing.

Appendix A Simulation parameters

We use the following parameters for a two-dimensional monolayer of $40\times 40$ simulated sites: the isotropic exchange constant between first neighbors is set to $J_{1}=1.6\times 10^{-20}$ J ($\sim 100$ meV) with lattice constant $a=1$ nm and saturation magnetization $M_{s}=1.1a^{3}$ MA.m2 Thiaville et al. (2012). The gyromagnetic ratio is that of the free electron , $\gamma=1.76\times 10^{11}$ rad s*-1* T*-1*, and the dimensionless damping factor is $\alpha=0.5$ Thiaville et al. (2012) (Intermediate-to-high damping regime of Langer’s theory).

Appendix B Langevin dynamics

The dynamics of the magnetic spin system $\{\mathbf{m}_{i}\}$, $i=1\dots N$, is governed by the set of coupled, dimensionless, stochastic Landau-Lifshitz-Gilbert (LLG) equations García-Palacios and Lázaro (1998):

[TABLE]

All the quantities are in units of the isotropic exchange coupling constant between first neighbors $J_{1}$, such that in the above expression, $\mathbf{b}_{\text{eff}}=-\dfrac{1}{J_{1}}\dfrac{\partial E}{\partial\mathbf{m}_{i}}$ is the local reduced effective field, and $\mathbf{b}_{\text{fl}}$ is a stochastic fluctuating field in the form of white noise, which accounts for fluctuations of the orientation of $\mathbf{m}_{i}$ caused by interactions with microscopic degrees of freedom of the environment. It is assumed to be Gaussian distributed with the following statistical properties, in agreement with the fluctuation-dissipation theorem García-Palacios and Lázaro (1998):

[TABLE]

where $j,k$ are Cartesian indices, $\langle\rangle$ denotes an average over many realizations of the fluctuating field, $\delta(\bar{t}-\bar{t}^{\prime})$ is the Dirac distribution, and $\delta_{jk}$ is the Kronecker symbol. The reduced time $\bar{t}$ is linked to physical time $t$ via $t=\dfrac{M_{s}}{\alpha\gamma J_{1}}\bar{t}$, and we maintain $d\bar{t}=0.001$. We use the stochastic Heun scheme García-Palacios and Lázaro (1998), which converges to the solution of (4) when the multiplicative noise is interpreted in the sense of Stratonovich. Additionally, periodic boundary conditions are applied to prevent the skyrmions from escaping through the boundaries.