Lifetime of coexisting sub-10 nm zero-field skyrmions and antiskyrmions

Abstract

Magnetic skyrmions have raised high hopes for future spintronic devices. For many applications, it would be of great advantage to have more than one metastable particle-like texture available. The coexistence of skyrmions and antiskyrmions has been proposed in inversion-symmetric magnets with exchange frustration. However, so far only model systems have been studied and the lifetime of coexisting metastable topological spin structures has not been obtained. Here, we predict that skyrmions and antiskyrmions with diameters below 10 nm can coexist at zero magnetic field in a Rh/Co bilayer on the Ir(111) surface—an experimentally feasible system. We show that the lifetimes of metastable skyrmions and antiskyrmions in the ferromagnetic ground state are above one hour for temperatures up to 75 and 48 K, respectively. The entropic contribution to the nucleation and annihilation rates differs for skyrmions and antiskyrmions. This opens the route to the thermally activated creation of coexisting skyrmions and antiskyrmions in frustrated magnets with Dzyaloshinskii–Moriya interaction.

Introduction

The observation of magnetic skyrmions^1,2,3,4 has triggered the vision to use these magnetic entities in applications^5,6 ranging from data storage⁷ and logic devices to neuromorphic^8,9,10 or quantum computing¹¹. The Dzyaloshinskii–Moriya interaction (DMI) which occurs due to spin–orbit coupling in materials with broken inversion symmetry permits the stabilization of skyrmions at transition–metal interfaces such as ultrathin films at surfaces^4,12,13,14 or magnetic multilayers^15,16,17,18. In these systems, the variation of chemical composition, stacking sequence, and interface structure opens the possibility to tune the magnetic interactions and thereby skyrmion properties. This approach led to the realization of skyrmions with diameters down to 30 nm at room temperature and the demonstration of current-induced skyrmion motion^15,16,17,19.

A key limitation of using only skyrmions as information carriers, e.g., in racetrack memory devices^20,21, is that the classical bit can only be stored by the presence (“1”) or absence (“0”) of a skyrmion, which could be impractical since precise control over the skyrmion position is required in this case. It is highly desirable to have another kind of localized, metastable spin structure which can act as a bit in a binary logic^22,23,24. The use of distinct magnetic solitons such as degenerate skyrmions (Sk) and antiskyrmions (Ask) has been also proposed for quantum computing¹¹. Therefore, systems hosting several coexisting metastable spin structures of different type are not only fundamentally interesting but also relevant for applications²⁵.

A zoo of topological spin structures beyond skyrmions has been theoretically predicted^{22,26,27,28,29,30} and some of them have already been observed in experiments (see e.g., the review by Göbel et al.²³). In tetragonal inverse Heusler compounds antiskyrmions with diameters of about 150 nm can be stabilized in an external magnetic field due to anisotropic DMI and were observed at room temperature by Lorentz transmission electron microscopy³¹. The long lifetime of antiskyrmions in these bulk materials was reproduced based on transition-state theory³². In such materials with anisotropic DMI, which can also be realized at interfaces²² and in two-dimensional magnets^24,33,34, antiskyrmions are energetically preferred over skyrmions. By use of in-plane magnetic fields, it was experimentally possible to transform antiskyrmions into elliptical skyrmions in this class of materials^35,36,37. The observed coexistence phase of two magnetic objects with opposite topological charge Q arises in these systems from the competition of anisotropic Dzyaloshinskii–Moriya and dipole–dipole interactions³⁸. In ferrimagnetic multilayers, even purely dipolar-interaction stabilized skyrmions and antiskyrmions with diameters of about 200 nm have been observed³⁹. It has also been shown that stable skyrmion and antiskyrmion solutions of the continuum model coexist in a limited range of magnetic field and/or anisotropy⁴⁰ even without magnetostatic interaction. Recently, the creation and annihilation of such skyrmion-antiskyrmion pairs with diameters of about 150 nm have been observed in thin films of the cubic chiral magnet FeGe⁴¹. However, coexisting nanoscale skyrmions and antiskyrmions have not been reported so far.

In order to realize coexisting sub-10 nm skyrmions and antiskyrmions which are ideally suited to outperform competing future technologies⁶ the interplay of other magnetic interactions can be exploited. It has been shown that frustrated exchange^{26,42,43,44,45} as well as higher-order exchange^46,47,48 interactions can stabilize isolated Sk and Ask in the field-polarized phase of a material, i.e., at finite magnetic field, even without DMI. However, so far the lifetime of coexisting skyrmions and antiskyrmions has neither been obtained in experiments nor in theoretical work. Therefore, it is not clear whether coexisting topological spin structures in a real material system exhibit a sufficient stability to be observed experimentally and can be used for applications. In addition, the relative probability of coexisting metastable topological spin structures has not been reported yet.

Here, we predict the coexistence of isolated skyrmions and antiskyrmions with diameters below 10 nm at zero magnetic field in an atomic Rh/Co bilayer on the Ir(111) surface—an experimentally feasible ultrathin film system^14,49. Due to strong exchange frustration, both types of topological spin structures are stable down to zero magnetic field, i.e., in the ferromagnetic ground state of the system. We calculate the lifetimes of metastable skyrmions and antiskyrmions using transition-state theory based on an atomistic spin model with all interactions parameterized from density functional theory (DFT). Lifetimes above one hour are obtained for skyrmions and antiskyrmions in Rh/Co/Ir(111) at temperatures of up to 75 K and 48 K, respectively. The lifetimes of antiskyrmions are even longer than those of skyrmions in the prototypical skyrmion system Pd/Fe/Ir(111)^{4,50,51,52,53}.

Further we study the effect of DMI on skyrmions and antiskyrmions in a frustrated magnet, in particular, on the energy barriers and pre-exponential factors in the Arrhenius law describing kinetics of skyrmion and antiskyrmion annihilation and nucleation. We demonstrate that the difference in the energy spectra for skyrmions and antiskyrmions in the presence of DMI leads to distinct entropic contributions to their transition rates. Surprisingly, even in the regime of significant DMI the creation of coexisting skyrmions and antiskyrmions becomes possible via local heating.

Results

Zero-field skyrmions and antiskyrmions

We study skyrmions and antiskyrmions in an atomic Rh/Co bilayer on the Ir(111) surface using the corresponding DFT-parameterized atomistic spin model¹⁴ given by

$$\begin{array}{l}E=-\mathop{\sum}\limits_{i,j}{J}_{ij}({{{{\bf{m}}}}}_{i}\cdot {{{{\bf{m}}}}}_{j})-\mathop{\sum}\limits_{i,j}{{{{\bf{D}}}}}_{ij}\cdot ({{{{\bf{m}}}}}_{i}\times {{{{\bf{m}}}}}_{j})\\ \qquad\,\,-\mathop{\sum}\limits_{i}K{({{{{\bf{m}}}}}_{i}\cdot {\hat{{{{\bf{e}}}}}}_{\perp })}^{2}-\mathop{\sum}\limits_{i}M({{{{\bf{m}}}}}_{i}\cdot {{{\mathbf{{\cal B}}}}}),\end{array}$$

