1 Introduction

Magnetic skyrmions are topologically non-trivial configurations which occur in certain magnetic materials. It was first observed in [1] that particular examples of such configurations are minimisers of natural energy expressions for the magnetisation vector. They have since then become the subject of intense study, both experimentally and theoretically, not least because of their potential use as information carriers in magnetic storage devices, see [2] for a review.

In this paper we use tools from gauge theory, complex analysis and differential geometry to introduce models for magnetic skyrmions which can be solved explicitly. The models are grounded in the physics of magnetic skyrmions and belong to the general class already considered in [1], with an energy expression consisting of a Dirichlet term, a Dzyaloshinskii–Moriya (DM) interaction energy [3, 4] and a potential combining Zeeman and easy plane anisotropy terms. However, our formulation reveals that, for critical values of the coupling constants, the models are of Bogomol’nyi type, which means that static solutions can be obtained by solving a first order partial differential equation. Moreover, this equation can be solved in terms of an arbitrary holomorphic function. Both the Bogomol’nyi property and the existence of an infinite family of explicit solutions are novel in the context of magnetic skyrmions.

Models of Bogomol’nyi type have historically played an important role in the study of topological solitons [5]. They require a particular choice of coupling constants, but their mathematical properties allow for a far more detailed and explicit study of the solitons and their dynamics than would be possible in the generic case. Intricate and surprising features of soliton dynamics such as shapes and symmetries of multi-soliton configurations or scattering behaviour were first observed in models of Bogomol’nyi type and later found in generic soliton models.

The focus of this paper is the mathematical structure of the critically coupled models for magnetic skyrmions, and we only begin to explore the properties of our solutions. However, even at this stage it is clear that our infinite family of solutions contains several of the ground state configurations associated with the various phases of generic models [6,7,8], and that it illustrates elliptical deformations [9, 10] and the recently studied skyrmion bags [11] or ‘sacks’ [12].

The simplest topological soliton theory of Bogomol’nyi type is the O(3)-sigma model in the plane [13]. The basic field is a map from the plane to the sphere, and finiteness of the Dirichlet energy requires the field to tend to a constant at spatial infinity, and to extend to a map from sphere to sphere. Such maps have a topological and integer degree, which is physically interpreted as the soliton number. The energy is bounded below by a multiple of the absolute value of the degree, and this bound is attained by configurations which satisfy the first order Bogomol’nyi equation. In this particular case, the Bogomol’nyi equation requires the configuration to be a holomorphic or anti-holomorphic map to the Riemann sphere.

In analogy with the baby skyrme model [14], one would expect the inclusion of DM interactions, Zeeman potential and anisotropy terms inevitably to destroy the Bogomol’nyi property of the pure O(3) sigma model. However, here we shall show that, with a careful choice of potential and for a one-parameter family of DM interaction terms, our models preserve the Bogomol’nyi property and can be solved in terms of a fixed anti-holomorphic and an arbitrary holomorphic map to the Riemann sphere. The fixed anti-holomorphic part turns out to be an analytical version of the usual Bloch or Néel magnetic skyrmions, but the holomorphic part is new.

Our family of models is introduced in Sect. 2, and the solutions are studied in Sect. 5. Readers primarily interested in the models and their solutions are invited to skip directly from Sect. 2 to Sect. 5. In the intervening sections we derive the Bogomol’nyi equation in two different ways. In Sect. 3, we write the theory as a non-abelian gauge theory with a fixed non-abelian gauge field and apply a trick for constructing gauged sigma models of Bogomol’nyi type introduced in [15]. In Sect. 4, we derive the Bogomol’nyi equation in complex stereographic coordinates and give the general solution in terms of fixed anti-holomorphic and an arbitrary holomorphic function. We show that, when that holomorphic function is rational, the energy is generically positive and quantised in multiples of \(4\pi \). We also point out that the energy is not well-defined when the leading holomorphic term at spatial infinity is linear, and propose a regularisation which preserves the generic formula. We study detailed properties of rational solutions in Sect. 5. The final Sect. 6 contains our conclusion and a brief discussion of open questions.

2 The Model

2.1 Energy and symmetry

The basic field in any mathematical model for magnetic skyrmions is the magnetisation vector, which in the planar and static case is a map \({\varvec{n}}:{\mathbb {R}}^2 \rightarrow S^2\). Here we consider models where the energy of a configuration is measured by the sum of three terms: the Dirichlet or Heisenberg energy (quadratic in derivatives), a generalised DM interaction energy (linear in derivatives), and a potential term which may involve linear or quadratic terms in the Cartesian components of \({\varvec{n}}=(n_1,n_2,n_3)^t\) (with the constraint \(n_1^2+n_2^2+n_3^2=1\) always understood). General models of this sort were considered in the seminal paper [1] in which the possibility of topologically stable configurations now known as magnetic skyrmions was first pointed out, and have been widely studied since then.

Our one-parameter family of models belongs to the general family considered in [1] but requires a particular choice of coupling constants, which we call critical. The parameter in the family is an angle \(\alpha \in [0,2\pi )\), and describes a generalised DM interaction term. In order to define this term, we need some additional notation.

We use Cartesian coordinates \(x_1\) and \(x_2\) in the plane and write \(\partial _1\) and \(\partial _2\) for partial derivatives with respect to them. Thinking of \({\mathbb {R}}^2\) as embedded in Euclidean \({\mathbb {R}}^3\) and using three-dimensional notation, we also write \({\varvec{e}}_1, {\varvec{e}}_2,{\varvec{e}}_3\) for the canonical basis of \({\mathbb {R}}^3\), with \({\varvec{e}}_3={\varvec{e}}_1\times {\varvec{e}}_2\). In terms of the rotation \( R(\alpha )\) about the 3-axis by \(\alpha \in [0,2\pi )\), we define

$$\begin{aligned} {\varvec{e}}_1^\alpha = R(\alpha ) {\varvec{e}}_1 = \begin{pmatrix} \cos \alpha \\ \sin \alpha \\ 0 \end{pmatrix}, \qquad {\varvec{e}}_2^\alpha = R(\alpha ) {\varvec{e}}_2 = \begin{pmatrix} -\sin \alpha \\ \cos \alpha \\ 0 \end{pmatrix}, \qquad {\varvec{e}}_3= \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}. \end{aligned}$$
(2.1)

Extending the usual definition of the gradient \(\nabla = \sum _{i=1}^2 {\varvec{e}}_i\partial _i\) to write

$$\begin{aligned} \nabla _\alpha =\sum _{i=1}^2 {\varvec{e}}_i^\alpha \partial _i, \end{aligned}$$
(2.2)

and defining

$$\begin{aligned} {\varvec{n}}^\alpha = R(\alpha ) {\varvec{n}}, \end{aligned}$$
(2.3)

the family of DM interaction terms we are interested in is

$$\begin{aligned} {\varvec{n}}^\alpha \cdot \nabla \times {\varvec{n}}^\alpha = {\varvec{n}}\cdot \nabla _{-\alpha } \times {\varvec{n}}. \end{aligned}$$
(2.4)

In components, it consist of two familiar parts:

$$\begin{aligned} {\varvec{n}}\cdot \nabla _{-\alpha } \times {\varvec{n}}= \cos \alpha \, w_B +\sin \alpha \, w_N, \end{aligned}$$
(2.5)

where \(w_b\) and \(w_N\) are the following contractions of the chirality tensor \({\varvec{n}}\times \partial _i {\varvec{n}}\):

$$\begin{aligned} w_B&= n_1\partial _2n_3 -n_2\partial _1 n_3 + n_3(\partial _1n_2 -\partial _2n_1),\nonumber \\ w_N&= - n_1\partial _1n_3 -n_2\partial _2 n_3 + n_3(\partial _1n_1 +\partial _2n_2). \end{aligned}$$
(2.6)

Choosing energy units so that the coefficient of the Dirichlet energy term is unity, the family of energy functionals we want to consider is

$$\begin{aligned} E[{\varvec{n}}]=\int _{{\mathbb {R}}^2} \frac{1}{2}(\nabla {\varvec{n}})^2 + \kappa {\varvec{n}}\cdot \nabla _{-\alpha } \times {\varvec{n}}+ \frac{\kappa ^2}{2}(1-n_3)^2 \, \mathrm {d}x_1\mathrm {d}x_2. \end{aligned}$$
(2.7)

The coupling constant \(\kappa \) could also be set to unity by a choice of length unit, but we find it convenient to keep it in our calculations. We assume \(\kappa >0\) for the remainder of this paper. Since

$$\begin{aligned} \frac{1}{2} (1-n_3)^2 = (1-n_3) - \frac{1}{2} (1-n_3^2), \end{aligned}$$
(2.8)

our potential term combines a Zeeman and easy plane anisotropy potential, with coefficients chosen in such a way that the sum is a perfect square.

