Single Chiral Skyrmions in Ultrathin Magnetic Films

The stability and sizes of chiral skyrmions in ultrathin magnetic films are calculated accounting for the isotropic exchange, Dzyaloshinskii–Moriya exchange interaction (DMI), and out-of-plane magnetic anisotropy within micromagnetic approach. Bloch skyrmions in ultrathin magnetic films with B20 cubic crystal structure (MnSi, FeGe) and Neel skyrmions in ultrathin films and multilayers Co/X (X = Ir, Pd, Pt) are considered. The generalized DeBonte ansatz is used to describe the inhomogeneous skyrmion magnetization. The single skyrmion metastability/instability area, skyrmion radius, and skyrmion width are found analytically as a function of DMI strength d. It is shown that the single chiral skyrmions are metastable in infinite magnetic films below a critical value of DMI dc, and do not exist at d>dc. The calculated skyrmion radius increases as d increases and diverges at d→dc−0, whereas the skyrmion width increases monotonically as d increases up to dc without any singularities. The calculated skyrmion width is essentially smaller than the one calculated within the generalized domain wall model.


Introduction
The individual (single) magnetic skyrmions have attracted considerable attention from researchers assuming potential applications in spintronic and information processing devices [1]. To achieve efficient manipulation of the skyrmion spin textures and to realize skyrmion-based low energy consumption devices, it is essential to understand the magnetic skyrmion stability and dynamics, for instance, in ultrathin ferromagnetic films.
The chiral magnetic skyrmions are a kind of magnetic topological soliton [2] in 2D spin systems characterized by a non-zero skyrmion number (topological charge, degree of mapping) defined as N = d 2 ρm · (∂ x m × ∂ y m)/4π, where m(ρ) = M(ρ)/M s is the unit magnetization vector, M s is the material saturation magnetization, and ρ = (x, y) are in-plane spatial coordinates. The number N = ±1, ±2, . . . is an integer for an infinite film. This topological charge can be interpreted as a quantized flux of the emergent magnetic field [3] through the film surface, Φ = |N|Φ 0 , where Φ 0 = h/e is the flux quantum.
The relativistic Dzyaloshinskii-Moriya exchange interaction (DMI) leads to the stabilization of chiral Neel or Bloch skyrmions with a given sense of the magnetization rotation within their internal configuration [1]. The role of the DMI in skyrmion stabilization was discussed in Refs. [4][5][6][7]. Following the ideas of Dzyaloshinskii [4], in Ref. [5] it was found that adding the term D[m · (∇ × m)] (linear in spatial derivatives of magnetization) to the magnetic energy density of an infinite cubic ferromagnet leads to the stabilization of an inhomogeneous magnetization texture for any finite value of the DMI parameter D. Such terms are allowed in magnetic crystals whose symmetry group lacks the space inversion symmetry operation (e.g., in the B20 cubic crystals MnSi, FeGe, etc. [1]). Then, it was shown [6] that accounting for DMI in the form of the Lifshitz invariants in a bulk uniaxial ferromagnet results in the instability of the uniform ferromagnetic state at D > D c = (4/π) √ AK, where A is the exchange stiffness and K is a uniaxial anisotropy constant. The 1D spin spiral becomes the ground state at D > D c . Therefore, DMI can stabilize 2D vortices (Bloch skyrmions, in modern terminology) for moderate values of D. Ivanov et al. [7] showed that the Bloch skyrmions in infinite films with easy axis anisotropy can be stabilized either by DMI or a high-order exchange interaction. Another kind of single chiral skyrmion (Neel skyrmions) was recently observed at room temperature by Boulle where the unit vector z is normal to the interface. The DMI lowers the skyrmion energy for the proper skyrmion chirality.
In the case of an infinite ferromagnetic film, the critical D value presumably remains the same as for bulk crystals, D c = (4/π) √ AK, although the effective anisotropy constant K is different. The isolated skyrmions are metastable at D < D c at zero external magnetic field, and other configurations (e.g., spin spirals, skyrmion lattices, stripe domains) are stabilized at D > D c [12,13].
In this article, we calculate the magnetic energy of a single chiral skyrmion in ultrathin magnetic film and determine the area of the skyrmion metastability, skyrmion magnetization profiles, and the equilibrium skyrmion radius and width. The case of an effective out-of-plane magnetic anisotropy is analyzed.