(1)

which takes into account the Heisenberg exchange interaction and the DMI between normalized magnetic moments m_i and m_j at lattice sites i and j of the Co layer as well as the uniaxial magnetocrystalline anisotropy and the Zeemann energy due to an external magnetic field ${{{\bf{{{{\mathcal{B}}}}}}}}$ applied perpendicular to the film. The exchange constants, J_ij, the strength of the DMI, D_ij, the magnetocrystalline anisotropy constant, K, and the size of the magnetic moment per site, M, have been obtained based on DFT for both fcc- and hcp-stacking of the Rh layer¹⁴ (see “Methods” for details and Supplementary Note 1 and Supplementary Fig. 1). We do not explicitly include dipole–dipole interactions. However, for ultrathin films, this energy term is very small—on the order of 0.1 meV/atom—and it can be effectively included into the magnetocrystalline anisotropy energy^54,55.

Isolated skyrmions and antiskyrmions (Fig. 1a, b) are initialized within the ferromagnetic ground state (for simulation details, such as initial skyrmion profiles and optimizations, see “Methods”). The radii of the skyrmions and antiskyrmions are obtained from the relaxed profiles using the definition of Bogdanov et al.⁵⁶. Even in the absence of magnetic field both types of spin structures are stable in our simulations and exhibit diameters well below 10 nm (Fig. 1e). Note that nanoscale skyrmions of this size have already been experimentally observed in Rh/Co/Ir(111) at zero magnetic field^14,49.

Fig. 1: Zero-field skyrmions and antiskyrmions in fcc-Rh/Co/Ir(111).

a Isolated skyrmion (Sk) and b isolated antiskyrmion (Ask) at ${{{\mathcal{B}}}}=0$ T. c Chimera (skyrmion) saddle point (SP) and d antiskyrmion SP. a–d show equally sized parts of the much larger respective simulation boxes. e Radius of skyrmions (Sk, red) and antiskyrmions (Ask, blue) as a function of the magnetic field applied perpendicular to the film. f Energy barriers ΔE for skyrmion (Sk → FM, red) and antiskyrmion (Ask → FM, blue) annihilation. For skyrmions, there is a crossover from the Chimera collapse (orange) at low fields to the radial collapse (red) mechanism at larger magnetic fields. Minimum energy path (MEP) for g the collapse of a skyrmion (initial state) via the Chimera mechanism into the FM (final) state and h the antiskyrmion collapse via the radially symmetric mechanism at ${{{\mathcal{B}}}}=0$ T. The total energy (tot), with respect to the FM ground state, is decomposed into the contributions from the exchange interaction (exc), the DMI, and the magnetocrystalline anisotropy energy (MAE).

Antiskyrmions are by about 20–50% smaller than skyrmions in the entire range of magnetic fields due to the effect of DMI. In order to study the stability of both states, we calculate the energy barriers protecting them against collapse into the ferromagnetic ground state using the geodesic nudged elastic band (GNEB) method⁵⁷ (see “Methods”). The pairwise annihilation of skyrmions and antiskyrmions does not only depend on intrinsic but also dynamical properties of states, such as their spatial separation⁴², and cannot be sufficiently captured by the static methods used in this study. Therefore, we restrict our simulations to isolated states and neglect combined skyrmion-antiskyrmion processes.

As presented in Fig. 1f, we find that the skyrmion annihilation energy barrier ΔE^Sk→FM is significantly larger than that for antiskyrmions, ΔE^Ask→FM, for all magnetic fields. Surprisingly, the antiskyrmion barrier of about 150 meV at ${{{\mathcal{B}}}}=0$ T is of the same order of magnitude as the barrier for the skyrmion collapse reported for Pd/Fe/Ir(111)⁴³, a well-studied system for which experiment^{4,12,51,53,58} and first-principles based theory^43,50,53,59 are in good agreement. This suggests that antiskyrmions in Rh/Co/Ir(111) could possess a similar stability as that observed for skyrmions in Pd/Fe/Ir(111) at low temperatures^{4,12,51,53,58}.

The GNEB calculations also provide information about the mechanisms of magnetic transitions. For skyrmions, we find the recently discovered Chimera collapse mechanism^14,53 at ${{{\mathcal{B}}}}\, < \,1.26$ T denoted as Sk (chim.) in Fig. 1f. At larger fields, the radially symmetric transition, denoted as Sk (rad.), is energetically favored over the Chimera transition. From the minimum energy path for the Chimera collapse at ${{{\mathcal{B}}}}=0$ T (Fig. 1g), we see that the energy of the skyrmion state (Fig. 1a) itself is lowered by the DMI while the energy costs of the saddle point (Fig. 1c), which determines the energy barrier, stems mostly from the frustrated exchange interaction.

The situation is different for antiskyrmions (Fig. 1b). There, we find only the radial shrinking annihilation mechanism in the entire investigated range of magnetic fields in which the initial state (Fig. 1b) collapses into the ferromagnetic ground state via a saddle point (Fig. 1d) exhibiting a core with in-plane pointing spins. The minimum energy path at ${{{\mathcal{B}}}}=0$ T (Fig. 1h) shows that the energy barrier is dominated by the frustrated exchange interaction. In contrast to skyrmions, the influence of the DMI on the antiskyrmion collapse is very small.

The stability of metastable topological spin structures can be quantified by their mean lifetime, τ. Here we use transition-state theory in harmonic approximation to the energy of states⁶⁰ where the expression for the lifetime takes the form of an Arrhenius law

$$\tau ={\tau }_{0}\exp (\Delta E/{k}_{{{{\rm{B}}}}}T)$$

(2)

where ΔE is the energy barrier of the transition and τ₀ a pre-exponential factor which can be explicitly calculated (see “Methods”).

The obtained lifetimes of skyrmions and antiskyrmions in Rh/Co/Ir(111) at zero magnetic field are displayed in Fig. 2 for both fcc- and hcp-stacking of the Rh layer in a logarithmic plot as functions of inverse temperature. Considering the logarithmic scaling of the ordinate

$$\ln \tau =\ln {\tau }_{0}+\Delta E/{k}_{B}T,$$

(3)

the slope of the curves thus gives the energy barrier, while the offset at T⁻¹ → 0 yields the pre-exponential factor. As expected from the energy barriers, the lifetime of skyrmions is much larger than that of antiskyrmions at low temperatures. Nevertheless, the lifetime of antiskyrmions is also above one hour below about 48 K and 43 K in hcp-Rh/Co and fcc-Rh/Co on Ir(111), respectively. For skyrmions, the one-hour lifetime is already exceeded below about 75 K and 66 K for the hcp and fcc stacking of Rh, respectively. Note, that at large temperatures there is a crossing of skyrmion and antiskyrmion lifetime due to different pre-exponential factors.

