A theory on skyrmion size

A magnetic skyrmion is a topological object consisting of a skyrmion core, an outer domain, and a wall that separates the skyrmion core from the outer domain. The skyrmion size and wall width are two fundamental quantities of a skyrmion that depend sensitively on material parameters such as exchange energy, magnetic anisotropy, Dzyaloshinskii–Moriya interaction, and magnetic field. However, quantitative understanding of the two quantities is still very poor. Here we present a general theory on skyrmion size and wall width. The two formulas we obtained agree almost perfectly with simulations and experiments for a wide range of parameters, including most of the existing materials that support skyrmions. Skyrmions are magnetic topological features which are expected to play an important role in future data storage and information processing devices. The authors outline a theoretical method to calculate the size and wall width of an isolated skyrmion.

Although much knowledge about the DMI-supported magnetic skyrmions has been accumulated after intensive studies including skyrmion generation 10,13,21,22 and manipulation 18,23,24 , the dependence of skyrmion size (R) on material parameters such as exchange energy, magnetic anisotropy energy, and DMI strength is still poorly understood at a quantitative, or even qualitative level 18,[25][26][27][28][29] . Even worse, the physical pictures behind most previous studies are not clear. In this paper, we confirm by micromagnetic simulations that the skyrmion profiles agree well with the 360°domain wall formula 30,31 in a wide extent. By using this formula as an ansatz and minimizing the energy with respect to the skyrmion size R and the skyrmion wall width w as two independent quantities, we obtain the analytical expressions of R and w as functions of A, D, K, and B. Interestingly and in contrast to previous theories 25,26,29 , R always decreases with A. In general, both R and w increase with D. These results agree perfectly with micromagnetic simulations and are consistent with experiments.