The energy expression (2.7) is invariant under translations in the plane and under the group O(2) of rotations and reflection of the plane, combined with simultaneous rotations and reflections of the target sphere. On fields, the rotations act as

$$\begin{aligned} {\varvec{n}}(x_1,x_2) \mapsto R(\sigma ) {\varvec{n}}(\cos \sigma \, x_1 - \sin \sigma \, x_2, \sin \sigma \,x_1 + \cos \sigma \, x_2), \quad \sigma \in [0,2\pi ), \end{aligned}$$
(2.9)

and the generator of reflections acts as

$$\begin{aligned} {\varvec{n}}(x_1,x_2) \mapsto R(2\gamma ) {\bar{{\varvec{n}}}}( x_1, {-}x_2), \quad \text {with} \quad {\bar{{\varvec{n}}}} {=} \begin{pmatrix} n_1 \\ {-}n_2 \\ n_3 \end{pmatrix}\quad \text {and} \quad \gamma {=} \frac{\pi }{2} {-}\alpha . \end{aligned}$$
(2.10)

The symmetry group is smaller than that of generic baby skyrme models [14] because the DM interaction breaks the product of the orthogonal groups in space and target space to a diagonal subgroup. It is worth noting that it would be mathematically natural to consider an alternative version of the DM interaction with the opposite chirality:

$$\begin{aligned} {\bar{{\varvec{n}}}} \cdot \nabla _{-\alpha } \times {\bar{{\varvec{n}}}}. \end{aligned}$$
(2.11)

If this term was used instead of the standard DM interaction, the symmetry would consist of spatial rotations as in (2.9) but combined with rotations \(R(-\sigma )\) of \({\varvec{n}}\). Both the DM interaction term (2.4) and the flipped version (2.11) were considered in [16], together with a generic linear combination of the Zeeman and anisotropy potential \((1-n_3^2)\). We will work with the DM interaction (2.4) and primarily consider the critical linear combination (2.8), but we discuss the more general potential in Sect. 2.2 and comment on how our results would change if we had used the opposite chirality in the Conclusion.

At this point we should really specify which boundary conditions we impose on the field \({\varvec{n}}\) at spatial infinity. It is a priori not clear if one should, as in the discussion of the O(3) sigma model, consider only configurations which extend to continuous maps \(S^2\rightarrow S^2\). If they did, then

$$\begin{aligned} Q[{\varvec{n}}]=\frac{1}{4\pi }\int _{{\mathbb {R}}^2} {\varvec{n}}\cdot \partial _1{\varvec{n}}\times \partial _2{\varvec{n}}\, \mathrm {d}x_1 \mathrm {d}x_2 \end{aligned}$$
(2.12)

would automatically be an integer, giving the degree of the extended map.

As we shall see, we should in fact allow for configurations which do not have a continuous extension \(S^2\rightarrow S^2\). Such maps do not have a topological degree, but the integral expression for Q still plays an important role. We will refer to it as degree throughout this paper.

In our model, the degree occurs invariably in conjunction with a term which also depends on the boundary behaviour, namely the total vortex strength

$$\begin{aligned} \Omega [{\varvec{n}}]= \frac{1}{4\pi }\int _{{\mathbb {R}}^2} \omega \, \mathrm {d}x_1 \mathrm {d}x_2, \end{aligned}$$
(2.13)

where the integrand is the vorticity of the first two components of \({\varvec{n}}^\alpha \):

$$\begin{aligned} \omega = \kappa ( \partial _1n_2^\alpha - \partial _2n_1^\alpha ). \end{aligned}$$
(2.14)

The expressions we have given for the total energy, the degree and the total vortex strength should all be interpreted as functionals on the space of magnetisation fields. For some magnetisation fields \({\varvec{n}}\), the relevant integrals may not be well-defined. One might therefore want to restrict the following discussion to a class of configurations which have a well-defined total energy, total vortex strength and degree. However, as we shall see, it is impossible to do this without discarding some of the most interesting configurations which arise as solutions in our model. In order to keep the discussion general but also mathematically rigorous, we therefore need notation for the restrictions of the integrals (2.7), (2.12) and (2.13) to compact subsets \(D\subset {\mathbb {R}}^2\). We define

$$\begin{aligned} E_D[{\varvec{n}}]&=\int _{D } \frac{1}{2}(\nabla {\varvec{n}})^2 + \kappa {\varvec{n}}\cdot \nabla _{-\alpha } \times {\varvec{n}}+ \frac{\kappa ^2}{2}(1-n_3)^2 \, \mathrm {d}x_1\mathrm {d}x_2, \nonumber \\ Q_D[{\varvec{n}}]&=\frac{1}{4\pi }\int _{D} {\varvec{n}}\cdot \partial _1{\varvec{n}}\times \partial _2{\varvec{n}}\, \mathrm {d}x_1 \mathrm {d}x_2, \nonumber \\ \Omega _D[{\varvec{n}}]&= \frac{1}{4\pi }\int _D \omega \, \mathrm {d}x_1 \mathrm {d}x_2, \qquad \qquad \qquad D\subset {\mathbb {R}}^2. \end{aligned}$$
(2.15)

Clearly, \(\omega \, \mathrm {d}x_1\wedge \mathrm {d}x_2 =\mathrm {d}\Theta \), where

$$\begin{aligned} \Theta =\kappa (n_1^\alpha \mathrm {d}x_1+ n_2^\alpha \mathrm {d}x_2) \end{aligned}$$
(2.16)

is a differential one-form which plays a central role in this paper. It then follows that

$$\begin{aligned} \Omega _D[{\varvec{n}}]= \frac{1}{4\pi }\int _{\partial D} \Theta . \end{aligned}$$
(2.17)

In particular, one can take D to be a disk of radius R and centred at the origin, so that \(\partial D= C^R\) is the circle of radius R. For some of the configurations we consider, the limit

$$\begin{aligned} \Omega ^\circ [{\varvec{n}}] = \frac{1}{4\pi }\lim _{R\rightarrow \infty } \int _{C^R} \Theta \end{aligned}$$
(2.18)

exists even when the integral defining \(\Omega [{\varvec{n}}]\) does not. We will treat \(\Omega ^\circ [{\varvec{n}}]\) as a regularised total vortex strength in those cases. When \(\Omega \) is well-defined it necessarily coincides with \(\Omega ^\circ \).

For the solutions we construct in this paper, the total vortex strength is generically well-defined and finite, and, as a consequence, the total energy turns out to be given by a simple formula. The regularisation (2.18) is such that the resulting regularised energy for the non-generic solutions naturally fits into this general formula. However, we should stress that the regularisation procedure is not essential for our main results.Footnote 1

Postponing a more detailed discussion of allowed configurations and topological invariants to Sect. 4 and the Conclusion, we now derive the variational equation for (2.7), only assuming that \({\varvec{n}}\) is twice differentiable. By considering the variation \(\delta {\varvec{n}}= \varvec{\epsilon }\times {\varvec{n}}\) for an infinitesimal vector function \(\varvec{\epsilon }\) which vanishes rapidly at spatial infinity, we obtain

$$\begin{aligned} 2\kappa ({\varvec{n}}\cdot \nabla _{-\alpha }) {\varvec{n}}= \left( \Delta {\varvec{n}}+ \kappa ^2 (1-n_3){\varvec{e}}_3\right) \times {\varvec{n}}. \end{aligned}$$
(2.19)

We will show that this equation is in fact implied by a first order equation.

2.2 Hedgehog fields

In the skyrmion literature, the magnetisation \({\varvec{n}}\) is often described in terms of spherical polar coordinates, defined via

$$\begin{aligned} {\varvec{n}}=\begin{pmatrix} \sin \theta \cos \phi \\ \sin \theta \sin \phi \\ \cos \theta \end{pmatrix}, \end{aligned}$$
(2.20)

where \(\theta \) and \(\phi \) are functions on the plane. This parametrisation is particularly useful when considering hedgehog fields. By definition, and using polar coordinates \((r,\varphi )\) in the plane, hedgehog fields have a profile \(\theta \) which depends on r only and a longitudinal angle \(\phi \) which is related to \(\varphi \) according to

$$\begin{aligned} \phi = \varphi + \gamma , \end{aligned}$$
(2.21)

for a constant angle \(\gamma \). Such fields are invariant under the rotational symmetry (2.9). They are additionally invariant under the reflection symmetry (2.10) if and only if we choose \(\gamma \) to be the complementary angle of \(\alpha \) as in (2.10), and we now make this choice. With the boundary condition

$$\begin{aligned} \theta (0)=\pi , \qquad \theta (\infty ) = 0, \end{aligned}$$
(2.22)

one checks that hedgehog fields have degree \(Q=-1\).

Before developing the general machinery for generating solutions of the equation (2.19), we note some properties of the much simpler hedgehog solutions in our model. We do this in a slightly more general family of models, obtained from (2.7) by replacing