Fig. 2: Lifetime of skyrmions and antiskyrmions in Rh/Co/Ir(111).

Lifetime of skyrmions (Sk, orange lines) and antiskyrmions (Ask, blue lines) in hcp-Rh/Co/Ir(111) (solid lines) and fcc-Rh/Co/Ir(111) (dashed lines) at ${{{\mathcal{B}}}}=0\,T$ shown over inverse temperature. For comparison, the lifetimes are also given for skyrmions and antiskyrmions in Pd/Fe/Ir(111) (dashed-dotted lines) at ${{{\mathcal{B}}}}=3.9\,T$, i.e., just above the critical field of the field-polarized phase (see methods for details).

For comparison, we have also displayed in Fig. 2 the calculated lifetimes for isolated metastable skyrmions in the field-polarized phase of Pd/Fe/Ir(111)⁵² at ${{{\mathcal{B}}}}=3.9\,$T. This value has been chosen because at this field the annihilation barriers for Sk in Pd/Fe/Ir(111) and Ask in Rh/Co/Ir(111) are similar, despite the fact, that the Ask are significantly smaller (see Supplementary Note 2 and Supplementary Fig. 2). Remarkably, the lifetime of Ask in Rh/Co/Ir(111) is even larger than that of Sk in Pd/Fe/Ir(111). Note, that the calculated Sk lifetimes^52,61 are in good agreement with experiments for Pd/Fe/Ir(111)^51,53. We have also calculated the lifetime of antiskyrmions in Pd/Fe/Ir(111) at ${{{\mathcal{B}}}}=3.9\,$T. However, even at the lowest considered temperatures, their lifetime is very short (Fig. 2). This can be explained based on the much smaller energy barrier as apparent from the small slope of the curve. The origin of this difference between the two film systems is, in addition to the applied magnetic field, the significantly enhanced exchange frustration in Rh/Co/Ir(111)¹⁴ (see Supplementary Note 2 and Supplementary Fig. 2).

DMI vs. frustrated exchange

To understand the impact of the DMI on skyrmions and antiskyrmions in a two-dimensional system and its competition with frustrated exchange interactions, we vary the DMI strength in our atomistic spin simulations. The exchange energy is identical for skyrmions and antiskyrmions with the same radial profiles (see methods). In contrast, DMI provides different contributions to the energy of skyrmions and antiskyrmions. For skyrmions, it leads to a preference of Néel-type skyrmions with a unique rotational sense due to the symmetry at the interface. For undistorted antiskyrmions the DMI energy vanishes. The origin of this difference is that each radial cut through the skyrmion (Fig. 1a) exhibits a clockwise spin rotation favored by the DMI. The antiskyrmion (Fig. 1b) exhibits opposite rotational senses along its two major axes which leads to zero total DMI energy for an undistorted antiskyrmion. In our simulations, however, we observe a small elongation of the antiskyrmion along the favored axis, leading to a lowering in energy. The maximal distortion can be seen in Fig. 1b, where the favored axis points in the y-direction. Similar but larger elongations of Ask have also been observed in refs. ^40,62.

In an inversion-symmetric material, the DMI cancels out leading to the degeneracy of skyrmions and antiskyrmions. Model systems of frustrated magnets in the absence of DMI have been studied in the past^26,42. In the considered Rh/Co bilayer on Ir(111), there is a significant DMI due to the heavy metal substrate which energetically prefers skyrmions over antiskyrmions while the exchange frustration stabilizes both. Because of this exchange frustration we found comparatively large lifetimes for both types of topological spin structures (cf. Fig. 2), but the skyrmion lifetimes exceed those of antiskyrmions due the DMI.

Since the DMI strength may vary depending on the considered interface, it is interesting to study the impact of DMI on the transition rates and lifetimes systematically. We parameterize the energy term due to DMI in our atomistic spin model, Eq. (1), by η ∈ [0, 1]:

$${E}_{{{{\rm{DMI}}}}}(\eta )=-\eta \mathop{\sum}\limits_{i,j}{{{{\bf{D}}}}}_{ij}\cdot ({{{{\bf{m}}}}}_{i}\times {{{{\bf{m}}}}}_{j})$$

(4)

such that η = 1 describes the full atomistic spin model for fcc-Rh/Co/Ir(111) as obtained from DFT and η = 0 an inversion-symmetric film with frustrated exchange but without DMI (for similar results obtained for hcp-Rh/Co/Ir(111) see Supplementary Note 3 and Supplementary Fig. 3).

So far, we have considered the stability of skyrmions and antiskyrmions versus collapse into the ferromagnetic ground state, i.e., their mean lifetime τ as given by the Arrhenius law Eq. (2). In order to discuss the annihilation and nucleation of skyrmions or antiskyrmions on equal footing, it is convenient to resort to the rate for a transition between state A and B defined by

$${\Gamma }^{{\text{A}\to\text{B}}}={\Gamma }_{0}^{{\text{A}\to\text{B}}}\exp (-\Delta {E}^{{\text{A}\to\text{B}}}/{k}_{{{{\rm{B}}}}}T),$$

(5)

where the activation energy ΔE^A→B = E^SP − E^A is obtained as the energy difference between the initial state A and the transition-state SP, and the pre-exponential factor ${\Gamma }_{0}^{\,{\text{A}\to\text{B}}\,}$ takes into account the entropic and dynamic contributions to the transition^52,61,63,64. Therefore, in our notation, Γ^X→FM denotes the annihilation rate and Γ^FM→X the nucleation rate for a state X ∈ {Sk, Ask}. Note that the superposition of all transition rates for collapses of a skyrmion or antiskyrmion is the inverse of the respective lifetime τ.

From our simulations, we find that skyrmions and antiskyrmions are metastable in the full range of η (Fig. 3). As expected, skyrmions and antiskyrmions exhibit the same radius (inset of Fig. 3b) for vanishing DMI (η = 0). With increasing DMI, the skyrmion diameter rises, while there is only a very small change for antiskyrmions. If we imagine skyrmions for simplicity as domain wall rings, they possess the correct geometry and rotational sense to be lowered in energy by the DMI. An extension of the domain wall length and consequently an increase of skyrmion size are therefore favored. For antiskyrmions, on the other hand, only segments of the ring are favored by the DMI, while other parts are disfavored due to exhibiting the opposite rotational sense, which leads to an elongation of antiskyrmions rather than an increase in size.

Fig. 3: Energy barriers and pre-exponential factors vs. DMI strength.

a Energy barriers for skyrmion (Sk, red and orange) and antiskyrmion (Ask, blue) annihilation (X → FM with X ∈ {Sk, Ask}) as a function of DMI strength η. b Pre-exponential factors for skyrmion and antiskyrmion annihilation at T = 1K. c Activation energies for skyrmion and antiskyrmion nucleation (FM → X with X ∈ {Sk, Ask}) and d corresponding pre-exponential factors at T = 1K. Here the helicity partition functions in the pre-exponential factor for skyrmion annihilation and nucleation shown in (b, d) have been calculated numerically using the temperature T = 1 K (see “Methods”). This is essential in order to assure continuity between harmonic and zero mode approximation at η → 0. In all panels, the red curves correspond to the radially symmetric skyrmion transition mechanism while the orange curves correspond to the Chimera mechanism. The inset in (b) shows the radius of skyrmions (red) and antiskyrmions (blue) as a function of η.

