Chiral magnetic skyrmions across length scales

Chiral magnetic skyrmions are topological magnetic configurations that are stabilized in materials with the Dzyaloshinskii–Moriya (DM) interaction [1–3]. 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 [4–7] and in a skyrmion lattice [8–10]. 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 forms 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 [11, 12] or antiferromagnetic skyrmion excitations [13], 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 one-dimensional domain wall profiles [14] has been suggested and is widely used to examine structural and dynamic properties [4, 15–17]. The domain wall profile with an adjustable wall thickness [11, 18] and a related ansatz [6] 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 [19] and of large radius [20]. Similar methods were used for a model including the magnetostatic field [21, 22]. In the present paper, we give a comprehensive analysis of the skyrmion profile for the case of easy-axis anisotropy perpendicular to the film and external magnetic field. We use the methods that were developed for the case of easy-axis anisotropy only in [19, 20]. The above references focus on asymptotic formulae for skyrmions of small (close to zero) and large (tending to infinity) radius. Specifically, the results of the asymptotics for small radii R are valid for $R \lesssim 0.3$ and those for large radii are valid for $R \gtrsim 3$. The present paper gives a set of formulae for determining the skyrmion radius that are good approximations in the whole parameter range, i.e. including the intermediate radii, $R\sim O(1)$. The present approximation formulae are derived so as to be consistent with the asymptotic ones. We further present the results of a set of numerical simulations for the effect of the magnetostatic field on the axially-symmetric skyrmion. This is included for completeness and in order to be able to make a direct comparison between the profiles of chiral skyrmions (obtained in the bulk of the text) and bubbles. Beyond basic differences, the two entities share many features. A comprehensive study of the effect of the non-local magnetostatic field has been given in [17, 23] where skyrmions and bubbles are compared and furthermore bi-stability phenomena have been revealed.

We, finally, give a study of the skyrmion energy, made possible by the availability of a detailed analytical description of the skyrmion profile. We extend previous asymptotic analysis results for skyrmions of small radius and we give fitting formulae covering almost the entire range of radii. The analytical calculations and formulae will help focus the efforts for applications of skyrmions.

The paper is arranged as follows. Section 2 formulates the problem. Section 3 gives the far field for the skyrmion profile. Section 4 gives the formulae for the skyrmion profile and radius for the case of an external field. Section 5 gives the corresponding formulae for the case of easy-axis anisotropy. Section 6 studies the effect of the magnetostatic field on the skyrmion and gives comparisons with magnetic bubbles. Section 7 extends asymptotic analysis results and introduces fitting formulae for the skyrmion energy. Section 8 contains our concluding remarks. The appendix gives a complete review of the method of asymptotic matching and this is expanded to include the case of combined external field and anisotropy so that it can readily be applied for the derivation of formulae in the main text.

We consider a thin film of a ferromagnetic material on the xy-plane with symmetric exchange, 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 $\boldsymbol{m} = \boldsymbol{m}(x,y)$ with a fixed length normalized to unity, $\boldsymbol{m}^2 = 1$. The energy is

$$\begin{align} E = A \int (\partial_\mu\boldsymbol{m}\cdot\partial_\mu\boldsymbol{m}) {\mathrm{d}}^2 x + D \int e_{\textrm{DM}}\,{\mathrm{d}}^2x + K \int (1-m_3^2) {\mathrm{d}}^2x - \mu_0 M_s^2 \int \mathbf{h}\cdot\boldsymbol{m}\,{\mathrm{d}}^2x \end{align} \tag{ 1 }$$

where $A, D, K$ are the exchange, DM, and anisotropy parameters, respectively, and the external field $\mathbf{h} = \mathbf{H}/M_{\mathrm{s}}$ is normalized to the saturation magnetization $M_{\mathrm{s}}$. The DM energy density $e_\textrm{DM}$ may be chosen to have the bulk form $e_\textrm{DM } = \boldsymbol{\hat{e}}_\mu\cdot(\partial_\mu\boldsymbol{m}\times\boldsymbol{m})$ or the interfacial form $e_\textrm{DM} = \epsilon_{\mu\nu} \boldsymbol{\hat{e}}_\mu\cdot(\partial_\nu\boldsymbol{m}\times\boldsymbol{m})$, where $\mu,\nu = 1,2$ and the summation convention is invoked, while $\boldsymbol{\hat{e}}_\mu$ are the unit vectors for the magnetization in the respective directions.

It is useful to consider the length scales that arise naturally,

$$\begin{equation} \ell_{\textrm{ex}} = \sqrt{\frac{2A}{\mu_0 M\,_{\mathrm{s}}^2}},\quad \ell_{\textrm{w}} = \sqrt{\frac{A}{K}},\quad \ell_{\textrm{D}} = \frac{2A}{D},\quad \ell_{\mathrm{S}} = \frac{D}{2K}. \end{equation} \tag{ 2 }$$

$\ell_\textrm{ex}$ is the exchange length, $\ell_\textrm{w}$ is the domain wall width, $\ell_\textrm{D}$ gives the pitch of the spiral solution for strong DM interaction, and $\ell_{S}$ is related to the radius of small skyrmions. Using $\ell_\textrm{ex}$ as the unit of length, and assuming an external field $\mathbf{h} = h\boldsymbol{\hat{e}}_3$ perpendicular to the film, the scaled form of the energy is

$$\begin{align} E = \frac{1}{2} \int (\partial_\mu\boldsymbol{m}\cdot\partial_\mu\boldsymbol{m}) {\mathrm{d}}^2 x + \lambda \int e_{\textrm{DM}}\,{\mathrm{d}}^2x + \frac{\kappa}{2} \int (1-m_3^2) {\mathrm{d}}^2x + h \int (1-m_3)\,{\mathrm{d}}^2x \end{align} \tag{ 3 }$$