$$\begin{aligned} \frac{\kappa ^2}{2}(1-n_3)^2 \rightarrow \frac{\mu ^2}{2}(1-n_3)^2 \end{aligned}$$
(2.23)

for a further real constant \(\mu \). Minimisers of the resulting energy functional were studied in [17], and those of degree \(Q=-1\) were shown to have holomorphicity properties similar to the ones which we will demonstrate more generally for stationary points of (2.7). We will exhibit these properties in our discussion of example solutions in Sect. 5. Here we derive the profile of hedgehog solutions by a more pedestrian method.

For hedgehog fields and \(\gamma = \frac{\pi }{2}-\alpha \), the energy expression with the replacement (2.23) is

$$\begin{aligned} E=2\pi \int _0^\infty rdr \left( \frac{1}{2} \left( \frac{d \theta }{dr}\right) ^2 + \frac{\sin ^2\theta }{2 r^2} +\kappa \left( \frac{d\theta }{dr} +\frac{\sin (2\theta )}{2r}\right) +\frac{\mu ^2}{2}(1-\cos \theta )^2\right) .\nonumber \\ \end{aligned}$$
(2.24)

The Euler–Lagrange equation is

$$\begin{aligned} \frac{d^2\theta }{dr^2} =- \frac{1}{r} \frac{d\theta }{dr} + \frac{\sin (2\theta )}{2r^2} - 2 \kappa \frac{\sin ^2 \theta }{r} +\mu ^2 \sin \theta (1-\cos \theta ). \end{aligned}$$
(2.25)

With the boundary condition (2.22), this is solved by

$$\begin{aligned} \theta = 2\tan ^{-1}\left( \frac{2\kappa }{\mu ^2 r}\right) . \end{aligned}$$
(2.26)

As already advertised, we will recover this profile from the simplest solution of a first order Bogomol’nyi equation in the case \(\mu =\kappa \) in equation (5.2).

3 A Bogomol’nyi Equation for Magnetic Skyrmions

We will now show that the energy functional (2.7) can be written as the sum of a squared expression and a linear combination of the integral expression for the degree (2.12) and the total vortex strength (2.13). The vanishing of the squared expression gives a first order Bogomol’nyi equation which implies the variational second order equation (2.19).

Our derivation of the Bogomol’nyi equation is inspired by a similar treatment of gauged sigma models in [15] and [18]. As noticed in [19], the combination

$$\begin{aligned} \partial _i{\varvec{n}}-\kappa {\varvec{e}}_i\times {\varvec{n}}, \qquad i=1,2, \end{aligned}$$
(3.1)

which occurs in many calculations involving magnetic skyrmions and which is often called ‘helical derivative’ can be thought of as a covariant derivative with respect to a non-abelian gauge field. To see the benefits of this, we take a more general viewpoint and consider more general su(2) gauge fields.

To minimise notation, we identify the su(2) Lie algebra with \({\mathbb {R}}^3\) and the Lie algebra commutator with the vector product. Defining the covariant derivative of \({\varvec{n}}\) as

$$\begin{aligned} D_i{\varvec{n}}= \partial _i {\varvec{n}}+ {\varvec{A}}_i \times {\varvec{n}}, \end{aligned}$$
(3.2)

and the non-abelian field strength

$$\begin{aligned} {\varvec{F}}_{ij}=\partial _i {\varvec{A}}_j - \partial _i {\varvec{A}}_j + {\varvec{A}}_i\times {\varvec{A}}_j, \qquad i,j=1,2, \end{aligned}$$
(3.3)

we note

$$\begin{aligned} (D_1{\varvec{n}}+ {\varvec{n}}\times D_2{\varvec{n}})^2 = (D_1{\varvec{n}})^2 + (D_2{\varvec{n}})^2 - 2D_1{\varvec{n}}\times D_2{\varvec{n}}\cdot {\varvec{n}}. \end{aligned}$$
(3.4)

and also, as already observed by ’t Hooft [20],

$$\begin{aligned} {\varvec{n}}\cdot D_1{\varvec{n}}\times D_2{\varvec{n}}- {\varvec{n}}\cdot {\varvec{F}}_{12} = {\varvec{n}}\cdot \partial _1{\varvec{n}}\times \partial _2{\varvec{n}}+ \partial _2({\varvec{n}}\cdot {\varvec{A}}_1) - \partial _1({\varvec{n}}\cdot {\varvec{A}}_2). \end{aligned}$$
(3.5)

This equation shows that the particular combination of the degree density (the integrand of (2.12)) with a two-dimensional curl on the right hand side can be expressed in a manifestly gauge invariant way.

We can now state and prove the main result in this section.

Lemma 3.1

The energy for magnetic skyrmions at critical coupling associated with a compact subset \(D\subset {\mathbb {R}}^2\) can be written as

$$\begin{aligned} E_D[{\varvec{n}}] = 4\pi (Q_D[{\varvec{n}}] + \Omega _D[{\varvec{n}}])+ \int _{D} (D_1{\varvec{n}}+{\varvec{n}}\times D_2{\varvec{n}})^2 \, \mathrm {d}x_1\mathrm {d}x_2, \end{aligned}$$
(3.6)

where we used the covariant derivative

$$\begin{aligned} D_i{\varvec{n}}= \partial _i{\varvec{n}}-\kappa {\varvec{e}}_i^{-\alpha }\times {\varvec{n}}, \quad i=1,2, \end{aligned}$$
(3.7)

defined in terms of (2.1). In particular, the equality

$$\begin{aligned} E_D[{\varvec{n}}] = 4\pi (Q_D[{\varvec{n}}] +\Omega _D[{\varvec{n}}]) \end{aligned}$$
(3.8)

holds for all compact subsets \(D\subset {\mathbb {R}}\) iff the Bogomol’nyi equation

$$\begin{aligned} D_1{\varvec{n}}=-{\varvec{n}}\times D_2{\varvec{n}}\end{aligned}$$
(3.9)

is satisfied. This equation implies the variational equation (2.19).

Proof

Consider the gauge field given by the constant, Lie algebra-valued one-form with Cartesian components

$$\begin{aligned} {\varvec{A}}_i= -\kappa {\varvec{e}}_i^{-\alpha }, \qquad i=1,2. \end{aligned}$$
(3.10)

Then

$$\begin{aligned} {\varvec{F}}_{12}= \kappa ^2 {\varvec{e}}_3, \quad {\varvec{n}}\cdot {\varvec{A}}_i =-\kappa {\varvec{n}}_i^{\alpha }, \quad i=1,2, \end{aligned}$$
(3.11)

and therefore combining the results (3.4) and (3.5) for this gauge field gives

$$\begin{aligned} (D_1{\varvec{n}}+ {\varvec{n}}\times D_2{\varvec{n}})^2 {= } (D_1{\varvec{n}})^2 {+} (D_2{\varvec{n}})^2 {-} 2\left( {\varvec{n}}\cdot \partial _1{\varvec{n}}{\times } \partial _2{\varvec{n}}{+} \kappa (\partial _1 n_2^\alpha {-}\partial _2n_1^\alpha ) {+} \kappa ^2n_3\right) .\nonumber \\ \end{aligned}$$
(3.12)

One also checks that

$$\begin{aligned} \frac{1}{2}(D_1{\varvec{n}}^2 +D_2{\varvec{n}}^2)= \frac{1}{2}(\nabla {\varvec{n}})^2 + \kappa {\varvec{n}}\cdot \nabla _{-\alpha } \times {\varvec{n}}+\frac{1}{2} \kappa ^2(1+n_3^2), \end{aligned}$$
(3.13)

and so the energy density of (2.7) can be written as

$$\begin{aligned} \frac{1}{2}(\nabla {\varvec{n}})^2 +&\kappa {\varvec{n}}\cdot \nabla _{-\alpha } \times {\varvec{n}}+ \frac{\kappa ^2}{2}(1-n_3)^2 \nonumber \\&=\frac{1}{2} (D_1{\varvec{n}}+ {\varvec{n}}\times D_2{\varvec{n}})^2 +{\varvec{n}}\cdot \partial _1{\varvec{n}}\times \partial _2{\varvec{n}}+\kappa ( \partial _1n_2^\alpha -\partial _2n_1^\alpha ). \end{aligned}$$
(3.14)

Integrating and using the definitions (2.15), we deduce that the energy associated with a compact subset \(D\subset R^2\) can be written as claimed in (3.6). The equality (3.8) holds for all \(D\subset {\mathbb {R}}^2\) iff

$$\begin{aligned} D_1{\varvec{n}}=- {\varvec{n}}\times D_2{\varvec{n}}\Leftrightarrow D_2{\varvec{n}}= {\varvec{n}}\times D_1{\varvec{n}}, \end{aligned}$$
(3.15)

where the equivalence follows by applying \({\varvec{n}}\times \).

Showing that the Eq. (3.15) implies the variational equation is a lengthy but standard calculation. We indicate the main steps. Spelling out the Bogomol’nyi equation, we have