Similar to the radii, we observe a degeneracy of the annihilation barriers for skyrmions and antiskyrmions (Fig. 3a) at zero DMI (η = 0) and a similar collapse via the radial transition mechanism. For increased strength of the DMI (η > 0) the skyrmion collapse barrier rises. At η ≈ 0.78 the Chimera transition mechanism sets in and the slope of the rise changes. For η > 0.92 the Chimera state becomes metastable, which can be seen from the very shallow minimum in the MEP after the saddle point (Fig. 1g). However, it exhibits only a tiny energy barrier of ΔE < 8 meV against a collapse into the FM state barely visible on the energy scale of Fig. 1g. In relation to all other transitions this energy barrier is negligibly small and therefore it is ignored in the skyrmion annihilation/nucleation rate calculations.

There is only a minor influence of the DMI on the antiskyrmion annihilation energy barrier, as expected from the simplified micromagnetic model (see “Methods”). The small increase in the annihilation barrier with the DMI strength (Fig. 3a) is due to the antiskyrmions’ elliptical deformation that decreases their energy.

The DMI-dependent changes in activation energies for the nucleation of skyrmions (FM → Sk) and antiskyrmions (FM → ASk) (Fig. 3c), can be explained in a similar way. For skyrmions, we find a lowering of the nucleation barrier with rising DMI, while the barrier is nearly constant for antiskyrmions. This indicates that for low temperatures, the nucleation of skyrmions is favored over the nucleation of antiskymions in the presence of DMI.

The obtained pre-exponential factors ${\Gamma }_{0}^{\,{\text{A}\to\text{B}}\,}$, essential for the computation of transition rates (Eq. (5)), are shown in Fig. 3b, d for the annihilation and nucleation, respectively. For the annihilation of skyrmions, the pre-exponential factor rises monotonically with η for the radial collapse mechanism and changes by about two orders of magnitude. A rising pre-exponential factor increases the likelihood of the respective transition. Therefore, skyrmion annihilation occurs more often with rising η, effectively decreasing the state’s average lifetime. Thus DMI causes entropic destablization of skyrmions, opposite to the energetic stabilization. In the regime of the Chimera collapse, there is a slight drop with the DMI strength. For skyrmion nucleation (Fig. 3d), the pre-exponential factor displays a local minimum at η ≈ 0.5.

In contrast, for antiskyrmions, the pre-exponential factors change only slightly with varying DMI strength (Fig. 3b, d). To understand this behavior, we analyze the skyrmion and antiskyrmion eigenmodes.

Skyrmion and antiskyrmion modes

Within harmonic transition-state theory (HTST), the entropic contribution to the pre-exponential factor is given by the ratio of partition functions for the initial state and the transition state defined by the saddle point^60,64, where the partition functions are inversely proportional to the square root of the eigenvalues of the Hessian at the given state (see “Methods”). Therefore, we analyze the eigenvalue spectra of the skyrmion and its saddle point (Fig. 4a, b) as well as the antiskyrmion and its saddle point (Fig. 4c, d). We focus on eigenvalues, Ω_n, describing the curvature of the energy surface along localized modes below the magnon continuum⁶⁵, i.e., 0 meV ≤ Ω_n ≤ 2K ≈ 2.33 meV, where K is the MAE constant. We identify the associated modes by analyzing the eigenvectors of the Hessian matrix as in refs. ^52,64,66.

Fig. 4: Energy modes of skyrmions and antiskyrmions vs. DMI strength.

Eigenvalue spectra as functions of η: a for skyrmion initial state, b for skyrmion saddle point state, radial (red) and chimera (orange), c for antiskyrmion initial state, and d antiskyrmion radial saddle point state. Numbers in brackets give the amount of degenerate modes. The shaded blue area on top of each panel marks the magnon continuum with the lower boundary at Ω_mag = 2K ≈ 2.33 meV.

As expected, the energy spectra for skyrmions and antiskyrmions coincide for η = 0, i.e., in the absence of DMI. The same holds for the corresponding saddle point states. With increasing η, key differences arise:

Twofold mode

The twofold mode^52,63 describes a unidirectional distortion of the state. At the saddle point state, it can be seen as the movement of the Bloch point-like defect along a lattice direction. For the skyrmion, as a radial symmetric state, both distortions along different lattice directions increase the energy of the state by the same amount. Therefore, the twofold eigenmodes are degenerate (Fig. 4a, b). For the antiskyrmion, in contrast, the distortions along the directions, where spins rotate clockwise or counterclockwise, are only degenerate if the DMI is zero. With rising DMI, the distortion along one direction is favored over the other, which consequently leads to a splitting of the eigenvalues as observed in Fig. 4c, d.

Fourfold mode

This mode^52,63 quenches the state along a specific direction and elongates it perpendicular to this direction. As for the twofold mode we find two fourfold modes with interchanged directions for the quenching/elongation, which are degenerate for both skyrmions and antiskyrmions.

Helicity mode

The helicity mode^42,67 transforms between the Néel-type and the Bloch-type configuration of a skyrmion as illustrated in Fig. 5a–d. Due to DMI the helicity γ only affects the energy of skyrmions and not of undistorted antiskyrmions (see methods). In our film system, the DMI favors the Néel-type configuration with a right-rotating sense which corresponds to γ = π (Fig. 5d). For η = 0, i.e., in the absence of DMI, Néel- and Bloch-type skyrmions are degenerate. Therefore, the mode which transforms one into the other is a zero (Goldstone) mode, i.e., its eigenvalue is zero (see Fig. 4a, b). For finite values of DMI strength, the eigenvalue corresponding to the helicity mode becomes positive for skyrmions. This is contrary to antiskyrmions, where the helicity mode stays a zero mode for all values of η (Fig. 4c, d).

Fig. 5: Illustration of helicity modes for skyrmions and antiskyrmions.

a–d Transformation of a Néel-type skyrmion with a counterclockwise rotational sense (a, γ = 0) via mixed Bloch-Néel skyrmions (b, $\gamma =\frac{\pi }{3}$ and c, $\gamma =\frac{2\pi }{3}$) into a Néel-type skyrmion with the opposite rotational sense (d, γ = π). e–h The same sequence of helicities starting from a Néel-type antiskyrmion (e) leads to a 90^∘ rotation of the antiskyrmion in real space.

The different behavior of skyrmions and antiskyrmions can be understood by looking at their symmetry. Changing the skyrmions helicity results in a twist of the spin structure (Fig. 5a–d), which costs energy in the presence of DMI. For the antiskyrmion, on the other hand, a change in helicity by Δγ results in a real space rotation of the spin structure by Δγ/2 (Fig. 5e–h), leaving the energy of the state unchanged.