Results
Energy functional and approximate profile of isolated skyrmions. We consider an ultra-thin ferromagnetic film in xy plane with an exchange constant A, an interfacial DMI coefficient D, a perpendicular easy-axis anisotropy K, and a perpendicular magnetic field B. The total energy E of the system consists of the exchange energy E ex , the DMI energy E DM , the anisotropy energy E an , and the Zeeman energy E Ze , where ÞdS. m is the unit vector of magnetization of a constant saturation magnetization M s and the integration is over the whole film. The energy reference is chosen in such a way that the energy of single domain state of m z = 1 is E = 0. We have assumed m is uniform in thickness direction. The demagnetization effect is included in E an through the effective anisotropy K ¼ K u À μ 0 M 2 s =2 corrected by the shape anisotropy, where K u is the perpendicular magnetocrystalline anisotropy, which is a good approximation when the film thickness d is much smaller than the length w over which magnetization varies substantially. It is convenient to use polar coordinates so that a point r in the plane is denoted by r and ϕ. Magnetization at r is described by polar and azimuthal angles Θ (r, ϕ) and Φ(r, ϕ) so that m = (sin Θ cos Φ, sin Θ sin Φ, cosΘ). A skyrmion centered at r = 0 can be described 19 by, with boundary conditions of Θ(0) = 0(π) and Θ(∞) = π(0). v is the vorticity (v = 1 for a skyrmion and v = −1 for an antiskyrmion), and γ is a constant classifying type of skyrmions (γ = 0 or π for Néel skyrmions and γ = ± π/2 for Bloch skyrmions). Similar to the magnetic bubbles 14,16,17 that refer normally to local metastable structures from dipole interactions, a skyrmion consists of a skyrmion core, an outer domain, and a wall separating the core and the outer domain. We define the skyrmion size R as the radius of the m z = 0 contour. The wall width w is another fundamental skyrmion quantity often ignored or not treated properly in previous studies [25][26][27][28] . One can also define the skyrmion polarity as p = [m z (r = ∞) − m z (r = 0)]/2 so that p = 1 (− 1) corresponds to spins in the core and outer domain pointing respectively to the − z(+ z) and + z(− z) directions.
In terms of Θ, four energy terms for v = 1 are For a skyrmion of p = 1 as shown in Fig. 1a, Θdr is small because sin Θcos Θ = 0 far from the skyrmion wall and sin Θ cos Θ changes its sign from negative near the skyrmion core to positive near the outer domain. Thus, to lower the total energy, one needs γ = 0(π) for D > 0 (< 0). This corresponds to a Néel skyrmion. Along a radial direction, the magnetization profile looks like a 360°Néel domain wall as illustrated in Fig. 1b. This leads us to model a skyrmion profile by a 360°domain wall originally proposed by Braun 31 (see Supplementary Note 1), It is worth mentioning that 360°domain wall profile has already been used to describe the skyrmion profile obtained in experiments 30,32 .
To test how good ansatz (7) is for a skyrmion, we use MuMax3 33 to simulate various magnetic stable states (see Methods). A = 15 pJm −1 , M s = 580 kAm −1 , and perpendicular easy-axis anisotropy K u = 0.8 MJm −320 are used to mimic Co layer in Pt/Co/MgO system. The initial state is m z = 1 for r > 10 nm and m z = − 1 for r ≤ 10 nm. The final stable state depends on the values of D and B. The lower left inset of Fig. 1 is three typical stable states. Figure 1c is a skyrmion for D = 3.7 mJm −2 and B = 0. Figure 1d is a singledomain state of m z = 1 (or m z = − 1) for D = 0 and B = 0. Figure 1e is a stripe domains state for D = 5 mJm −2 and B = 0. Figure 1b shows the spatial distribution of m z of the skyrmion in Fig. 1c along three radial directions, ϕ = 0 (crosses), 45°(triangles), and 90°(circles). All three sets of data are on the same smooth curve, showing m z is a function of r, but not ϕ. The curve can fit perfectly to Eq. (7) with R = 25.77 nm, w = 4.94 nm. We plotted also Φ(ϕ) at randomly picked spins from the simulated skyrmion. All numerical data (red circles) are perfectly on Φ = ϕ. These results not only confirm the validity of Eq. (2), but also suggest that m z (r) ≡ cos[Θ dw (r)] follows the 360°DW profile (7), although it is not exact that can be easily verified by substituting Eq. (7) into the nonlinear partial differential equations for skyrmions 29 .
Skyrmion size. The energy of a skyrmion can then be obtained from Eq. (3) by using the 360°domain wall profile Θ dw (r). Substituting Θ dw (r) into Eqs. (1) and (3), the total energy is, in general, a function of R and w [instead of a functional of Θ(r) and Φ(ϕ)] as where f i (x) (i = 1~6) are non-elementary functions defined in Methods. The skyrmion size and wall width R and w are the values that minimize E(R, w). In ref. 25 , w was assumed to be a fixed value such as w ¼ ffiffiffiffiffiffiffiffiffi A=K p that eventually leads to a wrong skyrmion size. Figure 2 are D− (a), A− (b), K− (c), and B− (d) dependences of the skyrmion size R (left y axis) and the skyrmion wall width w (right y axis), with other parameters fixed to the values for Co mentioned earlier. The symbols are the micromagnetic simulation data [R is the size of m z = 0 contour and w is the fit of the skyrmion profile to Θ dw (r)]. The sample size used in the simulations is large enough to avoid the boundary effects (see Supplementary Fig. 1). Solid lines are numerical results from Eq. (8). The simulation results agree almost perfectly (except slight deviation in the D-dependence of w for smaller D) with our analytic results of Eq. (8). Both micromagnetic simulations and analytical results clearly show that the skyrmion can exist for D < 3.8 mJm −2 in the current case. Above the upper limit, the stable state is not a skyrmion, but stripe domains as shown in Fig. 1e for D = 5 mJm −2 . E(R, w) of (8)
[subfigures a-c] and solution of Eqs (11) and (12) (11) and (12) for non-zero B] (red lines), and previous studies ref. 25  has a minimum as long as |D| < 3.8 mJm −2 that indicates existence of the skyrmion. However, micromagnetic simulations show that the skyrmion can only exist in the window of 1.2 mJm −2 < D < 3.8 mJm −2 when the skyrmion size is larger than 1 nm in the current case. Below 1.2 mJm −2 , the stable state is a single domain with all spins pointing up or down as shown in Fig. 1d. This discrepancy may be due to the discretization of continuous LLG equation in micromagnetic simulations. In principle, the mesh size should be much smaller than the skyrmion size. For a skyrmion of 1 nm, our mesh size is 0.01 nm. Due to the limited precision of the MuMax3 package, the mesh size cannot be further decreased. There exists also a minimal A of around 14 pJm −1 as shown in Fig. 2b and a minimal K of around 0.56 MJm −3 as shown in Fig. 2c, below which the skyrmion does not exist, and the stable state is stripe domains as shown in Fig. 1e. The skyrmion size decreases with B, which is consistent with the experimental observations 7,8,30 . It is still unclear how R and w depend on A, D, and K although Eq. (8) agrees almost perfectly with the simulation results. Thus, it is highly desirable to have simple approximate expressions for R and w in terms of material parameters. The exchange and DMI energies come from the spatial magnetization variation rate. For a skyrmion, the magnetization variation rates in the radial and tangent directions scale respectively as 1/w and 1/R. The exchange energy is then proportional to the skyrmion wall area of πRw multiplying the square of the magnetization variation rates 1/R 2 + 1/w 2 , i.e., E ex ∝ (R/w + w/R). The integrand of the second term (tangential term) of E DM is proportional to sin2Θ due to the triple product of the interfacial DMI vector (along the radial direction), the tangential variation of m (along the tangential direction), and m itself. As sin2Θ is almost asymmetric about r = R (positive when r > R and negative when r < R) and is non-zero only near r = R, after integral over r, the second term in E DM is vanishingly small. The DMI energy is then mainly contributed by the first term (radial term), which is proportional to wall area (Rw) multiplying the magnetization variation rate along radial direction (1/w), i.e., E DM ∝ R. The isotropic energy is mainly from the skyrmion wall area. Thus, E an ∝ Rw. The Zeeman energy of the skyrmion comes from the skyrmion core proportional to its area of π(R − cw) 2 , where c is a coefficient depending on the magnetization profile and from the wall area proportional to its area of πRw. To obtain the proportional coefficients, one needs to find approximate expressions for f i (R/w) (i = 1, …, 6) in Eq. (8).
In the case of R ) w (or x R=w ) 1), sinh(x) ≈ cosh(x) ≈ e x . Thus, function gðt; xÞ is positive and significantly non-zero only near t = x, reflecting the fact that E ex , E DM , and E an are mainly from the skyrmion wall region that is assumed to be very thin. Furthermore, the area bounded by g(t, x)-curve and t-axis is 1 so that g(t, x) ≈ δ(t − x) resembles the properties of a delta function. We can evaluate f i 's under this approximation (see Methods and Supplementary  Fig. 2). For example, f 1 (x) is The total energy is then Due to the specific form of the magnetization profile of Θ dw (r), Rw-term in E Ze vanishes and E Ze % 4πM s B R 2 2 þ π 2 24 w 2 À Á . The skyrmion size and wall width are then the values that make E(R, w) minimal, or For B = 0, Eqs. (11) and (12) can be analytically solved. The results are The dashed lines in Fig. 2a-c are the approximate formulas that compare quite well with simulation results too. For B ≠ 0, we do not have a closed-form solution of Eqs. (11) and (12) because they are a sextic equation for R or w (see Supplementary Note 2), but their numerical solutions are easily obtained that are plotted as dashed lines in Fig. 2d. In summary, our approximate formula agrees very well with the simulations for R ) w as expected from our approximation. For smaller skyrmions, the approximation is still not bad, and qualitatively gives correct parameter dependence. We can also determine the upper limit of D and lower limits of A, K, and B from the approximate formula. Since R must be real and finite, we have for B = 0, so that the upper limit of D and lower limits of A and K are obtained as D< 4 π ffiffiffiffiffiffiffi AK p , A> π 2 D 2 16K , and K> π 2 D 2 16A . The limit is also the critical value separating the uniform state from helical state 34 . These critical values are plotted in Fig. 2a-c as vertical dashed lines that agrees also well with simulations.