$$\begin{aligned} \partial _1{\varvec{n}}&= -{\varvec{n}}\times \partial _2 {\varvec{n}}+ \kappa ({\varvec{e}}_1^{-\alpha }\times {\varvec{n}}+ {\varvec{n}}\times ({\varvec{e}}^{-\alpha }_2\times {\varvec{n}})), \nonumber \\ \partial _2{\varvec{n}}&= {\varvec{n}}\times \partial _1 {\varvec{n}}+ \kappa ({\varvec{e}}^{-\alpha }_2\times {\varvec{n}}- {\varvec{n}}\times ({\varvec{e}}_1^{-\alpha }\times {\varvec{n}})). \end{aligned}$$
(3.16)

Therefore

$$\begin{aligned}&\partial _1^2 {\varvec{n}}+ \partial _2^2 {\varvec{n}}= 2 \partial _2{\varvec{n}}\times \partial _1{\varvec{n}}+\kappa ({\varvec{e}}_1^{-\alpha }\times \partial _1{\varvec{n}}+ {\varvec{e}}_2^{-\alpha }\times \partial _2{\varvec{n}}) \nonumber \\&+ \kappa (\partial _1{\varvec{n}}\times ({\varvec{e}}_2^{-\alpha }\times {\varvec{n}})+ {\varvec{n}}\times ({\varvec{e}}_2^{-\alpha }\times \partial _1 {\varvec{n}})- \partial _2{\varvec{n}}\times ({\varvec{e}}_1^{-\alpha }\times {\varvec{n}}))- {\varvec{n}}\times ({\varvec{e}}_1^{-\alpha } \times \partial _2{\varvec{n}})). \end{aligned}$$
(3.17)

Taking a cross product with \({\varvec{n}}\) and noting

$$\begin{aligned} {\varvec{n}}\times ({\varvec{e}}_i^{-\alpha }\times \partial _j{\varvec{n}}) = -n_i^\alpha \partial _j{\varvec{n}}, \end{aligned}$$
(3.18)

we arrive at

$$\begin{aligned} {\varvec{n}}\times \Delta {\varvec{n}}= -\kappa (n_1^\alpha \partial _1 + n_2^\alpha \partial _2) {\varvec{n}}+ \kappa ( n_1^\alpha {\varvec{n}}\times \partial _2{\varvec{n}}- n_2^\alpha {\varvec{n}}\times \partial _1{\varvec{n}}). \end{aligned}$$
(3.19)

Now we use the Bogomol’nyi equation again in the last term to conclude that

$$\begin{aligned} \kappa ( n_1^\alpha {\varvec{n}}\times \partial _2{\varvec{n}}- n_2^\alpha {\varvec{n}}\times \partial _1{\varvec{n}}) = -\kappa (n_1^\alpha \partial _1 + n_2^\alpha \partial _2) {\varvec{n}}+\kappa ^2(1-n_3){\varvec{e}}_3\times {\varvec{n}}. \end{aligned}$$
(3.20)

Therefore

$$\begin{aligned} {\varvec{n}}\times \Delta {\varvec{n}}= -2\kappa (n_1^\alpha \partial _1 + n_2^\alpha \partial _2) {\varvec{n}}+ \kappa ^2(1-n_3){\varvec{e}}_3\times {\varvec{n}}, \end{aligned}$$
(3.21)

which is the Eq. (2.19) obtained by variation. \(\square \)

As often in the O(3) sigma model or its gauged versions, the Bogomol’nyi equations are best studied in complex, stereographic coordinates. We do this in the next section.

4 Magnetic Skyrmions in Complex Coordinates

4.1 The Bogomol’nyi equation in stereographic coordinates

We use stereographic coordinates on the sphere defined by projection from the south pole. With the abbreviation

$$\begin{aligned} \nu = n_1+in_2, \end{aligned}$$
(4.1)

our stereographic coordinate is

$$\begin{aligned} w= \frac{\nu }{1+n_3}. \end{aligned}$$
(4.2)

When the magnetisation tends to the minimum of the potential term, \({\varvec{n}}\rightarrow (0,0,1)^t\), then \(w\rightarrow 0\). This makes w a natural choice of coordinate, but in describing our solutions we also need

$$\begin{aligned} v =\frac{1}{w}. \end{aligned}$$
(4.3)

For later use we also note the inverse relation

$$\begin{aligned} \nu = \frac{2w }{1+|w|^2} ,\quad n_3= \frac{1-|w|^2 }{1+|w|^2}. \end{aligned}$$
(4.4)

We introduce the complex coordinate \(z=x_1+ix_2\) in the plane, and use the standard holomorphic and anti-holomorphic derivatives

$$\begin{aligned} \partial _z= \frac{1}{2} (\partial _1-i\partial _2), \qquad \partial _{{\bar{z}}} = \frac{1}{2} (\partial _1+i\partial _2). \end{aligned}$$
(4.5)

Observing that, with the notation (2.3), \( e^{i\alpha } \nu = n_1^\alpha + in_2^\alpha \), the DM interaction term (2.4) can be written in stereographic coordinates as

$$\begin{aligned} \kappa {\varvec{n}}\cdot \nabla _{-\alpha } \times {\varvec{n}}= 2\kappa \text {Im}(e^{i\alpha }( n_3\partial _z \nu - \nu \partial _z n_3))= 4\kappa \text {Im} \left( e^{i\alpha } \frac{\partial _z w + w^2 \partial _z {\bar{w}}}{(1+|w|^2)^2}\right) . \end{aligned}$$
(4.6)

The other terms in the energy functional have standard expressions in stereographic coordinates, and so the energy (2.7) is

$$\begin{aligned} E[w]= \int _{{\mathbb {R}}^2} \frac{2|\nabla w|^2 + 4 \kappa \text {Im} (e^{i\alpha }(\partial _z w + w^2 \partial _z {\bar{w}})) + 2\kappa ^2| w|^4}{(1+|w|^2)^2}\, \mathrm {d}x_1 \mathrm {d}x_2. \end{aligned}$$
(4.7)

The integral (2.12) defining the degree is

$$\begin{aligned} Q[w]= \frac{ i}{2\pi } \int _{{\mathbb {R}}^2} \frac{\partial _1w \partial _2{\bar{w}} - \partial _2 w\partial _1{\bar{w}}}{(1+|w|^2)^2} \, \mathrm {d}x_1 dx_2, \end{aligned}$$
(4.8)

while the vorticity (2.14) is

$$\begin{aligned} \omega = 2\kappa \text {Im}(e^{i\alpha } \partial _z \nu ) = 4\kappa \text {Im}\left( e^{i\alpha } \frac{\partial _z w-w^2 \partial _z {\bar{w}}}{(1+|w|^2)^2}\right) . \end{aligned}$$
(4.9)

The one-form \(\Theta \) can be written as

$$\begin{aligned} \Theta = \kappa \text {Re} (e^{-i\alpha } {\bar{\nu }} \mathrm {d}z) =\kappa \frac{2\text {Re} (e^{-i\alpha } {\bar{w}} \mathrm {d}z)}{1+|w|^2} =\kappa \frac{2\text {Re} (e^{-i\alpha } v \mathrm {d}z)}{1+|v|^2}. \end{aligned}$$
(4.10)

In the following we write \(E_D[w], Q_D[w]\) and \(\Omega _D[w]\) for the integrals (2.15) over \(D\subset {\mathbb {R}}^2\) with the integrands expressed in terms of the complex field w. We can now state the main result of this paper.

Theorem 4.1

The energy (2.15) associated with a compact subset \(D\subset {\mathbb {R}}^2\) can be written as

$$\begin{aligned} E_D [w] = 4\pi (Q_D[w] +\Omega _D[w]) + \int _{D}8 \frac{( \partial _{{\bar{z}}} w - \frac{i}{2} \kappa e^{i\alpha } w^2 )(\partial _z \bar{w } + \frac{i}{2} \kappa e^{-i\alpha } {\bar{w}}^2 )}{(1+|w|^2)^2} \, \mathrm {d}x_1 \mathrm {d}x_2. \end{aligned}$$
(4.11)

The equality

$$\begin{aligned} E_D[w] = 4\pi (Q_D[w] +\Omega _D[w]) \end{aligned}$$
(4.12)

holds for all compact \(D\subset {\mathbb {R}}^2\) iff the field v defined in (4.3) satisfies the Bogomol’nyi equation

$$\begin{aligned} \partial _{{\bar{z}}} v = -\frac{i}{2} \kappa e^{i\alpha }. \end{aligned}$$
(4.13)

The general solution is

$$\begin{aligned} v = -\frac{i}{2} \kappa e^{i\alpha }{\bar{z}} + f(z), \end{aligned}$$
(4.14)

where f is an arbitrary holomorphic map from the plane to the Riemann sphere.

Proof

Using the standard identity