To our knowledge, the treatment of the helicity zero mode in HTST or Langer’s theory has not been reported in the literature so far. For antiskyrmions at arbitrary η and skyrmions at η = 0 we present an analytic expression for the zero mode partition function (see methods). For η > 0 we find temperature-dependent deviations in the zero mode- and harmonic approximation of the skyrmion helicity compared with numerically computed results (see Supplementary Note 4 and Supplementary Fig. 4). For a correct treatment of this mode we therefore use only numerically computed partition functions for the skyrmions and its radial SP’s helicity in the regime η > 0 (see “Methods”). Other zero modes are the translation of skyrmions and antiskyrmions treated as in ref. ⁶⁸ and the rotation of the chimera saddle point configuration which is determined by numerical integration along the mode, similar to skyrmions helicity for η > 0. With these methods we are able to compute the pre-exponetial factor of skyrmions and antiskyrmions in the whole range of η.

The eigenvalue spectra (Fig. 4) allow us to understand the dependence of the pre-exponential factors on the DMI strength η (Fig. 3b, d). The nearly constant values for pre-exponential factors of annihilation and nucleation of antiskyrmions can hardly be understood by looking at individual modes. The DMI is influencing modes all over the spectrum and consequently affecting the corresponding eigenvalues. Its negligible influence on the pre-exponential factor must therefore be seen as a global phenomenon.

For skyrmions, the diverging behavior of the annihilation pre-exponential factor with η → 1 can be explained by the behavior of the saddle point twofold modes. Considering the dependence ${\Gamma }_{0}^{\,{\text{Sk}\to\text{FM}}}\propto {\Omega }_{{\text{tf}}\,}^{-2}$ with Ω_tf being the eigenvalue of one of the two degenerated twofold modes, the pre-exponential factor diverges for Ω_tf → 0, which is the case for high η. The divergence of the skyrmions nucleation pre-exponential factor (Fig. 3d) can be explained similarly.

Transition rates

Based on the calculated activation energies and pre-exponential factors (Fig. 3), we can compute the annihilation and nucleation rates of skyrmions and antiskyrmions according to Eq. (5). The temperature dependence of the rates stems not only from the exponential decay over the energy barrier, but also from the pre-exponential factors, which inherit the temperature dependence from the single modes Boltzmann partition functions (see “Methods”).

The annihilation rates (Fig. 6a) range over many orders of magnitude depending on the temperature. We notice that the skyrmion annihilation rates decrease tremendously at low temperatures with rising DMI strength, i.e., with η. This effect is due to rising of the annihilation barrier with the DMI strength (Fig. 3a) and is reflected in Fig. 6a by the larger slopes of the annihilation rates for the larger η values. The effect of the DMI strength on the antiskyrmion annihilation (Fig. 6b) is quite small in agreement with its impact on the energy barrier (cf. Figs. 1h and 3a). The slope of the nucleation rate plots (Fig. 6c, d) displays a similar dependence on the DMI strength as the annihilation rates. Considering the logarithmic scaling of the ordinate (cf. Eq. (3)), the linear character of the curves also indicates that the activation energies ΔE are mostly responsible for the scaling of the rates over the shown range of temperatures. Nevertheless, the pre-exponential factor acts not only as an offset at T⁻¹ → 0 but also as small correction to the curves for increasing temperatures due to its intrinsic temperature dependence listed in Table 2.

Fig. 6: Annihilation and nucleation rates of skyrmions and antiskyrmions.

Transition rates of a skyrmion and b antiskyrmion annihilation as well as c skyrmion and d antiskyrmion nucleation. All rates are displayed vs. inverse temperature and for 11 different values of η evenly spaced in the interval [0, 1]. Radial and chimera transition as well as the superposition of both (skyrmion rates for η ∈ [0.78, 0.84]) are taken into account at the respective values of η.

Due to the increase of the skyrmions annihilation pre-exponential factor with η (cf. Fig. 3c), this leads to a crossing of annihilation rates for high and low η at high temperatures (cf. Fig. 6a). Considering that the antiskyrmion annihilation rates stay mostly unchanged with η, this is equal to a crossing of skyrmion and antiskyrmion annihilation rates for high η. Therefore, we can expect skyrmions to be destroyed more often than antiskyrmions at high temperatures in the realistic system with high DMI. In the next section, we point out how this entropic effect influences the equilibrium between skyrmions and antiskyrmions.

Skyrmion and antiskyrmion probabilities

We can use the transition rates from Fig. 6 to compute the relative probabilities of zero-field skyrmions and antiskyrmions at a given temperature and DMI strength. To do so, we model the problem as a Markov chain as illustrated in Fig. 7a and solve a three-state master equation for a system that can be either in the skyrmion, in the antiskyrmion, or in the FM state (see “Methods” for details). Note, that the transitions Sk → ASk and ASk → Sk have not been found in our GNEB calculations of minimum energy paths. These always end in a transition Sk → FM → ASk or vice versa which in turn is included in our description by the master equation.

Fig. 7: Ratio of probability distributions of skyrmions and antiskyrmions in thermal equilibrium.

a Schematic representation of Markov chain containing all transitions considered in the master equations. b The relative probabilities of skyrmions and antiskyrmions in thermal equilibrium are shown for different values of the DMI strength η. The insets show the time-dependence of the probabilities for the skyrmion, antiskyrmion, and FM state obtained by solving the master equation for η = 0.4 at c T = 500 K and d T = 200 K.

From the numerical solution of the master equation, we obtain the time-dependent probabilities p^Sk(t), p^ASk(t), and p^FM(t) to find the system in the Sk, Ask, or FM state, respectively. Two examples for the time-dependent probabilities are given for a DMI strength of η = 0.4, in Fig. 7c, d, one for T = 200 K, the other for T = 500 K. Note, that a rule of thumb in HTST, which states the theories validity when ΔE ≳ 5k_BT is fulfilled for all energy barriers, gives us the upper limit T ≲ 325 K considering the lowest energy barrier of 140 meV for the antiskyrmion annihilation at η = 0. For orientation, this limit is denoted as a tick on the upper x-axis in Fig. 7b. In both insets (Fig. 7c, d), the FM state is initialized at t = 0 in order to simulate the nucleation of skyrmions and antiskyrmions. In these examples, we find that antiskyrmions and skyrmions are nucleated on a similar time scale. However, for higher temperatures skyrmions decay faster due to the crossing of annihilation rates observed in Fig. 6a.