where we have assumed that $\boldsymbol{m} = \boldsymbol{\hat{e}}_3$ in the uniform magnetization of the film and we have used the scaled DM and anisotropy parameters

$$\begin{equation} \lambda = \frac{\ell_{\textrm{ex}}}{\ell_{\textrm{D}}} = \frac{D}{\sqrt{2A\mu_0 M\,_{\mathrm{s}}^2}},\qquad \kappa = \frac{\ell_{\textrm{ex}}^2}{\ell_{\textrm{w}}^2} = \frac{2K}{\mu_0 M\,_{\mathrm{s}}^2}. \end{equation} \tag{ 4 }$$

Static magnetization configurations satisfy the time-independent Landau–Lifshitz equation

$$\begin{equation} \boldsymbol{m} \times \left( \partial_\mu\partial_\mu\boldsymbol{m} - 2\lambda\, \mathbf{h}_{\textrm{DM}} + \kappa\,m_3 \boldsymbol{\hat{e}}_3 + h\boldsymbol{\hat{e}}_3 \right) = 0. \end{equation} \tag{ 5 }$$

The bulk DM form is $\mathbf{h}_\textrm{DM} = \boldsymbol{\hat{e}}_\mu\times\partial_\mu\boldsymbol{m}$ and the interfacial one is $\mathbf{h}_\textrm{DM} = \epsilon_{\mu\nu} \boldsymbol{\hat{e}}_\mu\times\partial_\nu\boldsymbol{m}$.

Skyrmion solutions with axial symmetry can be described in polar coordinates $(r,\phi)$ using the spherical angles ($\Theta, \Phi)$ for the magnetization. We choose $\Phi = \phi$ for the case of interfacial DM interaction and $\Phi = \phi+\pi/2$ for the case of bulk DM interaction. In both cases, the angle $\Theta = \Theta(r)$ satisfies

$$\begin{equation} \Theta^{^{\prime\prime}} + \frac{\Theta^{^{\prime}}}{r} - \frac{\sin(2\Theta)}{2 r^2} + 2\lambda \frac{\sin^2\Theta}{r}- \frac{\kappa}{2} \sin(2\Theta) - h \sin\Theta = 0. \end{equation} \tag{ 6 }$$

Existence, uniqueness, stability and minimality of axially symmetric skyrmions in the regime $\lambda^2/ \kappa\ll 1$ and $\lambda^2/h \ll 1$, respectively, has been examined rigorously in [24, 25]. In the following, we will use the convention that $\Theta(r = 0) = \pi$ in the skyrmion center and $\Theta = 0$ at spatial infinity $r\to\infty$. The skyrmion radius R will be defined to be at the radial distance where $\Theta = \pi/2$.

In the far field, where Θ is small, equation (6) reduces by linearization to

$$\begin{equation} r\,^2\Theta^{^{\prime\prime}}+r\Theta^{^{\prime}} - [1 + (\kappa+h) r\,^2]\Theta = 0. \end{equation} \tag{ 7 }$$

Under a scaling transformation $\sqrt{\kappa+h}\,r \to r$ this gives the modified Bessel equation. The appropriate solution, decaying for $r\to\infty$, is the modified Bessel function of the second kind $K_1(r)$ [26], that also appeared in a similar context in [27]. We have

$$\begin{equation} \begin{split} \Theta(r) = C\,K_1(r) = \frac{C}{r} \left[ 1 + (2\gamma-1) \left( \frac{r}{2} \right)^2 \right. \left. + 2 \sum_{n=1}^\infty \frac{(\ln\left( \frac{r}{2} \right) - \ln n) \left( \frac{r}{2} \right)^{2n}}{(n-1)!\, n!} - 2 \sum_{n=2}^\infty \frac{\xi_n \left( \frac{r}{2} \right)^{2n}}{(n-1)!\, n!} \right] \end{split} \end{equation} \tag{ 8 }$$

where C is an arbitrary constant, γ ≈ 0.57721 is the Euler–Mascheroni constant, and

$$\begin{equation} \begin{split} \xi_n & = \frac{1}{2} (\psi(n)+\psi(n+1))-\ln n, \\ \psi(n) & = -\gamma+\sum_{k=1}^{n-1}\frac{1}{k}, \quad n=2,\ldots. \end{split} \end{equation} \tag{ 9 }$$

The asymptotic form of $K_1(r)$ for large values of r gives [26]

$$\begin{equation} \Theta \sim C \sqrt{\frac{\pi}{2r}}\;e^{-r},\quad r\gg 1. \end{equation} \tag{ 10 }$$

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.

We consider that we only have an external field and no anisotropy, κ = 0. Then, we may choose $\ell_\textrm{ex}\sqrt{M_{\mathrm{s}}/H}$ as the unit of length which leads to setting h → 1 in equation (6) and the DM parameter is

$$\begin{equation} \lambda = \frac{\ell_{\textrm{ex}}}{\ell_{\textrm{D}}} \sqrt{\frac{M_{\mathrm{s}}}{H}} = \frac{D}{\sqrt{2A\mu_0 M_{\mathrm{s}}H}}. \end{equation} \tag{ 11 }$$

We initially focus on skyrmions of small radius. In this limit, it is convenient to work with the variable