$$\begin{aligned} \partial _z \bar{ w} \partial _{{\bar{z}}} w = \frac{1}{4} \left( |\partial _1w|^2 + |\partial _2 w|^2 - i (\partial _1w \partial _2{\bar{w}} - \partial _2 w\partial _1{\bar{w}})\right) . \end{aligned}$$
(4.15)

and the expression for the vorticity (4.9) in complex coordinates, we have the following identities for the energy density

$$\begin{aligned}&\frac{2|\nabla w|^2 + 4 \kappa \text {Im} (e^{i\alpha }(\partial _z w + w^2 \partial _z {\bar{w}})) + 2\kappa ^2| w|^4}{(1+|w|^2)^2} \\&= \omega + \frac{2|\nabla w|^2 + 8 \kappa \text {Im} (e^{i\alpha } w^2 \partial _z {\bar{w}}) + 2\kappa ^2| w|^4}{(1+|w|^2)^2} \\&= \omega + 2i \frac{\partial _1w \partial _2{\bar{w}} - \partial _2 w\partial _1{\bar{w}}}{(1+|w|^2)^2} +\frac{8 \partial _z \bar{ w} \bar{\partial } w+ 8 \kappa \text {Im} (e^{i\alpha } w^2 \partial _z {\bar{w}}) + 2\kappa ^2| w|^4}{(1+|w|^2)^2} \\&=\omega + 2i \frac{\partial _1w \partial _2{\bar{w}} - \partial _2 w\partial _1{\bar{w}}}{(1+|w|^2)^2}+ 8 \frac{( \partial _{{\bar{z}}} w - \frac{i}{2} \kappa e^{i\alpha } w^2 )(\partial _z \bar{w } + \frac{i}{2} \kappa e^{-i\alpha } {\bar{w}}^2 )}{(1+|w|^2)^2} . \end{aligned}$$

Integrating and using the expression (4.8) yields the claimed expression (4.11) for the energy.

It follows immediately that the equality (4.12) holds for all compact \(D\subset {\mathbb {R}}^2\) iff

$$\begin{aligned} \partial _{{\bar{z}}} w =\frac{i}{2} \kappa e^{i\alpha } w^2, \end{aligned}$$
(4.16)

which is the Bogomol’nyi equation (3.9) in complex coordinates. With v as defined, this is equivalent to

$$\begin{aligned} \partial _{{\bar{z}}} v =-\frac{i}{2}\kappa e^{i\alpha }, \end{aligned}$$
(4.17)

whose general solution is \( v = -\frac{i}{2}\kappa e^{i\alpha } {\bar{z}} + f(z)\), where f is an arbitrary holomorphic function, as claimed. Since f takes values in the Riemann sphere, it is allowed to take the value \(\infty \). \(\square \)

One checks the equation (4.13) is equivalent to the Bogomol’nyi equation (3.9), and that it therefore implies the variational equation (2.19). For later use we note that the energy density for configurations which satisfy the Bogomol’nyi equation is the vorticity plus \(4\pi \) times the integrand of the degree. This sum can be written as

$$\begin{aligned} \epsilon (x_1,x_2)= 4 \frac{ \partial _z w\partial _{{\bar{z}}} {{\bar{w}}} - \partial _{{\bar{z}}} w \partial _z {\bar{w}} +\kappa \text {Im}\left( e^{i\alpha } (\partial _z w -w^2\partial _z{{\bar{w}}})\right) }{(1+|w|^2)^2}. \end{aligned}$$
(4.18)

4.2 Degree and vorticity of magnetic skyrmions at critical coupling

Before discussing the topology of the magnetic skyrmions defined by (4.14), it is worth revisiting the simpler case of the standard O(3) sigma model in the plane, defined by the Dirichlet energy functional [5, 13]. The requirement of finite energy in that model leads to the condition that fields tend to a constant at spatial infinity and may be extended to smooth maps \(S^2\rightarrow S^2\). The Bogomol’nyi equations are then equivalent to the map being either holomorphic or anti-holomorphic. Considering the holomorphic case for definiteness, the energy is proportional to the degree and for this to be finite, the configuration has to be a rational map, i.e., of the form p(z)/q(z), where p and q are polynomials of degree m and n. The topological degree of the map is simply max(m, n) in that case.

The DM term, which is a crucial feature of all models of magnetic skyrmions, is not positive definite, and therefore the energy expression for magnetic skyrmions may be finite even for configurations which do not tend to a constant value at spatial infinity. As a result, even finite energy configurations do not necessarily extend to smooth maps \(S^2\rightarrow S^2\) and do not necessarily have a well-defined topological degree. Moreover, it is not clear a priori if solutions of the Bogomol’nyi equation in our models have well-defined total vortex strength and total energy.

We shall now illustrate these issues for our infinite family of solutions (4.14), and show that, for rational holomorphic functions, the combination \(4\pi (Q+\Omega ^0) \) of degree and the regularised total vortex strength (2.18) nonetheless always yields a positive integer multiple of \(4\pi \).

Before we enter a general discussion, it is illuminating to consider linear examples of the form

$$\begin{aligned} v= -\frac{i}{2} \kappa e^{i\alpha } ({\bar{z}} + A e^{i\chi }z), \qquad A\in {\mathbb {R}}^{\ge 0}. \end{aligned}$$
(4.19)

As we shall see, this family captures the essence of the problems one encounters when defining the degree and the total vortex strength.

The evaluation of the integral defining Q is elementary. Switching to polar coordinates according to \(z=re^{i\varphi }\), we find, after completing the radial integration,

$$\begin{aligned} Q[w]= \frac{1}{2\pi } \int _0^{2\pi } \frac{A^2 - 1}{1+ A^2 +2A\cos (2\varphi +\chi )}\mathrm {d}\varphi . \end{aligned}$$
(4.20)

The evaluation of the total vortex strength is more subtle. The one-form \(\Theta \) for the field (4.19) is