Discussion
In our model, we took into account the shape anisotropy, but did not include the demagnetizing field due to the inhomogeneous part of the magnetization distribution. For very large and thin films, where the film thickness d is much smaller than R and w, the neglected demagnetizing field is vanishing small (see Supplementary Note 3 for a detailed proof). This is the case for those inversion-symmetry breaking systems with a single ferromagnetic layer, such as Pd/Fe/Ir 8,30 and Pt/Co 10,20 , but this approximation is not good for thick films with bulk DMI, such as MnSi 19 and FeGe 7 , or those repeated multilayer systems 12 . Since most device proposals are based on thin-film systems, our model has a wide applicability. For thick films, the magnetization is nonuniform in the third dimension, so that the skyrmion is no longer a twodimensional (2D) magnetic texture. Twisted skyrmion tubes or bobbers may be energetically preferred 35 . Our model as well as all the 2D theories at present is not applicable, and more complicated three-dimensional theories are needed for this case. For the exchange interaction, we only consider a nearest-neighbor exchange interaction characterized by a single parameter A. For those systems with significant exchange interactions beyond nearest neighbor 36,37 , our model may not be directly applicable. However, if the extra exchange interaction (for instance, a nextnearest neighbor interaction) is small compared to A, its effect can be approximately considered by modifying A to an effective value A eff 38 . Our model is a classical one, in which quantum fluctuations and finite-size effects are not considered. However, the comparison with ref. 30 (see Supplementary Fig. 3) shows that even for skyrmions as small as 1 nm, our model still gives good results. The quantum fluctuations induce zero-point energy that helps stablizing the skyrmions, and may be included in a phenomenological way by replacing the model parameters by effective ones 39 . It is also worth mentioning that we consider the isolated skyrmions in infinite media. The edge confinement effect and skyrmion lattice are not the concern of this study.
We compare our theoretical results of the skyrmion size with the experimental results for Pd/Fe/Ir 8,30,40  It is natural to extend our approach to Bloch skyrmions in the systems with bulk inversion symmetry broken. The bulk DMI energy ÞdS can be rewritten as where γ = ± π/2 gives minimal energy. Since all other discussions are the same as those for Néel skyrmions, the formulas of R and w are applicable to the Bloch skyrmions. There are already a lot of studies on the size of skyrmions. In comparison with previous studies 18,25-29 , our results show not only the correct behavior of the skyrmion size as a function of D, A, and K, but also the correct limit at D = 0. The correct result can be obtained only by treating R and w as independent variables. In Ref. 18 , the skyrmion size in a skyrmion lattice with a spiral structure was found to be determined by D/A to some extent. This result is not true for an isolated skyrmion because the isolated skyrmion size depends sensitively on the anisotropy K 20,43 and the external magnetic field B 8,30 . Another well-cited formula is R ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi  14)), and same critical behavior near the critical condition. However, the behavior at large A is incorrect: micromagnetic simulations show R always decreases with A, but according to this formula R increases with A at large A. Ref. 26 gives another formula R ¼ Dπ 2 2Kπþ8M s B=π that does not depend on A. However, this behavior is inconsistent with simulations. There are also many other papers for the skyrmion size based on different ansatz about the skyrmion profile 27,28,44 , but in those studies the skyrmion energy is assumed to be a function of the skyrmion size R only, while the skyrmion wall width w is either ignored or treated as a constant.
In contrast, we treated R and w as independent variables and obtained more satisfactory results. We compare micromagnetic simulation data against our approximate formula of the skyrmion size R (Eq. (13) or the numerical solution of Eqs. (11) and (12) (for B ≠ 0), as well as the analytic formulars in refs. 25,26 in Fig. 2e-h to demonstrate the advance of our results. During the reviewing process of our paper, there are several new studies on the same problem. In ref. 45 , similar approximate skyrmion energy in terms of R and w (Eq. (10)) was also obtained although two 180 domain walls were used instead of a 360°d omain wall, the external magnetic field was not included, and no detailed derivation was given. In ref. 46 , similar approach was used and the magnetostatic energy due to both surface and volume magnetic charges was included. A bistability region in parameter space was found, in which two kinds of skyrmions with different sizes are predicted. The one with a larger size is stabilized mainly by demagnetizing field and the other is stabilized mainly by DMI (which is also discussed in ref. 47 ). In our model, because of the ultra-thin film is assumed, the demagnetizing field is dominated by the the shape anisotropy of uniform magnetization, and the field due to the inhomogeneity of magnetization is vanishingly small. Thus, there is no bistable skyrmion structures that is also consistent with ref. 46 and all our skyrmions are stabilized by DMI.
In conclusion, a single skyrmion can be well described by a 360°domain wall profile parametrized by two fundamental quantities, skyrmions size and wall width. Through the minimization of total energy with respect to the skyrmion size and wall width, we obtain analytic formulas for skyrmion size and wall width as a function of exchange stiffness, anisotropy coefficient, DMI strength, and external field. The formulas agree very well with simulations and experiments.

Methods
Numerical simulations. The numerical simulations are carried out by Mumax package. We consider a magnetic disk of diameter from 128 nm to 1024 nm and thickness 0.4 nm (see Supplementary Fig. 1). The mesh size of 1 nm × 1 nm × 0.4 nm is used in our simulations for those skyrmions larger that 5 nm in radius. For skyrmions of R < 5 nm, the mesh size is 0.1 nm × 0.1 nm × 0.4 nm. The skyrmion size R is obtained by finding the radius of the m z = 0 contour line. The skyrmion wall with is obtained by fitting the simulated radial profile by the 360°domain wall profile with R obtained before.
Derivation of energy expressions. To derive the functions f i (x) in the energy expression, we substitute Eq. (7) into Eqs. (3)- (6). For the exchange energy, by defining x = R/w, t = r/w, we have Thus, we define where the coefficient I 1 is determined by 2e 2x e 2x þ 1 ð Þ 2 dx ¼ 1: where I 3 is determined by The integrand in f 4 is 0 at r = 0, r = ∞ and r = R. Furthermore, it has opposite signs for r < R and r > R. So f 4~0 after the integration.
For the anisotropy energy, we have Similar to the exchange energy and DM energy, the anisotropy energy is only non-zero near r = R, too. The approximate form of function f 5 is the same as f 1 , The Zeeman energy is non-zero for both the wall region and the skyrmion core. The function f 6 is where Li s (x) is the polylogorithm function defined by Li s x ð Þ ¼ P 1 n¼1 x n n s . The asymptotic form of À 1 4 Li 2 Àe 2x ð Þis So we have f 6 x ð Þ % x 2 2 þ π 2 24 . The functions f 1~f6 as well as their approximate forms are plotted in Supplementary Fig. 2.
Data availability. The data that support the plots within this paper and other findings of this study are available from the corresponding author on reasonable request.