Methods
Let us consider an infinite magnetic film with thickness L of about 1 nm, and parameterize the unit magnetization vector by the spherical angles, m = m(Θ, Φ). The spatial distribution of magnetization is assumed to be independent of the thickness coordinate z. The angles Θ, Φ are functions of the polar radius vector ρ = (ρ, φ) located in the film plane. For this kind of magnetization configuration, the total magnetic energy functional is E[m] = L d 2 ρε(m) [6,7], with the energy density where A is the material exchange stiffness, ε DMI is the DMI energy density, with D being the DMI parameter, K u > 0 is the out-of-plane uniaxial anisotropy constant, m z is the magnetization z-component, and ε m is the magnetostatic energy. The interface DMI density is ε DMI (m) = D[m z (∇ · m) − (m · ∇)m z ] for the Neel skyrmions, or ε DMI (m) = D[m · ∇ × m] for the Bloch skyrmions, in thin films of the B20 cubic crystals. The magnetostatic energy ε m (m) is non-local in a general case. The volume and surface magnetic charges contribute to the magnetostatic energy. However, within the limit of ultrathin film, the volume magnetic charges can be neglected, and only surface magnetic charges on the film top/bottom surfaces related to the out-of-plane magnetization component m z contribute to the magnetostatic energy. Then, the magnetostatic energy density can be essentially simplified and written in the local form ε m (m) = µ 0 M 2 s m 2 z /2 [2,7] for both kinds of skyrmion. Therefore, the energy is accounted via an effective uniaxial anisotropy constant K = K u − µ 0 M 2 s /2 > 0. We also define the characteristic magnetic material length l = √ A/K, and the reduced dimensionless DMI strength d = Dl/A. We search for axially symmetric inhomogeneous magnetization configurations (m depends only on the radial coordinate ρ), that is, the magnetization angles are Θ = Θ(ρ), Φ = ϕ + ϕ 0 (ϕ 0 = 0, π for the Neel skyrmions or ϕ 0 = ±π/2 for the Bloch skyrmions). The total skyrmion magnetic energy as a functional of the skyrmion magnetization is represented by the polar magnetization angle Θ(ρ), E = E[Θ(ρ)]. The DMI energy depends on the skyrmion chirality C = ±1, which is defined as C = sin ϕ 0 for the Bloch skyrmions and C = cos ϕ 0 for the Neel skyrmions. The sign of DMI strength D depends on the particular ferromagnetic material. Appropriate choice of the sign of chirality at given D ensures that the product DC corresponds to negative DMI energy. We use the total reduced energy of the radially symmetric Bloch or Neel skyrmion (in units of 2π AL) which depends only on one material parameter: the reduced DMI strength, d.
The simplest magnetization distribution Θ(ρ) = 0 corresponds to the energy E[0] = 0 and describes the magnetic film ground state. However, there are metastable magnetization configurations with non-trivial dependence Θ(ρ), which can be found from the solution of the Lagrange-Euler equation corresponding to the energy functional given by Equation (2). The Lagrange-Euler equation is a non-linear differential equation and cannot be solved analytically. Therefore, we use the different approximate solutions below or trial functions for the skyrmion magnetization profile Θ(ρ). Introducing a trial function (skyrmion ansatz) to the energy functional (2), one can calculate the energy of the skyrmion configuration. The simplest trial function, sometimes used in the theory of domain walls and skyrmions [12], is a linear ansatz: is the skyrmion radius), and Θ(ρ) = 0 otherwise. The simplicity of this ansatz allows conduction of the integration in Equation (2) to get the energy E lin (r s ) = λ + r 2 s − πdr s , where d > 0, λ = 6.154. The skyrmion equilibrium radius r s = R s /l within the model is r s = πd/2, and the skyrmion energy is E lin = λ − π 2 d 2 /4. The linear model predicts that the skyrmion is in a metastable state at d < 2 √ λ/π ≈ 1.58 and that its energy is lower than the energy of the collinear out-of-plane magnetization state at d > 2 √ λ/π. We can write the Lagrange-Euler equation for the function Θ(ρ) to minimize the skyrmion energy (2) using the substitution tan (Θ(r)/2) = exp (− f (r)) [14]. The boundary conditions for the function Θ(ρ) are Θ(0) = π and Θ(∞) = 0 [2,15] or f (0) = −∞ and f (∞) = ∞. We define the skyrmion radius by the equation m z (R s ) = 0, Θ(R s ) = π/2 or f (r s ) = 0, where the reduced radius is r s = R s /l.
The approximate solution of the Lagrange-Euler equation at r >> 1, far from the skyrmion center r = 0, and d = 0, is f (r) = (r − r s ). This is an often-used radial domain wall ansatz taken from the theory of bubble domains in infinite films [16]. This ansatz does not satisfy the boundary condition f (0) = −∞, resulting in singularity of the exchange energy at r = 0. Many authors, including Rohart et al. [17] and Buettner et al. [18], considered the skyrmion magnetization configuration as a circular domain wall (DW) located at the skyrmion radius position R s , described by the singular domain wall ansatz tan (Θ(ρ)/2) = exp (±(ρ − R s )/∆), where ∆ is the wall width. In the limit of large radius skyrmion with a sharp edge R s /∆ >> 1, the radial DW model becomes asymptotically exact. Recently, it was generalized by Kravchuk et al. [15] considering the domain wall width as a variable δ different from its nominal value ∆ = √ A/K. The generalized DW ansatz can be used with caution only within the limit r s /δ >> 1 (i.e., for the large radius skyrmions) if one conducts integration in Equation (2) in the interval r ∈ [r s − δ, r s + δ] near the skyrmion edge. To avoid singularity at the origin r = 0 and describe the whole range of the skyrmion radii r s , we use the trial function f (r) = ln (r/r s ) + (r − r s )/δ suggested by DeBonte [19] to describe the bubble domains in infinite films. Although such a function is not a solution of the Lagrange-Euler equation, it is evident that f (r) leads to finite exchange energy and satisfies the boundary conditions. Below, we use the generalized DeBonte ansatz f (r) = ln (r/r s ) + (r − r s )/δ, where the skyrmion radius r s and the skyrmion width δ are variable and depend strongly on the DMI strength, d. The equalities cos Θ(r) = tanh f (r), sin Θ(r) = 1/ cosh f (r) allow us to calculate the skyrmion energy (2) and find the areas of the skyrmion metastability/stability. We consider that a skyrmion's state is stable when it has the lowest energy (ground state) in comparison with other magnetization states. A skyrmion state is metastable when it corresponds to a minimum of the magnetic energy, however, its energy is higher than that of some other magnetization configurations (a local minimum of the Materials 2018, 11, 2238 4 of 9 energy). The Bloch (Neel) skyrmion energy E(r s , δ) within the generalized DeBonte model is a function of two parameters, r s and δ. Accounting Θ r = −(1/δ + 1/r) sin Θ we rewrite the skyrmion exchange energy in the form The exchange energy (3) was calculated by DeBonte, yielding the simple expression where ξ = r s /∆ ≥ 1 is the reduced skyrmion radius, and 1/∆ = 1/r s + 1/δ is the reduced inverse skyrmion width. In the limit of small radius skyrmion r s → 0 , when the exchange energy dominates over other contributions to the energy density, the exchange energy is reduced to the well-known Belavin-Polyakov soliton limit [20], E ex ( ξ → 1 ) = 4, which is determined solely by the skyrmion charge |N| (|N| = 1 for the skyrmions considered here). The magnetic anisotropy and DMI energy can be represented using DeBonte ansatz as where x = ξ − 1, and the functions F a (x), F(x) are defined as integrals The function F(x) > 0, therefore we chose the sign of dC > 0 and below use the substitution dC → d .