$$\begin{equation} u(r) = \tan\frac{\Theta(r)}{2} \end{equation} \tag{ 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 [28]

$$\begin{equation} u(r) = \frac{R}{r}, \end{equation} \tag{ 13 }$$

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

$$\begin{equation} \tilde{v}(r) = r u(r). \end{equation} \tag{ 14 }$$

For comparison purposes, we note that the variable $v$ in [19] is related to the present $\tilde{v}$ by $\tilde{v} = 2v$.

The equation for $\tilde{v}$ is given in appendix in equation (A1). In section 'Near field', we consider the case of small values of λ and find the near field at the center of the skyrmion, in equation (A6), and well beyond its radius, in equation (A8). This generalises the results of [19] for the case where both an external field and an anisotropy are present. The analysis is based on the fact that $\tilde{v}(r)$ is small in the entire space, when R is small. In section 'Asymptotic matching', by matching the near field of equation (A8) with the far field in equation (8) in a region where both are valid, a relation between the system parameters and the skyrmion radius is obtained in equation (A13).

For h = 1 and κ = 0, equation (A13) reduces to

$$\begin{equation} \lambda = -R \left[ \gamma + \frac{1}{2} + \ln(R/2) \right] \end{equation} \tag{ 15 }$$

or

$$\begin{equation} \lambda = -R\, \ln\left( \frac{R}{\alpha_h} \right), \qquad \alpha_h = 2 {\mathrm{e}}^{-(\gamma+\frac{1}{2})} \approx 0.6811. \end{equation} \tag{ 16 }$$

We solve equation (A1) (or equation (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 figure 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.

Zoom In Zoom Out Reset image size

Isolated skyrmions exist in this model for $\lambda \leqslant \lambda_0 \approx 1.12$ while at $\lambda = \lambda_0$ a transition to a skyrmion lattice occurs [29]. Inverting relation (16) gives the leading order of approximation

$$\begin{equation} R = -\frac{\lambda}{\ln\lambda}. \end{equation} \tag{ 17 }$$

Following [30], a modification of the denominator of equation (17),

$$\begin{equation} R = - \frac{\lambda}{\ln(\lambda/1.8)} \end{equation} \tag{ 18 }$$

leads to a very good fit of the numerical data over a wide range as shown in figure 1. Equation (18) is consistent with equation (16) at λ → 0, but the factor 1.8 is found only by fitting the numerical data. For larger λ, the linear relation

$$\begin{equation} R = 1.6(\lambda-0.29) \end{equation} \tag{ 19 }$$

fits the numerical data. The combination of formulae (18) for $\lambda \lesssim 0.6$ and (19) for $\lambda \gtrsim 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. Equation (A6) gives the profile at the skyrmion inner core and equation (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 equation (A6) close to r = 0 (applied for $h = 1, \kappa = 0$) while the parameter λ is substituted from equation (15). The calculations are given in section 'Field at the skyrmion center' and lead to equation (A15). For the present case, this reduces to the approximation

$$\begin{equation} \tan\left(\frac{\Theta}{2}\right) = \frac{R}{r} \left[ 1 + R^2\,\left( \gamma + \frac{3}{4} + \ln(R/2) \right)\, \frac{r^2}{2} \right] \end{equation} \tag{ 20 }$$

and it is valid for r < R for small radius skyrmions.

In the far field, Θ and the field u in equation (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 = R K_1(r)$, where the factor R follows from the asymptotic matching conditions. We write

$$\begin{equation} \tan\left(\frac{\Theta}{2}\right) = R\,K_1(r),\quad r \gg R. \end{equation} \tag{ 21 }$$

Figure 2 shows the profile for a skyrmion of small radius, for λ = 0.4. 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, equation (20) is the correct approximation while equation (21) cannot be justified in this region. Note that the approximation becomes progressively less successful when applied for larger λ.

Zoom In Zoom Out Reset image size

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 $\lambda r^2 \ll \tilde{v}_0$ [19]. Using equation (16) and $\tilde{v}_0\approx R$, we have that equation (20) is valid for

$$\begin{equation} r\,^2\ll -\frac{1}{\ln R}. \end{equation} \tag{ 22 }$$

The form of equation (20) indicates that, in the regime where condition (22) is valid, the skyrmion profile is close to the BP profile. The exponentially decaying behavior (10) is valid for $r\gg 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.

The results of this section are relevant for materials presenting a small anisotropy, e.g. as a result of compensation between perpendicular magnetocrystalline anisotropy and the magnetostatic field.

We consider the case of easy-axis anisotropy and no external field, h = 0. Then, we may choose $\ell_\textrm{w}$ as the unit of length which leads to setting κ → 1 in equation (6) and

$$\begin{equation} \lambda = \frac{\ell_{\mathrm{S}}}{\ell_{\textrm{w}}} = \frac{D}{2\sqrt{AK}}. \end{equation} \tag{ 23 }$$

This section is based on the work presented in [19, 20].

5.1. Small radius

For skyrmions of small radius, it is convenient to work in the variables u(r) of equation (12) and $\tilde{v}(r)$ of equation (14) (see appendix). For κ = 1 and h = 0, equation (A13) reduces to

$$\begin{equation} \lambda = -R\, [\gamma + 1 + \ln(R/2)] \end{equation} \tag{ 24 }$$

or

$$\begin{equation} \lambda = -R\, \ln\left( \frac{R}{\alpha_\kappa} \right), \qquad \alpha_\kappa = 2 {\mathrm{e}}^{-(\gamma+1)} \approx 0.4131. \end{equation} \tag{ 25 }$$

This asymptotic formula was derived in [19]. It is shown in figure 3(a) that it matches excellently the numerically found skyrmion radius for small R in the range $0 \lt \lambda \lesssim 0.1$.

Zoom In Zoom Out Reset image size

Inverting formula (25) gives, to the lowest order of approximation, equation (17). A modification of the denominator of equation (17),

$$\begin{equation} R = - \frac{\lambda}{\ln(\lambda/0.72)}, \end{equation} \tag{ 26 }$$

as obtained in [30], leads to a very good fit of the numerical data over the very wide range $0 \lt \lambda \lesssim 0.6$, as shown in figure 3(a).

Using asymptotic results from appendix, we can give approximation formulae for the skyrmion profile. The field close to the skyrmion center is given in equation (A15). In the present case ($\kappa = 1, h = 0$), this reduces to a formula identical to equation (20) in section 4. The far field is given in section 3 and it is identical to equation (21) in section 4.

Figure 4 shows formulae (20) and (21) compared to numerically calculated skyrmion profiles for two values of the parameter.

Zoom In Zoom Out Reset image size

We finally mention that including the magnetostatic field is expected to be modeled as an easy-plane anisotropy for skyrmions of small radius and should thus not change qualitatively the present results [22]. This is supported also by the results presented in section 6.

5.2. Large radius

At the parameter value $\lambda = 2/\pi$, 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 powers of R [19],

$$\begin{equation} \Theta = \Theta_0 + \tilde{\Theta},\qquad \tilde{\Theta} = \frac{\Theta_1}{R}+\frac{\Theta_2}{R^2}+\frac{\Theta_3}{R^3}+\cdots \end{equation} \tag{ 27 }$$

where $\Theta_0, \Theta_1, \Theta_2, \Theta_3,\ldots$ are functions of r − R. The leading order contribution satisfies

$$\begin{equation} \tan \left(\frac{\Theta_0}{2}\right) = {\mathrm{e}}^{-(r-R)}. \end{equation} \tag{ 28 }$$

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 [20]

$$\begin{equation} \Theta(r) = \pi - {\mathrm{e}}^{-R}\sqrt{8\pi R}\, I_1(r) \end{equation} \tag{ 29 }$$

where $I_1(r)$ is the modified Bessel function of the first kind. For the far field, we have

$$\begin{equation} \Theta(r) = {\mathrm{e}}^R \sqrt{\frac{8R}{\pi}}\, K_1(r) \end{equation} \tag{ 30 }$$

that simplifies to

$$\begin{equation} \Theta(r) = 2 \;\frac{{\mathrm{e}}^{-(r-R)}}{\sqrt{r/R}},\qquad r - R \gg 1 \end{equation} \tag{ 31 }$$

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 [20].

Figure 5 shows a numerically calculated skyrmion profile and this is compared to the leading order formulae (28)–(30).

Zoom In Zoom Out Reset image size

The relation between the parameter λ and the radius R is derived during the asymptotic process. It is given as the series [19],

$$\begin{equation} \lambda = \frac{2}{\pi} + \frac{\lambda_2}{R^2} + \frac{\lambda_4}{R^4} + \cdots,\qquad R \gg 1 \end{equation} \tag{ 32 }$$

where

$$\begin{equation} \lambda_2 \approx -0.3057,\quad\lambda_4 \approx -0.8792. \end{equation} \tag{ 33 }$$

Retaining terms up to $O(1/R^2)$ and solving for R, we obtain

$$\begin{equation} R = \frac{|\lambda_2|^{1/2}}{\left(\frac{2}{\pi}-\lambda \right)^{1/2}}. \end{equation} \tag{ 34 }$$

This is shown in figure 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.

One could go on to consider the case when both perpendicular anisotropy and an external field are present. Section 3 and the appendix already contain the basic asymptotic results for the far field and the near field respectively, in the general case, and, therefore, a study for the case of combined anisotropy and external field can be done following the steps presented in this section and in section 4.

The magnetostatic field is known to be crucial for stabilizing magnetic bubbles [31] which are cylindrical domains that share topological features with skyrmions. Magnetic bubbles, even ones with high skyrmion numbers, as well as topologically trivial bubbles, with skyrmion number zero, have been observed and studied [17, 23]. When the DM interaction is present, the magnetostatic field is usually considered to be of secondary importance for skyrmion generation and stability. However, it has been shown that the study of the combined effect of DM interaction and magnetostatic field is necessary for a proper understanding of the relation between chiral skyrmions and magnetic bubbles, and it has revealed unusual phenomena such as bi-stability [17, 23]. A detailed calculation of the skyrmion energetics including the magnetostatic field [17] leads to the understanding of skyrmion and bubble stability. In this section, we present the results of numerical simulations for the profile of skyrmions when the magnetostatic field is included. We give two examples aiming to make a comparison between bubble and skyrmion profiles in specific cases.

We first summarize the results on the stability of cylindrical domains [32] 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. The bubble domain wall is primarily of Bloch type as this is favored 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 [33, 34].

The static Landau–Lifshitz equation has the form

$$\begin{equation} \boldsymbol{m} \times \left( \partial_\mu\partial_\mu\boldsymbol{m} - 2\lambda\, \boldsymbol{\hat{e}}_\mu\times\partial_\mu\boldsymbol{m} + m_3 \boldsymbol{\hat{e}}_3 + \frac{h}{\kappa}\boldsymbol{\hat{e}}_3 + \frac{\mathbf{h}_m}{\kappa} \right) = 0 \end{equation} \tag{ 35 }$$

where the bulk DM interaction term has been chosen, h_m is the magnetostatic field normalized to the saturation magnetization, and $\ell_\textrm{w}$ is used as the unit of length. The parameters are given here again for convenience

$$\begin{equation} \lambda = \frac{D}{2\sqrt{AK}},\qquad \kappa = \frac{2K}{\mu_0 M\,_{\mathrm{s}}^2}. \end{equation} \tag{ 36 }$$

We often consider the substitution $K \to K - \frac{1}{2}\mu_0 M_{\mathrm{s}}^2$ 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 setting $\kappa\to\kappa-1$ and $\lambda \to \kappa/(\kappa-1)\lambda$ in equation (35). After these substitutions, the field h_m would represent only the deviation of the magnetostatic field from its approximation as an easy-plane anisotropy term.

We find solutions of equation (35) using an energy relaxation algorithm [34]. 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 \leqslant z \leqslant d_f/2$. We find that solutions for the magnetization $\boldsymbol{m} = (m_r,\ m_\phi, m_z)$ satisfy the parity relations $m_r(r,-z) = -m_r(r,z),\, m_\phi(r,-z) = m_\phi(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\ell_\textrm{w}$ and parameter values $\lambda = 0.33,\;\kappa = 4$. We discretize space with lattice spacings $\Delta r = 0.1$ and $\Delta z = 0.5$ in the r and z directions respectively. The numerical mesh extends to r = 20. Figure 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.

Zoom In Zoom Out Reset image size

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 $\lambda = 2/\pi$ 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 destabilized by the magnetostatic field.

It is interesting to compare the chiral skyrmion profiles with magnetic bubble ones that are stabilized without the DM interaction. We choose a film with thickness $d_f = 11\ell_\textrm{w}$, extending over $-d_f/2 \leqslant z \leqslant d_f/2$, and the parameter values $\kappa = 2,\;h = 0.2$ (and set λ = 0). The lattice spacings are $\Delta r = 0.2,\;\Delta z = 1.0$ and the numerical mesh extends to r = 40. The energy minimization algorithm converges to a static bubble with radius R ≈ 8. Figure 7(a) 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 toward Neél type closer to the boundaries. The hybrid wall structure has been shown to persist when DM interaction is included [35, 36]. 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 [33, 34].

Zoom In Zoom Out Reset image size

Figure 7(b) 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. The issue of the differences between a chiral skyrmion and a bubbles was earlier debated in [37, 38]. We finally note that our comparison between a chiral skyrmion and a bubble is only possible for large radii as magnetic bubbles are not stabilized for small radii [32].

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 [21, 25, 39].

We consider the case of easy-axis anisotropy and no external field so that the energy is given in equation (3) for $\kappa = 1, h = 0$, and it is measured in units of $\mu_0 M_{\mathrm{s}}^2 \ell_\textrm{ex}^2 = 2A$. A useful result for the energy components of a skyrmion is obtained by employing a standard scaling argument [29],

$$\begin{equation} 2E_{\textrm{a}} + E_{\textrm{DM}} = 0. \end{equation} \tag{ 37 }$$

7.1. Small radius

For small skyrmion radius, asymptotic analysis gives [19]

$$\begin{equation} E_{\textrm{ex}} = 4\pi + o\left(\frac{\lambda^2}{\ln\lambda} \right), \quad E_{\textrm{DM}} = 8\pi \frac{\lambda^2}{\ln\lambda} + o\left(\frac{\lambda^2}{\ln\lambda} \right) \end{equation} \tag{ 38 }$$

while $E_\textrm{a}$ can be obtained from equation (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 λ₀ and the energy $E = E_{\lambda_0}(\Theta)$ for any configuration Θ. This is written as

$$\begin{equation} E_{\lambda_0}(\Theta) = E_{\textrm{ex}}(\Theta) + \lambda_0 \tilde{E}_{\textrm{DM}}(\Theta) + E_{\textrm{a}}(\Theta) \end{equation} \tag{ 39 }$$

where we have set

$$\begin{equation} E_{\textrm{DM}} = \lambda_0 \tilde{E}_{\textrm{DM}}, \end{equation} \tag{ 40 }$$

so that all terms $E_\textrm{ex}, \tilde{E}_\textrm{DM}, E_\textrm{a}$ are integrals over Θ that do not contain parameters.

We confine ourselves to the configurations $\Theta = \Theta_\lambda$ that are the profiles corresponding to static skyrmions for any parameter value λ (that may be different from λ₀). The virial relation (37) holds,

$$\begin{equation} E_{\textrm{a}} = -\frac{\lambda}{2}\tilde{E}_{\textrm{DM}}, \end{equation} \tag{ 41 }$$

and it is used in equation (39) to obtain

$$\begin{equation} E_{\lambda_0}(\lambda) = E_{\textrm{ex}} + \left(\lambda_0-\frac{\lambda}{2} \right)\tilde{E}_{\textrm{DM}} \end{equation} \tag{ 42 }$$

for the energy of profile $\Theta_\lambda$ in a system with parameter value λ₀.

We take the derivative of the energy expression (42) with respect to λ,

$$\begin{equation} E_{\lambda_0}^{^{\prime}} = E_{\textrm{ex}}^{^{\prime}} + \left( \lambda_0 - \frac{\lambda}{2} \right) \tilde{E}_{\textrm{DM}}^{^{\prime}} - \frac{1}{2}\tilde{E}_{\textrm{DM}} \end{equation} \tag{ 43 }$$

and require that $E_\lambda^{^{\prime}}(\lambda = \lambda_0) = 0$. This gives the condition

$$\begin{equation} E_{\textrm{ex}}^{^{\prime}} + \frac{1}{2} \left(\lambda \tilde{E}_{\textrm{DM}}^{^{\prime}} - \tilde{E}_{\textrm{DM}} \right) = 0\qquad \textrm{at}\,\; \lambda=\lambda_0. \end{equation} \tag{ 44 }$$

An explicit calculation of the terms of equation (44) utilizing equation (38), shows that terms of order $O(\lambda^2/(\ln\lambda)^m)$ with $m\geqslant 2$ should necessarily be included in the energy asymptotic series for the condition (44) to be satisfied. Specifically,

$$\begin{equation} \begin{split} E_{\textrm{ex}} & = 4\pi + 4\pi \sum_{m=2}^\infty a_m\,\frac{\lambda^2}{(\ln\lambda)^m},\\ \tilde{E}_{\textrm{DM}} & = 4\pi \sum_{m=1}^\infty b_m\,\frac{\lambda}{(\ln\lambda)^m} \end{split} \end{equation} \tag{ 45 }$$

where $a_m, b_m$ are constants (with $b_1 = 2$ from equation (38)). Substituting these in equation (44) obtains

$$\begin{equation} \sum_{m=1}^\infty \frac{2 a_{m+1} - m a_m - \frac{m}{2} b_m)}{(\ln\lambda)^{m+1}} = 0,\quad m=1,2,\ldots. \end{equation} \tag{ 46 }$$

In order that equation (46) be satisfied to all orders in $1/(\ln(\lambda))^m$, it should hold

$$\begin{equation} a_{m+1} = \frac{m}{4} \left( 2a_m + b_m \right)\quad m=1,2,\ldots. \end{equation} \tag{ 47 }$$

For m = 1, this gives (since $a_1 = 0$)

$$\begin{equation} a_2 = \frac{b_1}{4} = \frac{1}{2}. \end{equation} \tag{ 48 }$$

We thus have

$$\begin{equation} E_{\textrm{ex}} = 4\pi \left( 1 + \frac{1}{2}\,\frac{\lambda^2}{(\ln\lambda)^2} \right) + O\left( \frac{\lambda^2}{(\ln\lambda)^3} \right). \end{equation} \tag{ 49 }$$

More information would be required to obtain explicit values for $a_{m+1}, b_m$ with $m\geqslant 2$. Comparing the second of equation (45) with numerical data, we propose the following formula,

$$\begin{equation} E_{\textrm{DM}} = 4\pi \left( 2\frac{\lambda^2}{\ln\lambda} - \frac{1}{2} \frac{\lambda^2}{(\ln\lambda)^2} \right) + O\left( \frac{\lambda^2}{(\ln\lambda)^3} \right). \end{equation} \tag{ 50 }$$

Using formulae (49) and (50), the energy is

$$\begin{equation} E = 4\pi \left( 1 + \frac{\lambda^2}{\ln\lambda} + \frac{1}{4} \frac{\lambda^2}{(\ln\lambda)^2} \right) + O\left( \frac{\lambda^2}{(\ln\lambda)^3} \right). \end{equation} \tag{ 51 }$$

This is in agreement with rigorous results in [25] where an error term $O(\lambda^2/(\ln\lambda)^2)$ is given. The numerical data for the energy are fitted well with formulae (49)–(51), for parameter values $\lambda \lesssim 0.25$.

Figure 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,

$$\begin{equation} E = 4\pi \left( 1 + \frac{\lambda^2}{\ln(\lambda/0.75)} (1-\lambda) \right). \end{equation} \tag{ 52 }$$

This formula is plotted in figure 8. For λ close to the critical value $2/\pi$ formula (55) (of the next subsection) should be used instead.

Zoom In Zoom Out Reset image size

Equation (51) gives

$$\begin{equation} \begin{split} E \approx 4\pi \left( 1 + R^2\ln R \right),\quad R \ll 1. \end{split} \end{equation} \tag{ 53 }$$

Equation (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 equation (53),

$$\begin{equation} E = 4\pi \left[ 1 + 0.12\,R^2\ln(R/5.2) \right], \end{equation} \tag{ 54 }$$

that is plotted in figure 9.

Zoom In Zoom Out Reset image size

We finally note that the techniques used in this subsection can be readily applied to the case of skyrmions of small radius in applied external field.

7.2. Large radius

For large skyrmion radius, asymptotic analysis gives [20]

$$\begin{equation} E \sim 4\pi^2 |\lambda|^{1/2}\,\sqrt{\frac{2}{\pi}-\lambda},\quad \left|\frac{2}{\pi}-\lambda\right| \ll 1 \end{equation} \tag{ 55 }$$

or (cf equation (34))

$$\begin{equation} E = \frac{4\pi^2|\lambda_2|}{R} +O\left( \frac{1}{R^3} \right),\quad R \gg 1. \end{equation} \tag{ 56 }$$

In the range $0.6 \lesssim \lambda \lt 2/\pi$, the numerical data of figure 8 are fitted well with formula (55).

Figure 9 shows the skyrmion energy as a function of its radius. For large R, the data are fitted very well with the formula

$$\begin{equation} E = 4\pi^2|\lambda_2| \left( \frac{1}{R} + \frac{5}{2}\frac{1}{R^3} \right) \end{equation} \tag{ 57 }$$

which is a sharpening of equation (56) and it is plotted in the figure. We further note that each one of the energy components is fitted very well, for $R \gtrsim 3$, by the forms

$$\begin{equation} \begin{split} E_{\textrm{ex}} & = 2\pi R + \pi^2|\lambda_2| \left(\frac{3}{R} + \frac{5}{R^3} \right) \\ E_{\textrm{a}} & = 2\pi R - \pi^2|\lambda_2| \left(\frac{1}{R} + \frac{5}{R^3} \right). \end{split} \end{equation} \tag{ 58 }$$

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 [19, 20]. The analysis is extended to include the effect 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 in two specific cases. It tends to increase the skyrmion radius while it appears to destabilize skyrmions of large radius. This indicates that attention should be focused on skyrmions of small radius while larger radius skyrmions should be relevant for very thin films where the magnetostatic field is equivalent to an anisotropy term. We compared the profile of magnetic bubbles stabilized without the chiral DM 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. The present study of the effect of the magnetostatic field is confined to specific cases, and it is consistent with more extensive studies that study extensive details of this fairly complicated issue [17, 23].

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 results of section 3 for the far field, the formulae of appendix for the near field, and the formulation for skyrmions of large radius apply for every axially symmetric skyrmion configuration, such as the skyrmionium or kπ skyrmions [40, 41]. Although not a straightforward task, it may be possible to use the methods employed here in order to derive further detailed results for kπ skyrmions, such as the relation between radius and DM parameter.

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. We also hope to motivate work that may prove that some of the fitting formulae, or variations of them, can be obtained by accurate calculations.

The publication of the article in OA mode was financially supported by HEAL-Link.

The data that support the findings of this study are available upon reasonable request from the authors.

Appendix. Near field

We will mainly work with the field $\tilde{v}(r) = ru(r)$ defined in equation (14). It satisfies the equation

$$\begin{equation} \tilde{v}^{^{\prime\prime}} + \frac{\tilde{v}^{^{\prime}}}{r}\, \frac{3\tilde{v}^2-r^2}{\tilde{v}^2+r^2} - \frac{2 \tilde{v} \tilde{v}^{^{\prime} 2}}{\tilde{v}^2 + r^2} + 4\lambda \frac{\tilde{v}^2}{\tilde{v}^2+r^2} + \kappa\frac{\tilde{v}^2-r^2}{\tilde{v}^2+r^2} \tilde{v} - h \tilde{v} = 0. \end{equation} \tag{ A1 }$$

The ensuing analysis will follow the techniques developed in [19].

An assumption which is found to be consistent with our results is to neglect the term containing $\tilde{v}^{^{\prime} 2}$ in equation (A1),

$$\begin{equation} \tilde{v}^{^{\prime\prime}} + \frac{\tilde{v}^{^{\prime}}}{r}\, \frac{3\tilde{v}^2-r^2}{\tilde{v}^2+r^2} + 4\lambda \frac{\tilde{v}^2}{\tilde{v}^2+r^2} + \kappa\frac{\tilde{v}^2-r^2}{\tilde{v}^2+r^2} \tilde{v} - h \tilde{v} = 0. \end{equation} \tag{ A2 }$$

In the near field, we replace $\tilde{v}$ in equation (A2) with its initial value $\tilde{v}(0) = \tilde{v}_0$. The approximation is valid for the range of r over which $\tilde{v}_0-\tilde{v}\ll \tilde{v}_0$. The equation obtained in this way,

$$\begin{equation} \tilde{v}^{^{\prime\prime}} + \frac{\tilde{v}^{^{\prime}}}{r}\, \frac{3\tilde{v}_0^2-r^2}{\tilde{v}_0^2+r^2} = -4\lambda \frac{\tilde{v}_0^2}{\tilde{v}_0^2+r^2} - \kappa\frac{\tilde{v}_0^2-r^2}{\tilde{v}_0^2+r^2} \tilde{v}_0 + h \tilde{v}_0 \end{equation} \tag{ A3 }$$

is integrable by the integrating factor

$$\begin{equation*} \mu(r) = \frac{r^3}{(\tilde{v}_0^2+r^2)^2}. \end{equation*}$$

We obtain

$$\begin{equation*} \frac{{\mathrm{d}}}{{\mathrm{d}}r} \left[ \frac{r^3}{(\tilde{v}_0^2+r^2)^2} v^{^{\prime}} \right] = -4\lambda \frac{\tilde{v}_0^2 r^3}{(\tilde{v}_0^2+r^2)^3} -\kappa\frac{\tilde{v}_0^2 - r^2}{(\tilde{v}_0^2+r^2)^3} \tilde{v}_0 r^3 + h \frac{\tilde{v}_0 r^3}{(\tilde{v}_0^2+r^2)^2}, \end{equation*}$$

which integrates to

$$\begin{equation*} \frac{r^4}{(\tilde{v}_0^2+r^2)^2}\frac{\tilde{v}^{^{\prime}}}{r} = -\lambda \frac{r^4}{(\tilde{v}_0^2+r^2)^2} + \frac{\tilde{v}_0}{2} \left[ (\kappa+h) \ln\left( 1 + \frac{r^2}{\tilde{v}_0^2} \right) - \kappa \frac{\tilde{v}_0^2+2r^2}{(\tilde{v}_0^2+r^2)^2} r^2 - h \frac{r^2}{\tilde{v}_0^2+r^2} \right]. \end{equation*}$$

We finally have

$$\begin{equation} \frac{\tilde{v}^{^{\prime}}}{r} = -\lambda + \frac{\tilde{v}_0}{2} \left[ (\kappa+h) \frac{(\tilde{v}_0^2+r^2)^2}{r^4} \ln\left( 1 + \frac{r^2}{\tilde{v}_0^2} \right) - \kappa \frac{\tilde{v}_0^2+2r^2}{r^2} - h \frac{\tilde{v}_0^2+r^2}{r^2}\right]. \end{equation} \tag{ A4 }$$

The constant of integration has been judiciously chosen to eliminate the fourth-order singularity r⁻⁴ on the right.

We make the change of variables $w(\tau) = \tilde{v}/\tilde{v}_0$ and $\tau = r^2/\tilde{v}_0^2$. Then equation (A4) becomes

$$\begin{equation} \frac{{\mathrm{d}}w}{{\mathrm{d}}\tau} = -\lambda \frac{\tilde{v}_0}{2} +(\kappa + h) \frac{\tilde{v}_0^2}{4} \left[ \ln(1+\tau) + 2 \frac{\ln(1+\tau)}{\tau} + \frac{\ln(1+\tau)}{\tau^2} - \frac{1}{\tau} \right] - (2\kappa+h) \frac{\tilde{v}_0^2}{4}. \end{equation} \tag{ A5 }$$

This is integrated to give

$$\begin{equation*} w(\tau) = 1 - \lambda\frac{\tilde{v}_0}{2} \tau + (\kappa + h) \frac{\tilde{v}_0^2}{4} \left[ \tau \ln(1+\tau) - \frac{\ln(1+\tau)}{\tau} + 1 - 2{\textrm{Li}}_2(-\tau) - \tau \right] - (2\kappa+h) \frac{\tilde{v}_0^2}{4} \tau \end{equation*}$$

or

$$\begin{equation} w(\tau) = 1 - \left( \lambda + \frac{3}{2}\kappa \tilde{v}_0 + h \tilde{v}_0 \right) \frac{\tilde{v}_0}{2} \tau + (\kappa + h) \frac{\tilde{v}_0^2}{4} \left[ \tau \ln(1+\tau) - \frac{\ln(1+\tau)}{\tau} + 1 - 2{\textrm{Li}}_2(-\tau) \right] \end{equation} \tag{ A6 }$$

where $\textrm{Li}_2(x)$ is the dilogarithm function and the constant of integration has been chosen so that $w(0) = 1$.

Keeping the dominant terms for large τ we have

$$\begin{equation} w(\tau) = 1 - \left(\lambda + \frac{3}{2}\kappa \tilde{v}_0 + h \tilde{v}_0 \right) \frac{\tilde{v}_0}{2}\tau + (\kappa + h) \frac{\tilde{v}_0^2}{4} \tau\,\, \ln\tau \end{equation} \tag{ A7 }$$

which, in the original variables, gives the near field

$$\begin{equation} \tilde{v}_N(r) = \tilde{v}_0 - \left[ \lambda + \kappa \tilde{v}_0 \left(\frac{3}{2} + \ln \tilde{v}_0 \right) + h \tilde{v}_0 ( 1 + \ln\tilde{v}_0 ) \right] \frac{r^2}{2} + (\kappa + h) \tilde{v}_0 \frac{r^2}{2} \ln r. \end{equation} \tag{ A8 }$$

The skyrmion radius is found by the condition $\tilde{v}(r = R) = R$. In the case that the skyrmion radius R is small, this is given by

$$\begin{equation} R\approx \tilde{v}_0. \end{equation} \tag{ A9 }$$

Appendix. Asymptotic matching

The solution of the equation (7) for the far field is

$$\begin{equation} \Theta \propto \frac{1}{\tilde{s}} + (2\gamma-1) \tilde{s} + 2 \sum_{n=1}^\infty \frac{(\ln \tilde{s} - \ln n) \tilde{s}^{2n-1}}{(n-1)!\, n!} - 2 \sum_{n=2}^\infty \frac{\xi_n \tilde{s}^{2n-1}}{(n-1)!\, n!},\qquad \tilde{s}=\sqrt{\kappa+h}\,\frac{r}{2} \end{equation} \tag{ A10 }$$

where ξ_n is given in equation (9) and the solution is valid for any $\kappa, h$. For small values of the radial coordinate, we keep only the first three terms in the series,

$$\begin{equation*} \Theta \propto \frac{1}{\tilde{s}} \left[ 1 + (2\gamma-1) \tilde{s}^2 + 2 \tilde{s}^2 \ln \tilde{s} \right]. \end{equation*}$$

The far field for the field $\tilde{v}$ is derived from the latter formula noting that $\tilde{v} \approx \frac{r}{2}\Theta$ for small values of $\tilde{v}, \Theta$, as seen from equations (14) and (12). We find

$$\begin{equation*} \tilde{v}_F(r) = C \left[ 1 + (\kappa + h) (2\gamma-1) \frac{r^2}{4} + (\kappa + h) \frac{r^2}{2} \ln \left(\frac{\sqrt{\kappa+h}}{2}\,r \right) \right] \end{equation*}$$

or

$$\begin{equation} \tilde{v}_F(r) = C \left[ 1 + (\kappa + h) \left( \gamma-\frac{1}{2} + \ln \left(\frac{\sqrt{\kappa+h}}{2} \right) \right) \frac{r^2}{2} + (\kappa + h) \frac{r^2}{2} \ln r \right]. \end{equation} \tag{ A11 }$$

Matching the near field in equation (A8) with the far field in equation (A11) order by order we have $C = \tilde{v}_0$ and

$$\begin{equation*} \lambda + \kappa\tilde{v}_0 \left( \frac{3}{2} + \ln\tilde{v}_0 \right) + h\tilde{v}_0 (1+\ln\tilde{v}_0) = - (\kappa + h) \left[ \gamma-\frac{1}{2} + \ln \left(\frac{\sqrt{\kappa+h}}{2} \right) \right]\,\tilde{v}_0 \end{equation*}$$

or

$$\begin{equation} \lambda = -\tilde{v}_0 \left[ \kappa (\gamma + 1) + h \left( \gamma + \frac{1}{2} \right) + (\kappa+h) \ln\left( \frac{\sqrt{\kappa+h}}{2}\, \tilde{v}_0 \right) \right]. \end{equation} \tag{ A12 }$$

Combined with equation (A9), relation (A12) gives implicitly the skyrmion radius as a function of the parameters,

$$\begin{equation} \lambda = -R \left[ \kappa (\gamma + 1) + h \left( \gamma + \frac{1}{2} \right) + (\kappa+h) \ln\left( \frac{\sqrt{\kappa+h}}{2}\, R \right) \right]. \end{equation} \tag{ A13 }$$

Appendix. Field at the skyrmion center

The expression in the square bracket in equation (A6) has the expansion

$$\begin{equation*} \tau \ln(1+\tau) - \frac{\ln(1+\tau)}{\tau} + 1 - 2{\textrm{Li}}_2(-\tau) = \frac{5}{2} \tau + \frac{\tau^2}{6} + \ldots,\qquad \tau \ll 1. \end{equation*}$$

Thus, the near field, close to the skyrmion center is

$$\begin{equation*} w(\tau) \approx 1 - \left( \lambda + \frac{3}{2}\kappa \tilde{v}_0 + h \tilde{v}_0 \right) \frac{\tilde{v}_0}{2} \tau + (\kappa + h) \frac{\tilde{v}_0^2}{4} \left( \frac{5}{2} \tau + \frac{\tau^2}{6} \right). \end{equation*}$$

Using equation (A6), this gives

$$\begin{equation} w(\tau) \approx 1 + (\kappa+h) \left[ \gamma + \frac{3}{4}+ \ln\left( \frac{\sqrt{\kappa+h}}{2}\, \tilde{v}_0 \right) \right] \frac{\tilde{v}_0^2}{2}\tau + (\kappa+h) \frac{\tilde{v}_0^2}{24}\tau^2. \end{equation} \tag{ A14 }$$

For the field $\tilde{v}$, the latter gives

$$\begin{equation} \frac{\tilde{v}(r)}{\tilde{v}_0} \approx 1 + (\kappa+h) \left[ \gamma + \frac{3}{4}+ \ln\left( \frac{\sqrt{\kappa+h}}{2}\, \tilde{v}_0 \right) \right] \frac{r^2}{2} + (\kappa+h) \frac{r^4}{24\tilde{v}_0^2},\qquad r \ll \tilde{v}_0. \end{equation} \tag{ A15 }$$