Chiral magnetic skyrmions across length scales

The profile and energy of chiral skyrmions, found in magnetic materials with the Dzyaloshinskii-Moriya interaction, can be approximated by formulae obtained through asymptotic analysis in the limits of small and large skyrmion radius. Using fitting techniques, we modify these formulae so that their validity extends to almost the entire range of skyrmion radii. Such formulae are obtained for skyrmions that are stabilised in the presence of an external field or easy-axis anisotropy or a combination of these. We further study the effect of the magnetostatic field on the skyrmion profile. We compare the profile of magnetic bubbles, stabilized without the chiral Dzyaloshinskii-Moriya interaction to that of a chiral skyrmion.


INTRODUCTION
Chiral magnetic skyrmions are topological magnetic configurations that are stabilised in materials with the Dzyaloshinskii-Moriya interaction [1]. They have been observed in a number of experiments and the detailed features of individual skyrmions have been resolved experimentally to an impressive degree for isolated skyrmions [2][3][4][5][6][7] and in a skyrmion lattice [8]. The experimental works have mapped the profile of the skyrmion, i.e., the magnetization as a function of the distance from its center. The skyrmion profile is the foundation for the study of the statics and of dynamical behaviors of skyrmions and the subsequent derivation of quantitative results. The profile enters in formulae for dynamical phenomena, for example, skyrmion translation and oscillation modes [9,10] or antiferromagnetic skyrmion excitations [11], and it is crucial for quantitative calculations. Furthermore, it quantifies the localization of the entity which in turn establishes its particle-like nature.
Approximation formulae that capture basic features of the skyrmion profile have been proposed. For skyrmions of large radius, an ad-hoc ansatz based on explicit onedimensional domain wall profiles [12] has been suggested and is widely used to examine structural and dynamic properties [2,[13][14][15]. The domain wall profile with an adjustable wall thickness [9,16] and a related ansatz [4] have been employed as more flexible forms in order to capture the skyrmion features over a wider parameter range. For the case of a model with easy-axis anisotropy, accurate skyrmion profiles in the asymptotic sense have been obtained for skyrmions of small [17] and of large radius [18]. Similar methods were used for a model including the magnetostatic field [19,20]. In the present paper, we give a comprehensive analysis of the skyrmion profile for the case of easy-axis anisotropy and external magnetic field. We use the methods that were developed for the case of easy-axis anisotropy only in Refs. [17,18]. While the work in the above references focuses on asymptotic formulae for skyrmions of small and large radius, the present paper gives a set of formulae for determining the skyrmion radius that are good approximations in the whole parameter range, i.e., also for intermediate radii.
The approximation formulae are derived so as to be consistent with the asymptotic ones. We further present the result of a set of numerical simulations for the effect of the magnetostatic field on the axially-symmetric skyrmion. We make a direct comparison between the profiles of chiral skyrmions and bubbles which reveals that, beyond basic differences, the two entities share many features.
We, finally, make a detailed study of the skyrmion energy, made possible by the availability of a detailed analytical description of the skyrmion profile. The analytical calculations and formulae will help focus the efforts for applications of skyrmions.
The paper is arranged as follows. Sec. 2 formulates the problem. Sec. 3 gives the far field for the skyrmion profile. Sec. 4 gives the formulae for the skyrmion profile and radius for the case of an external field. Sec. 5 gives the corresponding formulae for the case of easy-axis anisotropy. Sec. 6 studies the effect of the magnetostatic field on the skyrmion and gives comparisons with magnetic bubbles. Sec. 8 contains our concluding remarks. In the Appendix, we give a complete review of the method of asymptotic matching and we expand it so that it can readily be applied for the derivation of formulae in the main text.

FORMULATION
We consider a thin film of a ferromagnetic material on the xy-plane with symmetric exchange, Dzyaloshinskii-Moriya (DM) interaction, anisotropy of the easy-axis type perpendicular to the film, and an external field. The micromagnetic structure is described via the magnetization vector m = m(x, y) with a fixed length normalized to unity, m 2 = 1. The energy is where A, D, K are the exchange, DM, and anisotropy parameters, respectively, and the external field h = H/M s is normalized to the saturation magnetization M s . The DM energy density e DM may be chosen to have the bulk form e DM =ê µ · (∂ µ m × m) or the interfacial form e DM = µνêµ · (∂ ν m × m), where µ, ν = 1, 2 and the summation convention is invoked, whileê µ are the unit vectors for the magnetization in the respective directions. It is useful to consider the length scales that arise naturally, (2) ex is the exchange length, w is the domain wall width, D gives the pitch of the spiral solution for strong DM interaction, and S is related to the radius of small skyrmions. Using ex as the unit of length, and assuming an external field h = hê 3 perpendicular to the film, the scaled form of the energy is where we have assumed that m =ê 3 in the uniform magnetization of the film and we have used the scaled DM and anisotropy parameters Static magnetization configurations satisfy the timeindependent Landau-Lifshitz equation The bulk DM form is h DM =ê µ ×∂ µ m and the interfacial one is h DM = µνêµ × ∂ ν m. Skyrmion solutions with axial symmetry can be described in polar coordinates (r, φ) using the spherical angles (Θ, Φ) for the magnetization. We choose Φ = φ for the case of interfacial DM interaction and Φ = φ + π/2 for the case of bulk DM interaction. In both cases, the angle Θ = Θ(r) satisfies Existence, uniqueness, stability and minimality of axially symmetric skyrmions in the regime λ 2 /κ 1 and λ 2 /h 1, respectively, has been examined rigorously in [21,22]. In the following, we will use the convention that Θ(r = 0) = π in the skyrmion center and Θ = 0 at spatial infinity r → ∞. The skyrmion radius R will be defined to be at the radial distance where Θ = π/2.

FAR-FIELD
In the far field, where Θ is small, Eq. (6) reduces by linearization to Under a scaling transformation √ κ + h r → r this gives the modified Bessel equation. The appropriate solution, decaying for r → ∞, is the modified Bessel function of the second kind K 1 (r) [23], that also appeared in a similar context in Ref. [24]. We have where C is an arbitrary constant, γ ≈ 0.57721 is the Euler-Mascheroni constant, and ξ n = 1 2 (ψ(n) + ψ(n + 1)) − ln n, The asymptotic form of K 1 (r) for large values of r gives [23] Θ ∼ C π 2r e −r , r 1.
It is interesting to note that the result for the far field does not depend on the DM interaction. The exponential decay (10) is due to the anisotropy and the external field.

EXTERNAL FIELD
We consider that we only have an external field and no anisotropy, κ = 0. Then, we may choose ex M s /H as the unit of length which leads to setting h → 1 in Eq. (6) and the DM parameter is We initially focus on skyrmions of small radius. In this limit, it is convenient to work with the variable u(r) = tan Θ(r) 2 (12) which is the modulus of the stereographic projection of the magnetization vector. In the case, of exchange interaction only, we have the well-known axially-symmetric Belavin-Polyakov (BP) solution [25] u(r) = R r , that represents a skyrmion of unit degree with a radius R. Inspired by the BP solution (13), we remove the singularity at the origin by defining the field v(r) = ru(r).
For comparison purposes, we note that the variable v in Ref. [17] is related to the presentṽ byṽ = 2v.
The equation forṽ is given in Appendix A in Eq. (A1). In Appendix A 1, we consider the case of small values of λ and find the near field at the center of the skyrmion, in Eqs. (A6), and well beyond its radius, in Eq. (A8). This generalises the results of Ref. [17] for the case where both an external field and an anisotropy are present. The analysis is based on the fact thatṽ(r) is small in the entire space, when R is small. In Appendix A 2, by matching the near field of Eq. (A8) with the far field in Eq. (8) in a region where both are valid, a relation between the system parameters and the skyrmion radius is obtained in Eq. (A13).
For h = 1 and κ = 0, Eq. (A13) reduces to We solve Eq. (A1) (or Eq. (6)) numerically using a shooting method and we find the skyrmion profiles for various values of the parameter λ. Given the profile, we detect the skyrmion radius. In Fig. 1 we show the skyrmion radius found numerically as a function of the parameter. Formula (16) is given in entry (b) of the figure and it fits excellently the numerically found skyrmion radius for small R.
Isolated skyrmions exist in this model for λ ≤ λ 0 ≈ 1.12 while at λ = λ 0 a transition to a skyrmion lattice occurs [26]. Inverting relation (16) gives the leading order of approximation Following Ref. [27], a modification of the denominator of Eq. (17), leads to a very good fit of the numerical data over a wide range as shown in Fig with Eq. (16) at λ → 0, but the factor 1.8 is found only by fitting the numerical data. For larger λ, the linear relation fits the numerical data. The combination of formulae (18) for λ 0.6 and (19) for λ 0.6 give a very good fit of the skyrmion radius for the entire range of the parameter

space.
We now turn to the skyrmion profiles. Eq. (A6) gives the profile at the skyrmion inner core and Eq. (8) gives the profile at large distances. Based on the above asymptotic results, we will give simple approximation formulae for the skyrmion profile.
The field around the skyrmion center is found by a Taylor expansion of the solution in Eq. (A6) close to r = 0 (applied for h = 1, κ = 0) while the parameter λ is substituted from Eq. (15). The calculations are given in Appendix A 3 and lead to Eq. (A15). For the present case, this reduces to the approximation and it is valid for r < R for small radius skyrmions.
In the far field, Θ and the field u in Eq. (12) are proportional to each other and both satisfy the modified Bessel equation (7). We will continue the discussion using the field u, instead of Θ, because it eventually gives a good fit for the skyrmion profile over a wider interval in r. For small R, we have u = RK 1 (r), where the factor R follows from the asymptotic matching conditions. We write Fig. 2 shows the profile for a skyrmion of small radius. Surprisingly, the far field formula (21), shown by the blue dashed line, gives a good approximation almost in the entire space. Close to r = 0, Eq. (20) is the correct approximation. Eq. (21) cannot be justified in this region.
We conclude the section by giving a description of the complete skyrmion profile over the full spatial range. A validity condition for the near field is λr 2 ṽ 0 [17]. Using Eq. (16) andṽ 0 ≈ R, we have that Eq. (20) is valid for The form of Eq. (20) indicates that in the regime where condition (22) is valid the skyrmion profile is close to the BP profile. The exponentially decaying behaviour (10) is valid for r 1. Between the regimes of validity of the BP profile and the exponentially decaying profile, the modified Bessel function (21) is still a good approximation.

EASY-AXIS ANISOTROPY
We consider the case of easy-axis anisotropy and no external field, h = 0. Then, we may choose w as the unit of length which leads to setting κ → 1 in Eq. (6) and This section is based on the work presented in Refs. [17,18].

A. Small radius
For skyrmions of small radius, it is convenient to work in the variables u(r) of Eq. (12) or This asymptotic formula was derived in Ref. [17]. It is shown in Fig. 3a that it matches excellently the numerically found skyrmion radius for small R in the range 0 < λ 0.1.
Inverting formula (25) gives, to the lowest order of approximation, Eq. (17). A modification of the denominator of Eq. (17), as obtained in Ref. [27], leads to a very good fit of the numerical data over the very wide range 0 < λ 0.6, as shown in Fig. 3a. Using asymptotic results from Appendix A, we can give approximation formulae for the skyrmion profile. The field close to the skyrmion center is given in Eq. (A15). In the present case (κ = 1, h = 0), this reduces to a formula identical to Eq. (20) in Sec. 4. The far field is given in Sec. 3 and it is identical to Eq. (21) in Sec. 4. Fig. 4 shows formulae (20) and (21) compared to numerically calculated skyrmion profiles for two values of the parameter.

B. Large radius
At the parameter value λ = 2/π, the skyrmion radius diverges to infinity and a phase transition occurs from the uniform to the spiral state. When the skyrmion is large, an asymptotic series for the profile is given in negative where Θ 0 , Θ 1 , Θ 2 , Θ 3 , · · · are functions of r − R. The leading order contribution satisfies This is the standard one-dimensional wall profile. We now give approximation leading order formulae for the core and the far field of the skyrmion. The field at the core is [18] where I 1 (r) is the modified Bessel function of the first kind. For the far field, we have that simplifies to in the region far from the domain wall. The profile at the skyrmion domain wall region (28) is matched to the skyrmion core and to the far field to all asymptotic orders O(R −n ) in Ref. [18]. Fig. 5 shows a numerically calculated skyrmion profile and this is compared to the leading order formulae (28), (29), and (30).
The relation between the parameter λ and the radius R is derived during the asymptotic process. It is given as the series [17], where Retaining terms up to O(1/R 2 ) and solving for R, we obtain This is shown in Fig. 3 by the green line. The fitting formula (26) and the asymptotic result (34) for large radius give a very good approximation of the numerically found radii for the entire parameter space.

MAGNETOSTATIC FIELD
The magnetostatic field is known to be crucial for stabilising magnetic bubbles [28] which are cylindrical domains that share topological features with skyrmions. The relation between skyrmions and bubbles has been already explored [15,29]. Given the role of the magnetostatic field in stabilizing magnetic bubbles, we expect that it will often (but not always) have a stabilizing effect also on chiral skyrmions. When the DM interaction is present, the magnetostatic field is usually considered to be of secondary importance for skyrmion generation and stability. In this section, we present the results of numerical simulations for the profile of skyrmions when the magnetostatic field is included. We also give a comparison between bubble and skyrmion profiles.
We first summarize the results on the stability of cylindrical domains [30] that will be useful in the interpretation of the numerical results in this section. Easy-axis anisotropy perpendicular to the film is necessary in order to have stable magnetic bubbles in a material. The magnetostatic field is favoring the expansion of the bubble and it thus has a demagnetization effect. The bubble domain wall tends to shrink in order to minimize the bubble size, however, the demagnetizing effect is typically stronger and thus no energy balance can be achieved. The stabilization of bubbles is actually obtained by the addition of an external (bias) field perpendicular to the film.
Magnetic bubbles, even ones with high skyrmion numbers, have been observed and studied [28]. Topologically trivial magnetic bubbles, with skyrmion number zero, have also been observed. The bubble domain wall is primarily of Bloch type as this is favoured by the magnetostatic interaction. Both Bloch chiralities are energetically equivalent. The structure of the bubble wall is though complicated by the effects of the film boundaries [31,32].
The static Landau-Lifshitz equation has the form where the bulk DM interaction term has been chosen, h m is the magnetostatic field normalized to the saturation magnetization, and w is used as the unit of length. The parameters are given here again for convenience We often consider the substitution K → K − 1 2 µ 0 M 2 s in order to take into account the fact that the magnetostatic field is equivalent to easy plane anisotropy for ultra-thin films. This substitution is equivalent to seting κ → κ − 1 and λ → κ/(κ−1)λ in Eq. (35). After these substitutions, the field h m would represent only the deviation of the magnetostatic field from its approximation as an easyplane anisotropy term.
We find solutions of Eq. (35) using an energy relaxation algorithm [32]. Our method works in cylindrical coordinates (r, z) and it is confined to axially symmetric configurations only. The magnetostatic field is calculated using a conjugate gradient method. We consider an infinite film achieved by employing a long enough numerical mesh in the r direction and assuming no magnetic charges on the lateral boundary. We place the film of thickness d f so that its central plane is at z = 0, i.e., the film extends over −d f /2 ≤ z ≤ d f /2. We find that solutions for the magnetization m = (m r , m φ , m z ) satisfy the parity relations m r (r, −z) = −m r (r, z), m φ (r, −z) = m φ (r, z), m z (r, −z) = m z (r, z).
We first explore the case of easy axis anisotropy only and no bias field. We consider a thin film with thickness d f = 3.5 w and parameter values λ = 0.33, κ = 4. We discretize space with lattice spacings ∆r = 0.1 and ∆z = 0.5 in the r and z directions respectively. The numerical mesh extends to r = 20. Fig. 6 shows the cylindrical components of the magnetization for a skyrmion solution as functions of r at various levels of z. The skyrmion radius is R ≈ 1.5, a significant increase in comparison to the skyrmion radius R = 0.42 found for λ = 0.33 when the magnetostatic field is neglected. Furthermore, the profile acquires a nonzero m r component, which is an odd function in the direction perpendicular to the film.
When the parameter λ is increased, the skyrmion radius increases. In the presence of the magnetostatic field, we find no static skyrmions for λ > 0. 35. This should be compared to the critical value λ = 2/π when the magnetostatic field is neglected. This finding suggests that skyrmions with a large radius cannot be obtained (without a bias field) because they are destabilised by the magnetostatic field.
It is interesting to compare the chiral skyrmion profiles with magnetic bubble ones that are stabilised without the DM interaction. We choose a film with thickness d f = 11 w , extending over −d f /2 ≤ z ≤ d f /2, and the parameter values κ = 2, h = 0.2 (and set λ = 0). The lattice spacings are ∆r = 0.2, ∆z = 1.0 and the numerical mesh extends to r = 40. The energy minimization algorithm converges to a static bubble with radius R ≈ 8. Fig. 7a shows the three components of the magnetization as functions of r at various levels of z. The stray field at the bubble domain wall induces large values for m r close to the film boundaries. Correspondingly, the component m φ depends strongly on z. As a result, the bubble wall is of a hybrid character; Bloch at the center and tilted towards Neél type closer to the boundaries. The hybrid wall structure has been shown to persist when DM interaction is included [33,34]. The component m z also has a dependence on z that results in the bubble domain wall width to depend on z and to the bulging shape of the bubble [31,32]. Fig. 7b shows the spherical angle Θ for the bubble by blue solid lines at various z-levels. For comparison, a red dashed line shows the chiral skyrmion profile that is obtained for a DM parameter λ = 0.6316 and no magnetostatic field. For that value of λ, the chiral skyrmion has a similar radius as the bubble in the figure. We finally note that our comparison between a chiral skyrmion and a bubble is only possible for large radii as magnetic bubbles are not stabilised for small radii [30].

ENERGY OF A SKYRMION
The skyrmion radius, studied in the previous sections, is a length scale that determines the balance of the individual interaction terms and thus provides information about the skyrmion energy. The energetics of the model is important not only as a means to understand the skyrmions as energy minima but also in order to build arguments concerning their stability and excitations. In the regime of small radii, the energy asymptotics has been examined rigorously in [19,22,35].
We consider the case of easy-axis anisotropy and no external field so that the energy is given in Eq. (37)

A. Small radius
For small skyrmion radius, asymptotic analysis gives [17] while E a can be obtained from Eq. (37). Result (38) can be used to produce asymptotic series for the energy terms as we will now demonstrate. We consider a specific DM parameter λ 0 and the energy E = E λ0 (Θ) for any configuration Θ. This is written as where we have set so that all terms E ex ,Ẽ DM , E a are integrals over Θ that do not contain parameters. We confine ourselves to the specific configurations Θ = Θ λ that are the profiles corresponding to static skyrmions for any parameter value λ (that may be different from λ 0 ). The virial relation (37) holds, and it is used in Eq. (39) to obtain for the energy of profile Θ λ in a system with parameter value λ 0 . We take the derivative of the energy expression (42) with respect to λ, and require that E λ (λ = λ 0 ) = 0. This gives the condition An explicit calculation of the terms of Eq. (44) utilizing Eq. (38), shows that terms of order O(λ 2 /(ln λ) m ) with m ≥ 2 should necessarily be included in the energy asymptotic series for the condition (44) to be satisfied. Specifically, where a m , b m are constants (with b 1 = 2 from Eq. (38)). Substituting these in Eq. (44) obtains In order that Eq. (46) be satisfied to all orders in 1/(ln(λ)) m , it should hold For m = 1, this gives (since a 1 = 0) We thus have More information would be required to obtain explicit values for a m+1 , b m with m ≥ 2. Comparing the second of Eqs. (45) with numerical data, we propose the following formula, Using formulae (49) and (50), the energy is This is in agreement with rigorous results in Ref. [22] where an error term O(λ 2 /(ln λ) 2 ) is given. The numerical data for the energy are fitted well with formulae (49), (50), (51), for parameter values λ 0.25. Fig. 8 shows the numerically calculated skyrmion energy as a function of λ. A modification of formula (51) is found to fit the numerical data very well almost in the entire parameter range, This formula is plotted in Fig. 8. For λ close to the critical value 2/π formula (55) (of the next subsection) should be used instead.
Eq. (51) gives Eq. (53) gives a good approximation only in the range R < 0.5. It is remarkable that the numerical data are fitted well over a large range of R by a quite significant modification of Eq. (53), that is plotted in Fig. 9. We finally note that the techniques used in this section can be readily applied to the case of applied external field.

B. Large radius
For large skyrmion radius, asymptotic analysis gives [18] E ∼ 4π 2 |λ| 1/2 2 π − λ, 2 π − λ 1 (55) or (cf. Eq. (34)) In the range 0.6 λ < 2/π, the numerical data of Fig. 8 are fitted well with formula (55). Fig. 9 shows the skyrmion energy as a function of its radius. For large R, the data are fitted very well with the formula which is a sharpening of Eq. (56) and it is plotted in the figure. We further note that each one of the energy components is fitted very well, for R 3, by the forms

CONCLUDING REMARKS
We have presented fitting formulae of skyrmion profiles and energy that cover almost the entire range of the parameter values and of skyrmion radii. These are based on formulae obtained through asymptotic analysis [17,18]. The analysis is extended to the case of an external magnetic field in the Appendix. We were thus able to study both mechanisms that give stable skyrmions, i.e., the case of an external field and the case of easy-axis anisotropy, and we have given asymptotic formulae when both terms are present in the model.
We have further studied the effect of the magnetostatic field on the skyrmion profile. It tends to increase the skyrmion radius while it appears to destabilise skyrmions of large radius. We compared the profile of magnetic bubbles stabilized without the chiral Dzyaloshinskii-Moriya interaction to that of a chiral skyrmion. We found that the profiles are similar in the two cases, the main difference between the two being the variation of the bubble profile across the film thickness.
Finally, we have produced asymptotic series expressions for all skyrmion energy terms using an argument for the minimization of the energy. This has resulted in numerical expressions that sharpen previous asymptotic results. We also give approximate formulae that cover almost the entire range of the parameter values and of skyrmion radii.
The obtained formulae can be used for comparisons with experimental results and for a variety of investigations where the details of skyrmion features are important, notably, in skyrmion dynamics.