Full phase diagram of isolated skyrmions in a ferromagnet

Magnetic skyrmions are topological quasi particles of great interest for data storage applications because of their small size, high stability, and ease of manipulation via electric current. Theoretically, however, skyrmions are poorly understood since existing theories are not applicable to small skyrmion sizes and finite material thicknesses. Here, we present a complete theoretical framework to determine the energy of any skyrmion in any material, assuming only a circular symmetric 360$^\circ$ domain wall profile and a homogeneous magnetization profile in the out-of-plane direction. Our model precisely agrees with existing experimental data and micromagnetic simulations. Surprisingly, we can prove that there is no topological protection of skyrmions. We discover and confirm new phases, such as bi-stability, a phenomenon unknown in magnetism so far. The outstanding computational performance and precision of our model allow us to obtain the complete phase diagram of static skyrmions and to tackle the inverse problem of finding materials corresponding to given skyrmion properties, a milestone of skyrmion engineering.

While homochiral skyrmions were first observed in bulk materials with broken inversion symmetry, multilayers stacks, such as such as Pt/Co/Ir [17], Ta/CoFeB/TaO x [21,22], Pt/Co/Ta [19], and Pt/CoFeB/MgO [19,23] with arbitrary repetitions of these layers, have seen increasing popularity recently. In such structures, inversion symmetry breaking and spin-orbit coupling at the ferromagnet/heavy-metal interface can lead to a strong DMI that promotes magnetic skyrmions with well-defined chirality. Multilayers are particularly attractive since the relevant energy terms (e.g., anisotropy, DMI, magnetostatic) can be engineered by controlling the interfaces and the volume fraction of magnetic material in the stack. However, with this flexibility comes a considerable engineering challenge in that the parameter space for engineering is overwhelmingly large: It has six dimensions (five material parameters plus external magnetic field) and thanks to the interfacial origin of many magnetic properties and to the effective medium scaling of these properties with the thickness of non-magnetic spacer layers [19,24], all of these dimensions can be tuned individually. Blind experimental investigation of this parameter space is hence prohibitive. Existing theories are of limited help, since they either involve crude approximations that fail to reproduce the wealth of skyrmion states even qualitatively [12,[25][26][27][28] or contain unsolved complicated integrals and differential equations, which renders the evaluation of the theory extremely computationally expensive and slow (for instance, micromagnetic simulations and Refs. [29][30][31]). In particular, most theories ignore the non-local nature of stray field interactions [28,[30][31][32], which is responsible for many interesting features in the intermediate film thickness regime. A single coherent theory that quickly and accurately predicts the existence and the properties of isolated skyrmions for any given point in the six-dimensional parameter space remains elusive. Here, we provide a theoretical model for all energy terms of an isolated skyrmion in a given material with infinite in-plane extent. The model is fully analytical and accurate to 1 % in the entire parameter space. Thanks to the analytical nature, we can find energy minima extremely quickly by searching for roots of the partial derivatives. The resulting equilibrium states show excellent agreement with simulations and experiments. Using our theory, we find exotic new states, such as multi-stabilities (co-existence of skyrmions with different properties in the same sample under the same conditions), zero stiffness skyrmions, and zero field skyrmions. These new states can have many novel applications, some of which we suggest here. We obtained the minimum energy skyrmion states for millions of material parameter and field combinations in less than a week, yielding the full phase diagram of skyrmion states and demonstrating that our model is suitable to solve the inverse problem of engineering skyrmion properties through material selection.
Our theory takes as input the uniaxial anisotropy constant K u , saturation magnetization M s , exchange constant A, interface and bulk DMI strengths D i and D b , magnetic layer thickness d, and applied out-of-plane field H z . Given these parameters, we derive the energy function that determines the spin structure of skyrmions in any material. In general, stable skyrmions correspond to sufficiently deep minima of the energy functional E[m] of all possible spin structures m(r), where m is the unit magnetization vector and r is the position vector. In practice, minimizing the energy functional in its raw form, as done in micromagnetic simulations, is prohibitively slow for systematic skyrmion engineering. Our simple and efficient analytical model is enabled by the recent experimental confirmation [15,18] of an analytic and universal 360°DW model [33] for the spin structure m(r) of all skyrmions. The full analytic model for the total energy function is provided in the supplemental information, along with a detailed discussion of how to solve the integrals of the individual interactions. Here, we focus on its implications.

Effective energy contributions
The skyrmion spin structure (Fig. 1a) can be accurately described by four parameters: radius R, DW width ∆, DW angle ψ, and topological charge N. R and ∆ are independent parameters that determine the magnetization profile m z (x, y), whereas ψ specifies whether the in-plane component of the DW spins is radial (Néel, ψ = 0, π), azimuthal (Bloch, ψ = π/2, 3π/2), or intermediate (transient). For large ρ = R/∆, skyrmions consist of an extended out-of-plane magnetized domain bounded by a narrow circular DW, while for ρ ∼ 1 the inner domain is reduced to a point-like core resembling a magnetic vortex. We refer to these limiting cases as bubble skyrmions and vortex skyrmions, respectively, consistent with the literature, but note that many skyrmions observed recently [8,17,19,34] showed intermediate values of ρ and cannot be classified distinctly.
The case of bubble skyrmions is readily treated analytically through the so-called wall-energy model known for a long time [35]. In this limit, the skyrmion energy simplifies to with 2πdσ DW R being the DW energy (σ DW is the energy density of an isolated DW), aR + bR ln(R/d) the Zeeman-like surface stray field energy, and cR 2 the Zeeman energy. Here a, b, c are material-dependent parameters that include corrections to the original model [35] to account for DMI, volume charges, and finite ∆ (see supplemental information). The crucial assumption of Eq. (1) is that σ DW , ∆, and ψ do not depend on R and that hence ∆ and ψ can be obtained from minimizing σ DW while R is the minimum of E assuming a constant σ DW . The simplicity of this wall-energy model is reason for its popularity [28,36]. However, our accurate model shows that ∆ can change with R even up to unexpectedly large radii of R = 100∆. Indeed, we find that the wall-energy model fails quantitatively and qualitatively for almost all skyrmions that are of interest today. Our model fully accounts for the correlation between radius, DW width and DW angle, thereby providing a single theoretical framework that accurately describes any skyrmion. Importantly, our model reveals unexpected and qualitatively new exotic behaviors that are precluded by the approximations inherent in prior treatments.
Our analytic equations for the total energy function can be numerically minimized with respect to R, ∆, and ψ, for a given set of material parameters and external field H z , to obtain the equilibrium skyrmion configuration. The predictions of our model agree precisely with micromagnetic simulations and with the experimental data of Romming et al. [15], see Fig. 1b. Note that fields are negative in our convention, antiparallel to the skyrmion core. Heuristically, we find that the function fits the field dependence of the equilibrium radius for a wide range of parameters.
Our model yields the energy of a given skyrmion configuration in less than a millisecond on a regular personal computer, thereby providing dramatic improvement over micromagnetic simulations in terms of  [15], and the solid grey line is a fit with the simple model of Eq. (2). The insets show the relaxed spin structures obtained from micromagnetic simulations corresponding to the large solid data points. c, Total energy E, normalized to Ad, as a function of skyrmion radius at a field of µ 0 H z = −2 T. At each point of the solid red line the energy has been minimized with respect to ∆ and ψ while keeping R fixed. R eq is the equilibrium radius as plotted in b. d-f, Decomposition of the total energy in c into individual components. d, DW energies: DMI energy (inverted), exchange energy, effective anisotropy energy, and volume stray field energy (multiplied by 10). e, Bulk energies: effective Zeeman energy and remaining surface stray field energies (multiplied by 10). f, Sum of all DW and sum of all bulk energies, together with the total energy. In this particular example, the sum of DW energies has a negative slope at the minimum of the total energy, which makes this skyrmion DMI stabilized according to our definition. The parameters for this data set are experimental values from Ref. [15] and our predictions are in excellent agreement with the experimental observations. computation speed. Moreover, it gives access to information that cannot be obtained by simulations. For example, by virtue of taking partial derivatives, the model allows to minimize ∆ and ψ for any non-equilibrium skyrmion radius and, therefore, to obtain the energy as a function of radius E(R). Micromagnetic simulations can only yield the equilibrium R, ∆, ψ. The E(R) curve directly relates to the skyrmion stability (by quantifying the energy barriers) and to its rigidity (by the curvature near the energy minimum). Figure 1c shows E(R) calculated for the skyrmion described in Fig. 1b at µ 0 H z =-2T, which exhibits a single minimum corresponding to an isolated skyrmion. The only other stable state is the ferromagnetic ground state at E = 0. Despite the different topology of the skyrmion and the ferromagnetic state, there is a continuous path from one to the other, which goes through the singular R = 0 state. The singluar R = 0 state does not have a topology, which is why topological quantization is lifted here. Remarkably, the energy along the path towards R = 0 remains finite, in contrast to earlier believes [28]. At R = 0, the skyrmion energy takes a universal value, The zero radius energy depends only on A and d and not on DMI. By finding a topologically valid and energetically possible path to annihilation we prove that skyrmions are not protected by their topology, even in continuum models and in the presence of strong DMI. Note that a finite value for the energy barrier has been found before [32,37], but the role of DMI and implications for the topological stability were not discussed.
The skyrmion of Fig. 1c exhibits a finite annihilation energy E a = E 0 − E(R) and a nucleation energy barrier E n = E 0 . Note that E n and E a overestimate the energy barriers because skyrmions can deform in a way that is not covered by the 360°DW model underlying our calculations and therefore reduce the nominal energy barrier [38]. However, previous studies [38] and our own micromagnetic simulations indicate that the reduction of the energy barrier due to deformations is smaller than 2Ad (even though sometimes extremely small cell sizes are required). This is in excellent agreement with the observation of Belavin and Polyakov [37] that the minimum energy of a skyrmion state in a model that includes only exchange energy is 8π Ad, which is approximately 2Ad smaller than the zero radius exchange energy of our 360°DW theory. Including thermal effects, we therefore consider minima of the total energy to be stable in our discussions below if All energy contributions to the total energy can be classified into two categories, DW and bulk energies, and inspection of the individual energy terms allows one to identify the mechanism responsible for skyrmion stability. At large radii, DW energies are linear in R, whereas all non-linear terms are bulk energies. Exchange, anisotropy, DMI, and volume stray field energies are DW energies. The Zeeman energy is a bulk energy.
Surface stray fields contribute to both categories: The DW contribution leads to an effective reduction of the anisotropy and the bulk contribution effectively reduces the external field.
The decomposition of the total energy into DW and bulk terms is depicted in Figs. 1d- can exist only if energy terms with positive slope are compensated by terms with negative slope, and the latter can only arise through DMI and surface stray fields. One can therefore classify skyrmions as DMI (stray field) stabilized if the sum of DW (bulk) energies has a negative slope at the equilibrium radius R eq . Our model hence provides the first mathematical basis for a terminology commonly used in the literature without rigorous justification [1].

New phases
We now apply our model to gain further insight into skyrmion properties and to analyze the phase diagram of static isolated skyrmions. Features of particular interest are highlighted in Fig. 2. First, Fig. 2a (Figs. 2b, c). We find this pocket in the stray-field stabilized part of the phase diagram where the phase boundary of the instability region has a cusp. The two types of skyrmions in the bi-stability region have very different properties (Fig 2d), confirmed by micromagnetic simulations: Their radii differ by more than one order of magnitude and their spin structure is Néel-like for the small skyrmion and transient for the large skyrmion. The transient value of ψ for the larger skyrmion originates in a wider domain wall, which increases the importance of the volume stray field energy that favors a Bloch-like spin orientation. The different size and domain wall angle can be used to move the skyrmions in non-collinear directions by spin orbit torques, as detailed in the supplemental information.
The unexpected emergence of multiple minima in E(R) is a consequence of introducing ∆ and ψ as free parameters, resulting in ∆(R) being nonlinear and sometimes nonmonotonic. The existence of degenerated isolated domain states is new in the entire field of magnetism. Decomposition of E(R) into DW and bulk energies (Fig. 2e) reveals that the origin of this phenomenon is that each term individually exhibits a minimum.
The minimum in E(R) at small radius derives from the minimum in the DW energy terms, shifted towards larger radii by the negative-sloped bulk energies and vice versa for the large radius minimum. This observation helps explain why the phase boundaries between stray field stabilized and DMI stabilized skyrmions in Fig. 2b are vertical and horizontal (see also supplemental figure S1). The horizontal line marks the critical DMI value above which σ DW is negative, i.e., where the DW energies have a negative slope everywhere and all minima are DMI stabilized. The vertical line indicates the critical field value above which the applied field fully compensates the Zeeman-like surface stray field, meaning that the bulk energies are always positive with a positive slope beyond that point and again all minima are DMI stabilized. In any of these cases, either the DW or the bulk energies cease to have a minimum, which finally also explains why we find the bi-stable phase pocket in the stray field stabilized phase.
The last peculiar phenomenon we uncover in our analysis is the existence of zero stiffness skyrmions. Figure 2f shows E(R) for a system that manifests three energy minima, where the maxima between these minima can easily be overcome by thermal energies at room temperature. In this particular example, the skyrmion radius can thermally fluctuate between 2 nm and 11 nm, such that it exhibits effectively zero stiffness with respect to variations of radius within this range. We expect that such skyrmions have a very low resonance frequency associated with their breathing mode, which could be exploited in skyrmion resonators [40] and should have impact on their inertia [8] and on skyrmion Hall angle [22,23]. a and E ∞ a minus the estimated internal energy of 2Ad that can be drawn from deformations, at given material parameters. The light grey area indicates that no minima exist at zero fields or that E eff a is smaller than 2Ad + k B T. The dark grey area illustrates where spontaneous domain nucleation is expected at zero field, which is excluded from being useful for isolated skyrmion applications. b, The inset shows a slice of panel a at Q = 1.8, plotting the annihilation energy barrier and the equilibrium skyrmion radius. The maximum energy barrier and the corresponding skyrmion radius are highlighted by a red and grey circle, respectively. This maximum energy barrier and corresponding radius are plotted in the main panel as a function of Q (effectively showing a slice of panel a when moving along the maxima of E eff a ). The solid lines correspond to the regular material parameters used throughout this paper, the dotted lines correspond to a reduced M s and the dashed lines represent a larger exchange constant. All lines end when the maximum E a is found at D i > 4 mJ/m 2 , which to the best of our knowledge is the maximum DMI that can be engineered. If larger DMI values were allowed, all energy lines would continue to go up and all radii lines would continue to go down.

Applications
We now consider the design of skyrmions suitable for applications, such as racetrack-type memory devices in which bit sequences are encoded by the presence and absence of skyrmions that can be shifted by electric current [9][10][11]41]. Three key attributes for such applications are (i) small bit sizes, (ii) long term thermal stability, and (iii) skyrmion stability in zero applied field. Indeed, we find a section in the phase diagram that meets all these requirements, as illustrated in Fig. 3. As sketched in the inset of Fig. 3a, zero field skyrmions are local energy minima bounded by two annihilation energy barriers E 0 a and E ∞ a that prevent shrinking to zero size and infinite growth, respectively. By subtracting the internal deformation energy 2Ad from the minimum of E 0 a and E ∞ a , we obtain the effective annihilation energy barrier E eff a that can be used to estimate long term thermal stability. For properly chosen material parameters, E eff a can exceed the threshold 40k B T ( Fig. 3b) required for commercial storage devices. The corresponding radius is < 20 nm. Note that all energy contributions generally become larger at larger radii. Therefore, fluctuations of any material parameter or of the field affect the energy barrier E ∞ a much stronger than E 0 a , see inset of Fig. 3a. This is why E eff a increases slowly for increasing DMI values before it drops very rapidly after passing the maximum of E eff a , see inset of Fig. 3b. It is therefore advantageous to have a DMI value slightly below the maximum E eff a to ensure that E eff a is robust against moderate external fields. Finally, all zero field skyrmions have E > 0, consistent with earlier assessments excluding ground state zero field skyrmions below the Curie temperature [16,42]. It is reasonably easy to obtain E eff a = 10Ad, which explains why the largest energy barriers can be obtained when increasing A (or d).

The full phase diagram
To demonstrate the power of our model and to understand the effect of material parameters on skyrmion properties, we derived and analyzed the properties of skyrmions as a function of more than one million different material parameters and magnetic fields, a task that is impossible with existing theoretical tools. Fig. 4 illustrates some of the most interesting features of the derived full phase diagram. Specifically, the figure analyzes a magnetic multilayer in which the non-magnetic spacer layers (e.g., Pt and Ta) are three times as thick as the magnetic layers. All room-temperature skyrmion systems today are based on such multilayers [8,17,19,[21][22][23]. We employ the effective medium approach [19,24] to treat these multilayers. The total magnetic material thickness in such films is between 1 nm and 100 nm, the interfacial DMI strength is below 4 mJ/m 2 and the anisotropy quality factor Q = 2K u µ 0 M 2 s is typically between 1 and 2. For this parameter range, we show in the first two rows of Fig. 4 the radius R and the DW angle ψ of the smallest possible skyrmions, i.e., under the maximum field just before collapse. Below, we present the means of stabilization (DMI or stray fields). Similar diagrams with variable A, M s , and non-magnetic spacer layer thickness are provided in the supplemental information.
The most striking common feature of all panels in Fig. 4   Phase diagram as a function of total magnetic layer thickness d, interfacial DMI strength D i , and anisotropy quality factor Q for a multilayer with three layers of non-magnetic material for every layer of magnetic material. The panel rows show, from top to bottom: the radius R and the DW angle ψ of the smallest skyrmion just before collapse, and the possible means of stabilization (where "DMI" indicates that all skyrmions in this material are DMI stabilized, "SF" indicates that all skyrmions are stray field stabilized and "both" means that both types of stabilizations can be found, depending on the external field). The solid line in the ψ panels depicts the minimum DMI strength D SW cψ required to obtain fully Néel skyrmions (ψ = 0) in all skyrmions, independent of their size. but quantitatively unrelated trend: Above a critical DMI value of D cψ , skyrmions are of Néel type when they collapse. The transition in ψ is not as sharp as for R and, importantly, D cψ is consistently smaller than D cr , implying that extremely small skyrmions are always of Néel type. Note also that small skyrmions are more likely to be of Néel type than straight walls in the same material. In other words, the critical DMI value for finding isolated Néel walls D SW cψ is much larger than D cψ . The region between D cψ and D SW cψ is where the dependence of ψ on the skyrmion size is most pronounced and bi-stable states are most likely to have different spin orientations.
Sub-10 µm skyrmions exist in almost the entire phase diagram. Materials with purely DMI stabilized skyrmions exist, but almost exclusively at very small values of Q. Already at Q = 1.4, which is a typical value for cobalt based multilayers [19,23], purely DMI stabilized skyrmions exist only for DMI values larger than 4 mJ/m 2 , well beyond experimentally-reported values. Hence we conclude that most skyrmions investigated experimentally so are best described as stray field stabilized.
Finally, comparing the additional phase diagrams in the supplemental information, we can qualitatively note that small skyrmions are favored by a low anisotropy, a low M s , a small exchange constant, a large DMI value, and sizable non-magnetic spacer layers. Low M s , low A and thick spacer layers also lead to more abundant bi-stability regions, but note that here larger Q values are beneficial.

Conclusions and outlook
In summary, we presented an analytical model that allows exploration of the entire static phase diagram of isolated magnetic skyrmions via rapid, systematic calculations. We expect many new applications to arise from the exotic states found here, beyond what we already suggested. In principle, our model assumes infinite films and the behaviour in finite sized elements can be different [28]. However, in most cases confinement increases the stability of skyrmions, as long as skyrmions still fit into the element. Therefore, our predictions can be considered conservative and applicable to most nanostructures as well. Also, skyrmions in antiferromagnets [43,44] are covered by our theory by setting M s to zero. Still, some open challenges remain.
For instance, the dynamics of skyrmions, and the effects of in-plane fields, are not yet covered by our model. But we believe that the concepts presented here to solve the integral equations pave the path to tackle those issues as well.