$$\begin{aligned} \Theta = 4 \frac{ (1+ A\cos (2\varphi + \chi )) \mathrm {d}\varphi + A\sin (2\varphi +\chi ) \mathrm {d}\ln r}{r^{-2} + ( 1 + A^2 +2A \cos (2\varphi +\chi )}. \end{aligned}$$
(4.21)

When computing the total vortex strength, we need to integrate this form over a curve along which r is large. The leading term is

$$\begin{aligned} \Theta \sim 4 \frac{ (1+ A\cos (2\varphi + \chi )) }{ 1 + A^2 +2A \cos (2\varphi +\chi ) } \mathrm {d}\varphi +4\frac{A\sin (2\varphi +\chi )}{ 1 + A^2 +2A \cos (2\varphi +\chi )} \mathrm {d}\ln r. \end{aligned}$$
(4.22)

Clearly, the integral of the term proportional to \(\mathrm {d}\varphi \) gives the same answer for any simple curve enclosing the origin. However, the integral of the term proportional to \( A \, \mathrm {d}\ln r\) depends on the curve we choose, even in the limit of large radius. One can use the \(\varphi \)-dependence to introduce arbitrary contributions by deforming the contour with an outward bulge starting at some angle \(\varphi _1\) and ending at \(\varphi _2 > \varphi _1\). We conclude that the total vortex strength is not well-defined for configurations defined by (4.19) when \(A\ne 0\).

However, the integral of \(\Theta \) along a large circle \(C^R\) centred at the origin has a well-defined limit as the radius tends to infinity, precisely because \(\mathrm {d}\ln r\) does not contribute along such a circle. Thus, with the definition (2.18)

$$\begin{aligned} \Omega ^0[w] =\frac{1}{2\pi } \int _0^{2\pi }\frac{2+2 A\cos (2\varphi +\chi )}{1 + A^2 +2A \cos (2\varphi +\chi )} \mathrm {d}\varphi . \end{aligned}$$
(4.23)

It is immediate that

$$\begin{aligned} Q[w] + \Omega ^\circ [w] =1, \end{aligned}$$
(4.24)

regardless of the value of A. However, the contribution from the degree and the vortex strength depends crucially on A. Since

$$\begin{aligned} Q[w]= \frac{1}{2\pi } \int _0^{2\pi } \frac{A^2 - 1}{1+ A^2 + 2 A\cos (2\varphi +\chi )}\mathrm {d}\varphi ={\left\{ \begin{array}{ll} 1 \quad &{} \text {if} \quad A>1 \\ 0 \quad &{} \text {if} \quad A=1 \\ -1 \quad &{}\text {if} \quad A<1, \end{array}\right. } \end{aligned}$$
(4.25)

and

$$\begin{aligned} \Omega ^\circ [w]= \frac{1}{2\pi } \int _0^{2\pi } \frac{2+ 2A\cos (2\varphi +\chi )}{1 + A^2 +2 A\cos (2\varphi +\chi )}\mathrm {d}\varphi ={\left\{ \begin{array}{ll} 0 \quad &{} \text {if} \quad A>1 \\ 1\quad &{} \text {if} \quad A=1 \\ 2 \quad &{}\text {if} \quad A< 1, \end{array}\right. } \end{aligned}$$
(4.26)

we see that a configuration dominated by the holomorphic part (\(A>1\)) has degree 1 and vanishing vortex strength. A configuration dominated by the anti-holomorphic part (\(A<1\)) has degree -1 but vortex strength 2. In the intermediate case \(A=1\), the degree comes out as 0 and the vortex strength contributes 1.

The deeper reason behind the ‘jumping’ of the degree of the map defined by (4.19) lies in the extendibility of the map to one between spheres. An overall factor is irrelevant for this discussion, so we consider

$$\begin{aligned} v={\bar{z}} +Ae^{i\chi } z. \end{aligned}$$
(4.27)

Then \(w=1/v\) has a pole when \(re^{-2i\varphi }= -Are^{i\chi }\). This has no solution when \(A\ne 1\), but is solved by the entire line

$$\begin{aligned} \varphi = -\frac{\chi }{2} \pm \frac{\pi }{2} \end{aligned}$$
(4.28)

when \(A=1\). In particular, w therefore does not have a good limit for \(z\rightarrow \infty \) when \(A=1\): the result is infinity along the direction (4.28) but zero otherwise. It therefore does not extend to a smooth map between spheres in that case. When \(A\ne 1\) one checks, by considering the map in terms of \(\zeta =1/z\), that w does extend to a smooth map between spheres. Our integrations confirm this analysis for \(A\ne 1\), but also show that the combination of degree and vortex strength gives a stable result even when \(A=1\).

Our observations about the examples (4.19) generalise. In order to formulate this generalisation we define the regularised energy of a solution of the Bogomol’nyi equation as

$$\begin{aligned} E^\circ [w]= 4\pi (Q[w] + \Omega ^0[w]). \end{aligned}$$
(4.29)

Lemma 4.2

If p and q are polynomials in z of degree m and n and without common factor, the integral defining the total energy of the magnetic skyrmion solution determined via

$$\begin{aligned} v = -\frac{i}{2}\kappa e^{i\alpha }{\bar{z}} + \frac{p(z)}{q(z)}, \end{aligned}$$
(4.30)

is well-defined provided \(m\ne n+1\), i.e., provided p/q does not grow linearly for large z. The total energy is

$$\begin{aligned} E[w]=4\pi \, \text {max} (m,n+1)\qquad \text {if} \quad m\ne n+1. \end{aligned}$$
(4.31)

When \(m=n+1\), the total energy is not well-defined but the regularised total energy is

$$\begin{aligned} E^\circ [w]=4\pi \, m \qquad \text {if} \quad m=n+1. \end{aligned}$$
(4.32)

Proof

It is clear that

$$\begin{aligned} f(z)= \frac{p(z)}{q(z)} \end{aligned}$$
(4.33)

is a holomorphic map to the Riemann sphere, so (4.30) defines a magnetic skyrmion satisfying the Bogomol’nyi equation. Written in terms of v, the energy density (4.18) is

$$\begin{aligned} \epsilon (x_1,x_2)= 4 \frac{ \partial _z v\partial _{{\bar{z}}} {{\bar{v}}} - \partial _{{\bar{z}}} v \partial _z {\bar{v}} +\kappa \text {Im}\left( e^{i\alpha } (\partial _z {\bar{v}} -\bar{v}^2\partial _z v)\right) }{(1+|v|^2)^2}. \end{aligned}$$
(4.34)

This expression shows in particular that the energy density is smooth at the poles of w (zeros of v), so that any divergence in the energy integral must come from the behaviour at infinity. We first show that there is no divergence when \(m\ne n+1\).

For solutions of the form (4.30) and \(m-n> 1\), the leading term in the energy density for large r comes from

$$\begin{aligned} 4 \kappa \frac{ \text {Im}\left( -e^{i\alpha } \bar{v}^2\partial _z v\right) }{(1+|v|^2)^2}, \end{aligned}$$
(4.35)

leading to the asymptotic formula

$$\begin{aligned} |\epsilon (r,\varphi )|= C r^{-(m-n+1)} + {\mathcal {O}} \left( r^{-(m-n+2)}\right) , \quad \text {for some} \, C\in {\mathbb {R}}. \end{aligned}$$
(4.36)

This is integrable with respect to the integration measure \(r \mathrm {d}r \, \mathrm {d}\varphi \) for \(m-n> 1\).

When \(m -n < 1\), the leading large-r behaviour in the energy density is determined by the linear anti-holomorphic term in v, so the the energy density behaves asymptotically as

$$\begin{aligned} |\epsilon (r,\varphi )|= C r^{-4} + {\mathcal {O}} \left( r^{-5}\right) , \quad \text {for some} \, C\in {\mathbb {R}}. \end{aligned}$$
(4.37)

This is again integrable with respect to the integration measure \(r \mathrm {d}r \, \mathrm {d}\varphi \).

Next we turn to the evaluation of the energy integral. Our calculations will also show that, for \(m=n+1\), the integral requires regularisation. Our strategy is as follows. For any solution of the Bogomol’nyi equation we have \( E_D[w]=4\pi (Q_D[w] +\Omega _D[w])\) for any compact subset \(D\subset {\mathbb {R}}^2\). In order to evaluate the energy integral, we turn the integral defining the degree into boundary integrals and evaluate them together with the boundary integral defining the total vortex strength. In the cases where the total energy is well-defined, we will find that the boundary contribution from infinity is independent of the choice of contour at infinity. In the case where p/q grows linearly at infinity, we evaluate the boundary contribution on a circle at infinity, leading to our result for the regularised energy.

We recall that the expression (2.12) for the degree is the integral of the pull-back of the area form on \(S^2\), and that it can be written in terms of w and v as

$$\begin{aligned} 4\pi Q[w] =2i\int _{{\mathbb {R}}^2} \frac{\mathrm {d}w\wedge \mathrm {d}{\bar{w}} }{(1+|w|^2)^2} =2i\int _{{\mathbb {R}}^2} \frac{\mathrm {d}v\wedge \mathrm {d}{\bar{v}} }{(1+|v|^2)^2}. \end{aligned}$$
(4.38)

Moreover, the integrand can be written as an exact form

$$\begin{aligned} 2i\frac{\mathrm {d}v\wedge \mathrm {d}{\bar{v}} }{(1+|v|^2)^2} =\mathrm {d}\left( i\frac{v\mathrm {d}{\bar{v}} - {\bar{v}} \mathrm {d}v}{1+|v|^2} \right) , \end{aligned}$$
(4.39)

but the one-form in brackets is singular at the poles of v. We can make these singularities explicit as follows.

With f of the given form, we note

$$\begin{aligned} v=\frac{g}{q}, \quad \text {with} \quad g= -\frac{i}{2}\kappa e^{i\alpha }{\bar{z}} q(z) +p(z). \end{aligned}$$
(4.40)

For any v of this form, one checks that

$$\begin{aligned} i\frac{v \mathrm {d}{\bar{v}} - {\bar{v}} \mathrm {d}v}{1+|v|^2} = i\frac{g \mathrm {d}{\bar{g}} -{\bar{g}} \mathrm {d}g + q \mathrm {d}{\bar{q}} -{\bar{q}}\mathrm {d}q}{|g|^2 +|q|^2} +i( \mathrm {d}\ln q-\mathrm {d}\ln {\bar{q}}). \end{aligned}$$
(4.41)

The first term on the right hand side is manifestly smooth, but \( i( \mathrm {d}\ln q-\mathrm {d}\ln {\bar{q}}) \) is singular at the zeros of q. Thus, picking a compact region \(D\subset {\mathbb {R}}^2\) which contains open neighbourhoods of the zeros of q, and denoting negatively oriented circles of radius \(\epsilon \) around each of the zeros of q (possibly repeated) by \(C_j^\epsilon \), \(j=1,\ldots ,n\), we can write

$$\begin{aligned} 2i\int _{D} \frac{\mathrm {d}v\wedge \mathrm {d}{\bar{v}} }{(1+|v|^2)^2}&= i \lim _{\epsilon \rightarrow 0} \sum _{i=1}^n\int _{C^\epsilon _j}(\mathrm {d}\ln q-\mathrm {d}\ln {\bar{q}}) + i \int _{\partial D} \frac{v\mathrm {d}{\bar{v}} - {\bar{v}} \mathrm {d}v}{1+|v|^2} \nonumber \\&= 4\pi n + i \int _{\partial D} \frac{v\mathrm {d}{\bar{v}} - {\bar{v}} \mathrm {d}v}{1+|v|^2}. \end{aligned}$$
(4.42)

Therefore, the degree and the total vortex strength associated with the region D can be combined into

$$\begin{aligned} 4\pi (Q_D[w] + \Omega _D[w])= 4\pi n+ \int _{\partial D} \beta , \end{aligned}$$
(4.43)

where we introduced the one-form

$$\begin{aligned} \beta = \frac{iv\mathrm {d}{\bar{v}} - i {\bar{v}} \mathrm {d}v+ \kappa e^{-i\alpha } v \mathrm {d}z + \kappa e^{i\alpha } {\bar{v}} \mathrm {d}{\bar{z}} }{1+|v|^2}, \end{aligned}$$
(4.44)

which combines the one-form whose exterior derivative is the degree density (4.39) with the form \(\Theta \) (4.10) used in the definition of the total vortex strength.

Now, for solutions of the Bogomol’nyi equation,

$$\begin{aligned} \mathrm {d}v = -\frac{i}{2}\kappa e^{i\alpha } \mathrm {d}{\bar{z}} + \mathrm {d}f. \end{aligned}$$
(4.45)

It follows that

$$\begin{aligned} \beta = \frac{iv\mathrm {d}\bar{f} - i {\bar{v}} \mathrm {d}f+ \frac{\kappa }{2}e^{-i\alpha } v \mathrm {d}z + \frac{\kappa }{2} e^{i\alpha } {\bar{v}} \mathrm {d}{\bar{z}} }{1+|v|^2}. \end{aligned}$$
(4.46)

In order to evaluate the integral in (4.43) and its limit, we distinguish cases.

(i)   \(m > n+1\).  In this case the leading term in v for large r is \(az^{m-n}\) for some complex number a, and the leading term for f is also \(az^{m-n}\). Inserting these, we find

$$\begin{aligned} \beta \sim i(m-n) (\mathrm {d}\ln {\bar{z}} -\mathrm {d}\ln z) = 2(m-n)\mathrm {d}\varphi . \end{aligned}$$
(4.47)

The integral of the asymptotic form of \(\beta \) around any simple curve enclosing the origin is \(4\pi (m-n)\), and we conclude

$$\begin{aligned} 4\pi (Q[w] + \Omega [w]) = 4\pi (n + (m-n)) = 4\pi m \quad \text {if} \,\, m>n+1. \end{aligned}$$
(4.48)

(ii)   \(m = n+1\).  In this case we write the leading term in f as \(-\frac{i}{2}\kappa e^{i\alpha } Ae^{i\chi } z\) for some complex number \(Ae^{i\chi }\), and so that the leading terms in v for large r are \(-\frac{i}{2}\kappa e^{i\alpha }( {{\bar{z}}} + Ae^{i\chi } z)\). Then one checks, using essentially the calculation leading to (4.22), that the leading terms in \(\beta \) are

$$\begin{aligned} \beta \sim 2 \mathrm {d}\varphi +4\frac{A\sin (2\varphi +\chi )}{ 1 + A^2 +2A \cos (2\varphi +\chi )} \mathrm {d}\ln r. \end{aligned}$$
(4.49)

As already discussed in the paragraph following (4.22), the presence of the term proportional to \( \mathrm {d}\ln r\) means that, for \(A\ne 0\), the integral of \(\beta \) cannot be given a meaning independently of the curve, even in the limit of large radius. However, regularising by insisting on circular integration paths we observe

$$\begin{aligned} \lim _{R\rightarrow \infty } \int _{C^R}\beta = 4\pi , \end{aligned}$$
(4.50)

and hence

$$\begin{aligned} 4\pi (Q[w] + \Omega ^\circ [w]) = 4\pi (n+1) \quad \text {if} \, \, m=n+1. \end{aligned}$$
(4.51)

(iii)   \(m < n+1\).  Now the leading term in v is \(-\frac{i}{2}\kappa e^{i\alpha }{{\bar{z}}}\) and the holomorphic term is subleading for large r. The integration of \(\beta \) in the large r limit is therefore a special case of the calculation in (ii), obtained by setting \(A=0\). This eliminates the term proportional to \(\mathrm {d}\ln r\) in (4.49) and produces a limit independent of the chosen curve. We obtain

$$\begin{aligned} 4\pi (Q[w] + \Omega [w]) = 4\pi (n+1) \quad \text {if} \,\, m<n+1, \end{aligned}$$
(4.52)

which completes the proof. \(\square \)

For the remainder of the paper we focus on rational solutions of the form (4.30) and set

$$\begin{aligned} N= \text {max}(m,n+1). \end{aligned}$$
(4.53)

A simple counting argument shows that there is a 4N (real-)dimensional moduli space of rational maps of the form (4.33). For \(N=1\) (regularised energy \(4\pi \)), the 4-dimensional family of solutions is conveniently written as

$$\begin{aligned} v_1(z)=-\frac{i}{2}\kappa e^{i\alpha }\left( {\bar{z}} + az\right) +b \quad a,b \in {\mathbb {C}}. \end{aligned}$$
(4.54)

This includes the family (4.19) discussed in detail in the previous section for \(b=0\). For \(N=2\) (energy or regularised energy \(8\pi \)), the 8-dimensional family of solutions can be written as

$$\begin{aligned} v_2(z)=-\frac{i}{2}\kappa e^{i\alpha }{\bar{z}} + \frac{az^2 +bz+c}{dz+e} , \quad a,b,c,d,e \in {\mathbb {C}}, (a,b,c,d,e)\sim \lambda (a,b,c,d,e), \, \lambda \in {\mathbb {C}}^*, \end{aligned}$$
(4.55)

where the equivalence relation removes the redundant simultaneous rescaling of numerator and denominator by the same non-zero complex number, and we need to require that (i) a and d do not vanish simultaneously, (ii) d and e do not vanish simultaneously and (iii) the resultant of (a, b, c, d, e) is non-vanishing to ensure that numerator and denominator do no have common factors [5]. We will discuss the form and energy distribution of some of these solutions in the next section.

5 Solutions and Their Energy Density

The results obtained thus far add up to a simple recipe for constructing magnetisation fields \({\varvec{n}}\) which solve the Bogomol’nyi equation (and hence the variational equation (2.19)) out of two complex polynomials p and q. Inserting these polynomials into the expression (4.30) for the complex field v, and combining (4.3) and (4.4) to write the magnetisation \({\varvec{n}}\) as

$$\begin{aligned} n_1+in_2= \frac{2{\bar{v}} }{|v|^2+1}, \qquad n_3 = \frac{|v|^2-1}{|v|^2+1}, \end{aligned}$$
(5.1)

we obtain solutions of the Bogomol’nyi equation (3.9).

The simplest solutions are obtained when the holomorphic contribution to v vanishes, i.e., when \(f=0\). We use them to illustrate the translation from our coordinates into the ones conventionally used in the discussion of magnetic skyrmions in the literature. Translating \(v=-\frac{i}{2}\kappa e^{i\alpha }{\bar{z}}\) into the magnetisation field via (5.1), and comparing with the hedgehog parametrisation (2.20), we deduce

$$\begin{aligned} \theta = 2\tan ^{-1}\left( \frac{2}{\kappa r}\right) . \end{aligned}$$
(5.2)

This yields the Bloch skyrmion for \(\gamma =\frac{\pi }{2}\) (so \(\alpha =0\)) and the Néel skyrmion for \(\gamma = 0\) (so \(\alpha =\frac{\pi }{2}\)) in their standard form [2], but with a particularly simple profile function interpolating between \(\theta =\pi \) at \(r=0\) and \(\theta =0\) at \(r=\infty \). The solutions (5.2) agree with the hedgehog solutions (2.26) of the variational equations when \(\mu =\kappa \).

For the remainder of this section we set \( \kappa =1\), and study some example solutions of the Bogomol’nyi equation (3.9) in some detail. We adopt the convention, widely used in the magnetic skyrmion literature, to refer to configurations of negative degree (like the Bloch and Néel solutions) as skyrmions, and to the configurations of positive degree as anti-skyrmions.

We have organised our discussion according to the integer N defined in (4.53), and begin with the \(N=1\) family (4.54). Of the four real parameters in the two complex numbers a and b, three can be understood in terms of the symmetry group of translations and rotations (2.9). Rotations leave the basic skyrmion (\(a=b=0\)) invariant, but translations generate a shift in b. For the general configuration (4.54), rotations by \(\sigma \) act by mapping

$$\begin{aligned} -\frac{i}{2}\kappa e^{i\alpha }\left( {\bar{z}} + az\right) + b \mapsto -\frac{i}{2}\kappa e^{i\alpha } \left( {\bar{z}} + e^{-2i\sigma } a z\right) + e^{-i\sigma } b, \end{aligned}$$
(5.3)

so can be used to adjust the phase of a.

This action has an elementary but interesting consequence when \(b=0\). Since a rotation by some angle \(\sigma \) leaves the anti-holomorphic term invariant but rotates the phase of the linear holomorphic term by \(-2\sigma \), configurations with \(a\ne 0 \) and \(b=0\) are mapped to themselves after a rotation by \(\pi \). This clearly generalises to configurations with a homogeneous holomorphic part proportional to \(z^n\), which are mapped to themselves after a rotation by \(2\pi /n\).

Fig. 1
figure 1

Top from left to right: Bloch skyrmion \(v=-\frac{i}{2}{\bar{z}}\) in the theory with \(\alpha =0\) and Néel skyrmion \(v=\frac{1}{2}{\bar{z}}\) in the theory with \(\alpha =\frac{\pi }{2}\). Bottom: anti-skyrmions in the theory with \(\alpha =0\), with \(v=-\frac{i}{2}\left( {\bar{z}}+ 4z\right) \) shown on the left and \(v=-\frac{i}{2}\left( {\bar{z}}+6iz \right) \) on the right. Note that the magnetisation of anti-skyrmions rotates oppositely to that of the skyrmions when one traverses a positively oriented circle around the centre, where \(v=0\) or \({\varvec{n}}=(0,0,-1)^t\). Note also the different scale of the anti-skyrmion plots. Size and orientation of anti-skyrmions can be adjusted via the coefficient of z.

Full size image

The only parameter in the \(N=1\) family (4.54) which cannot be adjusted by a symmetry transformation is the magnitude |a| of the complex coordinate a. Varying it leads to the most interesting deformation in this family. We already know that for \(a=0\) we have the basic hedgehog skyrmion. According to our formula (4.25), we have degree \( -1\) configurations, i.e. skyrmions, for \(|a|<1\) and degree 1 configurations, i.e. anti-skyrmions, for \(|a| > 1 \). The interpolation between the two necessarily involves the case \(|a|= 1 \) where the degree is zero. As observed in the discussion after (4.25), the stereographic coordinate w has a pole along an entire line in this case. The magnetisation takes the value \({\varvec{n}}=(0,0,-1)^t\) and both the potential \(\frac{1}{2} (1-n_3)^2\) and the energy density (4.34) are maximal along this line. We therefore call \(N=1\) configurations with \(|a|=1\) line defects.

To sum up, as |a| increases from zero we deform a hedgehog skyrmion through a line defect into an anti-skyrmion. In this process, the axisymmetric energy distribution of a hedgehog skyrmion is first squeezed and stretched into an elliptical shape and then into a line at \(|a|=1\). As |a| is increased further the energy distribution contracts again into an elliptical shape and shrinks; it never recovers the axisymmetry of the original skyrmion. In Fig. 1 we illustrate our discussion by showing the magnetisation of the basic Bloch and Néel skyrmions in our model, and also the magnetisation of two anti-skyrmions in the theory with \(\alpha =0\).

Next we turn our attention to the family of solutions (4.55) with \(N=2\). We have not fully explored the eight parameters in this family, of which three can again be accounted for by the symmetry operations of translation and rotation. Here, we only exhibit two interesting phenomena and fix \(\alpha =0\) for definiteness. The first is the nonlinear superposition of skyrmions and anti-skyrmions in this model. Configurations of the form

$$\begin{aligned} v=-\frac{i}{2}\left( {\bar{z}}- \frac{ z^{n}}{R^{n-1}}\right) , \qquad R\in {\mathbb {R}}^{>0}, \quad n \in {\mathbb {Z}}^{>1} , \end{aligned}$$
(5.4)

have degree n, but the functions v have \((n+2)\) zeros: one at the origin and \((n+1)\) zeros at

$$\begin{aligned} z_k=R e^{\frac{2\pi k i}{n+1}} , \quad k=0,\ldots , n. \end{aligned}$$
(5.5)

The magnetisation takes the value \({\varvec{n}}=(0,0,-1)^t\) at the zeros, but inspection of the winding shows that the configuration consists of a skyrmion at the origin surrounded by \((n+1)\) anti-skyrmions at the locations (5.5). The energy density (4.34) is peaked at the anti-skyrmion locations. Such superpositions of solitons and anti-solitons do not solve the Bogomol’nyi equations of the pure O(3) sigma model, which require maps to be either holomorphic or antiholomorphic. In our model they are possible, essentially because the number of zeros for functions of the form (4.30) can exceed the absolute value of the degree. In fact, determining the number of zeros of (4.30) (and hence the number of skyrmions and anti-skyrmions counted without sign) is an interesting mathematical problem, see [21] for rigorous bounds and also [22, 23] for further results and the application of this problem to gravitational lensing.

Continuing with \(\alpha =0\), we observe that rational solutions of the form

$$\begin{aligned} v=-\frac{i}{2}\left( {\bar{z}} - \frac{R^2}{z}\right) , \quad R\in {\mathbb {R}}^{>0}, \end{aligned}$$
(5.6)

look like skyrmion bags [11] or sacks [12]. These particular bags have degree \(Q=0\), take the vacuum value \({\varvec{n}}=(0,0,1)^t\) at the origin and the value \({\varvec{n}}=(0,0,-1)^t\) on a circle of radius R centred at the origin. The energy density (4.34) is maximal on this circle. The solutions (5.6) are axisymmetric and an example of what is sometimes called skyrmionium in the literature. Clearly one can obtain more general bags by acting with translations on (5.6). We note that the circles of zeros of functions like (5.6) also play a special role in the context of gravitational lensing, where they are called Einstein rings [23].

6 Conclusion

In this paper we introduced models for magnetic skyrmions in the plane for which an infinite family of analytical solutions can be given explicitly. In our study we concentrated on the family of magnetic skyrmions determined by a rational holomorphic function according to (4.30). We showed that, with a suitable regularisation in the case of linear growth in the holomorphic function at infinity, the total energy takes the quantised values \(4\pi N\), where N is a positive integer which combines the degree with the (possibly regularised) total vorticity of a configuration. This integer does not appear to have been studied in the literature on magnetic skyrmions, but our study suggests that it plays an important role.

We also determined the collective coordinates or moduli of the solutions for given N, and studied example configurations for low values of N. They include remarkable deformations of Bloch and Néel skyrmions into line defects and anti-skyrmions. Finally, we exhibited some of the bag and multi-anti-skyrmion configurations included in the family (4.30).

Magnetic skyrmions are sometimes also called chiral skyrmions because the DM interaction breaks reflection symmetry. This is evident in our solutions through the mandatory and fixed anti-holomorphic part (of negative degree) but the optional holomorphic part (of positive degree). It is interesting and somewhat unexpected that our models allows for nonlinear superpositions of skyrmions and anti-skyrmions in static configurations. However, their roles are not symmetric. To obtain perfect mirrors of our models one would need to replace the DM interaction with the one of opposite chirality (2.11). One checks that this would lead to a Bogomol’nyi equation which enforces a fixed holomorphic part, but allows for an arbitrary anti-holomorphic part. This is consistent with the criterion derived in [16] for the preference of skyrmions over anti-skyrmions depending on the chirality.

The explicit family of solutions in our models and their chiral twins should be studied further in order to obtain a systematic understanding of the types of defects they capture. One would also like to understand how these defects relate to those observed in the various phases of generic models for magnetic skyrmions as discussed for example in [8].

Mathematically, one would like to understand more precisely the class of maps from the plane to the sphere for which the integrals defining the degree, total energy and total vortex strength are well-defined. It would also be important to ascertain if the energy functional (2.7) is bounded below for a suitably defined class of configurations (not just our solutions). In [19] it was shown that in a closely related model with a pure Zeeman potential the total energy is bounded below by a multiple of the absolute value of the degree. In our model, the energy is potentially unbounded below unless one imposes suitable behaviour at spatial infinity.

The second order variational equation (2.19) and the Bogomol’nyi equation (4.13) deserve further study. One would like to know, for example, if there are other finite-energy solutions (possibly after suitable regularisation) of the variational equation, not included in our rational family (4.30).

Finally, future work should include the study of time evolution and the effect of external fields on the solutions in our model.

Note added   While this paper was under review, research was reported in the literature which has implications for the definition of energy in our models. A modification of the energy expression for magnetic skyrmions by a boundary term already proposed for analytical reasons in [19] was generalised in [24] in the framework of gauged sigma models. Applied to the models discussed here, this modification amounts to subtracting the total vorticity \(\Omega _{{\mathbb {R}}^2}\) (2.15) from our energy (2.7) or, equivalently, to replacing \({\varvec{n}}\cdot \nabla _{-\alpha } \times {\varvec{n}}\) by \(({\varvec{n}}-{\varvec{e}}_3) \cdot \nabla _{-\alpha } \times {\varvec{n}}\) in (2.7). The modification does not affect the variational equations or the Bogomol’nyi equations studied here. The modified energy has the advantage of being finite without regularisation for all rational solutions of the form (4.30). In fact, for any such solution the modified energy is equal to \(4\pi Q\), where Q is the degree, instead of the value \(4\pi N\) for our energy (after regularisation if required). On the other hand, the geometrical investigation reported in [25] shows that the integer N (4.53) can be interpreted as the equivariant degree of the rational solutions (4.30), suggesting that the energy we used here is geometrically natural. It thus appears that there is a certain tension between the energy expressions preferred from an analytical and a geometrical point of view. Further work is required to resolve this.