Towards room temperature nanoscale skyrmions in ultrathin films

Breaking the dilemma between small size and room temperature stability is a necessary prerequisite for skyrmion-based information technology. Here we demonstrate by means of rate theory and an atomistic spin Hamiltonian that the stability of isolated skyrmions in ultrathin ferromagnetic films can be enhanced by the concerted variation of magnetic interactions while keeping the skyrmion size unchanged. We predict film systems where the lifetime of sub-10 nm skyrmions can reach years at ambient conditions. The long lifetime of such small skyrmions is due to exceptionally large Arrhenius pre-exponential and the stabilizing effect of the energy barrier is insignificant at room temperature. A dramatic increase in the pre-exponential is achieved thanks to softening of magnon modes of the skyrmion, thereby lowering the entropy of the skyrmion with respect to the transition state for collapse. Increasing the number of skyrmion deformation modes should be a guiding principle for the realization of nanoscale, room-temperature stable skyrmions.

Although potential feasibility of this concept has been demonstrated by experimental detection of skyrmions in various magnetic systems [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23], realization of viable skyrmion-based digital devices is an ongoing, challenging problem. Ultimately, the size of the skyrmionic bits should not exceed 10 nm so as to improve on the data density level already achieved in conventional technology. For such small skyrmions, thermal stability becomes an issue as thermal fluctuations can induce spontaneous collapse of the skyrmion state and, therefore, corrupt the stored data. Indeed, sub-10 nm skyrmions have so far been detected only at very low temperatures [7-10, 17, 23]. On the other hand, room temperature stability has been reported for larger skyrmions [11][12][13][14][15][16][18][19][20][21][22]. The challenge is to design materials where skyrmions are both small and long-lived at ambient conditions. Theoretical calculations can help engineer skyrmions towards applications by evaluating skyrmion lifetime as a function of material parameters, applied stimuli, and temperature. The lifetime τ , i.e. the mean time it takes the magnetic system coupled to the heat reservoir of temperature T to escape from the energy well corresponding to the metastable skyrmion state, can be described by the Arrhenius expression: (1) * bessarab@hi.is and characterized by the annihilation energy barrier ∆E and the pre-exponential factor τ 0 . Previous theoretical studies have focused on calculations of the energy barrier with respect to two mechanisms of skyrmion decay into a polarized, ferromagnetic configuration: Radial collapse of a skyrmion [24][25][26] and escape of the skyrmion through the system's boundary [27][28][29]. Asymmetric collapse of a skyrmion in systems with frustrated magnetic exchange interaction has also been reported in recent studies [23,30,31]. Minimum energy path calculations applied to an atomistic spin Hamiltonian have revealed dependencies of the skyrmion annihilation barrier on relevant material parameters [27,32], applied magnetic field [33,34], presence of impurities [28,35], and frustration in the magnetic exchange interaction [26].
Recently, Büttner et al. [36] have undertaken an extensive analysis of the phase diagram for isolated skyrmions within the continuous magnetization framework, where the annihilation energy barrier was approximated by the universal energy of the zero-diameter skyrmion relative to the skyrmion energy minimum. The study concluded that ultrasmall skyrmions in materials exhibiting local ferromagnetic order of spins could not be stable at ambient conditions [22].
However, knowledge about the energy barrier alone does not provide a complete picture of the thermal stability. Entropic and dynamical effects must also be evaluated for the prediction of the lifetime. These contributions are incorporated in the pre-exponential factor τ 0 , inverse of which is often referred to as the attempt frequency. Definite evaluation of the prefactor appears to be particularly important for skyrmionic systems, where τ 0 assumes unusual values [37,38], depends strongly on the mechanism of skyrmion collapse [33,39] and demonstrates extreme sensitivity to applied mag-netic field [33,40,41]. Clearly, assuming τ 0 to have some fixed value, e.g. associated with the Larmor precession, can lead to incorrect conclusions about thermal stability of magnetic skyrmions. Recent studies have revealed the key role of the skyrmion's internal modes in the unusual behaviour of the prefactor [39,40]. In particular, the modes can cause large entropy barriers, lowered attempt frequencies and enhanced thermal stability of skyrmions [38][39][40][41].
In this article, it is explored to what extent the lifetime of nanoscale skyrmions in ultrathin ferromagnetic films can be enhanced under ambient conditions. The focus on the ultrathin systems is motivated not only by fundamental importance, but also by technological relevance due to improved control over the spin structure and reduced power dissipation associated with currentdriven skyrmion motion (see e.g. Ref. [42] and references therein). The thermal stability diagram of an isolated skyrmion is systematically analyzed using the harmonic transition state theory for magnetic systems [43,44] and atomistic spin Hamiltonian. Both the collapse energy barrier ∆E and the pre-exponential factor τ 0 are calculated and the skyrmion lifetime is evaluated as a function of magnetic interaction parameters and temperature.
The results provide information about systems hosting long-lived magnetic skyrmions at room temperature and zero applied magnetic field. It is problematic to reach energy barriers exceeding thermal energy by a factor of 40-50 at room temperature -a commonly used criterion for reliable information storage -with realistic material parameter values while keeping the skyrmion size at nanoscale. This conclusion is in agreement with previous studies [36]. However, sub-10 nm skyrmions in ultrathin ferromagnetic films can indeed be stable on a timescale of years at ambient conditions thanks to the remarkably large value of the Arrhenius prefactor τ 0 .