By taking the values in the limit t → ∞, we can obtain the equilibrium distribution of the states. These values are independent of the initial state chosen. In Fig. 7b, the ratio of the equilibrium probabilities for skyrmions and antiskyrmions is shown. In the absence of DMI (η = 0), skyrmions and antiskyrmions are degenerate, and thus the ratio of probabilities is unity. Increasing the DMI strength η stabilizes skyrmions versus antiskyrmions at low temperatures. However, for high temperatures the antiskyrmion becomes more likely than the skyrmion due to the effect discussed for Fig. 7c. This entropic effect becomes even more significant for high values of η but the crossover (${p}_{\infty }^{\,{\text{Ask}}}/{p}_{\infty }^{{\text{Sk}}\,}=1$) also shifts to higher temperatures. Nevertheless, it is remarkable that there exists a regime of temperatures below the rule-of-thumb-limit of T ≲ 325 K, where antiskyrmions become more likely to be found than skyrmions.

This general effect, which is based on the simultaneous stabilization of skyrmions and antiskyrmions due to exchange frustration and entropic destablization of skyrmions due to the interplay of DMI and temperature, opens the door for thermal manipulation of the relative probability of skyrmions and antiskyrmions. Although their lifetimes drop exponentially with temperatures (cf. Fig. 6), short and local heating pulses, the heat of which dissipates faster than the characteristic time for the equilibration between skyrmions and antiskyrmions (cf. insets of Fig. 7b), can nucleate antiskyrmions away from the equilibrium distribution due to the Kibble–Zurek mechanism^69,70, which has recently also been proposed based on Monte–Carlo simulations for the creation of skyrmions and antiskyrmions in inversion-symmetric materials⁴² and obeys a memory effect, which has already been observed in experiments⁷¹. According to our simulations, this effect permits for the thermal nucleation of antiskyrmions even for low global temperatures in frustrated magnets with broken inversion symmetry and significant DMI. This greatly extends the type of material systems in which coexisting metastable skyrmions and antiskyrmions can be created.

Discussion

Based on an atomistic spin model parametrized from DFT we have demonstrated that sub-10 nm skyrmions and antiskyrmions can coexist at zero magnetic field in the ferromagnetic ground state of the Rh/Co bilayer on the Ir(111) surface and that they exhibit large lifetimes of longer than an hour up to temperatures of 50 K. Since skyrmions have already been observed experimentally in Rh/Co/Ir(111)¹⁴, our prediction may enable the discovery of nanoscale antiskyrmions and their coexistence with skyrmions. It has been shown previously by simulations that skyrmions and antiskyrmions are distinguishable by means of spin-polarized STM^27,72. Our work further revealed that frustrated magnets even with significant DMI can be a promising class of materials to study the coexistence of spin structures with different topological charges.

While we have studied a specific experimentally feasible system, there is a variety of material systems in which the interplay of frustrated exchange and DMI can be achieved that may lead to coexisting Sk and Ask with similar properties as in our work. These include transition–metal films and multilayers^14,73,74 as well as 2D van der Waals magnets^47,48,75,76. In such systems, magnetic interactions can be engineered by varying the interface composition or structure that may further enhance the lifetimes of coexisting topological spin structures.

By systematically studying the influence of the DMI strength on the transition rates, we found that only skyrmions undergo entropic destablization with increasing DMI. Due to this effect, the relative probability of skyrmions and antiskyrmions can be tuned by thermal manipulation in frustrated magnets with broken inversion symmetry, e.g., by local heating using lasers or electrical currents. Together with the Kibble–Zurek mechanism, this can be exploited for the probabilistic nucleation of antiskyrmions away from equilibrium as needed for neuromorphic or Bayesian computing⁹.

Methods

Atomistic spin model of Rh/Co/Ir(111)

For all atomistic spin simulations on Rh/Co/Ir(111) presented in this work, we used the interaction constants given in Table 1 to parameterize the spin model, Eq. (1). All interaction constants used in our simulations have been obtained previously by means of DFT calculations¹⁴. We have performed simulations using the constants for fcc- and hcp-stacking of the Rh overlayer. Note, that the DMI strength has been reduced to 50% of the DFT value (η = 0.5) in case of hcp-Rh/Co/Ir(111) which leads to a very good agreement with experimental data on skyrmions in this system¹⁴. The values of the magnetocrystalline anisotropy energy are K = −1.166 meV and K = −1.635 meV for fcc-Rh/Co/Ir(111) and hcp-Rh/Co/Ir(111), respectively, where the negative sign indicates an easy out-of-plane magnetization axis and the magnetic moment is M = 2.5μ_B per atom. For all calculations, we used system sizes of either 50 × 50 or 70 × 70 lattice points with periodic boundary conditions representing the magnetic moments of the Co layer.

Table 1 DFT interaction constants of Rh/Co/Ir(111).

Note, that we have used a classical atomistic spin model in our work which is often an excellent approximation of quantum magnets in two and three dimensions⁷⁷. Nevertheless, quantum fluctuations may become relevant for atomic scale skyrmions⁷⁸. However, we are dealing with skyrmions and antiskyrmions of much larger size, on the order of a few hundred atoms. Previously, the same first-principles-based atomistic spin model as in our study has been successfully applied to Pd/Fe/Ir(111), a similar ultrathin film system, which has been extensively studied by spin-polarized STM experiments^4,12,51,58. Good quantitative agreement has been obtained with experimental work on skyrmion properties in Pd/Fe/Ir(111) such as diameters and lifetime^43,53,61 as well as for the temperature-dependent phase diagram⁵⁹.

Skyrmion and antiskyrmion profile

For initializing a skyrmion or antiskyrmion state in our simulation box, we use the following parametrization⁵

$${{{\bf{m}}}}(\phi ,\rho )=\left(\begin{array}{c}\cos [\Phi (\phi )]\sin [\Theta (\rho )]\\ \sin [\Phi (\phi )]\sin [\Theta (\rho )]\\ \cos [\Theta (\rho )]\end{array}\right)$$

(6)

where ϕ, ρ are the polar coordinates with respect to the skyrmion center. The angular functions are given by

$$\Phi (\phi )=k\phi +\gamma$$

(7)

$$\Theta (\rho )=\pi +\mathop{\sum}\limits_{\alpha =\pm 1}\arcsin \left[\tanh \left(\frac{2(\alpha c+\rho )}{w}\right)\right]$$

(8)

with vorticity $k\in {\mathbb{Z}}$, helicity γ ∈ [0, 2π] and profile parameters c, w > 0. Note that in our setups, the topological charge Q of states generated in this way is given by⁵

$$Q=\frac{1}{4\pi }\int\nolimits_{-\infty }^{\infty }{{{\bf{m}}}}\cdot ({\partial }_{x}{{{\bf{m}}}}\times {\partial }_{y}{{{\bf{m}}}})\,dxdy=-k\,,$$

(9)

which is therefore directly linked to the skyrmion (antiskyrmion) vorticity of k = 1 (k = − 1). After initialization of the skyrmions (antiskyrmions), the structure is optimized with respect to energy by the velocity projection optimization algorithm⁵⁷. From these configurations, the radius is estimated according to ref. ⁵⁶.

Micromagnetic energies

Inserting the parametrization given by Eqs. (6)–(8) into the micromagnetic energy functional⁷⁹ yields the following expression for the exchange energy