Results and Discussion
The total skyrmion magnetic energy within the model can be represented as function of two variable parameters, ξ and ∆: The equation ∂E(ξ, r s )/∂r s = 0 leads to r s (ξ) = dF(ξ − 1)/2F a (ξ − 1) and allows us to exclude r s from the minimization procedure and write an analytical equation for the equilibrium skyrmion radius as an inverse function of the DMI parameter d(ξ) It immediately follows from Equation (6) that the reduced skyrmion radius ξ is a function of d 2 , ξ = φ(d 2 ), and for large radius skyrmions ξ >> 1, d(ξ >> 1) → d c = 4/π , or the equilibrium skyrmion reduced radius diverges, ξ(d) → ∞ , at d → d c − 0 . In the vicinity of d c , Equation (6) [15]. At ξ >> 1 F a (ξ − 1) = ξ −2 ln (1 + exp (2ξ)) ≈ 2/ξ, F(ξ − 1) ≈ π + O(e −ξ ), the skyrmion energy is essentially simplified, E(ξ, r s ) = 2ξ − πdr s + 2(1 + r 2 s )/ξ, and is reduced to one, accounted in the generalized DW model. We note that the DMI and anisotropy energies are proportional to r s , whereas the exchange energy is not: it contains the term 1/r s even within the simplified DW model. This is in disagreement with the statement by Bernand-Mantel et al. [21] that the exchange energy is linearly proportional to r s . Note that the critical value of d c = 2/π, two times smaller than d c = 4/π ≈ 1.273, was calculated for the isolated chiral skyrmions in infinite films in zero external magnetic field by Kiselev et al. [12], and later this value was corrected by Leonov et al. [13] to be d c = 1.224.
In the limit of small DMI strength d << 1, r s (d) cannot be directly determined from the equation r s (d) = dF(0)/2F a (0) because F(0) = 4 is finite, but F a ( x → 0 ) → ∞ is singular. The non-analytic behavior of the function F a (x) at x → 0 can be approximately presented as F a (x) = F a (1)/x α . To calculate ξ(d), we need to analyze the exchange energy. The approximate Equation (3b) has very good accuracy at ξ ≥ 2, but it predicts a wrong asymptotic behavior at ξ → 1 + 0 and the exact Equation (3a) should be used instead within this limit. We rewrite the exchange energy in the form The functions n (x) are not analytic at x = ξ − 1 → 0 , but it is possible to calculate 0 (0) = π and −1 (x) = −1 (0) + −1 (0)x, −1 (0) = 2, −1 (0) = −π. Therefore, the asymptotic behavior of the function E ex (x) is determined by the first term in Equation (7), Using this expression, we can solve Equation (6) in the limit x = ξ − 1 → 0 and get ξ(d) . Numerical calculation of the asymptote of F a ( x → 0 ) showed that the exponent α = 2/3, and F a (1) = 1.121. Therefore, the skyrmion radius calculated within the generalized DeBonte model at d << 1, r s (d) = ∆(d) = (4/F 3 a (1))d 3 , is essentially smaller than the radius predicted by the generalized DW model. The skyrmion energy is E(d) = 4 1 − F −3 a (1)d 4 . It is slightly higher than the DW model energy E DW (d) = 4 1 − (d/d c ) 2 /2 . This is not a surprise because the DW model containing integration in the vicinity of r s always underestimates the skyrmion energy at d < d c .
The equilibrium skyrmion radius, width, and energy vs. the DMI strength are shown in Figures 1-3 generalized domain wall (DW) ansatz [15] (dashed red line); (3) linear ansatz [12] (dotted blue line). The radius obtained from numerical minimization of the skyrmion energy (2) is shown by deep green squares.   Figure 1). The generalized DW model [15] predicts the skyrmion width for intermediate values of d, which is approximately two times larger than one calculated within the generalized DeBonte ansatz (see Figure 2). The skyrmion energies calculated within the DeBonte and DW models are very close for 0 c dd , whereas the linear model [12] overestimates the skyrmion energy up to 50% and predicts the wrong value of c d (see Figure 3). The skyrmion radius   s rd ( Figure 1) and skyrmion energy (Figure 3) calculated analytically using the DeBonte ansatz and numerically practically coincide.
Above, we calculated the stability of the chiral Bloch and Neel skyrmion magnetization The DW ansatz and linear skyrmion ansatz result in the incorrect dependence r s (d), especially at small d (d < 1) (see Figure 1). The generalized DW model [15] predicts the skyrmion width for intermediate values of d, which is approximately two times larger than one calculated within the generalized DeBonte ansatz (see Figure 2). The skyrmion energies calculated within the DeBonte and DW models are very close for 0 < d < d c , whereas the linear model [12] overestimates the skyrmion energy up to 50% and predicts the wrong value of d c (see Figure 3). The skyrmion radius r s (d) (Figure 1) and skyrmion energy (Figure 3) calculated analytically using the DeBonte ansatz and numerically practically coincide.
Above, we calculated the stability of the chiral Bloch and Neel skyrmion magnetization configurations in ultrathin films as a function of the DMI strength. The second derivative of the skyrmion energy (5) ∂ 2 E/∂r 2 s = 2F a (x) > 0. Therefore, the sufficient condition of existence of the skyrmion local energy minimum (∂ 2 E/∂r 2 s )(∂ 2 E/∂ξ 2 ) − (∂ 2 E/∂r s ∂ξ) 2 > 0 is satisfied for the skyrmion solution ξ(d), r s (d) within the interval 0 < d < d c . The isolated skyrmions are metastable within a range of the values of d satisfying the inequality d < d c and do not exist at d > d c (the skyrmion minimum transforms to an energy maximum at d = d c ).
To describe skyrmion magnetization analytically we used the DeBonte radial domain wall ansatz [19], the accuracy of which was numerically checked for circular dots in Ref. [22]. The calculated equilibrium skyrmion radius R s (d) and the skyrmion width ∆(d) increase with increasing DMI strength ( Figure 1). However, the continual model becomes inaccurate for the sizes below 1 nm.
The typical values of A= 10 pJ/m and K= 0.1 MJ/m 3 yield the magnetic length l = 10 nm for ultrathin films. The conditions R s (d) ≥ 1 nm, ∆(d) ≥1 nm mean that the continual model can be applied if the reduced DMI strength d ≥ 0.2 or |D| ≥ 0.2 mJ/m 2 in absolute units. Simulations [23] within a discrete model on a simple cubic lattice with period a showed that the skyrmion state collapses to the uniformly magnetized state at R s ≈ (4 ÷ 5)a or R s ≈1.0-1.3 nm for Co. We note that in restricted geometry (circular dots) the skyrmion radius dependence on the DMI strength R s (d) has an inflection point at d ≈ d c [17,24] and the skyrmion width ∆(d) reveals a broad maximum in the vicinity of d c [24]. The typical value of DMI strength, D, accessible in experiments with ultrathin films like X/Co (X = Pt, Ir, Pd) is 1-2 mJ/m 2 [8][9][10][11]. Therefore, all observed Neel skyrmions in these nanostructures are metastable  [11], we could calculate l = 6.8 nm and d = 1.231. The skyrmion radius, measured by Lorentz transmission electron microscopy [11], is R exp s = 45 nm, or r exp s = R exp s /l = 6.6, whereas the calculations yielded the value r cal s = 4.1. The agreement of the skyrmion sizes measured by X-ray imaging [8,9] and by our calculations is reasonably good. The skyrmion size measured by Lorentz transmission microscopy was larger than the calculated one. This can be explained by the different mechanisms of image formation in these experiments. The image contrast is proportional to the magnetization out-of-plane component m · z for the X-ray imaging [8-10], whereas the contrast is proportional to the out-of-plane component of the magnetization curl (∇ × m) · z, for Lorentz microscopy imaging. The parameters K and D can be extracted with reasonable accuracy from independent experiments. The exchange stiffness A is poorly defined for ultrathin films with ferromagnetic layer thickness 0.5-1 nm. The skyrmion sizes measured in Refs. [8,9,11] are quite large, 40-65 nm. This means that the DMI parameter d is also large and close to its critical value d c , and the value of the skyrmion radius is very sensitive to the exact value of d (see Figure 1). According to its definition d = |D|/ √ AK, the DMI parameter depends on A. This leads to an uncertainty in the interpretation of the experimental data [8, 9,11]. This uncertainty may lead to the case d > d c for Ir/Co/Pt multilayer films [9]. Decreasing the out-of-plane magnetic field can essentially increase the skyrmion sizes (see Figure 2 in ref [9]), indicating that the single skyrmion state is unstable in zero out-of-plane field. We note that the dependences of the skyrmion radius R s on the DMI strength D for different values of A, simulated in Ref. [11], can be reduced to the universal curve R s (d) presented in Figure 1 if one changes the variable D to dimensionless variable d.
The case of magnetic dots considered in Refs. [14,17,22,24,25] is more complicated because the skyrmion configuration can be the dot ground state. The Neel skyrmions in circular dots can be metastable or stable even at D > D c = (4/π) √ AK. The calculated value of D for a transition between the metastable and stable Neel skyrmions in ultrathin circular dots is 1.5-2 times larger than one for infinite films for weak effective magnetic anisotropy 2K/µ 0 M 2 s << 1 [14]. It was also shown that the Bloch skyrmions can be the dot ground state for in-plane magnetic anisotropy K < 0 and D = 0 [25].
In the investigated case of out-of-plane effective magnetic anisotropy K > 0, the large values of the Dzyaloshinskii-Moriya interaction strength D > D c cause the nucleation of more complicated magnetization configurations (nπ-skyrmions [17], spin spirals, labyrinth domain, etc.), that is, the individual Neel or Bloch magnetic skyrmion state with the topological charge |N| ≈ 1 is no longer metastable.

Conclusions
We found that the isolated Bloch and Neel skyrmions in ultrathin magnetic films are metastable within the range of the DMI strength 0 ≤ d < d c , where d c = 4/π or D c = 4A/πl in absolute units, A is the material exchange stiffness, and l = A/(K u − µ 0 M 2 s /2) is the material magnetic length. The calculated skyrmion radius R s increases as d increases and diverges at d → d c − 0 , whereas the skyrmion width ∆ increases monotonically as d increases without any singularities at d → d c − 0 . The calculated skyrmion width is essentially smaller than the one calculated within the generalized domain wall model. The generalized DeBonte ansatz is a very good approximation to calculate the skyrmion radius, width, and energy. The linear skyrmion model cannot be applied for quantitative analysis of the skyrmion energy and size.
Author Contributions: A.A. performed the numerical simulations; K.G. performed analytical calculations; A.A. and K.G. analyzed the data; K.G. supervised the research and wrote the paper.
Funding: K.G. acknowledges support by IKERBASQUE (the Basque Foundation for Science). This research was funded by the Spanish MINECO grant FIS2016-78591-C3-3-R and the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant No. 644348.

Conflicts of Interest:
The authors declare no conflict of interest.