RESULTS
Phase diagram. An extended ultrathin skyrmionic system is described here by a classical atomistic spin Hamiltonian on a monolayer hexagonal lattice. The total energy E of the system includes three contributions due to the Heisenberg exchange, Dzyaloshinskii-Moriya (DM) interaction and magnetic anisotropy, respectively, where each contribution is characterized by an effective interaction parameter: J, D, and K (see Methods section for the detailed description of the simulated system).
The technologically-relevant case of zero applied magnetic field is in focus of the present study, therefore the Zeeman term is not included in Eq.
(2). A hexagonal lattice is used here so as to directly mimic ultrathin systems where the presence of skyrmions has been confirmed experimentally [7][8][9][10]23], but the results of this work Figure 1. Magnetic phase diagram for a monolayer of magnetic moments on a hexagonal lattice. A spin spiral is the ground state of the system in the SS region (blue), while the ferromagnetic state is the only stable configuration in the FM domain (green). Metastable, isolated skyrmions exist in the Sk sector. The color scheme (see top bar) shows the calculated HTST lifetime of isolated skyrmions in terms of the intrinsic precession time τint = µ/Jγ, which is on the order of femtoseconds for a typical magnetic system. A thermal energy of 2.6J was assumed in the lifetime calculations. Some of the contours of equal skyrmion radius R and lifetime τ are shown with blue and black solid lines, respectively. Material parameter domain corresponding to skyrmions with R < 8a and τ > 10 6 τint is marked with hatching. Insets show skyrmion profiles and definition of the skyrmion radius for the two points, marked I and II, lying on the same R-isoline. Black and white diamonds indicate effective material parameters for Fe-and Co-based ultrathin skyrmionic systems, respectively, which are listed in the legend. For Ref. [23], the parameters were extracted using the method described in Supplementary Note 2.
can be interpreted in terms of parameters of the squarelattice model or continuous magnetization Hamiltonian using the material parameter transformations provided in Supplementary Note 1. Figure 1 shows the phase diagram of the system in zero applied magnetic field. The diagram was obtained by relaxing a trial single-skyrmion profile to a local energy minimum and examining how the resulting magnetic configuration depends on the reduced parameters of magnetic anisotropy, K/J, and DM interaction, D/J. Within a certain domain in the material parameter space, isolated Néel-type skyrmions emerge as metastable states in the ferromagnetic (FM) background. The skyrmion domain occupies a rather narrow region of the phase diagram [1, 4,32]. In fact, most of skyrmionic materials known so far exhibit spin-spiral states at zero external magnetic field and a magnetic field needs to be applied for skyrmions to start forming. This scenario is realized, for example, in atomic Pd/Fe bilayers on Ir(111) [9,26,38] and on Rh(111) [45]. Based on the values of the effective parameters of magnetic interactions one can place these systems on the phase diagram (see black diamonds labeled 1 − 6 in Fig. 1). Nevertheless, zero-field skyrmions have been reported in a few systems, including atomic Rh/Co bilayers on Ir(111) [23], for which the stacking of the Rh layer affects the material parameters (see white diamonds labeled 7 and 8 in Fig. 1, which correspond to hcp and fcc stacking, respectively).
Minimum energy path (MEP) calculations (see Methods section) revealed the mechanism of skyrmion annihilation, which turned out to be the same for all material parameter values within the chosen domain. The mechanism corresponds to usual radial collapse, where the skyrmion symmetrically shrinks and eventually disappears [24][25][26]32].
Based on the calculated MEPs for skyrmion annihilation, the lifetime τ characterizing the stability of the skyrmions against thermally activated decay into the FM state was evaluated using Eq. (1) where both the energy barrier and the pre-exponential factor were calculated according to the harmonic transition state theory (HTST) (see Methods section). The thermal energy of 2.6J was assumed in the lifetime calculations, which roughly corresponds to room temperature for J = 10 meV, a typical value of the Heisenberg exchange parameter in ultrathin magnetic films [9]. The calculated lifetime as a function of material parameters is superimposed on the phase diagram in the form of a contour plot. For the sake of generality, the lifetime is presented in Fig. 1 in units of the intrinsic precession time τ int = µ/Jγ, with µ, γ being the on-site magnetic moment and gyromagnetic ratio, respectively. τ int lies in the femtosecond range for usual magnets.
The radius R of relaxed, energy-minimum skyrmion states was also quantified in terms of the nearest-neighbor distance a using the definition of Bogdanov and Hubert [3]: where Θ(r) is the polar angle of the spin located at the distance r from the skyrmion center and the subscript 0 denotes the point of the steepest slope of Θ(r). Eight contours of constant R were carefully traced, with R ranging from 5a to 12a (see Methods Section). Two of the contours are shown in Fig. 1 as examples.
Distributions of the lifetime and radius can now be compared straightforwardly, giving a clear picture of the relationship between skyrmion size and skyrmion stability. Both R and τ increase with the DM interaction strength but decrease as the anisotropy parameter gets larger (see Fig. 1), i.e. both quantities behave alike as functions of relevant material parameters. This is consistent with the experimentally observed trend that large skyrmions tend to be more stable than small ones. A more detailed analysis shows, however, that R-and τisolines do intersect indicating that skyrmion lifetime at a given temperature is not uniquely defined by the skyrmion size. In other words, the stability of skyrmions can be tuned by a concerted material parameter transformation that preserves the skyrmion radius. This conclusion is further supported by a direct comparison of lifetimes of the skyrmions that belong to the same Rcontour. For example, parameter sets K I = 0.09J, D I = 0.16J and K II = 0.54J, D II = 0.52J (the sets are marked with crosses on the phase diagram, see Fig. 1) result in the same equilibrium skyrmion radius, R = 6a. However, the corresponding HTST-estimates of the lifetime differ by five orders of magnitude: τ I = 4.45 × 10 2 τ int , while τ II = 2.43 × 10 7 τ int . Due to crossing of R-and τ -isolines, it is in fact possible to indicate a material parameter domain which corresponds to skyrmions with R < R α and τ > τ α for a given time span τ α and radius R α . Such a domain is highlighted in Fig. 1 for τ α = 10 6 τ int and R α = 8a as an example.
Interestingly, the shape of the skyrmion profile is not conserved along the R-isolines either. Specifically, the profile has an arrow-like shape for parameter set I; For set II, the circular domain wall encompassing the skyrmion core becomes thinner and the profile changes to a form resembling a magnetic bubble where more spins at the skyrmion core point almost antiparallel to the magnetization of the FM background (see the insets in Fig. 1). In general, bubble-like skyrmions turn out to be more stable than arrow-like skyrmions of the same size. This observation could facilitate practical realization of long-lived, nanoscale skyrmions at room temperature.
It seems clear from Fig. 1 that the skyrmion lifetime can indeed be enhanced via the concerted increase in both K/J and D/J while keeping the skyrmion size unchanged. To elucidate the origin of the lifetime enhancement, it is informative to examine the collapse energy barrier ∆E and the Arrhenius pre-exponential factor τ 0 separately as functions of displacement along the contours of equal skyrmion radius R. This analysis is presented in the following.
Collapse energy barriers. The variation of ∆E along the contours of equal skyrmion radius is shown in Fig. 2. For all considered radii, the barrier grows monotonically as K/J and D/J increase in a concerted way, i.e. the shape of the skyrmion gradually changes from arrow to bubble. Decomposition of the energy barrier into individual interaction-resolved components (see Fig. 2b) demonstrates positive contribution from the DM interaction and negative contribution from the Heisenberg exchange and magnetic anisotropy. This result is in agreement with earlier studies of the skyrmion collapse [25,27,32]. Al- though all of the contributions increase steadily in absolute value along the R-isolines and largely compensate each other, the balance between them changes, which results in the barrier enhancement.
Enhancement of the energy barrier along the R-isolines is the consequence of the mechanism of skyrmion annihilation. Indeed, progressive shrinking of the skyrmion is produced by a symmetrical rotation of the spins towards the FM state [24,25,33]. Increase in the amplitude of this rotation as well as in the number of spins involved in the process favor enhancement of the corresponding energy barrier. This is consistent with the tendency of larger skyrmions to have a higher energy barrier for collapse. On the other hand, the overall rotation of spins involved in the radial collapse can be enlarged without changing the skyrmion size simply by making the skyrmion core thicker. This change in the skyrmion shape, leading to the energy barrier enhancement, is indeed realized along the skyrmion radius isolines studied here. As the size of the skyrmion core can not exceed the skyrmion diameter, the growth of the barrier is expected to stop at a certain level for a given skyrmion radius. Figure 2 shows that the barrier enhancement in-deed becomes weaker for large values of K/J and D/J corresponding to bubble-like skyrmions for all considered skyrmion radii.
Within the HTST, the bottleneck for the skyrmion collapse is represented by the relevant first order saddle point (SP) on the energy surface of the system and the energy barrier ∆E is defined as the energy difference between the SP and the local minimum corresponding to the skyrmion state: MEP calculations have shown that SPs for the radial collapse correspond to a very small, Bloch point-like defect in the FM background [24,25]. With this knowledge, it seems reasonable to approximate the bottleneck state by a zero-size skyrmion, whose energy E 0 is given by the universal energy of a Belavin-Polyakov soliton with topological charge |Q| = 1 [46,47], and define the barrier as E 0 − E min , as it was suggested in previous studies [36]. This method is, however, not justified by the rate theories, and care must be taken when applying it to the skyrmion stability problem. Indeed, E 0 , which amounts to 4 √ 3πJ for the atomistic spin model on a hexagonal lattice (see Supplementary Note 3), is significantly larger than the SP energy obtained from the MEP calculations for all skyrmions considered in this study (see the inset in Fig. 2). As a result, the approximate method, although correctly capturing the trend for the barrier, gives a large overestimate as compared to the more accurate HTST predictions (see Fig. 2). It is concluded here that the barrier estimates based on the energy of the zero-diameter skyrmion can be insufficient for the quantitative determination of the skyrmion stability and lifetime. This conclusion is in agreement with the recent study by Heil et al. [31] demonstrating that the Belavin-Polyakov energy needs to be properly corrected in order to obtain an accurate result for the SP energy. The increase in the energy barrier along the R-isolines partially explains why bubble-like skyrmions are more stable than arrow-like ones [32]. However, it is evident from Fig. 2 that the barrier never exceeds 10J even for largest skyrmions considered here (R = 12a). For realistic J on the order of 10 meV, this translates into roughly 4k B T at ambient conditions, which seems to be too little to ensure room-temperature stability on macroscopic time scale [36]. It will be shown below that the stabilizing effect of the barrier is indeed insufficient for the small skyrmions considered here to be stable at room temperature. However, a long lifetime of nanoscale skyrmions at ambient temperature can still be achieved thanks to extremely large prefactor τ 0 , as demonstrated in the following Sections.
Pre-exponential factor. Calculated results for the prefactor τ 0 (see Methods section for the details of HTST calculations) are shown in Fig. 3. Overall, τ 0 depends strongly on the material parameters and demonstrates a steady growth along the isolines of skyrmion radius in the direction of increasing K/J and D/J. This behavior of the pre-exponential factor in fact enhances the stabilizing effect of the energy barrier for bubble-like skyrmions, i.e. contributions of τ 0 and ∆E to the skyrmion lifetime are both in the same direction. Strong variation of τ 0 with the material parameter values confirms that the constant prefactor approximation can indeed be inapplicable for the skyrmion lifetime calculations, which is in agreement with the conclusions of previous studies [33,[38][39][40]. Remarkably, the calculated values of τ 0 can span twenty orders of magnitude within the chosen material parameter range. These vast changes in the prefactor are actually comparable to what has recently been observed experimentally for skyrmions in the Fe 1−x Co x Si system [41].
Inspection of the results obtained for various R-isolines also shows sensitivity of τ 0 to the skyrmion radius. For fixed D/J, the prefactor grows with R, contributing to the lifetime enhancement for large skyrmions (see the inset in Fig. 3). The prefactor depends on R stronger for larger values of the reduced DM interaction. Strong Rdependence of the prefactor was also reported in Ref. [40], but there the change in the skyrmion size was induced by an external magnetic field rather than by the variation int int Figure 3. Results of pre-exponential factor calculations within HTST. Variation of the Arrhenius pre-exponential factor τ0 along several contours of equal skyrmion radius R. Color codes distribution of the reduced anisotropy parameter. The inset shows R dependence of τ0 for the two values of the reduced DM interaction parameter. Filled squares indicate the calculated data points.
of the material parameters.
The dramatic change in the prefactor stems from the entropic effects [38][39][40][41], which are represented within HTST by the Hessian of the energy of the system at the SP and skyrmion state minimum (see Methods section): where the determinants can be computed as a product of the eigenvalues. The numerical results for the prefactor can be understood by analysing the eigenvalues of the Hessian representing the energy spectra of magnetic excitations at the skyrmion state as well as at the transition state. Figure 4 shows how the excitation spectra evolve along one of the R-isolines. At the transition state, the spectrum only slightly differs from that of the unperturbed FM system due to a rather small non-collinear region at the SP (see insets in Fig. 4a). In particular, only a few modes localized on the defect are present within the anisotropy-induced gap, ε gap = 2K, and the number of such modes does not change with the material parameter values. In contrast, the presence of a skyrmion can produce significant changes to the spectrum of the system [48][49][50][51]. In addition to the modes corresponding to the in-plane translation and uniform breathing of the skyrmion as well as rotation of the skyrmion core, progressively more localized modes describing various skyrmion deformations emerge in the gap as K/J and D/J increase in a concerted way while the equilibrium skyrmion radius remains unchanged (see Fig. 4b). The deformation modes can be described as periodic modulations of local skyrmion radius along the azimuthal direction. For fixed material parameters, the deformation with more modulation periods along the perimeter of the skyrmion corresponds to a larger eigenvalue of the Hessian. On the other hand, the energy of each deformation mode decreases monotonically along the R-isoline in the direction of increasing K/J and D/J. As a result of the difference between the spectra at the transition state and the skyrmion state, the ratio of the products of eigenvalues -i.e. the ratio of determinants -becomes large [see Eq. (10)]. The increase in the number of localized modes in the magnon gap and their softening, which are characteristic features of the skyrmion state, do not occur at the transition state. This increases the entropy difference between the skyrmion state and transition state which in turn leads to larger values of the pre-exponential factor τ 0 for the skyrmions with the bubble-like shape. Skyrmion lifetime. Figure 5 summarizes the calculated results for the lifetime τ of nanoscale isolated skyrmions in an ultrathin FM film at ambient conditions. Absolute values for τ were obtained by setting temperature T to 300 K and assuming typical values for the exchange parameter J and the magnitude of the magnetic moment µ: J = 10 meV and µ = 3µ B , with µ B being Bohr magneton [9]. Strong variation of the prefactor along the R-isolines (see Fig. 3) has a clear impact on the skyrmion lifetime, which can exceed the level of ten years for relatively large reduced DM interaction and anisotropy. Note that the diameter of the skyrmions considered here does not exceed 24a, which translates into 6.5 nm for the samples grown on the Ir(111) surface [7]. Therefore, it is still possible to obtain long skyrmion lifetime for more moderate values of the material parameters by increasing the skyrmion size while keeping it within the nanoscale range.
To emphasize the decisive role of the Arrhenius pre-exponential factor in the stabilization of nanoscale skyrmions at room temperature, the lifetime was also calculated using the commonly applied constant prefactor approximation. For each skyrmion radius, the value of τ 0 calculated at D/J = 0.12 was used to evaluate the lifetime along the entire R-isoline. For D/J = 0.12, the Variation of the HTST-predicted skyrmion lifetime τ (black solid lines) along several contours of equal skyrmion radius R. The following values of temperature T , Heisenberg exchange parameter J and magnitude of the magnatic moments µ were assumed in the calculations: T = 300 K, J = 10 meV, µ = 3µB . Color codes distribution of the reduced anisotropy parameter. Filled squares indicate the calculated data points. The inset shows the Arrhenius plot for the two values of the reduced DM interaction parameter on R = 12a isoline. The magnitude of the lifetime calculated in the constant Arrhenius prefactor approximation is within the pink area for all skyrmion radii from 5a to 12a. Blue area indicates lifetimes greater than ten years. magnitude of τ 0 changes from 30 ps for R = 5a to 1 ns for R = 12a, which is actually in the range of typical prefactor values assumed in studies of thermally activated magnetic transitions [27,[53][54][55][56]. The values of the skyrmion lifetime calculated with these values of τ 0 are enclosed in the pink area in Fig. 5. Since the prefactor is assumed to remain constant along each R-isoline in this case, the slight enhancement of the lifetime is due to the increase in the energy barrier. Clearly, this enhancement is insignificant compared to what HTST predicts. Relatively small increase in the energy barrier but a large change in the pre-exponential factor are further illustrated by the Arrhenius plots for the two points along R = 12a isoline (see the inset in Fig. 5): The plots have similar slopes but intercept the vertical axis at very different levels.

DISCUSSION
The qualitative difference in the behavior of the excitation spectra at the transition state and at the skyrmion state is the origin of large changes in the prefactor. This phenomenon is observed for all skyrmion sizes considered in the present study (see Supplementary Fig. 1). However, larger skyrmions can accommodate more modulation periods along their circumference. This leads to stronger variation of τ 0 along R-isolines corresponding to larger skyrmion radii (see Fig. 3). The possibility to introduce more localized modes in the magnon gap by increasing the skyrmion size also explains the sensitivity of the prefactor to the skyrmion radius.
The emergence of the modes localized on the skyrmion and their uncompensated softening are the primary reasons for the large variation of the prefactor along the Risolines. To gain a qualitative insight into these effects, it is convenient to treat an isolated skyrmion as a composition of a core, an outer FM domain, and a wall separating the core and the outer domain. Under this representation, transformation of the skyrmion shape along the Risoline can be interpreted in terms of variation of the domain wall (DW) width, with wider DWs corresponding to arrow-like skyrmions and thinner DWs corresponding to bubble-like skyrmions (see the insets in Fig. 1). On the other hand, the skyrmion deformations can be viewed as transverse fluctuations of the DW and parametrized by Fourier harmonics. Neglecting the curvature of the DW and disregarding the DM interaction as well as long-range dipolar interaction, the deformation energy ε n associated with the nth harmonic is given by [52]: where L is the DW length (see Supplementary Note 4). The number of harmonics N in the magnon gap is then defined by the largest n for which ε n < ε gap = 2K: where ∆ ∼ a J/K is the DW width. N gives a qualitative estimate of the number of the modes localized on the skyrmion with perimeter 2πR = L. This increases monotonically as the DW becomes thinner, which explains why progressively more modes localize on the skyrmion as it gradually changes its shape along the R-isoline to form a bubble-like structure. Equation (5) does not explain well all features of the excitation spectrum of the skyrmion because it is based on the asumption of straight, achiral DW. In particular, the deformation mode softening for the bubble-like skyrmions is not reproduced by Eq. (5). However, the mode softening can still be obtained using the straight DW model if one includes the DM interaction in the Hamiltonian and considers a concerted increase in both K and D. In particular, energy of transverse deformation modes of the DW does show a decrease when the variation of K and D is identical to that along the contours of equal skyrmion radius (see Supplementary Fig. 2).
It seems clear from the above results that room temperature stability of nanoscale skyrmions can be achieved almost entirely due to exceptionally large pre-exponential factor τ 0 rather than high energy barrier ∆E. Large values of the prefactor are an indication of large entropy barriers [38][39][40][41]. Therefore, the concept of obtaining long skyrmion lifetime by means of the enhanced entropy barriers, which was pointed out independently in Refs. [39,40], becomes critical for sub-10 nm skyrmions at room temperature. A key strategy for the realization of large entropy barriers stabilizing the skyrmion state is to establish as many localized modes corresponding to skyrmion deformations as possible, because these modes introduce uncompensated increase in the entropy of the skyrmion state. For nanoscale skyrmions at zero applied magnetic field, the sufficient number of localized modes is realized when the skyrmions resemble magnetic bubbles, i.e. when reduced parameters of DM interaction and anisotropy acquire relatively large values compared to the Heisenberg exchange parameter. This regime seems difficult to achieve in typical transition-metal systems where J assumes usual values on the order of 10 meV. However, desired values of D/J and K/J can actually be obtained by reducing the value of J (calculated skyrmion lifetime for several values of J and T is shown in Supplementary  Fig. 3). Note that the parameter J should be understood as an effective exchange constant for systems with frustrated exchange. There, the effective exchange parameter characterizes well low-energy excitations of the system, but may fail to describe the states far from the energy minimum [26], and, therefore, small value of J does not necessarily mean low ordering temperature [57], especially for the anisotropic systems [58].
Interestingly, the material parameters approach the desired values in ultrathin Fe and Co films (see Fig. 1) which have already been studied [9,23,26,38,45], making these systems promising for ultrasmall, roomtemperature skyrmions stabilized by entropy barriers. Systems combining layers of 3d, 4d and 5d transitionmetal elements appear to be particularly interesting since they make it possible to tune the strengths of the Heisenberg exchange, DM interaction and anisotropy independently and over a wide range via layer composition, alloying, and intermixing at the 4d/3d and 5d/3d interfaces. The feasibility of this approach has been proposed based on density functional theory calculations for [Pd 1−x Rh x /Fe/Ir] structures [59] reporting an order of magnitude variation of the effective exchange interaction parameter in the magnetic layer. Large variations of the DM interaction have been reported for [Rh/Co/Pt] systems [60] and it has been shown that tuning of exchange interactions is possible for Co-based films [23]. Therefore, it is anticipated that a variation of the composition of [4d/3d/5d] film structures -which are yet to be explored experimentally -could make it possible to reach the regime of isolated zero-field skyrmions which are stable at room temperature.
The results presented here will re-establish interest in ultrathin films as technologically-relevant skyrmionic systems. In skyrmionics, the focus has largely been shifted toward multilayered systems comprising heavymetal/ferromagnet layers, with a large number of repetitions. However, increase in the number of layers in such systems results in large stray fields which can induce twisted spin textures with a nonuniform chirality across the film thickness and, as a result, complicate control over the magnetic structure. Overall increase in the thickness of the system also increases power dissipation associated with current-driven skyrmion motion [42]. These unfavorable effects can be avoided in ultrathin film systems, where nanoscale skyrmions at ambient conditions can be stabilized thanks to ultralow attempt frequency.

METHODS
Simulated system. A two-dimensional skyrmionic system is modeled as a single monolayer of classical magnetic vectors localized on vertices of a hexagonal lattice. The energy of the system is given by Eq.
(2) where each contribution is defined as follows: Here, m i is the unit vector in the direction of the magnetic moment at lattice site i; J and D are the parameters of Heisenberg exchange and DM interaction between nearest neighbor spins, respectively. The unit DM vector d ij lies in the monolayer plane and points perpendicular to the bond connecting sites i and j. The model also includes out-of-plane anisotropy with the easy axis defined by the unit vector e K pointing perpendicular to the monolayer plane and effective parameter K incorporating both magnetocrystalline and magnetostatic contributions [61]. Only one single skyrmion is placed in the simulated system. The size of the computational domain is chosen to be 80×80 lattice sites, which is large enough for the isolated equilibrium skyrmion solution not to be affected by the boundaries. Periodic boundary conditions are applied so as to model extended two-dimensional systems.
Identification of isolines of skyrmion radius. Each isoline α along which the skyrmion radius R [see Eq. (3)] assumes a constant value of R α was obtained numerically using the following technique. At first, profiles of R as functions of the magnetic anisotropy K for several fixed values of the DM interaction, D = D l , are calculated. For each profile, the sought-for value of the anisotropy parameter K l corresponding to the predefined skyrmion radius R α is isolated by interpolation between the data points and then refined using the bisection method until the desired accuracy has been achieved. As a result, a set of P points (K l , D l ) α , l = 1, . . . , P , is obtained, where each point corresponds to a skyrmion with the radius R = R α ± 0.1a. This set of points in the material parameter space gives a discrete representation of the R-isoline. Eight isolines corresponding to R = 5a, 6a, . . . , 12a have been calculated in this manner. They are presented in Supplementary Table 1.
Evaluation of skyrmion lifetime.
The mean skyrmion lifetime which characterizes the stability of the skyrmion state with respect to thermal fluctuations is calculated using harmonic transition state theory (HTST) for magnetic systems [43,44] extended to include the presence of Goldstone modes [33]. HTST presents a rigorous foundation for the Arrhenius law [see Eq. (1)] and provides means for definite evaluation of both the energy barrier ∆E and the pre-exponential factor τ 0 for thermally activated decay processes [43,62]. Within HTST, analysis of skyrmion stability relies on the identification of relevant saddle points (SPs) on the energy surface [see Eqs.
(2), (7)-(9)] characterizing transition states for the skyrmion decay into the ferromagnetic state or other states available in the system: The SP energy relative to the skyrmion energy minimum gives the energy barrier ∆E [see Eq. (4)], but the prefactor τ 0 is defined by the curvature of the energy surface at the saddle point and at the skyrmion-state minimum. HTST predicts the following formula for the pre-exponential factor [33,43]: Here, λ is the dynamical factor defining the system's velocity orthogonal to the surface separating the skyrmion state from the FM state in the configuration space, det H Sk and det H SP denote the determinants of the Hessian matrices at the skyrmion state minimum and the SP, respectively. The Hessians can be obtained using a projection operator approach presented in Ref. [63]. The determinants are computed as a product of the eigenvalues and the prime means that the negative one is omitted. At elevated temperatures which are of interest here, in-plane translations of both the skyrmion and SP configuration should be treated as Goldstone modes with volumes V Sk and V SP , respectively, and corresponding eigenvalues not included in the determinants. Being a quasi-equilibrium theory, transition state theory assumes that Boltzmann distribution is established in the region of configuration space that corresponds to the skyrmion state before the system escapes due to thermal fluctuations. This assumption is usually justified when the energy barriers are large compared to the thermal energy. A more general prerequisite is that the escape events are rare on the intrinsic time scale of magnetization dynamics of the system. For long-lived skyrmions which are of interest here, this condition is expected to be met a priori. Note that the separation of time scales can arise from both an energy barrier and an entropy barrier. For skyrmions, stabilization by entropy barriers is particularly important.
Minimum energy path calculations. Calculation of minimum energy paths (MEPs) is needed for a definitive identification of transition state saddle points, which define the mean lifetime within the harmonic transition state theory. A MEP between two minima is the path in configuration space which lies lowermost on the energy surface. The saddle point corresponds to a maximum along the MEP. A Geodesic nudged elastic band (GNEB) method [24,64], an extension of the nudged elastic band method [65] to magnetic systems, is used to find MEPs of skyrmion annihilation. The GNEB method involves taking some initial guess of the path represented by a discrete chain of states of the system, and systematically bringing that to the nearest MEP by zeroing the transverse component of the gradient force at each point along the path. In order to distribute the states evenly along the path, virtual springs are introduced between adjacent states. At each state, a local tangent to the path is estimated and the GNEB force guiding the states towards the MEP is defined as the sum of the transverse component of the negative energy gradient and the component of the spring force along the tangent. The positions of states are then adjusted using some optimization algorithm, e.g. velocity projection optimization [24], so as to zero the GNEB force. Both the GNEB force and the path tangent are defined in the local tangent space to the curved configuration space, which is needed to satisfy constraints on the length of magnetic moments and to properly decouple the perpendicular component of the energy gradient from the spring force. The shortest path connecting the skyrmion state and the FM state is taken to be the initial guess for the GNEB calculations.   A spin spiral is the ground state of the system in the SS region (blue), while the ferromagnetic state is the only stable configuration in the FM domain (green). Metastable, isolated skyrmions exist in the Sk sector. The color scheme (see top bar) shows the calculated HTST lifetime of isolated skyrmions in terms of the intrinsic precession time τ int = µ/Jγ, which is on the order of femtoseconds for a typical magnetic system. A thermal energy of 2.6J was assumed in the lifetime calculations. Some of the contours of equal skyrmion radius R and lifetime τ are shown with blue and black solid lines, respectively. Material parameter domain corresponding to skyrmions with R < 8a and τ > 10 6 τ int is marked with hatching. Insets show skyrmion profiles and definition of the skyrmion radius for the two points, marked I and II, lying on the same R-isoline. Black and white diamonds indicate effective material parameters for Fe-and Co-based ultrathin skyrmionic systems, respectively, which are listed in the legend. For Ref. [23], the parameters were extracted using the method described in Supplementary Note 2. Color codes distribution of the reduced anisotropy parameter. Filled squares indicate the calculated data points. The magnitude of the energy difference between the Belavin-Polyakov soliton state and relaxed, energy-minimum skyrmion state is within the pink area for all skyrmion radii from 5a to 12a. The inset shows the variation of the saddle point energy (dark grey area) and universal energy of the Belavin-Polyakov soliton (red dashed line) relative to the FM state along the R-isolines. b, Variation of interaction-resolved contributions to the collapse energy barrier due to magnetic exchange ∆E ex , DM interaction ∆E DM and magnetic anisotropy ∆E ani along R = 8a isoline. Filled circles indicate the calculated data points.   The following values of temperature T , Heisenberg exchange parameter J and magnitude of the magnatic moments µ were assumed in the calculations: T = 300 K, J = 10 meV, µ = 3µ B . Color codes distribution of the reduced anisotropy parameter. Filled squares indicate the calculated data points. The inset shows the Arrhenius plot for the two values of the reduced DM interaction parameter on R = 12a isoline. The magnitude of the lifetime calculated in the constant Arrhenius prefactor approximation is within the pink area for all skyrmion radii from 5a to 12a. Blue area indicates lifetimes greater than ten years.