$$\begin{array}{l}{E}_{{\text{exc}}}\,=\,J{\displaystyle\int\nolimits_{-\infty }^{\infty }}\left[{\left({\nabla }_{x}{{{\bf{m}}}}\right)}^{2}+{\left({\nabla }_{y}{{{\bf{m}}}}\right)}^{2}\right]\,{d}^{2}{{{\bf{r}}}}\\ \qquad\quad=\,2\pi J\left[{\displaystyle\int\nolimits_{0}^{\infty }}\rho {\left(\frac{\partial \Theta }{\partial \rho }\right)}^{2}+\frac{{\sin }^{2}(\Theta )| Q{| }^{2}}{\rho }\,d\rho \right]\end{array}$$

(10)

and similarly for the DMI energy

$$\begin{array}{l}{E}_{{\text{DM}}}\,=\,D{\displaystyle\int\nolimits_{-\infty }^{\infty }}\left[{m}_{z}(\nabla \cdot {{{\bf{m}}}})-({{{\bf{m}}}}\cdot \nabla ){m}_{z}\right]\,{d}^{2}{{{\bf{r}}}}\\ \qquad\quad=\,2\pi D\left[{\displaystyle\int\nolimits_{0}^{\infty }}\rho \frac{\partial \Theta }{\partial \rho }+\frac{\sin (2\Theta )Q}{2}\,d\rho \right]{\delta }_{Q,-1}\cos \gamma \end{array}$$

(11)

where J and D are the micromagnetic exchange and DMI constant, respectively. From Eq. (10), it is evident that skyrmions with a topological charge of Q = −1 and antiskyrmions with Q = +1 are degenerate with respect to the exchange energy. In contrast, Eq. (11) demonstrates that only skyrmions are affected by the DMI and that E_DM vanishes for any other value of the topological charge Q due to the Kronecker delta term, δ_Q,−1. Note, that the DMI energy also depends on the helicity γ which prefers Néel-type skyrmions, i.e., γ = 0 (clockwise rotational sense) or γ = π (counterclockwise rotational sense) depending on the sign of D.

Computation of transition rates

The minimum energy path for a transition between two states A, B is computed using the geodesic nudged elastic band (GNEB) method⁵⁷. The highest energy point along the minimum energy path gives the SP, which is needed for the calculation of the transition rates Γ^A→B within the harmonic transition-state theory^61,80

$${\Gamma }^{{\text{A}\to\text{B}}}={\Gamma }_{0}^{{\text{A}\to\text{B}}}\exp (-\beta \Delta {E}^{{\text{A}\to\text{B}}}),$$

(12)

with $\beta ={({k}_{B}T)}^{-1}$. Here, the activation energy ΔE^A→B = E^SP − E^A is obtained as the energy difference between the initial state A and the transition-state SP, and the pre-exponential factor ${\Gamma }_{0}^{\,{\text{A}\to\text{B}}\,}$ takes into account the entropic contributions to the transition^52,63,64 and is given by^60,61

$${\Gamma }_{0}^{\,{\text{A}\to\text{B}}}=\frac{v}{\sqrt{2\pi \beta }}\frac{\mathop{\prod }\nolimits_{n=2}^{2N}{Z}_{n}^{{{\text{SP}}}}}{\mathop{\prod }\nolimits_{n=1}^{2N}{Z}_{n}^{{{\text{A}}}\,}}\,.$$

(13)

Here v denotes the dynamical factor and ${Z}_{n}^{X},X=\,{{\mbox{A, SP}}}\,$ are the partition functions of the initial and the transition states, which can either be treated in harmonic or zero mode approximation

$${Z}_{n}^{X}=\left\{\begin{array}{c}\sqrt{\frac{2\pi }{\beta {\Omega }_{n}^{X}}},\quad {\Omega }_{n}^{X}\, > \,0,\\ {L}_{n}^{X},\quad {\Omega }_{n}^{X}=0,\end{array}\right.$$

(14)

where ${\Omega }_{n}^{X}$ are Hessian eigenvalues and ${L}_{n}^{X}$ denotes the length of the zero mode in space of spin configurations. Here only the harmonic approximation depends on temperature. This leads to a temperature dependence of Γ₀ ∝ T^ΔG/2 for the pre-exponential factor where ΔG = G^X − G^SP is the difference in zero modes at the initial and saddle point state⁶⁰. However, the numerical treatment of the skyrmions helicity mode makes the temperature dependence even more complicated. The resulting temperature dependencies of the pre-exponential factors are summarized in Table 2. In this work, we deal with three different zero modes, which are addressed in detail in the following sections.

Table 2 Temperature dependence of pre-exponential factors.

Note, that we neglect quantum mechanical tunneling of the magnetization. It has been shown previously that this effect needs to be considered only at very low temperatures well below 10 K⁸¹, a regime in which the lifetime of skyrmions and antiskyrmions considered in our study is extremely long.

Partition function for zero modes

Zero modes^82,83 describe deformations of spin configurations that leave the energy of the state unchanged. Therefore, the curvature of the energy surface along this mode is zero, reflected by a vanishing Hessian eigenvalue of Ω ≈ 0 for the respective zero mode. A Boltzmann partition function Z^X for such a mode becomes the length of a curve ${{{\mathcal{C}}}}$,

$${Z}^{X}={\int}_{{{{\mathcal{C}}}}}{{{\mbox{e}}}}^{-\frac{\beta }{2}\Omega {q}^{2}}dq\approx {\int}_{{{{\mathcal{C}}}}}dq\,={L}^{X},$$

(15)

where ${{{\mathcal{C}}}}$ describes the modes trajectory in configuration space.

Translation zero mode

Translation zero modes describe in-plane displacements of a localized spin structure, such as the skyrmion. In our two-dimensional film we therefore have two translation modes. To calculate the partition function for one of these modes the spin structure M(r) = (m(r₁), . . . , m(r_N)) is translated by an in-plane vector a which gives M(r + a) The line integral from Eq. (15) can then be evaluated using the geodesic distance between starting and end point of the curve⁶⁴, reading

$${Z}_{{{\text{tr}}}}=\sqrt{\mathop{\sum }\limits_{n=1}^{N}| | {{{\bf{m}}}}({{{{\bf{r}}}}}_{n}+{{{\bf{a}}}})-{{{\bf{m}}}}({{{{\bf{r}}}}}_{n})| {| }^{2}}$$

(16)

where the geodesic norm ∣∣ ⋅ ∣∣ is the shortest distance between two points on the unit sphere.

Helicity zero mode

Let now ${{{\mathcal{C}}}}$ be a closed curve associated with the helicity change as illustrated in Fig. 5. In order to obtain an analytic form for the line integral from Eq. (15), we parameterize the curve using the helicity γ ∈ [0, 2π]

$${Z}_{{{\text{hl}}}}={\oint }_{{{{\mathcal{C}}}}}dq=\int\nolimits_{0}^{2\pi }| | \dot{{{{\mathcal{C}}}}}(\gamma )| | d\gamma$$

(17)

where the absolute of the gradient field is defined as $| | \dot{{{{\mathcal{C}}}}}(\gamma )| | =\sqrt{\mathop{\sum }\nolimits_{n = 1}^{N}{\left\vert \frac{d{{{\bf{m}}}}}{d\gamma }({{{{\bf{r}}}}}_{n},\gamma )\right\vert }^{2}}$. We now observe, that following the helicity mode, every spin m(r_n) performs a full rotation in the xy-plane. On the unit sphere it therefore travels on a ring with length ${L}_{n}=2\pi \sqrt{{\left({m}^{x}({{{{\bf{r}}}}}_{n})\right)}^{2}+{\left({m}^{y}({{{{\bf{r}}}}}_{n})\right)}^{2}}=2\pi \sqrt{1-{\left({m}^{z}({{{{\bf{r}}}}}_{n})\right)}^{2}}$. Since this rotation is homogeneous, the absolute value of the spins gradient along the curve becomes constant with $\left\vert \frac{d{{{\bf{m}}}}}{d\gamma }({{{{\bf{r}}}}}_{n},\gamma )\right\vert =\frac{{L}_{n}}{{\gamma }_{1}-{\gamma }_{0}}=\sqrt{1-{\left({m}^{z}({{{{\bf{r}}}}}_{n})\right)}^{2}}$ where γ₀ = 0 and γ₁ = 2π are the integration limits. Inserting this into Eq. (17), we obtain the analytic form for the helicity modes partition function in zero mode approximation

$${Z}_{\text{hl} }=2\pi \sqrt{\mathop{\sum }\limits_{n=1}^{N}1-{\left({m}^{z}({{{{\bf{r}}}}}_{n})\right)}^{2}}$$

(18)

In the case of skyrmions, the helicity mode cannot be treated as a zero mode anymore in the presence of DMI. We therefore model the homogeneous rotation of spins along the mode by a series of states M_k, where spins are iteratively rotated around the z axis by the change in helicity Δγ

$${{{{\bf{m}}}}}_{k}({{{\bf{r}}}})={{{{\bf{R}}}}}_{\Delta \gamma }^{z}{{{{\bf{m}}}}}_{k-1}\left({{{\bf{r}}}}\right)\,,$$

(19)

From this series of images, the Boltzmann partition function can be numerically computed

$$Z=\mathop{\sum }\limits_{k=1}^{K}{e}^{-\beta ({E}_{k}-{E}_{0})}\sqrt{\mathop{\sum }\limits_{n=1}^{N}| | {{{{\bf{m}}}}}_{k}({{{{\bf{r}}}}}_{n})-{{{{\bf{m}}}}}_{k-1}({{{{\bf{r}}}}}_{n})| {| }^{2}}$$

(20)

Chimera rotation zero mode

Another zero mode that occurs in our simulations is the rotation mode of the chimera SP, describing a full rotation of the structure around the Bloch point (cf. ref. ⁵³). To our knowledge, no analytic expression for the partition function of this particular mode has been presented yet. We compute the corresponding Boltzmann partition function numerically as in Eq. (20). For that, we parameterize the images along the mode by the angle φ so that spins at the same position r_n but in consecutive images k − 1 and k read⁸⁴

$${{{{\bf{m}}}}}_{k}({{{{\bf{r}}}}}_{n}-{{{\bf{c}}}})={{{{\bf{R}}}}}_{\Delta \varphi }^{z}{{{{\bf{m}}}}}_{k-1}\left({{{{\bf{R}}}}}_{-\Delta \varphi }^{z}[{{{{\bf{r}}}}}_{n}-{{{\bf{c}}}}]\right)\,.$$

(21)

Here c is the center of rotation on the lattice, which has been chosen to be the location of the saddle points Bloch point-like configuration. This position is computed as the point where the absolute value of the states interpolated magnetization vanishes, therefore ${{{\bf{c}}}}=\arg {\min }_{{{{\bf{r}}}}}| {{{\bf{m}}}}({{{\bf{r}}}})|$.

Master equations

Consider a system, that can be in three different states: the skyrmion state, the antiskyrmion state, and the FM state. Based on the calculated transition rates Γ^m→n between those states, we can simulate the time evolution of the probabilities of the states, as described in the following. Let p^A be the time-dependent probability for the system to be in a specific state A = FM, Sk, Ask. We construct the vector

$${{{\bf{p}}}}(t)={\left({p}^{{{\text{FM}}}}(t),{p}^{{{\text{Sk}}}}(t),{p}^{{{\text{Ask}}}}(t)\right)}^{T},\quad \mathop{\sum}\limits_{\,{{\text{A}}}}{p}^{{{\text{A}}}}(t)=1$$

(22)

The equation describing the evolution of p(t) is the master equation

$$\frac{d}{dt}{{{\bf{p}}}}(t)={{{\bf{U}}}}{{{\bf{p}}}}(t)$$

(23)

with transition matrix U

$${\left({{{\bf{U}}}}\right)}_{AB}=\left\{\begin{array}{c}{\Gamma }^{B\to A}\quad \,{{\mbox{if}}}\,\quad \,{{\mbox{A}}}\,\ne \,{{\mbox{B}}}\,\\ -\mathop{\sum}\limits_{A}{\Gamma }^{B\to A}\quad \,{{\mbox{else}}}\,\end{array}\right.$$

(24)

The general solution to the master equation can then be computed from the three eigenpairs (λ_n, v_n) of U

$${{{\bf{p}}}}(t)=\mathop{\sum }\limits_{n=1}^{3}{c}_{n}{{{\mbox{e}}}}^{{\lambda }_{n}t}{{{{\bf{v}}}}}_{n},\quad {{{\bf{U}}}}{{{{\bf{v}}}}}_{n}={\lambda }_{n}{{{{\bf{v}}}}}_{n}$$

(25)

In the last step, the constants c_n have to be calculated from the initial conditions which, in our case, correspond to the initialization of the system in the FM state. The equilibrium distribution p_∞ for large times t → ∞, which fulfills

$$0={{{\bf{U}}}}{{{{\bf{p}}}}}_{\infty }$$

(26)

is directly computed from the nullspace of U.

The temperature dependence of the solutions stems from the temperature dependence of the rates themselves (cf. Table 2). In the setup of the transition matrix, we also take into account that the probability for finding a skyrmion or antiskyrmion increases with the search area. This is done by scaling only the nucleation rates by the number of possible nucleation saddle points in a given area F, reading on the hexagonal lattice $2F={a}^{2}\sin (\pi /3){N}_{sp}$ with the Iridium lattice constant of a = 0.2715 nm. Physically, this reflects the fact that a nucleation can take place at any interstitial site of the lattice^53,61. In this regard, the probabilities in the inset of Fig. 7b show the probability for finding at least one skyrmion, antiskyrmion or none of both (FM) in an area F = 1 μm². Note that the ratio of skyrmion and antiskyrmion probability has been found to be independent of F.