Stability and asymptotic interactions of chiral magnetic skyrmions in a tilted
magnetic field

Bruno Barton-Singer 1$\star$ and Bernd J. Schroers2

1 Maxwell Institute of Mathematical Sciences and Department of Mathematics,
Heriot-Watt University, Edinburgh EH14 4AS, United Kingdom<EMAIL_ADDRESS>
2 Maxwell Institute of Mathematical Sciences and School of Mathematics,
University of Edinburgh, Edinburgh EH9 3FD, United Kingdom

## Abstract

Using a general framework, interaction potentials between chiral magnetic
solitons in a planar system with a tilted external magnetic field are
calculated analytically in the limit of large separation. The results are
compared to previous numerical results for solitons with topological charge
$\pm 1$. A key feature of the calculation is the interpretation of
Dzyaloshinskii-Moriya interaction (DMI) as a background $SO(3)$ gauge field.
In a tilted field, this leads to a $U(1)$-gauged version of the usual equation
for spin excitations, leading to a distinctive oscillating interaction
profile. We also obtain predictions for skyrmion stability in a tilted field
which closely match numerical observations.

###### Contents

  1. 1 Introduction
  2. 2 Linearised theory of chiral magnets
    1. 2.1 Energy functional of the chiral magnet
    2. 2.2 Asymptotic form of the chiral magnet energy functional
    3. 2.3 Soliton tails in the form of a multipole expansion
    4. 2.4 Soliton tails in terms of an effective source
    5. 2.5 Numerical calculation of multipole moments
    6. 2.6 Stability of magnetic solitons in a tilted magnetic field
  3. 3 Interaction potential for magnetic solitons
    1. 3.1 The interaction potential in terms of soliton tails
    2. 3.2 The interaction potential in terms of effective sources
    3. 3.3 Comparison to numerically calculated interaction potential in a tilted field
    4. 3.4 Soliton interactions in axisymmetric potential
  4. 4 Conclusion
  5. A Numerical method
  6. B Interaction potential for general energy functional

##  1 Introduction

Magnetic skyrmions are examples of topological solitons, topologically non-
trivial field configurations of finite energy that minimise an energy
functional [1, 2]. In general, these configurations are localised in a
particular region, and behave like particles in the sense that the lowest-
energy excitations of the system involve translating or rotating the localised
fields. It is therefore interesting to look at the effective dynamics of these
emergent particles, whether that is in the context of second-order Lorentz-
invariant dynamics, first-order gradient flow dynamics or first-order
Schrödinger dynamics, as it is an important part of the low-energy dynamics of
the system.

One may try to understand the dynamics of solitons by reducing the infinite-
dimensional dynamics of the field theory to a finite number of degrees of
freedom [3], often the positions and orientations of the solitons, and
calculating the energy’s dependence on these parameters. The assumption is
that the dynamics on this finite-dimensional space given by this restricted
energy closely approximate the full dynamics on the infinite-dimensional
space. Depending on the context, this is known as the moduli space or
collective co-ordinate approach. This assumption is also the essence of the
Thiele equation [4], where it is assumed that motion of a soliton can be
approximated by rigid translation of the statically stable configuration.
Generalisations of the Thiele method introduce a finite number of extra
parameters [5, 6] or allow the shape of the soliton to change depending on
position [7].

The accuracy of this approximation has been assessed rigorously for some
specific models. Some examples with both first-order Schrödinger dynamics and
second-order Lorentz-invariant dynamics are reviewed in [8]. In particular,
the moduli space approximation is shown to be justified in certain limits for
Chern-Simons vortices evolving according to Schrödinger dynamics [9], a
situation which is in some ways analogous to magnetic skyrmion dynamics.

Magnetic skyrmions are topological solitons in a magnetisation field for a
specific energy functional modelling chiral magnets [10], and the
magnetisation field evolves, in the simplest case, according to the Landau-
Lifshitz(-Gilbert) equation [11, 12]. Mathematically, this is a combination of
gradient flow dynamics, where the field changes so as to decrease the energy
functional as fast as possible, and Schrödinger dynamics, where the field
evolves in a way that preserves the energy functional. The Schrödinger
dynamics comes from the quantum-mechanical evolution of the magnetisation on a
microscopic level, which preserves energy, while the gradient flow component
represents damping forces that cause energy to leave the system.

An interaction potential is an example of energy restricted to a moduli space
in the specific case of two solitons. Conventionally, we subtract the energies
of the isolated solitons in order to separate out the potential arising from
interaction. In this paper we calculate this interaction potential in the
asymptotic limit where the distance between the two solitons becomes large.
There is a considerable literature on interaction potentials for solitons. In
many cases - for example in the study of nonabelian monopoles [13], nuclear
Skyrmions [14, 15, 16], baby Skyrmions [17] and abelian vortices [18] \- the
long-range asymptotics of the interaction potential can be derived in terms of
point sources interacting in a much simpler linear field theory. Such linear
approximations can be justified at various levels of rigor in each of those
models, but there is no general understanding of when and why linear point-
particle pictures capture the asymptotics of soliton interactions.

The question of estimating the interaction energy of magnetic skyrmions
analytically was first considered in [19]. More recently, [20, 21] calculated
the interaction energy in skyrmions supported by frustration and
Dzyaloshinskii-Moriya interaction (DMI) respectively. The latter of these drew
on methods used, for example, in [17, 18]. The interaction potential of
magnetic skyrmions in a tilted magnetic field has been numerically simulated
[22] and partially calculated [23], but an expression for the analytical
potential in terms of the separation of the skyrmions has not been given. The
advantages of an explicit formula are several: firstly, it gives new
understanding of the interaction potential in terms of multipoles in a $U(1)$
gauge theory which is, in principle, applicable to solitons of any degree.
Secondly, it gives an understanding of how the interaction potential depends
on the DMI, the external field tilt and the soliton separation, and thirdly it
gives us rigorous results, for example whether the interaction energy leads to
attraction between solitons for large separations.

The analytical calculation in this paper is based on the following general
observations. The field far from the core of any soliton (the ‘tail’) can be
approximated by the solution of the linearised Euler-Lagrange equations, in
the presence of a combination of multipole sources (the ‘effective source’)
positioned within the core of the soliton. It follows that to leading order in
inverse separation the effective source for the field far away from two well-
separated solitons is the sum of the two individual effective sources; this is
equivalent to saying that their soliton tails are approximately linearly
superposed. The leading contribution to the interaction between two solitons
can be expressed in terms of this far field so that, given certain assumptions
which are satisfied here, we can approximate the interaction potential in
terms of the individual soliton tails (see Appendix B). A surprising feature
of the interaction energy of solitons in general is that when it is fully
calculated in terms of the effective sources, we often find that it looks like
the interaction energy between the effective sources [17, 24]: in other words,
solitons do not only act like multipoles as sources of their tails, but also
respond like those same multipoles to the tail of the other soliton. Here we
extend this observation to chiral magnetic solitons in a tilted external
field, with the modification that the DMI in general introduces an effective
background $U(1)$ gauge field that interacts with the tails, and thus adds a
modulation on top of the familiar multipole-multipole interaction energy.

This picture explains the characteristic oscillating behaviour of the
interaction energy of magnetic solitons as a function of their separation.
What is more, because the interaction formula is derived generally, we can use
it to look at other cases of interest. We discuss the interaction of novel
textures that have recently been numerically observed in magnetic field
applied normal to the plane [25, 21, 26].

The paper is split into two main sections. In Section 2, we investigate the
tails of an individual soliton in a chiral magnet. We first review the
interpretation of DMI as a gauge field in Section 2.1, then discuss how this
gauge picture changes as magnetisation fields approach the background
magnetisation far from the soliton in Section 2.2. We solve the resulting
linearised Euler-Lagrange equation in general in Section 2.3, then see how the
soliton can be seen as an ‘effective source’ for its tail in Section 2.4. We
then test the accuracy of the linearised Euler-Lagrange equation against
numerics in Section 2.5 with details of the numerical methods given in
Appendix A. Finally, we note that the behaviour of solutions to the linearised
Euler-Lagrange equation gives new insight into elliptical instability of
magnetic skyrmions and antiskyrmions in Section 2.6, and compare this to
another novel calculation of elliptical instability based on the energy of an
isolated $2\pi$-domain wall.

In Section 3.1, we define the interaction potential and describe how it can be
approximated in terms of the tails of the individual solitons, under the
assumption that these tails fall off exponentially and that the solitons
perturb each other a small amount that goes to zero as they go to infinite
separation. The details of this calculation are covered in Appendix B. We then
substitute in the solutions we have from 2.4 to write the interaction energy
in terms of the effective sources of the two solitons in Section 3.2. We
consider two cases: in section 3.3, we consider the tilted-field case, where
we can compare to numerical observation [22]. This involves the numerical
results from Section 2.5. The numerical process is less expensive and more
accurate than directly calculating interaction energies, since we only need to
simulate an isolated soliton. In Section 3.4, we consider when the applied
field is not tilted, simplifying the potential but allowing for a greater
diversity of solitons whose interactions have not yet been considered. We
derive features of the interaction potential without any numerics. In
particular we argue that the presence of ‘chiral kinks’ [26] on either soliton
leads to an attractive force between the solitons.

##  2 Linearised theory of chiral magnets

###  2.1 Energy functional of the chiral magnet

We consider solitons in the magnetisation field $\boldsymbol{n}(\vec{x})$ in
the plane, satisfying the constraint
$\boldsymbol{n}(\vec{x})\cdot\boldsymbol{n}(\vec{x})=1$. We adopt the
convention that boldface vectors $\boldsymbol{v}$ always have three
components, and vectors with an arrow $\vec{v}$ always have two. Throughout
this paper we are concerned with the chiral magnet energy functional, which we
write as

$E(\boldsymbol{n})=\int\left(\frac{1}{2}\partial_{i}\boldsymbol{n}\cdot\partial_{i}\boldsymbol{n}+\boldsymbol{D}_{i}\cdot(\boldsymbol{n}\times\partial_{i}\boldsymbol{n})+V(\boldsymbol{n})\right)d^{2}x,$
(1)

where the term
$\boldsymbol{D}_{i}\cdot(\boldsymbol{n}\times\partial_{i}\boldsymbol{n})$ is
the generalized DMI [27, 28] introduced in different notation in [29] namely
$\boldsymbol{D}_{i}=-\frac{1}{J}\hat{D}\boldsymbol{e}_{i}$ in terms of the
standard orthonormal basis
$\boldsymbol{e}_{1},\boldsymbol{e}_{2},\boldsymbol{e}_{3}$ of $\mathbb{R}^{3}$
and the matrix $\hat{D}$, called the DM tensor. Summation over repeated
indices with $i=1,2$ is assumed throughout. The most commonly considered cases
of Bloch-type and Néel-type DMI are obtained by taking $\hat{D}$ to be the
identity and a 90 degree rotation about the 3-axis respectively.

We can equally view the chiral magnet energy, or rather a whole family of
chiral magnet energies, as an $SU(2)$ or equivalently $SO(3)$-gauged sigma
model [30, 31, 32]:

$E(\boldsymbol{n})=\int\left(\frac{1}{2}\lvert\partial_{i}\boldsymbol{n}+\boldsymbol{A}_{i}\times\boldsymbol{n}\rvert^{2}+\bar{V}(\boldsymbol{n})\right)d^{2}x,$
(2)

For later use we note that the gauge transformations are given by a spatially
dependent rotation matrix, in axis-angle co-ordinates
$R(\theta(\vec{x}),\boldsymbol{e}(\vec{x}))$, under which the fields and DMI
transform as follows:

$\displaystyle\boldsymbol{n}$
$\displaystyle\mapsto\tilde{\boldsymbol{n}}=R(\theta,\boldsymbol{e})\boldsymbol{n}$
$\displaystyle\boldsymbol{A}_{i}$
$\displaystyle\mapsto\tilde{\boldsymbol{A}_{i}}=R(\theta,\boldsymbol{e})\boldsymbol{A}_{i}-\partial_{i}\theta\boldsymbol{e}-\sin\theta\partial_{i}\boldsymbol{e}-(1-\cos\theta)\boldsymbol{e}\times\partial_{i}\boldsymbol{e}.$
(3)

We make some choices in formulating this gauge theory: in [31, 32], the
potential term is $-\boldsymbol{F}\cdot\boldsymbol{n}$ where $\boldsymbol{F}$
is the curvature of the gauge field, and thus explicitly gauge-invariant. Here
we are more general and let the vector parameters that appear in
$\bar{V}(\boldsymbol{n})$ also rotate under $R(\theta,\boldsymbol{e})$. We can
see that if we expand the energy functional out for a given value of the
fields $\boldsymbol{n}$, $\boldsymbol{A}_{i}$ then we recover a particular
chiral magnet model, with

$\displaystyle\boldsymbol{D}_{i}$ $\displaystyle=\boldsymbol{A}_{i}$
$\displaystyle V(\boldsymbol{n})$
$\displaystyle=\bar{V}(\boldsymbol{n})+\frac{1}{2}\lvert\boldsymbol{A}_{i}\times\boldsymbol{n}\rvert^{2}.$
(4)

Therefore different configurations of $\boldsymbol{A}_{i}$ correspond to
different material parameters $\boldsymbol{D}_{i}$ and $V(\boldsymbol{n})$,
and gauge transformations link a configuration $\boldsymbol{n}$ in a system
with DMI parameters $\boldsymbol{D}_{i}=\boldsymbol{A}_{i}$ to a different
configuration $\tilde{\boldsymbol{n}}$ in a different system with DMI
parameters $\boldsymbol{D}_{i}=\tilde{\boldsymbol{A}}_{i}$ and a different
(possibly spatially varying) potential. Although this differs from the usual
interpretation of gauge transformations as relating physically equivalent
configurations of fields, the language and technology of gauge theory has
proved useful in finding exact chiral skyrmion solutions [31], and we will see
here that it also provides a conceptually simple way of understanding
asymptotic skyrmion interactions.

Configurations of the magnetisation field are classified by their degree,
which under suitable assumptions on $\boldsymbol{n}$ is an integer:

$Q(\boldsymbol{n})=\frac{1}{4\pi}\int\boldsymbol{n}\cdot(\partial_{1}\boldsymbol{n}\times\partial_{2}\boldsymbol{n})d^{2}x.$
(5)

We call a non-trivial configuration that minimizes the energy (1) for a given
degree a magnetic soliton. Degree $-1$ magnetic solitons are called skyrmions,
while degree $+1$ magnetic solitons are called antiskyrmions.

The first term in the energy is the Heisenberg interaction, which favours
alignment of $\boldsymbol{n}(\vec{x})$ at a single constant value. Its
prefactor can be fixed to $\frac{1}{2}$ by picking appropriate units of
energy.

The second term is the DMI. We are most interested in the case where

$\boldsymbol{D}_{1}=-k\begin{pmatrix}\cos\beta\\\ \sin\beta\\\
0\end{pmatrix},\quad\boldsymbol{D}_{2}=-k\begin{pmatrix}-\sin\beta\\\
\phantom{-}\cos\beta\\\ 0\end{pmatrix},$ (6)

for real parameters $k$ and $\beta$. We call this axisymmetric DMI, because it
is invariant under the physical rotation:

$\boldsymbol{n}(\vec{x})\mapsto
R(\theta,\boldsymbol{e}_{3})\boldsymbol{n}(R(-\theta,\boldsymbol{e}_{3})\vec{x}).$
(7)

Then $\beta=0$ corresponds to the normal Bloch-type DMI,
$k\boldsymbol{n}\cdot(\nabla\times\boldsymbol{n})$, while
$\beta=\frac{\pi}{2}$ corresponds to Néel-type DMI. A general angle $\beta$
corresponds to a linear combination of both kinds of DMI.

The third term $V(\boldsymbol{n})$ is the potential function. It attains its
minimal value on a set of vacuum configurations. We choose one of them, or
take the unique vacuum if the set has only one element, call it
$\boldsymbol{n}_{0}$ in the following and impose it as the boundary value at
spatial infinity. Asymptotically we can therefore approximate
$\boldsymbol{n}(\vec{x})$ in terms of tangent vector fields to
$\boldsymbol{n}_{0}$. To do this we use the exponential map
$\exp_{\boldsymbol{n}(\vec{x})}:T_{\boldsymbol{n}}S^{2}\to S^{2}$ which takes
a tangent vector $\boldsymbol{\epsilon}(\vec{x})$ to the sphere at
$\boldsymbol{n}(\vec{x})$ (so
$\boldsymbol{\epsilon}(\vec{x})\cdot\boldsymbol{n}(\vec{x})=0$) to the point
on the sphere along the great circle in the direction of
$\boldsymbol{\epsilon}(\vec{x})$, at a distance
$\lvert\boldsymbol{\epsilon}(\vec{x})\rvert$:

$\exp_{\boldsymbol{n}(\vec{x})}(\boldsymbol{\epsilon}(\vec{x}))=\boldsymbol{n}(\vec{x})\cos\lvert\boldsymbol{\epsilon}(\vec{x})\rvert+\sin\lvert\boldsymbol{\epsilon}(\vec{x})\rvert\frac{\boldsymbol{\epsilon}(\vec{x})}{\lvert\boldsymbol{\epsilon}(\vec{x})\rvert}.$
(8)

The map allows us to describe fields $\boldsymbol{n}(\vec{x})$ near
$\boldsymbol{n}_{0}$ in terms of a linear tangent vector field
$\boldsymbol{\psi}_{\boldsymbol{n}}$ to $\boldsymbol{n}_{0}$, defined via

$\boldsymbol{\psi}_{\boldsymbol{n}}(\vec{x})=\exp^{-1}_{\boldsymbol{n}_{0}}(\boldsymbol{n}(\vec{x})),$
(9)

and this is essential for our formulation of the linearised theory. The
tangent plane $T_{\boldsymbol{n}_{0}}S^{2}$ is two-dimensional, so that
$\boldsymbol{\psi}_{\boldsymbol{n}}$ has two coordinates
$(\psi_{\boldsymbol{n}})_{1}$, $(\psi_{\boldsymbol{n}})_{2}$ with reference to
an orthonormal frame $\boldsymbol{e}^{n}_{1},\boldsymbol{e}^{n}_{2}$ oriented
so that
$\boldsymbol{e}^{n}_{1}\times\boldsymbol{e}^{n}_{2}=\boldsymbol{n}_{0}$.

We assume that $V(\boldsymbol{n})$, expanded around $\boldsymbol{n}_{0}$, is
rotationally symmetric around $\boldsymbol{n}_{0}$ to quadratic order.
Defining
$\tilde{V}(\boldsymbol{\psi}_{\boldsymbol{n}})=V(\exp_{\boldsymbol{n}_{0}}(\boldsymbol{\psi}_{\boldsymbol{n}}))$
this amounts to assuming that

$\frac{\partial^{2}\tilde{V}}{\partial(\psi_{\boldsymbol{n}})_{i}\partial(\psi_{\boldsymbol{n}})_{j}}\bigg{\rvert}_{\boldsymbol{n}_{0}}=\mu^{2}\delta_{ij},\quad\mu^{2}\geq
0,$ (10)

i.e. it costs equally to perturb in any direction away from the vacuum. Note
that this expression is independent of our choice of basis
$\boldsymbol{e}^{n}_{1},\boldsymbol{e}^{n}_{2}$.

One physically relevant case where (10) holds is when the potential is the
typical combination of a Zeeman interaction from a magnetic field applied
normal to the plane and an anisotropy term:

$V(\boldsymbol{n})=h_{z}(1-n_{3})+h_{a}(1-n_{3}^{2}),$ (11)

with $h_{z}+2h_{a}\geq 0$. Then $\boldsymbol{n}_{0}=\boldsymbol{e}_{3}$ and
$\mu^{2}=h_{z}+2h_{a}$. Both $V(\boldsymbol{n})$ and $\boldsymbol{n}_{0}$ are
symmetric under rotations around $\boldsymbol{e}_{3}$, so we call such a
potential axisymmetric. Another case of interest where (10) holds is when we
have just a Zeeman interaction, but from a tilted magnetic field with
direction $\boldsymbol{e}_{h}$:

$V(\boldsymbol{n})=h_{z}(1-\boldsymbol{e}_{h}\cdot\boldsymbol{n}).$ (12)

We take $h_{z}>0$. Then $\boldsymbol{n}_{0}=\boldsymbol{e}_{h}$ and
$\mu^{2}=h_{z}$. Specifically and without loss of generality in the case of
axisymmetric DMI (6), we consider a magnetic field tilted at angle
$\theta_{h}$ to $\boldsymbol{e}_{3}$ so that
$\boldsymbol{e}_{h}=\cos\theta_{h}\boldsymbol{e}_{3}+\sin\theta_{h}\boldsymbol{e}_{1}$
and we pick the orthonormal frame
$\boldsymbol{e}^{n}_{1}=\cos\theta_{h}\boldsymbol{e}_{1}-\sin\theta_{h}\boldsymbol{e}_{3},\;\boldsymbol{e}^{n}_{2}=\boldsymbol{e}_{2}$.
In the case of axisymmetric potential (11), we just take
$\boldsymbol{e}^{n}_{1}=\boldsymbol{e}_{1}$,
$\boldsymbol{e}^{n}_{2}=\boldsymbol{e}_{2}$.

We briefly comment on two cases where the condition (10) is not satisfied. One
is where we use the potential (11), but with $h_{z}+2h_{a}<0$. In that case
the minimum of $V(\boldsymbol{n})$ is attained on a circle and
$\boldsymbol{n}_{0}$ is some particular point on this circle which
spontaneously breaks the rotational symmetry of the potential. Thus we call
these potentials symmetry-breaking. We now have a zero-mode in one direction
away from $\boldsymbol{n}_{0}$, so that picking a basis
$(\boldsymbol{\psi}_{\boldsymbol{n}})_{1},(\boldsymbol{\psi}_{\boldsymbol{n}})_{2}$
that diagonalises
$\frac{\partial^{2}\tilde{V}}{\partial((\psi_{\boldsymbol{n}})_{i}\partial(\psi_{\boldsymbol{n}})_{j}}\bigg{\rvert}_{\boldsymbol{n}_{0}}$
we have diagonal entries

$\frac{\partial^{2}\tilde{V}}{\partial(\psi_{\boldsymbol{n}})_{1}\partial(\psi_{\boldsymbol{n}})_{1}}\bigg{\rvert}_{\boldsymbol{n}_{0}}=\mu^{2},\quad\frac{\partial^{2}\tilde{V}}{\partial(\psi_{\boldsymbol{n}})_{2}\partial(\psi_{\boldsymbol{n}})_{2}}\bigg{\rvert}_{\boldsymbol{n}_{0}}=0,$
(13)

where $\mu^{2}=-2h_{a}\left(1-\left(\frac{h_{z}}{2h_{a}}\right)^{2}\right)$.

The second case is a tilted Zeeman interaction combined with anisotropy, where
we have two different non-zero masses:

$\frac{\partial^{2}V}{\partial(\psi_{\boldsymbol{n}})_{1}\partial(\psi_{\boldsymbol{n}})_{1}}\bigg{\rvert}_{\boldsymbol{n}_{0}}=\mu_{1}^{2},\quad\frac{\partial^{2}V}{\partial(\psi_{\boldsymbol{n}})_{2}\partial(\psi_{\boldsymbol{n}})_{2}}\bigg{\rvert}_{\boldsymbol{n}_{0}}=\mu_{2}^{2}.$
(14)

Both of these cases could still be treated by the same methods we will use
below, but the calculations become more complicated because the linearised
Euler-Lagrange equations have less symmetry.

###  2.2 Asymptotic form of the chiral magnet energy functional

Soliton solutions can generically be split into two parts: the soliton tail,
where the fields are close to the vacuum $\boldsymbol{n}_{0}$, and the soliton
core, which is everywhere else. We expect this to be a small compact region in
order to minimize the energy. When we calculate the asymptotic interaction
potential, it is expressed in terms of these soliton tails. So we must
approximate these in order to approximate the interaction potential.

We can approximate the tails of an exact solution to the non-linear Euler-
Lagrange equations in terms of a corresponding solution to the Euler-Lagrange
equations linearised around $\boldsymbol{n}_{0}$. We will quantify the
accuracy of this approximation in terms of the distance from the core below.
To find these linearised Euler-Lagrange equations, we expand the energy (1) to
quadratic order in $\boldsymbol{\psi}_{\boldsymbol{n}}$ according to (9):

$E_{(2)}(\boldsymbol{\psi}_{\boldsymbol{n}})=\int\left(\frac{1}{2}\partial_{i}\boldsymbol{\psi}_{\boldsymbol{n}}\cdot\partial_{i}\boldsymbol{\psi}_{\boldsymbol{n}}+\boldsymbol{D}_{i}\cdot(\boldsymbol{n}_{0}\times\partial_{i}\boldsymbol{\psi}_{\boldsymbol{n}})+\boldsymbol{D}_{i}\cdot(\boldsymbol{\psi}_{\boldsymbol{n}}\times\partial_{i}\boldsymbol{\psi}_{\boldsymbol{n}})+\frac{1}{2}\mu^{2}\lvert\boldsymbol{\psi}_{\boldsymbol{n}}\rvert^{2}\right)d^{2}x,$
(15)

bearing in mind that this energy density and thus the resulting Euler-Lagrange
equations are only valid at large $r$, where
$\boldsymbol{\psi}_{\boldsymbol{n}}$ is small. This is equivalent to expanding
the Euler-Lagrange equations directly, but helps us see how the DMI becomes a
$U(1)$ background gauge field for the tails.

Since
$\boldsymbol{D}_{i}\cdot(\boldsymbol{n}_{0}\times\partial_{i}\boldsymbol{\psi}_{\boldsymbol{n}})$
is a divergence, it can be turned into a boundary term which does not affect
the Euler-Lagrange equations. It can be removed entirely by a redefinition of
the energy as in [30]. For the purpose of determining these equations we
therefore ignore this term, but will revisit it when calculating interaction
energies.

Further, since
$(\boldsymbol{\psi}_{\boldsymbol{n}}\times\partial_{i}\boldsymbol{\psi}_{\boldsymbol{n}})\parallel\boldsymbol{n}_{0}$,

$E_{(2)}(\boldsymbol{\psi})=\int\left(\frac{1}{2}(\partial\boldsymbol{\psi})^{2}+\boldsymbol{D}_{i}^{\parallel}\cdot(\boldsymbol{\psi}\times\partial_{i}\boldsymbol{\psi})+\frac{1}{2}\mu^{2}\lvert\boldsymbol{\psi}\rvert^{2}\right)d^{2}x+\text{boundary
term},$ (16)

where we split $\boldsymbol{D}_{i}$ into components parallel and perpendicular
to $\boldsymbol{n}_{0}$, so
$\boldsymbol{D}_{i}=\boldsymbol{D}_{i}^{\parallel}+\boldsymbol{D}_{i}^{\perp}$.
We also replace $\boldsymbol{\psi}_{\boldsymbol{n}}$ with $\boldsymbol{\psi}$,
as it is a dummy variable at this point.

Defining $a_{i}=\boldsymbol{D}_{i}\cdot\boldsymbol{n}_{0}$, and collecting the
components of $\boldsymbol{\psi}$ into complex scalar field,
$\psi=\psi_{1}+i\psi_{2}$, this energy can be rewritten as

$E_{(2)}(\boldsymbol{\psi})=E_{\text{lin}}(\psi)+\text{boundary term},$ (17)

where

$E_{\text{lin}}(\psi)=\int\left(\frac{1}{2}\lvert(\vec{\partial}+i\vec{a})\psi\rvert^{2}+\frac{1}{2}(\mu^{2}-\lvert\vec{a}\rvert^{2})\lvert\psi\rvert^{2}\right)d^{2}x$
(18)

is the energy functional of a linear field theory for a complex scalar in the
background of a fixed abelian gauge field $\vec{a}$. The first term in the
energy is a gauged Dirichlet energy for a complex scalar, and the second term
is a ‘mass term’ in the language of quantum field theory. However, unusually
and importantly for us, the effective mass $m$ of the scalar is a function of
both the potential and the DMI, given by

$m=\sqrt{\mu^{2}-\lvert\vec{a}\rvert^{2}}.$ (19)

Note that $E_{\text{lin}}$ is only bounded below if $m^{2}\geq 0$. If
$m^{2}<0$, it can be made arbitrarily negative. This shows that we would be
expanding the energy around the wrong vacuum in this case: the state where the
field is everywhere equal to the minimum of the potential,
$\boldsymbol{n}(\vec{x})=\boldsymbol{n}_{0}$, is no longer a minimum of the
energy functional due to the effect of the DMI. Thus $m^{2}=0$ corresponds to
some sort of phase transition, the details of which will depend on the full
nonlinear potential. At $m^{2}=0$, $E_{\text{lin}}$ is bounded below and thus
a sensible energy functional viewed in its own right. However we would need to
expand the nonlinear energy functional to higher order to see if we are really
expanding around the correct vacuum. For the following we will consider
parameter regimes where $m^{2}>0$ and thus $m>0$, but we will return to the
significance of the phase transition in Sec. 2.6.

Another way to view the asymptotic abelianisation of the theory is to note
that $E_{\text{lin}}$ is the asymptotic form of a nonlinear $U(1)$ gauge
theory. To see this, write the energy as

$E(\boldsymbol{n})=\int\left(\frac{1}{2}\lvert\partial_{i}\boldsymbol{n}+\boldsymbol{A}_{i}^{\parallel}\times\boldsymbol{n}\rvert^{2}+\boldsymbol{D}_{i}^{\perp}\cdot(\boldsymbol{n}\times\partial_{i}\boldsymbol{n})+\bar{V}^{\parallel}(\boldsymbol{n})\right)d^{2}x,$
(20)

where now we can only act with rotations around $\boldsymbol{n}_{0}$ for our
gauge transformations, and the formula for the material parameters of a
particular theory reached by our gauge transformations is

$\displaystyle\boldsymbol{D}_{i}$
$\displaystyle=\boldsymbol{D}_{i}^{\perp}+\boldsymbol{A}_{i}^{\parallel}$
$\displaystyle V(\boldsymbol{n})$
$\displaystyle=\bar{V}^{\parallel}(\boldsymbol{n})+\frac{1}{2}\lvert\boldsymbol{A}_{i}^{\parallel}\rvert^{2}$
(21)

The potential of this new gauge theory $\bar{V}^{\parallel}(\boldsymbol{n})$
is different is different from both the physical potential $V(\boldsymbol{n})$
and the potential $\bar{V}(\boldsymbol{n})$. Crucially, as long as $m^{2}>0$
it has the same minimum $\boldsymbol{n}_{0}$ as $V(\boldsymbol{n})$. The term
involving $\boldsymbol{D}_{i}^{\perp}$ does not contribute to
$E_{(2)}(\boldsymbol{\psi})$ and thus does not appear in the linearised Euler-
Lagrange equations. So the potential $\bar{V}^{\parallel}(\boldsymbol{n})$ is
the most meaningful in terms of understanding the soliton tails and the
$m^{2}=0$ phase transition.

This phenomenon of asymptotic abelianisation is familiar from the simplest
non-abelian Higgs model [33] and underlies the theory of ’t Hooft-Polyakov
monopoles [34, 35], although the context is quite different. Here the gauge
field $\boldsymbol{A}_{i}$ is non-dynamical, so the concept of mass does not
apply and the Higgs mechanism has no analogue. We are instead concerned with
the fact that the scalar complex field $\psi$, corresponding to the would-be
Goldstone bosons of that theory, is acted on by a $U(1)$ gauge field, the
massless photon of that theory. In the theory of magnetic skyrmions, most
often
$\boldsymbol{D}_{1}^{\parallel}=\boldsymbol{D}_{2}^{\parallel}=\boldsymbol{0}$,
giving a trivial $U(1)$ theory, so this viewpoint has not previously been
applied.

Since $\vec{a}$ is constant and a $U(1)$ connection, it is flat and can be
‘gauged away’ by defining

$\tilde{\psi}=e^{i\vec{a}\cdot\vec{x}}\psi.$ (22)

Then the linearised energy takes the standard form

$E_{\text{lin}}(\tilde{\psi})=\int\left(\frac{1}{2}\lvert\vec{\partial}\tilde{\psi}\rvert^{2}+\frac{1}{2}m^{2}\lvert\tilde{\psi}\rvert^{2}\right)d^{2}x,$
(23)

where the dependence on $\vec{a}$ is only in the effective mass $m$. This
expression for the linearised energy can also be seen as coming from the
nonlinear energy functional (20): applying the gauge transformation given by
$\boldsymbol{n}\mapsto\tilde{\boldsymbol{n}}=R(\vec{a}\cdot\vec{x},\boldsymbol{n}_{0})\boldsymbol{n}$
eliminates $\boldsymbol{A}_{i}^{\parallel}$ according to (3), and then
expanding the resulting energy to quadratic order around
$\tilde{\boldsymbol{n}}_{0}$ gives us the same final result.

Like the energy expression (18), the equivalent form (23) is only a good
approximation to the energy density of a soliton in the asymptotic region far
from the centre of the soliton, as represented by large $r$. So solutions of
the Euler-Lagrange equation

$(-\partial_{i}\partial_{i}+\mu^{2})\psi-2i\vec{a}\cdot\vec{\partial}\psi=0,\qquad
r\;\text{large},$ (24)

are expected to provide good approximations to the soliton tail in the
asymptotic region, but not to the soliton core. Using the gauge transformation
(22), this becomes the static Klein-Gordon equation

$(-\partial_{i}\partial_{i}+m^{2})\tilde{\psi}=0,\qquad r\;\text{large}.$ (25)

Thus the gauge transformation described above maps the linearised equation
into a standard linear problem whose solutions are well known. In the language
of quantum field field theory, it describes static excitations of a scalar
field with mass $m$ given by (19).

This is a result of independent interest: we see that spin excitations with
mass $\mu$ in the presence of DMI such that
$\boldsymbol{D}_{i}\cdot\boldsymbol{n}_{0}\neq 0$ behave like spin excitations
with a lower mass $m$ which twist along a direction picked by
$\boldsymbol{D}_{i}\cdot\boldsymbol{n}_{0}$. This reduction in mass and
twisting was observed in [22] for the specific case of Bloch-type DMI and a
potential of the form (12), but we see here that this is a general feature of
chiral magnets when the DMI and potential have the relation
$\boldsymbol{D}_{i}\cdot\boldsymbol{n}_{0}\neq 0$. As we shall see, the
simplicity of the linear problem (25) together with the transformation (22)
explains remarkably subtle and surprising features of the interaction of
magnetic solitons in a tilted field.

Note that asymptotically isotropic potential (10) is crucial for the
simplicity of the linearised problem: if we had different ‘masses’ for
perturbing away from the vacuum in different directions, as in (13), (14), the
asymptotic potential would not be invariant under rotation around
$\boldsymbol{n}_{0}$. Therefore this gauge twist would still give us an
equivalent model where DMI vanishes asymptotically, but the potential would be
spatially dependent. The Euler-Lagrange equations can still be asymptotically
approximated in this case [36], but the $U(1)$ symmetry which simplifies
subsequent calculations has been lost.

###  2.3 Soliton tails in the form of a multipole expansion

We now aim to make concrete the claim at the beginning of Sec. 2.2, that a
given nonlinear solution can be approximated by a specific solution of (25).
Firstly, given a solution $\boldsymbol{n}$ to the full nonlinear Euler-
Lagrange equations, $\boldsymbol{\psi}_{\boldsymbol{n}}$ is the corresponding
linear field satisfying an appropriate nonlinear equation according to (9),
leading us to define
$\psi_{\boldsymbol{n}}=(\boldsymbol{\psi}_{\boldsymbol{n}})_{1}+i(\boldsymbol{\psi}_{\boldsymbol{n}})_{2}$
and
$\tilde{\psi}_{\boldsymbol{n}}=e^{i\vec{a}\cdot\vec{x}}\psi_{\boldsymbol{n}}$.
Our claim is that there is a specific solution to (25), which we call
$\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}}$, that approximates
$\tilde{\psi}_{\boldsymbol{n}}$. In particular
$\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}}\to 0$ at spatial infinity to
reflect the fact that $\boldsymbol{n}\to\boldsymbol{n}_{0}$.

To proceed we look for the general solution to (25), satisfying the same
external boundary condition. This is well-known in different coordinate
systems. For us, solutions in polar coordinates are particularly interesting
because they give rise to an interpretation of a soliton ‘from afar’ as the
source of a linear multipole field. The choice of a centre for our polar
coordinate system $(r,\phi)$ is arbitrary at this point, but we will need to
discuss this below. Having picked coordinates, we can expand as follows:

$\tilde{\psi}(r,\phi)=\sum_{M=-\infty}^{\infty}c_{M}^{\text{lin}}(r)e^{iM\phi}$
(26)

with $c_{M}^{\text{lin}}(r)$ a complex function. Substitution into (25) yields
the Bessel equation for each $c_{M}^{\text{lin}}$,

$\sum_{M=-\infty}^{\infty}\left(c_{M}^{\text{lin}\prime\prime}(r)+\frac{1}{r}c_{M}^{\text{lin}\prime}(r)-\frac{M^{2}}{r^{2}}c_{M}^{\text{lin}}(r)+m^{2}c_{M}^{\text{lin}}(r)\right)e^{iM\phi}=0,\qquad
r\;\text{large},$ (27)

With the boundary condition that $c_{M}^{\text{lin}}$ vanishes as $r\to\infty$
and provided $m>0$, this equation is solved by the modified Bessel function of
the second kind of order $M$ [37], that is $c_{M}^{\text{lin}}(r)=K_{M}(mr)$
[17]. For $m=0$, we find $c_{0}^{\text{lin}}(r)=0$, $c_{M\neq
0}^{\text{lin}}(r)=r^{-M}$. We need exponential falloff for the derivation
below, so we do not consider this case. For large $r$, the leading terms are

$K_{M}(mr)=\sqrt{\frac{\pi}{2}}\frac{e^{-mr}}{\sqrt{mr}}+O\left(\frac{e^{-mr}}{(mr)^{\frac{3}{2}}}\right).$
(28)

We can then write the general solution, with $C_{M}$ arbitrary complex
numbers:

$\tilde{\psi}(\vec{x})=\sum_{M=-\infty}^{\infty}C_{M}K_{M}(mr)e^{iM\phi}.$
(29)

The missing information here that picks out the specific solution
$\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}}$ that we are looking for is the
unspecified internal boundary conditions on (25). These reflect the non-linear
core of the specific soliton whose tails we want to approximate. To find
$\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}}$, we specify $C_{M}$ such that if
we analogously expand the nonlinear solution

$\tilde{\psi}_{\boldsymbol{n}}(r,\phi)=\sum_{M=-\infty}^{\infty}c_{M}(r)e^{iM\phi},$
(30)

then $c_{M}(r)\to C_{M}K_{M}(mr)$ as $r\to\infty$, i.e. we set
$C_{M}=\lim_{r\to\infty}\frac{c_{M}(r)}{K_{M}(mr)}$, assuming such a limit
exists.

At this point we we can quantify the accuracy of approximating
$\tilde{\psi}_{\boldsymbol{n}}$ by
$\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}}$. Above we derived the linearised
Euler-Lagrange equations by expanding the energy to quadratic order in $\psi$.
Equally, we can take the nonlinear Euler-Lagrange equation

$\boldsymbol{n}\times\left(-\partial_{i}\partial_{i}\boldsymbol{n}+2\partial_{i}\boldsymbol{n}\times\boldsymbol{D}_{i}+\frac{\partial
V}{\partial\boldsymbol{n}}\right)=0$ (31)

and expand
$\boldsymbol{n}=\exp_{\boldsymbol{n}_{0}}(\boldsymbol{\psi}_{\boldsymbol{n}})$
for large radius. If we keep only linear terms in
$\boldsymbol{\psi}_{\boldsymbol{n}}$, then we find Equation (24). This is
solved by (29) in general, but only
$\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}}$ as defined above will have the
property that
$\delta\tilde{\psi}=\tilde{\psi}_{\boldsymbol{n}}-\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}}$
is potentially subleading. If we expand to the next order and rewrite in terms
of this $\delta\tilde{\psi}$, then we see that

$\displaystyle(-\partial_{i}\partial_{i}+m^{2})\delta\tilde{\psi}+O(\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}})^{2})=0\implies\tilde{\psi}_{\boldsymbol{n}}(\vec{x})=\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}}(\vec{x})+O\left(\frac{e^{-2mr}}{mr}\right).$
(32)

The difference between $\tilde{\psi}_{\boldsymbol{n}}$ and
$\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}}$ is therefore subleading when we
substitute it into the interaction energy below.

###  2.4 Soliton tails in terms of an effective source

In the derivation above we solved Equation (25), always bearing in mind that
it only applies at large $r$. For what follows it is useful to consider the
field $\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}}$ as defined on the whole
plane. This can be done at the cost of introducing a complex ‘effective
source’ $\tilde{\rho}_{\boldsymbol{n}}$:

$\displaystyle(-\partial_{i}\partial_{i}+m^{2})\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}}=\tilde{\rho}_{\boldsymbol{n}}(\vec{x}).$
(33)

This effective source is just a different way of writing the information of an
asymptotic solution, and contains no new information. The effective source is
introduced because the interaction potential turns out to be written simply in
terms of the effective source of one soliton interacting with the tail of the
other, which means that the interaction of two solitons can be approximated by
the interaction of their effective sources.

We can rewrite the information contained in the constants $C_{M}$ in terms of
the function $\tilde{\rho}_{\boldsymbol{n}}$. It is useful to define the
derivative:

$D_{M}=\begin{cases}(\partial_{1}+i\partial_{2})^{M}&M>0\\\ 1&M=0\\\
(\partial_{1}-i\partial_{2})^{M}&M<0.\end{cases}$ (34)

By using the recurrence relations satisfied by the modified Bessel functions
of the second kind [37], we can see that these derivatives move us along the
series of independent solutions to (25):

$\displaystyle
D_{1}\left(K_{M}(mr)e^{iM\phi}\right)=-mK_{M+1}(mr)e^{i(M+1)\phi}$
$\displaystyle
D_{-1}\left(K_{M}(mr)e^{iM\phi}\right)=-mK_{M-1}(mr)e^{i(M-1)\phi},$ (35)

so in particular $\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}}(\vec{x})$ is
equal to an infinite sum of derivatives acting on $K_{0}(mr)$. Finally, we
note that $\frac{1}{2\pi}K_{0}(mr)$ is the fundamental solution to this
equation, i.e. the solution in the presence of a delta-function source. Using
integration by parts, we can therefore write the general solution
$\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}}(r,\phi)$ as the solution to (33)
when $\tilde{\rho}_{\boldsymbol{n}}$ is equal to a combination of sources of
the form $D_{M}\delta(\vec{x})$. These have the interpretation of being
idealised multipole sources [38], from which the multipole expansion gets its
name.

To be precise, if we write the source as follows:

${}\tilde{\rho}_{\boldsymbol{n}}(\vec{x})=2\pi\sum_{M=-\infty}^{\infty}m^{-\lvert
M\rvert}q_{M}e^{i\gamma_{M}}D_{M}\delta(\vec{x})$ (36)

with $q_{M}>0$, $\gamma_{M}\in[0,2\pi)$ real numbers that we can interpret as
the strength and orientation of the multipoles respectively, then we can
consider the solution to (33) in the presence of this source:

$\displaystyle\tilde{\psi}_{\boldsymbol{n}}^{\text{lin}}(\vec{x})$
$\displaystyle=\int\tilde{\rho}_{\boldsymbol{n}}(\vec{x}^{\prime})\frac{1}{2\pi}K_{0}(m\lvert\vec{x}-\vec{x}^{\prime}\rvert)d^{2}x^{\prime}$
$\displaystyle=\sum_{M=-\infty}^{\infty}(-1)^{M}m^{-\lvert
M\rvert}q_{M}e^{i\gamma_{M}}D_{M}K_{0}(mr)$
$\displaystyle=\sum_{M=-\infty}^{\infty}q_{M}e^{i\gamma_{M}}e^{iM\phi}K_{\lvert
M\rvert}(mr)$ (37)

and we see that it is equal to (29), with

$C_{M}=q_{M}e^{i\gamma_{M}}.$ (38)

So we see that (36) contains all of the data of the general solution to (25).
The freedom we had in picking an origin for our polar co-ordinate system
manifests itself here as a freedom in choosing where the effective source is
located. This can be anywhere in the compact region containing the soliton
core. If we choose a different location then the numerical values $q_{M}$,
$\gamma_{M}$ will change, so multipole moments are always defined relative to
a chosen centre.

When solitons have reflection or rotation symmetries it is advantageous to
adapt the polar coordinate to these symmetries to simplify the multipole
expansion. We have seen one example of symmetry already, for the case of
axisymmetric DMI (6). If a soliton is symmetric under the rotation symmetry
(7) around a particular point then picking that point as the origin of our co-
ordinate system we find that $q_{M\neq 1}=0$.

Axisymmetric DMI has another symmetry, given by reflection of space and the
magnetisation at different angles. If we call the reflection about the line at
an angle $\phi_{0}$ relative to the x-axis $P_{\phi_{0}}$:

$P_{\phi_{0}}=\begin{pmatrix}\cos(2\phi_{0})&\sin(2\phi_{0})\\\
\sin(2\phi_{0})&-\cos(2\phi_{0})\end{pmatrix},$ (39)

then this symmetry is given by

$(\vec{x},(n_{1},n_{2}))\mapsto(P_{\phi_{0}}\vec{x},P_{\phi_{0}+\beta-\frac{\pi}{2}}(n_{1},n_{2})).$
(40)

Again, a soliton may satisfy this symmetry for a given $\phi_{0}$ or set of
different $\phi_{0}$. If we pick the origin of our polar co-ordinate system to
be on the line of the symmetry, we can impose this symmetry on the multipole
expansion (37) to constrain the multipoles, bearing in mind that our answer
will depend on our choice of basis vectors of the tangent space
$\boldsymbol{e}^{n}_{1}$, $\boldsymbol{e}^{n}_{2}$ that we made above. For
$\boldsymbol{n}_{0}=\boldsymbol{e}_{3}$ and thus
$\boldsymbol{e}^{n}_{1}=\boldsymbol{e}_{1}$,
$\boldsymbol{e}^{n}_{2}=\boldsymbol{e}_{2}$, we find

$q_{M}\neq
0\implies\gamma_{M}=-\frac{\pi}{2}+\beta-(M-1)\phi_{0}\quad\mod\pi.$ (41)

In the case of axisymmetric potential (11) and axisymmetric DMI (6), the
energy in total has $O(2)\ltimes\mathbb{R}^{2}$ symmetry. Individual solutions
may break it to a discrete subgroup such as
$\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ for reflections at two perpendicular
lines or $\mathbb{Z}_{2}$ for a single reflection [26]. In these cases,
$\phi_{0}$ is a free parameter describing the orientation of the soliton and a
zero-mode of the energy. In the case of a single $\mathbb{Z}_{2}$ invariance,
(41) tells us how all $\gamma_{M}$ must change as the orientation of the
soliton changes. In the case of two reflection symmetries at right angles,
(41) not only sets all $\gamma_{M}$ but also sets $q_{M}=0$ for even $M$. In
particular the antiskyrmion in an axisymmetric potential obeys this latter
constraint.

Alternatively, solitons may retain the full $O(2)$ symmetry of the energy, as
illustrated by skyrmion and skyrmionium solutions [39]. If we centre the polar
coordinate system at the fixed point of the spatial rotations, in addition to
constraining $q_{M\neq 1}=0$ it also requires that
$\gamma_{1}=-\frac{\pi}{2}+\beta$. So solitons like the skyrmion and
skyrmionium act as dipole source of fixed orientation for their linear tails.

We can also consider what happens when the energy explicitly breaks
axisymmetry, for example in the case of tilted field (12). In this case the
energy has only a single $\mathbb{Z}_{2}$ symmetry, determined by the
direction of tilt. For an external magnetic field in the direction
$\boldsymbol{e}_{h}=\sin\theta_{h}\cos\phi_{h}\boldsymbol{e}_{1}+\sin\theta_{h}\sin\phi_{h}\boldsymbol{e}_{2}+\cos\theta_{h}\boldsymbol{e}_{3}$,
one finds $\phi_{0}=\phi_{h}-\beta+\frac{\pi}{2}$. Numerically, skyrmions and
antiskyrmions in a tilted field are found to retain this symmetry. As
described in Sec. 2.1, if we take $\phi_{h}=0$ without loss of generality, we
can choose $\boldsymbol{e}^{n}_{1}$, $\boldsymbol{e}^{n}_{2}$ accordingly, and
in this case Equation (41) holds with $\phi_{0}=\frac{\pi}{2}-\beta$.

For completeness we note that the real fields $\tilde{\psi}_{1}$,
$\tilde{\psi}_{2}$ can be approximated by similar multipole expansions in
terms of $q_{M}^{1}$, $q_{M}^{2}$, $\gamma_{M}^{1}$, $\gamma_{M}^{2}$, subject
to the constraint that $q_{-M}^{1,2}=q_{M}^{1,2}$,
$\gamma_{-M}^{1,2}=-\gamma_{M}^{1,2}$. The terms $q_{1}^{1,2}$, $q_{2}^{1,2}$,
$q_{3}^{1,2}\ldots$ can be interpreted as the effective dipole, quadrupole,
octupole etc. sources for the corresponding real field $\tilde{\psi}_{1,2}$,
and $\gamma_{M}^{1,2}$ as the orientations of these multipoles [38]. The
complex sources that we work with are simply linear combinations of these
sources
$q_{M}e^{i}\gamma_{M}=q^{1}_{M}e^{i\gamma_{M}^{1}}+iq_{M}^{2}e^{i\gamma_{M}^{2}}$,
which have no constraint as $\tilde{\psi}$ is a complex field. Equation (41)
shows us that the orientations of these multipoles as described by
$\gamma^{1,2}_{M}$ or $\gamma_{M}$ cannot be directly interpreted as the
orientation of the soliton, but they are closely linked.

###  2.5 Numerical calculation of multipole moments

The definition of the multipole moments of a soliton in the previous section
relies on a division of the soliton field into a nonlinear core and linear
tail, and the results depend on the choice of origin for polar coordinates.
The numerical determination of the multipole moments demonstrates both of
these features. We illustrate this by considering a skyrmion and an
antiskyrmion in the model studied in [22], with a Bloch-type DMI with
parameters $\boldsymbol{D}_{i}=-2\pi\boldsymbol{e}_{i}$ ($\beta=0$) and a
tilted magnetic field
$\boldsymbol{e}_{h}=\boldsymbol{n}_{0}=\sin(\frac{\pi}{3})\boldsymbol{e}_{1}+\cos(\frac{\pi}{3})\boldsymbol{e}_{3}$
of strength $h_{z}=0.8\cdot(2\pi)^{2}$. This means $\phi_{0}=\frac{\pi}{2}$
and $m=2\pi\sqrt{0.05}\simeq 1.4$, so that the decay length of the soliton
tail is $\simeq 0.7$. We choose the origin of our polar coordinate system and
thus the location of our multipole sources to be at the point where
$\boldsymbol{n}(\vec{x})=-\boldsymbol{e}_{3}$. It would be more natural to
choose the point at which $\boldsymbol{n}=-\boldsymbol{n}_{0}$, but throughout
this paper we choose $-\boldsymbol{e}_{3}$ for consistency with the numerical
simulations of [22], with which we make quantitative comparison. For either
choice, the centre lies on the fixed line of the $\mathbb{Z}_{2}$ reflection
symmetry (40) of the skyrmion and antiskyrmion solution in this model, so we
expect to find $\gamma_{M}$ fixed according to (41). In the following we refer
to this point as the ‘soliton centre’.

We generate the isolated skyrmion and antiskyrmion by minimising the energy
from appropriate initial conditions, see Appendix A for details. To determine
the soliton centre, we approximate it by the point where $n_{3}$ is most
negative, since $\boldsymbol{n}$ is not exactly equal to $-\boldsymbol{e}_{3}$
anywhere on the lattice. Using the soliton centre as the origin for polar
coordinates $(r,\phi)$ we then invert (30) to get

$c_{M}(r)=\frac{1}{2\pi}\int_{0}^{2\pi}e^{-iM\phi}\tilde{\psi}_{\boldsymbol{n}}(r,\phi)d\phi.$
(42)

We know from Equation (32) that $c_{M}(r)\to
c_{M}^{\text{lin}}(r)=C_{M}K_{M}(mr)$ as $r\to\infty$. Thus at sufficiently
large $r$, $\lvert c_{M}(r)\rvert$ should approach a multiple of $K_{M}(mr)$,
and $\arg(c_{M}(r))$ should be independent of the radius. Moreover,
$\gamma_{M}$ and thus $\arg(c_{M}(r))$ should be fixed according to (41).
However, the uncertainty in the location of the soliton centre, which in
reality will not be at the exact point at which
$\boldsymbol{n}=-\boldsymbol{e}_{3}$, will lead to (41) being only
approximately satisfied. Theory and numerical results are compared in Figs. 1
and 2, with error bars coming from this uncertainty in the location of the
soliton centre, see Appendix A. We only plot comparisons for $-3\leq M\leq 3$,
since $\lvert c_{M}(r)\rvert$ is observed to fall off rapidly with increasing
$\lvert M\rvert$.

(a) skyrmion

(b) antiskyrmion

Figure 1: Magnitude of angular Fourier terms $c_{M}(r)$ of the tails of the
skyrmion and antiskyrmion in a tilted field relative to the Bessel functions
$K_{M}(mr)$, where the inverse decay lengthscale $m$ depends on the DMI
strength $k$ and the Zeeman interaction strength $h_{z}$ through
$m=\sqrt{h_{z}-k^{2}\sin^{2}\theta_{h}}$. The deviation for small $r$ is
expected, as this represents the nonlinear core of the soliton. The graph
should reach a constant value as $r$ becomes larger, giving a measurement of
the multipole strengths $q_{M}$ (36). Error bars come from the uncertainty in
the position of the centre of the polar co-ordinate system.

(a) skyrmion

(b) antiskyrmion

Figure 2: Argument of the angular Fourier terms $c_{M}(r)$ of the tails of the
skyrmion and antiskyrmion in a tilted field. Error bars come from the
uncertainty in the position of the centre of the polar co-ordinate system. We
exclude those $c_{M}(r)$ where the absolute value is so small that the error
in the argument becomes $O(2\pi)$. The lines plotted should approach a
constant at large $r$ where the linear approximation is justified, while at
small $r$ deviation is expected. If the soliton respects the reflection
symmetry (40), then (41) further tells us the value of this constant up to
$\pm\pi$. This graph thus acts as both a test of the linearised Euler-Lagrange
equations approximation and a check that numerically found skyrmions and
antiskyrmions in a tilted field have the reflection symmetry described.

We see that in terms of $\lvert c_{M}(r)\rvert$, the fit is good beyond
$r\simeq 0.8$ for the skyrmion, and $r\simeq 1.25$ for the antiskyrmion.
Meanwhile, $\arg(c_{M}(r))$ is constant all the way down to $r=0.25$ for some
$M$, and even within the core only deviates when the value $\lvert
c_{M}(r)\rvert/K_{M}(mr)$ is small. There are some $c_{M}(r)$ for which
$\lvert c_{M}(r)\rvert/K_{M}(mr)<0.01$ even for large $r$: these lead to
widely varying $\arg(c_{M}(r))$ and so these are not plotted in Fig. 2. In all
cases plotted, the theoretically predicted value of $\gamma_{M}$ according to
(41) is within the errorbars.

###  2.6 Stability of magnetic solitons in a tilted magnetic field

One of the remarkable features of chiral magnetic skyrmions is that, despite
their topological nature, they have a variety of instability modes, some of
which are studied numerically in [22, 26]. They include elliptic instabilties
where skyrmions or antiskrymions elongate indefinitely into domain walls for
certain values of the material parameters and the external magnetic field.

When $m^{2}<0$ in (33), it follows that the uniform state
$\boldsymbol{n}(\vec{x})=\boldsymbol{n}_{0}$ is linearly unstable, as is any
solution that approaches $\boldsymbol{n}_{0}$ in all directions. That is not
in itself a good enough reason to think that this region tells us anything
about elliptical instability. However, as we approach this region, the
lengthscale of decay diverges and so the whole nonlinear solution, while it
exists, will expand. At the same time, the oscillation along a given direction
that is favoured by the DMI as described by $\vec{a}$ will therefore become
more pronounced. These are suggestions that before or at the same time as this
phase transition, solitons become elliptically unstable.

We therefore plot the $m^{2}<0$ region and compare it to numerical results
[22] for the onset of elliptical instability. In the case of axisymmetric DMI
with strength $k$ (6) and a tilted applied field with strength $h_{z}$ and
tilt $\theta_{h}$ as in (12), $m^{2}<0$ is equivalent to
$h_{z}<k^{2}\sin^{2}\theta_{h}$. In Fig. 3 it is is plotted as a grey region
bounded by a black solid line, along with the numerical data for the onset of
elliptical instability for skyrmions and antiskyrmions in red and blue crosses
respectively.

In the numerics $k=1$ and the DMI is Bloch-type, i.e. $\beta=0$, but this
leads to no loss of generality as any axisymmetric DMI is equivalent to Bloch
under rotation of $\boldsymbol{n}$ and by defining suitable units of length
$k$ can be set to 1. Similarly, in the numerics only $Q=-1$ solutions are
considered. However, there is a transformation

$\displaystyle n_{3}$ $\displaystyle\mapsto-n_{3}$ $\displaystyle\vec{x}$
$\displaystyle\mapsto-\vec{x}$ $\displaystyle\theta_{h}$
$\displaystyle\mapsto\pi-\theta_{h}$ (43)

which leaves the energy of a solution invariant while changing the sign of its
charge, $Q\mapsto-Q$. Since this transformation changes $\theta_{h}$ it is not
a symmetry of the energy but a transformation that links equal-energy
solutions in different systems, like the gauge transformation introduced in
Sec. 2.1. This transformation interchanges the role of skyrmion and
antiskyrmion as $\theta_{h}$ is varied from $0$ to $\pi$, as noted in [22]. It
allows us to use the same numerical data to insert the points at which $Q=1$
solutions experience elliptical instability.

Figure 3: Comparison of theoretical predictions and numerical observations of
the onset of soliton instability in a tilted magnetic field, with strength
$h_{z}$, tilt $\theta_{h}$, DMI strength $k$. The blue shaded area inside the
dashed blue line is the region of the $(h_{z},\theta_{h})$ phase diagram where
domain walls have negative energy per unit length. The grey shaded area inside
the solid black line is the region where the vacuum is linearly unstable.
These give two theoretical estimates for the onset of instability. Blue and
red crosses show the numerically observed values of $(h_{z},\theta_{h})$ at
which $Q=\pm 1$ solutions (related by the transformation (43)) become
elliptically unstable as $h_{z}$ is decreased[22].

The lower bound for instability turns out to closely fit the numerical data
for applied fields tilted closer to the plane, while for fields close to the
perpendicular, the lower bound is not very useful; there is a large region of
the phase diagram for which elliptical instability happens but is not seen by
the calculation above. To understand the instability in this region, we use an
energy comparison with domain walls, generalising the method employed [40, 26]
for axisymmetric potential (11) to include tilted fields.

If, without loss of generality, we take Bloch-type DMI, that is axisymmetric
DMI (6) with $\beta=0$, and tilted potential (12), we can construct the domain
wall ansatz

$\boldsymbol{n}(x)=\begin{pmatrix}\sin(\theta_{h}+f(x_{2})\\\ 0\\\
\cos(\theta_{h}+f(x_{2}))\end{pmatrix},$ (44)

where $f(x_{2})$ goes from $0$ to $2\pi$ as we go from $-\infty$ to $+\infty$.
This orientation is chosen to give a negative contribution to the DMI. This
domain wall, if we do not specify this orientation, is the most general
solution approaching $\boldsymbol{n}_{0}$ at $x_{2}\to\pm\infty$ retaining the
symmetries of the energy $x_{1}\mapsto x_{1}+a$ and
$(x_{1},n_{2})\mapsto(-x_{1},-n_{2})$ (40), and thus by the principle of
symmetric criticality [41] if we minimise the energy of this ansatz with
respect to $f$ we will find a true stationary point of the energy. For tilted
field and axisymmetric DMI we can always use the combination of translation
and a reflection-like symmetry to have an ansatz depending on a single
function $f$, but as a result of the loss of the $O(2)$ symmetry such an
ansatz constrains us to consider domain walls lying parallel to some axis. For
this choice of DMI and tilt it is the $x_{1}$ axis, so we are specifically
investigating if solitons will extend along the $x_{1}$ direction by assuming
their cross-section in the $x_{2}$ direction is well-approximated by an
isolated domain wall.

Substituting (44) into the energy, we get the Sine-Gordon energy modified by a
boundary term coming from the DMI, just as in the case of normally applied
magnetic field without anisotropy:

$E_{\text{min}}(f)=\left(\int_{-\infty}^{\infty}dx_{1}\right)\left(-2\pi
k+8\sqrt{h_{z}}\right).$ (45)

Hence the energy per unit length of such a domain wall changes sign at
$h_{z}=\left(\frac{k\pi}{4}\right)^{2}$. For smaller $h_{z}$, the domain wall
has negative energy per unit length. This suggests an elliptical instability,
because any soliton necessarily has a cross-section that looks like the domain
wall above, and if domain walls have negative energy per unit length then we
expect them to grow in length.

The instability line $h_{z}=\left(\frac{k\pi}{4}\right)^{2}$ is also plotted
in Fig. 3 as a dashed blue line, with the shaded blue region corresponding to
$h_{z}<\left(\frac{k\pi}{4}\right)^{2}$. We see that the domain wall
instability calculation, by contrast with the $m^{2}=0$ curve, is most useful
as $\theta_{h}$ approaches $0$ or $\pi$. In the previous calculations for
axisymmetric potential [26], it was shown that the curve $E_{\text{min}}(f)=0$
is a particularly good fit to numerical observations of elliptical instability
for antiskyrmions, where the cross-section perpendicular to the ‘long’ axis of
the antiskyrmions is well approximated by an isolated domain wall. Following
the $Q=-1$ solution as we vary $\theta_{h}$ from $0$ to $\pi$ produces a
soliton that is effectively an antiskyrmion, in the sense that it is in
correspondence with the antiskyrmion solution at $\theta_{h}=0$ under the
transformation (43). This antiskyrmion-like solution is oriented so that its
cross-section along $x_{2}$ is a good approximation to an isolated domain
wall. Therefore we see this good fit as $\theta_{h}\to\pi$. Meanwhile, in
normal magnetic field the domain wall method slightly overestimates critical
$h_{z}$, as some energy barrier separates the axisymmetric skyrmion from an
arbitrarily extended $Q=-1$ solution. We see this overestimation emerge as
$\theta_{h}\to 0$. The reverse applies for the $Q=+1$ solutions.

##  3 Interaction potential for magnetic solitons

###  3.1 The interaction potential in terms of soliton tails

As described in the introduction, to define an interaction potential of
solitons we must first describe how we construct a suitable moduli space of
two interacting solitons. Given a field $\boldsymbol{n}^{A}$ describing
soliton $A$ and a field $\boldsymbol{n}^{B}$ describing soliton $B$, we must
construct a field $\boldsymbol{n}^{AB}$ that describes soliton $A$ and $B$
located at points $\vec{R}^{A}$, $\vec{R}^{B}$ in the plane. There is no
canonical way to do this, or even say what it means for a soliton to be
‘located’ at a particular point. In the following calculation we do not in
fact commit to a particular method, but for concreteness it is useful to have
a particular procedure in mind, so that we can see if our assumptions are
reasonable.

Here we take inspiration from the way interaction potentials are defined
numerically [42, 7, 22], which we here call pinning: in this procedure the
interaction energy of two skyrmions located at $\vec{R}^{A}$, $\vec{R}^{B}$ is
taken to be the minimum energy configuration reached by gradient descent from
some appropriate starting ansatz, subject to the constraint that the lattice
points at $\vec{R}^{A}$, $\vec{R}^{B}$ are fixed to a particular value within
the soliton core, generally $-\boldsymbol{e}_{3}$. We call the point within a
skyrmion where this value is attained the soliton centre in what follows. We
can think of this as an adiabatic approximation, where under gradient descent
the timescale of solitons moving together or apart is much larger than all
other timescales, so the solitons instantaneously adjust to minimise their
energy at a given separation.

We suppose that we can analytically define a moduli space in an analogous way,
with some modifications. Instead of gradient descent, which depends on our
choice of starting ansatz, we define $\boldsymbol{n}^{AB}$ as the absolute
minimiser over all configurations satisfying pinning constraints which
approximate solitons $A$ and $B$ having their centres at points $\vec{R}^{A}$,
$\vec{R}^{B}$. For skyrmions in axisymmetric potential, the point of
rotational symmetry makes a natural choice for soliton centre. For general
solitons, the choice is more arbitrary. Given some assumptions we will discuss
below, the interaction potential we calculate is independent of the choice of
centres, although different choices of centre will lead to the same potential
being described in terms of a different independent co-ordinate. Since we
compare to numerics where $\boldsymbol{n}=-\boldsymbol{e}_{3}$ is chosen as
the centre, we must pick that here as well. Again in contrast to skyrmions in
axisymmetric potential, general solitons may have internal degrees of freedom
such as orientation, represented here by (possibly multi-component) quantities
$\phi^{A},\phi^{B}$. In such cases, we must refine the procedure by adding
constraints to fix these parameters also. In theories with a translationally
invariant energy expression we can choose $\vec{R}^{A}=\vec{0}$ and write
simply $\vec{R}$ for $\vec{R}^{B}$, and we do this here.

With all this in mind, we define the configuration
$\boldsymbol{n}^{AB}[\vec{R};\phi^{A},\phi^{B}]$, which models the
interactions of solitons $A$ and $B$ described by fields $\boldsymbol{n}^{A}$,
$\boldsymbol{n}^{B}$, determined by internal parameters $\phi^{A}$,
$\phi^{B}$, and with charges $Q^{A}$, $Q^{B}$, as the absolute minimiser
within the space of $Q^{A}+Q^{B}$ configurations where
$\boldsymbol{n}(\vec{0})=\boldsymbol{n}(\vec{R})=-\boldsymbol{e}_{3}$, and any
further necessary constraints to determine $\phi^{A}$, $\phi^{B}$ and rule out
other pairs of interacting solitons with the same total charge. Our moduli
space modelling the interaction of these two solitons is then obtained by
allowing $\vec{R}$ to vary over a suitable open set typically of the form
$\\{\vec{R}\mid\lvert\vec{R}\rvert>R_{c}\\}$, and the internal degrees of
freedom to vary arbitrarily. The interaction potential on this moduli space is
defined as

$V^{AB}(\vec{R};\phi^{A},\phi^{B})=E(\boldsymbol{n}^{AB})-E(\boldsymbol{n}^{A})-E(\boldsymbol{n}^{B}).$
(46)

We now derive an asymptotic expression for the interaction potential in the
limit of large soliton separation, that is
$R=\lvert\vec{R}\rvert\gg\frac{1}{m}$. Motivated by the analogy to the
numerical procedure, we assume that as $R$ becomes large,
$\boldsymbol{n}^{AB}$ approaches the field $\boldsymbol{n}^{A}$ near
$\vec{0}$, and the field $\boldsymbol{n}^{B}$ near $\vec{R}$. We also assume
that as we go far from both $\vec{0}$ and $\vec{R}$, the field
$\boldsymbol{n}^{AB}$ approaches linear superposition of the tails of the two
solitons, as discussed in the introduction. Finally, we restrict ourselves to
consider local energy functionals with the property that $\boldsymbol{n}^{A}$,
$\boldsymbol{n}^{B}$ and their derivatives fall off exponentially towards the
vacuum away from their respective centres. Note that this last assumption is
satisfied in the case considered in Sec. 2.3, with the decay lengthscale given
by $\frac{1}{m}$ (29). These are the central assumptions to the results that
follow.

The formula we obtain is remarkably general and simple, and is the main result
of this paper. However, the derivation is rather technical and therefore
relegated to Appendix B. The calculation involves dividing the plane into two
infinite parts, $\sigma^{A}$ and $\sigma^{B}$, which respectively contain the
core of soliton $A$ and the core of soliton $B$. The dividing curve is thus
called $\partial\sigma^{A}$, see Fig. 4.

$\vec{0}$$\vec{R}$$\sigma^{A}$Soliton $A$ core$\sigma^{B}$Soliton $B$ core
Figure 4: Schematic representation of the setup for calculating the
interaction potential.

The result of the calculation is an approximation for the interaction
potential purely in terms of the soliton tails
$\boldsymbol{\psi}_{\boldsymbol{n}^{A}}$,
$\boldsymbol{\psi}_{\boldsymbol{n}^{B}}$:

$V^{AB}(\vec{R};\phi^{A},\phi^{B})=2\int_{\partial\sigma^{A}}\boldsymbol{\psi}_{\boldsymbol{n}^{B}}\cdot(\partial_{i}\boldsymbol{\psi}_{\boldsymbol{n}^{A}}+\boldsymbol{D}_{i}\times\boldsymbol{\psi}_{\boldsymbol{n}^{A}})dS^{i}+O(e^{-\frac{3}{2}mR}).$
(47)

where $dS^{i}$ represents the vector surface element of $\partial\sigma^{A}$.
Although $\partial\sigma^{A}$ appears in this expression, its exact form and
location do not matter, as we shall see below.

We can rewrite this in terms of the complex field
$\psi_{\boldsymbol{n}}=(\psi_{\boldsymbol{n}})_{1}+i(\psi_{\boldsymbol{n}})_{2}$
and $a_{i}=\boldsymbol{D}_{i}\cdot\boldsymbol{n}_{0}$ as defined in (18) to
get

$V^{AB}(\vec{R};\phi^{A},\phi^{B})=2\Re\int_{\partial\sigma^{A}}(\bar{\psi}_{\boldsymbol{n}^{B}}(\vec{\partial}+i\vec{a})\psi_{\boldsymbol{n}^{A}})\cdot
d\vec{S}+O(e^{-\frac{3}{2}mR}).$ (48)

Moreover, from (32), $\psi_{\boldsymbol{n}^{A}}$, $\psi_{\boldsymbol{n}^{B}}$
can be replaced by $\psi_{\boldsymbol{n}^{A}}^{\text{lin}}$,
$\psi_{\boldsymbol{n}^{B}}^{\text{lin}}$ to the level of approximation we are
already working at, meaning we can describe the leading behaviour of the
potential in terms of the linearised soliton tails:

$V^{AB}(\vec{R};\phi^{A},\phi^{B})=V^{AB}_{\text{lin}}(\vec{R};\phi^{A},\phi^{B})+O(e^{-\frac{3}{2}mR}),$
(49)

where

$V^{AB}_{\text{lin}}(\vec{R};\phi^{A},\phi^{B})=2\Re\int_{\partial\sigma^{A}}(\bar{\psi}_{\boldsymbol{n}^{B}}^{\text{lin}}(\vec{\partial}+i\vec{a})\psi_{\boldsymbol{n}^{A}}^{\text{lin}})\cdot
d\vec{S}.$ (50)

Equation (48) is the $U(1)$-gauged version of the interaction energy of
complex scalar fields in a linear theory. This is analogous to the interaction
potential of baby skyrmions [17], where the linear interaction potential is
equivalent to the interaction energy of complex scalar fields in a linear
theory without a gauge field, and thus equivalent to the independent
interaction of two real scalar fields.

Equation (50) tells us that if $\psi_{\boldsymbol{n}^{A}}^{\text{lin}}$ or
$\psi_{\boldsymbol{n}^{B}}^{\text{lin}}$ winds around the origin in the
complex plane as we vary either $\vec{R}$ or some of the $\phi^{A}$,
$\phi^{B}$, then $V^{AB}_{\text{lin}}$ must change sign. We will use this fact
below to determine when attraction can exist between magnetic solitons.

We also see that we have a natural expression given the interpretation of
$\vec{a}$ as a $U(1)$ gauge field:

$V^{AB}_{\text{lin}}(\vec{R};\phi^{A},\phi^{B})=2\Re\left(e^{-i\vec{a}\cdot\vec{R}}\int_{\partial\sigma^{A}}(\bar{\tilde{\psi}}_{\boldsymbol{n}^{B}}^{\text{lin}}\vec{\partial}\tilde{\psi}_{\boldsymbol{n}^{A}}^{\text{lin}})\cdot
d\vec{S}\right).$ (51)

Again, we see that this is just like the interaction energy without DMI, but
with a gauge twist applied to the fields, and with a modulation representing
the parallel transport between the two soliton centres with respect to the
gauge field $\vec{a}$. Integrating by parts in either the bulk $\sigma^{A}$ or
$\sigma^{B}$ and using the fact that
$\tilde{\psi}_{\boldsymbol{n}^{A}}^{\text{lin}}$,
$\tilde{\psi}_{\boldsymbol{n}^{B}}^{\text{lin}}$ solve (33) with corresponding
effective sources $\tilde{\rho}_{\boldsymbol{n}^{A}}$,
$\tilde{\rho}_{\boldsymbol{n}^{B}}$, we can write (51) also as bulk integrals:

$V^{AB}_{\text{lin}}(\vec{R};\phi^{A},\phi^{B})=-\Re\left(e^{-i\vec{a}\cdot\vec{R}}\int_{\sigma^{A}}\bar{\tilde{\psi}}_{\boldsymbol{n}^{B}}^{\text{lin}}\tilde{\rho}_{\boldsymbol{n}^{A}}d^{2}x\right)=\Re\left(e^{-i\vec{a}\cdot\vec{R}}\int_{\sigma^{B}}\tilde{\psi}_{\boldsymbol{n}^{A}}^{\text{lin}}\bar{\tilde{\rho}}_{\boldsymbol{n}^{B}}d^{2}x\right).$
(52)

This is like the interaction energy in [17], but before we take the real part
of this complex multipole interaction, we first multiply it by the parallel
transport factor $e^{-i\vec{a}\cdot\vec{R}}$, complicating a representation in
terms of real multipoles. Moreover, without the high symmetry of that
scenario, we cannot reduce it to solely an interaction of dipoles. However, we
do see the phenomenon discussed in the introduction: not only do solitons act
as a combination of multipole sources for their tails, which can be thought of
as just as a convenient mathematical representation, but they also respond
like those same multipole sources to the presence of another soliton. Note the
irrelevance of the exact location of the boundary $\partial\sigma^{A}$, given
that the sources $\tilde{\rho}_{\boldsymbol{n}^{A}}$,
$\tilde{\rho}_{\boldsymbol{n}^{B}}$ are located at $\vec{0}$, $\vec{R}$.

The advantage of this representation is that, by substituting in Equations
(36), (37) it will allow the interaction to be explicitly written in terms of
the separation of the two multipole sources. However, as discussed in Sec.
2.4, these multipoles are only meaningful given a particular choice of centre
of expansion. The final expression will be written in terms of the separation
between the centres that we choose, with multipoles that will be different for
different choices of centre. Therefore for different choices of centre,
$V^{AB}_{\text{lin}}$, expressed as a function of the separation between the
multipoles (which is itself a function of $\vec{R}$) could look quite
different. However, the expression (51) for $V^{AB}_{\text{lin}}(\vec{R})$ is
fundamentally in terms of $\tilde{\psi}_{\text{lin}}$ and thus independent of
the choice we make.

The expression (51) is valid without constraint on $\boldsymbol{D}_{i}$ and
$V(\boldsymbol{n})$, provided the tails fall off appropriately, so could be
used more generally. However, as discussed above, the solutions to the
linearised Euler-Lagrange equations for general $V(\boldsymbol{n})$ are more
complicated.

In [23], the authors substitute numerical solutions for the tails of specific
solitons at given couplings into (48), for a variety of energy functionals.
However, here we continue to an explicit expression for $V^{AB}_{\text{lin}}$,
which has certain advantages as discussed in the introduction. We also use
numerical simulation of an isolated soliton at specific couplings when we
compare to numerical results for the interaction potential, but even in this
case the approach is more general: a single numerical simulation of each
isolated soliton will give us enough information to approximate the
interaction potential between any two solitons at arbitrary $\vec{R}$.

###  3.2 The interaction potential in terms of effective sources

We can now explicitly calculate the interaction potential. We consider two
solitons, $A$ and $B$, described by corresponding sources
$\tilde{\rho}_{\boldsymbol{n}^{A}}$ and $\tilde{\rho}_{\boldsymbol{n}^{B}}$
located at $\vec{0}$ and $\vec{R}$ respectively:

$\displaystyle\tilde{\psi}_{\boldsymbol{n}^{A}}^{\text{lin}}(\vec{x};\phi^{A})$
$\displaystyle=\sum_{M}q^{A}_{M}e^{i(M\phi+\gamma^{A}_{M})}K_{\lvert
M\rvert}(mr)$
$\displaystyle\tilde{\psi}_{\boldsymbol{n}^{B}}^{\text{lin}}(\vec{x};\phi^{B})$
$\displaystyle=\sum_{M}q^{B}_{M}e^{i(M\phi(\vec{x}-\vec{R})+\gamma^{B}_{M})}K_{\lvert
M\rvert}(m\lvert\vec{x}-\vec{R}\rvert),$ (53)

where $\phi(\vec{x}-\vec{R})$ is the polar angle in co-ordinates centred on
$\vec{R}$, and $q_{M}^{A,B}$, $\gamma_{M}^{A,B}$ are functions of
$\phi^{A,B}$. This dependence can be constrained from symmetry, which we will
see below.

In the cases where $q_{M}^{A,B}$, $\gamma_{M}^{A,B}$ are not constrained
enough from symmetry, we find them numerically from the asymptotics of a
single soliton, using the results from Sec. 2.5. Note that this is why we
chose to make our multipole expansion around the point at which
$\boldsymbol{n}=-\boldsymbol{e}_{3}$ in that section: the interaction
potential is ultimately written in terms of the separation between the
multipole sources, while it is defined as a function of $\vec{R}$, the
separation between the soliton centres (as well as internal degrees of
freedom). To be able to write our potential explicitly, these two separations
must be the same, so the multipole sources must be chosen to lie at the
soliton centres.

We now substitute this expansion into the integral in (52), using the relation
between the expansion of a field and the expansion of its source expressed in
equations (36), (37) :

$\displaystyle\int_{\sigma^{A}}\bar{\tilde{\psi}}_{\boldsymbol{n}^{B}}^{\text{lin}}\tilde{\rho}_{\boldsymbol{n}^{A}}d^{2}x$
$\displaystyle=2\pi(-1)^{N}\sum_{M,N}q^{A}_{M}q^{B}_{N}e^{i(\gamma^{A}_{M}-\gamma^{B}_{N})}K_{\lvert
M-N\rvert}(mR)e^{i(M-N)\chi},$ (54)

and thus

$\displaystyle
V^{AB}_{\text{lin}}(\vec{R};\phi^{A},\phi^{B})=2\pi\sum_{M,N}(-1)^{N+1}q^{A}_{M}q^{B}_{N}\cos(\gamma^{A}_{M}-\gamma^{B}_{N}-\vec{a}\cdot\vec{R}+(M-N)\chi)K_{\lvert
M-N\rvert}(mR),$ (55)

where from $\vec{a}$ defined in (18) we define $a=\lvert\vec{a}\rvert$ and
$\chi$ as the angle between $\vec{R}$ and $\vec{a}$.

This is the interaction potential between solitons in a chiral magnet in large
generality: as described above, we have assumed that the potential is
asymptotically isotropic about the vacuum (10), general DMI, $m$ as defined in
(19) is real, and no other interactions. Because of the independence of which
soliton we take at $\vec{0}$, this potential has the symmetry
$\chi\to\chi+\pi$, $A\leftrightarrow B$. This means for interactions between
like solitons, $\chi\to\chi+\pi$ is a symmetry.

At sufficiently large $R$, we can use the expansion of $K_{M}(mR)$ in powers
of $\frac{1}{mR}$ given in Equation (28):

$V^{AB}_{\text{lin}}(\vec{R},\phi_{0}^{A},\phi^{0}_{B})=\sqrt{2\pi^{3}}f(\chi;\phi_{0}^{A},\phi_{0}^{B})\frac{e^{-mR}}{\sqrt{mR}}+O\left(\frac{e^{-mR}}{(mR)^{\frac{3}{2}}}\right),$
(56)

where

$f(\chi;\phi_{0}^{A},\phi_{0}^{B})=\sum_{M,N}(-1)^{N+1}q^{A}_{M}q^{B}_{N}\cos(\gamma_{M}^{A}(\phi_{0}^{A})-\gamma_{N}^{B}(\phi_{0}^{B})-aR\cos\chi+(M-N)\chi).$
(57)

At this radius, we can similarly expand the individual soliton tails:

$\tilde{\psi}_{\boldsymbol{n}^{A,B}}^{\text{lin}}(\vec{x};\phi^{A,B})=\left(\sum_{M}q^{A,B}_{M}e^{i(M\phi+\gamma^{A,B}_{M})}\right)\sqrt{\frac{\pi}{2}}\frac{e^{-mr}}{\sqrt{mr}}+O\left(\frac{e^{-mr}}{(mr)^{\frac{3}{2}}}\right),$
(58)

from which it follows that $V^{AB}_{\text{lin}}$ can be approximated in terms
of a product of the two soliton tails at the midpoint between them:

$\displaystyle
V^{AB}_{\text{lin}}(\vec{R};\phi^{A},\phi^{B})=-\boldsymbol{\psi}_{\boldsymbol{n}^{A}}\left(\frac{R}{2},\chi;\phi^{A}\right)\cdot\boldsymbol{\psi}_{\boldsymbol{n}^{B}}\left(\frac{R}{2},\chi+\pi;\phi^{B}\right)\sqrt{2\pi
mR}+O\left(\frac{e^{-mR}}{(mR)^{\frac{3}{2}}}\right).$ (59)

This is a remarkably simple form of the interaction energy, if one wants to
quickly estimate whether two solitons will attract or repel at a large
distance, but for the rest of the paper we will continue to use the expression
(55) for greater accuracy at smaller $R$.

###  3.3 Comparison to numerically calculated interaction potential in a
tilted field

In the case of tilted field (12), there are two known stable solitons, the
skyrmion and antiskyrmion [22]. These are known to both retain the single
$\mathbb{Z}_{2}$ symmetry of the energy that remains for tilted field, of the
form (40) with $\phi_{0}$ fixed by the parameters of the theory, and have no
internal degrees of freedom. Therefore to construct our moduli space, we do
not need to add any constraints besides fixing
$\boldsymbol{n}^{AB}(\vec{0})=\boldsymbol{n}^{AB}(\vec{R})=-\boldsymbol{e}_{3}$
and specifying the overall topological charge as $-2$, $0$ or $2$. According
to (41), $\gamma^{A}_{M}$ and $\gamma^{B}_{M}$ are fixed to the same value
such that $\gamma^{A}_{M}-\gamma^{B}_{N}=(N-M)\phi_{0}$. Without loss of
generality, we consider Bloch-type DMI (6), so $\beta=0$, and field tilt in
the direction of $\boldsymbol{e}_{1}$, as discussed in Sec. 2.1, so that
$\phi_{h}=0$, $\beta=0$ and thus $\phi_{0}=\frac{\pi}{2}$, so we have

$V^{AB}_{\text{lin}}(R,\chi)=2\pi\sum_{M,N}(-1)^{N+1}q^{A}_{M}q^{B}_{N}\cos\left((N-M)\frac{\pi}{2}-aR\cos\chi+(M-N)\chi\right)K_{\lvert
M-N\rvert}(mR).$ (60)

Once $\gamma_{M}^{A}$, $\gamma_{N}^{B}$ satisfy this, it imposes a
corresponding symmetry on the potential of $\chi\to\pi-\chi$. Meanwhile, we
calculate $q_{M}$ by using the results from Sec. 2.5, which only apply to the
specific material parameters considered: $h_{z}=0.8\cdot(2\pi)^{2}$, $k=2\pi$.
While this means that even in our analytical approach some numerics is
required, the cost of simulating a single soliton in isolation is much lower
than calculating the interaction directly, which requires pinning the soliton
at every possible distance from every other soliton it could potentially
interact with. Moreover, our interaction energy can be calculated for
arbitrarily large $R$, while numerical error becomes dominant in the direct
interaction calculation as the energy difference from isolated solitons
becomes exponentially small.

To calculate $q_{M}=\lvert C_{M}\rvert$, we fit our numerically found $\lvert
c_{M}(r)\rvert$ to $\lvert C_{M}\rvert K_{M}(mr)$. To find
$\gamma_{M}=\arg(C_{M})$, we take the angular mean of $\arg(c_{M}(r))$. Both
the fitting and the mean are taken between $r=1$ and $r=2$, so as to exclude
the nonlinear soliton core. The errors on $\arg(c_{M}(r))$ and $\lvert
c_{M}(r)\rvert$ are propagated through to errors in $\lvert C_{M}\rvert$,
$\arg(C_{M})$, bearing in mind that the errors are not independent at each
$r$. The resulting $C_{M}$ and errorbars are plotted in Fig. 5.

(a) skyrmion

(b) antiskyrmion

Figure 5: Plotting the numerically extracted multipole sources $C_{M}$ seen in
the general solution of the linearised Euler-Lagrange equations for solitons
in a tilted field,
$\tilde{\psi}_{\text{lin}}=\sum_{M=-\infty}^{\infty}C_{M}e^{iM\phi}K_{M}(mr)$,
where the inverse decay lengthscale $m=\sqrt{h_{z}-k^{2}\sin^{2}\theta_{h}}$
for DMI strength $k$, Zeeman interaction strength $h_{z}$.

When $\theta_{h}=\frac{\pi}{2}$ there is a symmetry of the energy
$n_{3}\mapsto-n_{3}$, $\vec{x}\mapsto-\vec{x}$ interchanging skyrmion and
antiskyrmion solutions [22], which can be seen as a special case of the
transformation discussed in Sec. 2.6, and as a result skyrmion and
antiskyrmion sources are related by a reflection around the imaginary axis. In
Fig. 5 we see that the symmetry still holds approximately, suggesting that the
values of the sources change continuously with respect to the material
parameters.

Having numerically found $q_{M}$, $\gamma_{M}$, for both skyrmion and
antiskyrmion, we can substitute them into (55) and plot
$V_{\text{lin}}(R,\chi)$, comparing it to the numerical results in [22].
Because that paper defined $\vec{R}$ as the distance between the two points
where $\boldsymbol{n}$ is fixed to equal $-\boldsymbol{e}_{3}$, we have
defined our multipole locations to be in the same place, as discussed above.

(a) skyrmion-skyrmion

(b) skyrmion-antiskyrmion

(c) antiskyrmion-antiskyrmion

(d) skyrmion-skyrmion

(e) skyrmion-antiskyrmion

(f) antiskyrmion-antiskyrmion

Figure 6: Analytical (a,b,c) and numerical[22] (d,e,f) results for the
interaction potential between solitons in a tilted field, as a function of
their separation. Dashed lines represent the limits of the numerical
antiskyrmion-antiskyrmion interaction data, which is the most restrictive.

Our result is only asymptotically valid, while the data is at small distances,
but there is nevertheless good agreement for $R>1.5$, see Fig. 6. We see that
because of the error in $\gamma_{M}$ discussed in Sec. 2.5, there is some
small violation of the $\chi\to\pi-\chi$ symmetry in the analytical
prediction, which could be removed by using the theoretically predicted values
of $\gamma_{M}$.

This analysis thus explains the observed oscillation of the inter-soliton
potential in a tilted field as a manifestation of the emergent $U(1)$ gauge
theory (18). In general it is fixed by the DMI lengthscale but with a non-
trivial factor $\frac{1}{\sin\theta_{h}}$. As we take the tilt angle to $0$,
we would see the period of the oscillation diverge to infinity.

We also have an explanation for the unusual asymmetric lobe structure of the
skyrmion-antiskyrmion interaction potential, by analogy to the Roget’s
palisade illusion [43], or rolling shutter effect [44], where rotational
motion of a some object at angular speed $\omega$ is combined with the lateral
motion of a camera shutter (at speed $v$) to produce an image where the object
is distorted. This distortion can be described in terms of a map on polar co-
ordinates [45]:

${}(r,\phi)\mapsto\left(r,\phi-\frac{\omega}{v}r\cos\phi\right).$ (61)

Figure 7: “Airplane Prop + CMOS Rolling Shutter = WTF", Soren Ragsdale,
licenced under CC BY 2.0

We can see that each term in (55) where $M\neq N$ looks like the distorted
image of the potential

$2\pi(-1)^{N+1}q^{A}_{M}q^{B}_{N}\cos((N-M)\phi_{0}+(M-N)\chi)K_{\lvert
M-N\rvert}(mR)$ (62)

if it were rotating at $\frac{\omega}{v}=\frac{a}{N-M}$. If one were to plot
this function, it would look like a propellor with $\lvert M-N\rvert$ blades,
hence the likeness between this interaction potential and the photo of rolling
shutter effect in Fig. 7. The combination of all terms with $M\neq N$, then,
is analogous to a photo of a combination of propellors with different numbers
of blades spinning at different speeds, but the dominant multipoles will have
the largest contribution. In our case, the effective propellor is dominated by
the $M=0$, $N=-1$ term, together with the $M=1$, $N=0$ term which is
proportional to it and ‘rotates’ at the same speed. Meanwhile, the terms with
$M=N$ all add together to give a term proportional to $\cos(-aR\cos\chi)$,
i.e. oscillation only along $\vec{a}$. The combination of these two parts
gives us the interaction potential above.

Between two identical solitons, equal contributions from oppositely rotating
propellors will add up and obscure the rolling shutter picture, as they must
to retain the $\chi\to\chi+\pi$ symmetry, but for interactions between unlike
solitons we can see this will be a general feature, assuming only a small
number of multipole sources contribute significantly as is seen here. This
rolling shutter analogy also gives us an idea of how the potential will change
as we take $\theta_{h}\to 0$: the lobes of positive and negative interaction
energy will move upwards as they become ‘less distorted’. This change is
illustrated qualitatively in Fig 8, where we take $a=k\sin\theta_{h}$ to zero
while leaving all other parameters constant.

(a) $a\mapsto\frac{a}{2}$

(b) $a\mapsto\frac{a}{5}$

(c) $a\mapsto 0$

Figure 8: The skyrmion-antiskyrmion interaction potential for different values
of $a=k\sin\theta_{h}$, where $k$ is the DMI strength and $\theta_{h}$ is the
applied field tilt. At the same time we leave other parameters $m$,
$q^{Sk}_{M}$, $q^{ASk}_{M}$,$\gamma^{Sk}_{M}$, $\gamma^{ASk}_{M}$ constant.
This models how the interaction potential will change as the field tilt is
reduced, but it is only schematic since in reality the parameters we are
fixing are also functions of $\theta_{h}$. In particular, as $\theta_{h}\to 0$
then $q^{Sk}_{M\neq 1}\to 0$ and $q^{ASk}_{M\text{even}}\to 0$, as both the
skyrmion and antiskyrmion regain symmetry. The interaction potential would
accordingly gain a second reflection symmetry about the x axis.

###  3.4 Soliton interactions in axisymmetric potential

In the case of axisymmetric potential and axisymmetric DMI, there are two
possibilities: either a soliton retains the $U(1)$ symmetry of the energy (7)
and thus has no internal degrees of freedom, or it breaks it and has one
internal degree of freedom coming from its orientation. We call these degrees
of freedom $\phi_{0}^{A}$, $\phi_{0}^{B}$ for the corresponding solitons. By
considering the action of this $U(1)$ symmetry on a soliton, we can see that
$q_{M}^{A,B}$ will be independent of $\phi_{0}^{A,B}$ while $\gamma_{M}^{A,B}$
will have a prescribed dependence.

We can generally define these orientations up to $\pm\pi$ by considering the
dispersion tensor
$\Gamma_{ij}=\int\partial_{i}\boldsymbol{n}\partial_{j}\boldsymbol{n}d^{2}x$
[46]: when this matrix is not proportional to the identity matrix, it has a
one-dimensional eigenspace corresponding to its largest eigenvalue, picking a
direction in the plane. The angle between this line and the $x$-axis defines
$\phi_{0}^{A,B}$. When the soliton has at least one reflection-like symmetry
(40), $\phi_{0}^{A,B}$ is also the angle of one of the spatial reflections,
and the $\gamma^{A}_{M},\gamma^{B}_{N}$ can be expressed in terms of
$\phi_{0}^{A,B}$ according to (41). Note that the
$\gamma^{A}_{M},\gamma^{B}_{N}$ will not see the $\pm\pi$ ambiguity.

Considering also that $a=0$ (this also means that
$\tilde{\psi}_{\text{lin}}=\psi_{\text{lin}}$), the general formula (55)
becomes

$V^{AB}_{\text{lin}}(R,\chi;\phi_{0}^{A},\phi_{0}^{B})=2\pi\sum_{M,N}(-1)^{N+1}q^{A}_{M}q^{B}_{N}\cos(\gamma_{M}^{A}(\phi_{0}^{A})-\gamma_{N}^{B}(\phi_{0}^{B})+(M-N)\chi)K_{\lvert
M-N\rvert}(mR).$ (63)

where we now make explicit that $\gamma_{M}^{A,B}$ depend on $\phi_{0}^{A,B}$
respectively, while $q_{M}^{A,B}$ do not.

According to (56), $V^{AB}_{\text{lin}}$’s dependence on $R$ and $\chi$ splits
at sufficiently large $R$:

$V^{AB}_{\text{lin}}(R,\chi;\phi_{0}^{A},\phi_{0}^{B})=\sqrt{2\pi^{3}}f(\chi;\phi_{0}^{A},\phi_{0}^{B})\frac{e^{-mR}}{\sqrt{mR}}+O\left(\frac{e^{-mR}}{(mR)^{\frac{3}{2}}}\right),$
(64)

where

$f(\chi;\phi_{0}^{A},\phi_{0}^{B})=\left(\sum_{M,N}(-1)^{N+1}q^{A}_{M}q^{B}_{N}\cos(\gamma_{M}^{A}(\phi_{0}^{A})-\gamma_{N}^{B}(\phi_{0}^{B})+(M-N)\chi)\right).$
(65)

Therefore for sufficiently large $R$, $V^{AB}_{\text{lin}}$ decreases in
absolute value as a function of $R$, so for given $\chi$ the attraction is
either repulsive, if positive, or attractive, if negative.

In axisymmetric potential and DMI, the $O(2)$ symmetry of the skyrmion around
the soliton centre sets all sources except $q_{1}$ to zero, and sets
$\gamma_{1}=\frac{\pi}{2}-\beta$, as discussed in Sec. 2.4. Thus setting
$q^{A}_{1}=q_{1}^{B}=:q_{1}^{Sk}$ in this case, the skyrmion-skyrmion
potential is

${}V^{SkSk}_{\text{lin}}(R,\chi)=2\pi(q_{1}^{Sk})^{2}K_{0}(mR).$ (66)

With the approximation of
$K_{0}(mr)\simeq\sqrt{\frac{\pi}{2}}\frac{e^{-mr}}{\sqrt{mr}}$, this is the
potential derived in [21]. It is also in principle applicable to the
interaction of skyrmionium with a skyrmion, or two skyrmioniums, as this
soliton has the same symmetry, e.g.

${}V^{SkSkm}_{\text{lin}}(R,\chi)=2\pi q_{1}^{Sk}q_{1}^{Skm}K_{0}(mR)$ (67)

where $q_{1}^{Skm}$ is the corresponding dipole source of the far field of the
skyrmionium with its centre defined as the point at which
$\boldsymbol{n}=\boldsymbol{e}_{3}$. As discussed above, this requires us to
define an unambiguous way to define the combined field of a skyrmion and
skyrmionium at separation $\vec{R}$. In all these cases, the interaction is
independent of $\chi$, and repulsive for any $R$.

However, other solitons are supported in axisymmetric DMI and potential [26].
For example, an antiskyrmion can be supported for a small range of coupling
parameters. We label its effective sources $q_{M}^{ASk}$, $\gamma_{M}^{ASk}$.
As discussed in Sec. 2.4, because it breaks the $O(2)$ symmetry of the energy
to a $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ subgroup, it has an orientation
which is free to vary, which we here call $\phi_{0}^{ASk}$. The multipole
orientations $\gamma_{M}^{ASk}$ are fixed in terms of this orientation:

$\gamma_{M}^{ASk}=-\frac{\pi}{2}+\beta-(M-1)\phi_{0}^{ASk}+n_{M}\pi,\quad
n_{M}\in\mathbb{Z}.$ (68)

while $q_{M}^{ASk}=0$ for $M$ even.

If we then consider the interaction of a skyrmion with an antiskyrmion, we
find

${}V^{SkASk}_{\text{lin}}(R,\chi;\phi_{0}^{ASk})=2\pi
q^{Sk}\sum_{N\text{odd}}q^{ASk}_{N}\cos((N-1)(\phi_{0}^{ASk}-\chi)-n_{M}\pi)K_{\lvert
N-1\rvert}(mR).$ (69)

Here we see the dependence of $V^{AB}$ on an internal parameter, as discussed
in Sec. 3.1. This potential is stationary with respect to variations of
$\phi_{0}^{ASk}$ for $\phi_{0}^{ASk}=\chi,\chi+\frac{\pi}{2}$, but we cannot
guarantee that these are the only critical values of $\phi_{0}^{ASk}$, nor say
whether they are maxima or minima. However we can show that the interaction
potential can be made negative for any $R$, $\chi$ by a suitable choice of
$\phi_{0}^{ASk}$. To show this we rely on the link between winding of an
individual tail and the sign of the overall interaction potential (48). If we
assume $\psi^{ASk}_{\text{lin}}(r,\phi;\phi_{0}^{ASk})$ is not equal to zero
for all values of $r$ larger than the soliton core size, it must wind once
clockwise around the origin in the complex plane as we vary $\phi$ from $0$ to
$2\pi$ for topological reasons, meaning the product of the fields above winds
twice clockwise around the origin as we vary $\chi$ from $0$ to $2\pi$.
Equivalently, the product above winds twice anticlockwise around the origin as
$\phi_{0}^{ASk}$ varies from $0$ to $2\pi$. So we know that the interaction
can be made negative by varying $\phi_{0}^{ASk}$. Therefore according to
Equation (64), the antiskyrmion can be oriented so as to attract the skyrmion.

This observation can be extended more generally to the variety of soliton
solutions in axisymmetric potential, provided we can extend the pinning
procedure so as to unambiguously find an interaction between any two solitons.
The general result confirms the observation made in [26], namely that any
solution with chiral kinks on the outer domain wall, which necessarily has a
zero-mode $\phi_{0}$ associated to the $U(1)$ symmetry of the energy (7), can
be rotated so that it attracts another soliton, provided the tail does not
reach the vacuum at any finite distance from the core. The argument
generalises what goes before: as we vary $\phi$ from $0$ to $2\pi$,
$\psi_{\text{lin}}(r,\phi;\phi_{0})$ will wind $-N_{\text{kink}}+1$ times
anticlockwise around the origin in the complex plane, where $N_{\text{kink}}$
is the chiral kink number, as defined in the above paper, associated to the
outermost domain wall of the soliton. Now we use the fact that since
$\phi_{0}$ describes the $U(1)$ symmetry of the energy, we can link it to the
variation of $\phi$:
$\psi_{\text{lin}}(r,\phi;\phi_{0}+\alpha)=e^{-i\alpha}\psi_{\text{lin}}(r,\phi+\alpha;\phi_{0})$.
This means that as we increase $\alpha$ from $0$ to $2\pi$,
$\psi_{\text{lin}}(r,\phi;\phi_{0}+\alpha)$ will wind $-N_{\text{kink}}$
times. This means that as we vary $\phi_{0}$ in any potential involving this
soliton, there will be $\lvert N_{\text{kink}}\rvert$ regions of negative
interaction potential and thus at sufficiently large $R$, the two solitons can
be oriented to attract each other. If two solitons attract at arbitrary
distance, it implies a bound state of the two, although we cannot guarantee
stability of the combined configuration against collapse. Nevertheless, this
general proof of attraction between solitons with chiral kinks can be seen as
a partial explanation of why antiskyrmions are not found alone but with a
large number of other magnetic solitons [26].

##  4 Conclusion

In this paper we presented a framework for calculating soliton interactions in
general and applied it to get explicit formulae in two new cases: chiral
magnetic skyrmions and antiskyrmions interacting in a tilted applied magnetic
field, and a variety of magnetic solitons interacting in a chiral magnet with
normally applied magnetic field and anisotropy. The treatment is general
enough to include Bloch and Néel-type DMI. In the case of tilted field, we
found close agreement between the analytical formula and previous numerical
observations [22]. In the case of normal magnetic field combined with
anisotropy, we found that solitons with chiral kinks [26] can always be
oriented so as to attract another soliton.

This calculation generalised previous ones in the topological soliton
literature by incorporating the effect of a $U(1)$ background gauge connection
on the effective scalar field theory that describes the tails of solitons.
This arises when considering tilted field applied to a chiral magnet, as the
non-zero overlap between the minimum of the potential $\boldsymbol{n}_{0}$ and
the DMI vectors $\boldsymbol{D}_{i}$ leads to a new term in the linearised
Euler-Lagrange equations. This $U(1)$ gauge connection naturally descends from
understanding the DMI as an $SO(3)$ gauge connection. This then gives an
explanation of the oscillating interaction potential generically observed in a
tilted field.

We also generalised previous calculations in the magnetic skyrmion literature
specifically, by understanding the far field of a general soliton as being
effectively sourced by an infinite set of multipoles, rather than just a
single one. This is required for any soliton that does not have the especially
high symmetry of, for example, skyrmion and skyrmionium solutions. To predict
interaction potentials, we therefore had to extract the strengths of these
multipoles from looking at the isolated soliton. This sort of calculation has
been done before in the context of nuclear skyrmions [24], although our
methods here are different. This allowed us to consider interactions involving
non-axisymmetric magnetic solitons like the antiskyrmion, which had not
previously been calculated.

The discussion in this paper can be generalised in various directions.
Firstly, the interaction formula (55) could be applied directly to any energy
functional that satisfies the restrictions that the potential is isotropic
close to the vacuum (10), and with a sufficiently large effective mass, (19).
In particular, a case we have not considered in this paper is the same range
of potentials (11), (12) but with general DMI that is not axisymmetric.
Secondly, we could extend to cases where the potential does not satisfy the
isotropy condition, as in (13), (14). The equations (48) and (51) would still
be true, but even the linearised Euler-Lagrange equations cannot be solved
exactly, only approximated [36]. Finally, the methods used can be applied to
any case where exponentially localised configurations interact, in a chiral
magnet or any similar medium. This includes skyrmion strings, Bloch points and
more. For this the formula (80) could be used with the necessary
modifications.

Separately, we used the study of linearised Euler-Lagrange equations of the
chiral magnet, which were necessary to calculate the interaction potential, to
investigate the magnetic soliton elliptical instability. The divergence of
decay lengthscale of solutions to these equations coincides with the uniform
state becoming linearly unstable, and we predict elliptical instability as we
approach this region. The predictions closely match numerically observed
elliptical instability when the applied field is tilted close to the plane.
Remarkably, this instability calculation combines with a completely different
method of estimating soliton instability, by finding the point at which domain
walls become energetically favoured, to overall produce a good fit to to
numerical observations at both small and large tilts of the applied field. It
is interesting that a seemingly single phenomenon can be explained only by a
combination of two such different methods. One can ask whether this reflects a
real distinction between two forms of elliptical instability, or alternatively
how these two methods are related, and if they can be combined.

## Acknowledgements

We thank Vlad Kuchkin for providing numerical data from [22].

#### Funding information

B B-S acknowledges an EPSRC-funded PhD studentship.

##  Appendix A Numerical method

Skyrmion and antiskyrmion solutions were generated by arrested Newton flow
[47] on a 400x400 lattice with periodic boundary conditions, with grid size
$dx=0.02$. This means that the DMI lengthscale $\frac{1}{k}$ was approximately
equal to $8dx$, and the size of the domain $400dx$ was approximately ten times
$\frac{1}{m}$, the decay length. Exact solutions from the critically coupled
model [31] with the appropriate charge were used as initial configurations.

To Fourier transform the angular dependence of $\tilde{\psi}$, we first must
approximate the location of the soliton centre. This was done by looking for
the lattice point where $\boldsymbol{n}$ was closest to $-\boldsymbol{e}_{3}$.
This point then defines the centre of a polar co-ordinate system $(r,\phi)$.
We then took a series of circles at successive values of $r$ and a range of
$\phi$ and found $\tilde{\psi}$ by interpolation, then calculated the integral
(42) for each $r$. The maximum radius at which this was done was half the
distance to the edge of the domain. It was found that beyond this the effects
of the finite size of the domain caused $\lvert c_{M}(r)\rvert/K_{M}(mr)$ to
again deviate from a constant value.

This method of finding the soliton centre is fairly crude, and by using more
sophisticated methods we can reduce the discrepancy between our measurement
for $\gamma_{M}$ and the value predicted by (41). So to account for this
discrepancy in general, we consider the error in the location of the soliton
centre, $\delta x_{0}$, in the calculations below. This leads to errors in the
polar co-ordinates $(\delta r=\delta x_{0},\delta\phi=\frac{\delta
x_{0}}{r})$. Note that these are not independent at each point. This leads to
an error in the argument of $\tilde{\psi}$, $\delta\arg(\tilde{\psi})=a\delta
x_{0}$, and thus an error in $\arg(c_{M}(r))$. Then in the integral (42), the
error in $\phi$ leads to an error in $\arg(c_{M}(r))$ equal to $M\delta
x_{0}/r$, while the error in $r$ leads to an error in $\lvert c_{M}(r)\rvert$
equal to $c_{M}^{\prime}(r)\delta x_{0}$. These are the errorbars plotted in
Figs. 1, 2 and 5.

To test the accuracy of this method, we simulated the skyrmion in normal
magnetic field, $\boldsymbol{D}_{i}=-2\pi\boldsymbol{e}_{i}$,
$h_{z}=0.8(2\pi)^{2}$, $h_{a}=0$, $\theta_{h}=0$. Because of the axisymmetry
of the skyrmion, we expect only $C_{1}\neq 0$. The largest $\lvert C_{M\neq
1}\rvert$ we found numerically was $0.07$, compared to a value of $\sim 7.3$
for $\lvert C_{1}\rvert$. We can also compare the value of $\lvert
C_{1}\rvert$ found to that for the hedgehog solution of the Euler-Lagrange
equations as found by shooting. In Fig. 9, both solutions are shown to be
well-fitted by the appropriate Bessel function $K_{1}(mr)$ beyond a certain
radius, but there is a discrepancy between them that leads to a $1.9\%$ error
in the calculated value of $\lvert C_{1}\rvert$, which we see is well
contained within the errorbars described above.

Figure 9: Magnitude of angular Fourier term $c_{1}(r)$ of the axisymmetric
skyrmion tail relative to the Bessel function $K_{1}(mr)$, for the solution
found by arrested Newton flow (PDE) and for the solution of the hedgehog
Euler-Lagrange equations found by shooting (ODE), for the chiral magnet energy
functional (1) with axisymmetric DMI (6) of strength $k=2\pi$ and axisymmetric
potential (11) $h_{z}=0.8(2\pi)^{2}$, $h_{a}=0$.

##  Appendix B Interaction potential for general energy functional

Here we discuss a general framework for calculating the interaction energy of
solitons in a magnetisation field (i.e. target space $S^{2}$) with an
arbitrary local energy functional, provided that solutions of the Euler-
Lagrange equations fall off exponentially. The argument is in fact more
general still. It applies for target spaces that are any Riemannian manifold,
meaning that we can for instance consider the target space to be $SO(3)$,
which can be relevant for certain classes of antiferromagnet. It also holds
for any number of spatial dimensions, which would for instance be useful if we
want to calculate the interaction strength of skyrmion strings. However, for
the following we keep the number of spatial dimensions to two and the target
space to $S^{2}$. The steps of this argument are not new [17, 23], but the
purpose of this discussion is to show its generality and that it does not
require a superposition ansatz. Throughout we think in terms of the infinite-
dimensional configuration space consisting of all possible configurations of
the magnetisation field.

As described in Sec. 3.1, we define the interaction potential $V^{AB}$ in
terms of the difference between the energy of a combined configuration
$\boldsymbol{n}^{AB}$ and the energy of the two isolated soliton
configurations $\boldsymbol{n}^{A}$, $\boldsymbol{n}^{B}$. In the simplest
case, we construct $\boldsymbol{n}^{AB}$ as a function of $\vec{R}$ by finding
the absolute minimum of the energy, within configurations of an appropriate
topological degree, subject to the constraint that
$\boldsymbol{n}^{AB}=-\boldsymbol{e}_{3}$ at $\vec{0}$ and $\vec{R}$. This
creates a 2-dimensional moduli space within the infinite-dimensional
configuration space, parametrised by $\vec{R}$. This is motivated by numerical
pinning procedures to approximate interaction energies, with which we make
direct comparison [22]. This moduli space also has a significance when the
Landau-Lifshitz equation is dominated by damping, so that in mathematical
terms fields undergo gradient flow on the configuration space. If we assume
that under this dynamics the soliton centres will move together or apart while
remaining unique, then this pinning procedure consists of foliating the
configuration space transversally to the path of the true dynamics and finding
the minimum on each sheet. In particular, this means that gradient flows that
start within our moduli space will remain on it, so that the path of gradient
flow on the moduli space is the path of the dynamics on the configuration
space and so the moduli space approach is exact, not an approximation.

There are other choices that we could make to construct $\boldsymbol{n}^{AB}$.
Firstly, our choice of soliton centre is somewhat arbitrary and we only choose
the point at which $\boldsymbol{n}=-\boldsymbol{e}_{3}$ for comparison to
numerics, as discussed in the main paper. A more natural choice might be the
point $\vec{x}_{0}$ at which
$\boldsymbol{n}^{AB}(\vec{x}_{0})=-\boldsymbol{n}_{0}$, as this continues to
be unique as we rotate the groundstate from $\boldsymbol{e}_{3}$ to
$-\boldsymbol{e}_{3}$. The choice of soliton centre can also be on a case-by-
case basis for different solitons. As discussed in the main paper, we assume
that in general we can modify the constraints to specify between all the
soliton configurations under consideration.

A separate approach, if we are interested in the property that the moduli
space contains gradient flows of the full energy, is to explicitly construct
it as such a space. To do this we extend our space of configurations to
include the solitons at infinite separation, and look at the (un)stable
manifold of this configuration if the attraction is repulsive (attractive)
[48]. Another approach would be to construct $\boldsymbol{n}^{AB}$ by ansatz,
namely some pointwise superposition of the two fields [17, 23], that satisfies
the properties we require below. This would give an upper bound on the
interaction energy as defined by the pinning method, if it were calculated
exactly, as in [49].

For the purposes of this appendix, the composite field $\boldsymbol{n}^{AB}$
only has to satisfy a few basic features in terms of the fields
$\boldsymbol{n}^{A}$, $\boldsymbol{n}^{B}$ and $\vec{R}$. Firstly, it should
attain $-\boldsymbol{e}_{3}$ at $\vec{0}$ and $\vec{R}$. Secondly, as
$R\to\infty$, $\boldsymbol{n}^{AB}\to\boldsymbol{n}^{A}$ around soliton $A$,
and $\boldsymbol{n}^{AB}\to\boldsymbol{n}^{B}$ around soliton $B$. Thirdly,
the field and its derivatives fall off exponentially away from the soliton
core. Fourthly, the fields approach linear superposition far from both soliton
cores, as discussed in the introduction. We can make these assumptions more
precise below once we have introduced some notation.

We define the region closer to $\vec{0}$ than $\vec{R}$ as $\sigma^{A}$, and
its complement as $\sigma^{B}$. Since the energy functional is local, we can
write the energy functional $E$ as the sum of the integrals over the two
regions, $E=E_{\sigma^{A}}+E_{\sigma^{B}}$. We can now split the calculation
for $V^{AB}$:

$V^{AB}=E_{\sigma^{A}}(\boldsymbol{n}^{AB})-E_{\sigma^{A}}(\boldsymbol{n}^{A})-E_{\sigma^{A}}(\boldsymbol{n}^{B})+E_{\sigma^{B}}(\boldsymbol{n}^{AB})-E_{\sigma^{B}}(\boldsymbol{n}^{A})-E_{\sigma^{B}}(\boldsymbol{n}^{B})$
(70)

To continue we call a map
$\boldsymbol{\epsilon}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{3}$ a tangent
vector to a given magnetisation field $\boldsymbol{n}$ in configuration space
if $\boldsymbol{\epsilon}(\vec{x})\cdot\boldsymbol{n}(\vec{x})=0$ for all
$\vec{x}$. Using the exponential map on the sphere defined in the main text
(8), we define the exponential map in configuration space [50]
$\exp_{\boldsymbol{n}}$ which takes a tangent vector $\boldsymbol{\epsilon}$
to $\boldsymbol{n}$, and turns it into a magnetisation field that is pointwise
the exponential map $\exp_{\boldsymbol{n}(x)}(\boldsymbol{\epsilon}(x))$:

$(\exp_{\boldsymbol{n}}(\boldsymbol{\epsilon}))(\vec{x})=\exp_{\boldsymbol{n}(\vec{x})}(\boldsymbol{\epsilon}(\vec{x}))$
(71)

In addition to $\boldsymbol{\psi}_{\boldsymbol{n}^{B}}$ as defined in the main
text (9), we can also define
$\boldsymbol{\epsilon}^{B}=\exp^{-1}_{\boldsymbol{n}^{A}}(\boldsymbol{n}^{AB})$.
The first describes the ‘tail’ of the isolated soliton $B$, while the second
describes the perturbation of the tail of soliton $B$ on soliton $A$. As shown
in 2.3, $\boldsymbol{\psi}_{\boldsymbol{n}^{B}}$ can be approximated by
solving the linearised Euler-Lagrange equations for the particular energy
functional we are working with. Meanwhile, we have no explicit way of finding
$\boldsymbol{\epsilon}^{B}$ but we can link it to
$\boldsymbol{\psi}_{\boldsymbol{n}^{B}}$ using the assumptions on
$\boldsymbol{n}^{AB}$.

Note that according to (9),
$\exp^{-1}_{\boldsymbol{n}_{0}}(-\boldsymbol{n}_{0})$ is not defined, so we
can only define $\boldsymbol{\psi}_{\boldsymbol{n}^{A,B}}(\vec{x})$,
$\boldsymbol{\epsilon}^{A,B}(\vec{x})$ when $\boldsymbol{n}(\vec{x})$ is never
$-\boldsymbol{n}_{0}$. However, $\boldsymbol{n}^{AB}(\vec{x})$ does indeed
reach $-\boldsymbol{n}_{0}$ near or at the soliton centre. This is not a
problem in practice: note that we always define
$\boldsymbol{\psi}_{\boldsymbol{n}^{A,B}}(\vec{x})$ and
$\boldsymbol{\epsilon}^{A,B}(\vec{x})$ below in integrals over regions away
from the corresponding soliton centre. In actuality we are then using the fact
that the integral is independent of the value of
$\boldsymbol{n}^{AB}(\vec{x})$ outside its domain to replace
$\boldsymbol{n}^{AB}(\vec{x})$ with a field that is identical within the
domain of integration but outside the integral never reaches
$-\boldsymbol{n}_{0}$, and then defining $\boldsymbol{\epsilon}^{B}$,
$\boldsymbol{\psi}_{\boldsymbol{n}^{B}}$ in terms of the inverse exponential
function on that field. We just skip this cumbersome notation.

We can now state our second assumption on $\boldsymbol{n}^{AB}$ more
explicitly: as we go far from the centre of soliton $B$,
$\lvert\boldsymbol{\epsilon}^{B}(\vec{x})\rvert\sim\lvert\boldsymbol{\psi}_{\boldsymbol{n}^{B}}(\vec{x})\rvert$,
and thus $\lvert\boldsymbol{\epsilon}^{B}(\vec{x})\rvert\to 0$ as $R\to\infty$
for all $\vec{x}$ in region $\sigma^{A}$, and vice versa. We can also
formalise our assumption of linear superposition: in the region far from both
solitons,
$\boldsymbol{\epsilon}^{A,B}\to\boldsymbol{\psi}_{\boldsymbol{n}^{A,B}}$.
Finally, we have the assumption that the tails and their derivatives are
bounded by exponential decay:
$\boldsymbol{\psi}_{\boldsymbol{n}^{A,B}}(\vec{x})=O(e^{-mr})$,
$\partial_{i}\boldsymbol{\psi}_{\boldsymbol{n}^{A,B}}(\vec{x})=O(e^{-mr})$,
and so on. This last assumption is satisfied by solitons in the chiral magnet,
see Equations (29), (32), but in this appendix the discussion is more general
and $m$ can be taken to generally describe the inverse decay lengthscale of
whatever soliton is under consideration.

We define the variation of the energy functional:

$\delta_{\boldsymbol{\epsilon}}E(\boldsymbol{n})=\frac{d}{dt}\biggr{\rvert}_{t=0}(E(\exp_{\boldsymbol{n}}(t\boldsymbol{\epsilon})),$
(72)

and the second variation:

$\delta^{2}_{\boldsymbol{\epsilon}^{\prime},\boldsymbol{\epsilon}}E(\boldsymbol{n})=\frac{d}{dt}\biggr{\rvert}_{t=0}\delta_{\boldsymbol{\epsilon}}E(\exp_{\boldsymbol{n}}(t\boldsymbol{\epsilon}^{\prime})),$
(73)

then we can Taylor expand all three terms:

$E_{\sigma^{A}}(\boldsymbol{n}^{AB})-E_{\sigma^{A}}(\boldsymbol{n}^{A})=\delta_{\boldsymbol{\epsilon}^{B}}E_{\sigma^{A}}(\boldsymbol{n}^{A})+\frac{1}{2}\delta^{2}_{\boldsymbol{\epsilon}^{B},\boldsymbol{\epsilon}^{B}}E_{\sigma^{A}}(\boldsymbol{n}^{A})+O(\lvert\boldsymbol{\epsilon}^{B}\rvert^{3})$
(74) $\displaystyle E_{\sigma^{A}}(\boldsymbol{n}^{B})$
$\displaystyle=E_{\sigma^{A}}(\boldsymbol{n}_{0})+\delta_{\boldsymbol{\psi}_{\boldsymbol{n}^{B}}}E_{\sigma^{A}}(\boldsymbol{n}_{0})+\frac{1}{2}\delta^{2}_{\boldsymbol{\psi}_{\boldsymbol{n}^{B}},\boldsymbol{\psi}_{\boldsymbol{n}^{B}}}E_{\sigma^{A}}(\boldsymbol{n}_{0}))+O(\lvert\boldsymbol{\psi}_{\boldsymbol{n}^{B}}\rvert^{3})$
$\displaystyle=\delta_{\boldsymbol{\psi}_{\boldsymbol{n}^{B}}}E_{\sigma^{A}}(\boldsymbol{n}_{0})+\frac{1}{2}\delta^{2}_{\boldsymbol{\psi}_{\boldsymbol{n}^{B}},\boldsymbol{\psi}_{\boldsymbol{n}^{B}}}E_{\sigma^{A}}(\boldsymbol{n}_{0})+O(\lvert\boldsymbol{\psi}_{\boldsymbol{n}^{B}}\rvert^{3})$
(75)

where $\lvert\boldsymbol{\epsilon}\rvert$ is the maximum value of
$\lvert\boldsymbol{\epsilon}(x)\rvert$ over the whole region $\sigma^{A}$, in
practice its value at the midpoint $\vec{x}=\frac{1}{2}\vec{R}$. We use our
assumption that
$\lvert\boldsymbol{\epsilon}\rvert\sim\lvert\boldsymbol{\psi}\rvert$ to
replace $O(\lvert\boldsymbol{\epsilon}^{B}\rvert^{3})$ with
$O(\lvert\boldsymbol{\psi}_{\boldsymbol{n}^{B}}\rvert^{3})$, and so on.

This is where the exponential decay assumption becomes useful: since we are
expanding in terms of the maximum value of
$\lvert\boldsymbol{\psi}_{\boldsymbol{n}^{B}}(x)\rvert$ over the domain of
integration, this expansion only makes sense if in general the integral of the
function is bounded by its maximum value, and whatever derivative operators
may act within the integrals that we throw away do not change the order of
that maximum value. Both of these things are true for exponentially decaying
functions, and thus true for functions bounded below exponential decay with
all derivatives bounded below exponential decay. We will use this several
times more below to bound integrals in terms of the maximum value of a field
over the domain of integration. At this point we can bound the term we are
neglecting in terms of $\max_{\sigma^{A}}{e^{-3mr}}=e^{-\frac{3}{2}mR}$.

For what follows it is useful to separate the first variation into two
contributions. By integration by parts, a variation can always be separated
into a bulk integral where the variation field appears without derivatives,
and an integral along $\partial\sigma$ which we call
$\partial\delta_{\boldsymbol{\epsilon}}E_{\sigma}(\boldsymbol{n})$:

$\delta_{\boldsymbol{\epsilon}}E_{\sigma}(\boldsymbol{n})=\int_{\sigma}\boldsymbol{\epsilon}\cdot\mathbf{f}(\boldsymbol{n},\partial\boldsymbol{n},\ldots)d^{2}x+\partial\delta_{\boldsymbol{\epsilon}}E_{\sigma}(\boldsymbol{n}).$
(76)

Because $\boldsymbol{n}^{A}$ is a minimizer of the full energy
$E(\boldsymbol{n})$, $\boldsymbol{n}^{A}$ solves the Euler-Lagrange equations,
$\boldsymbol{n}\times\mathbf{f}(\boldsymbol{n},\partial\boldsymbol{n},\ldots)=\mathbf{0}$,
and therefore the term
$\delta_{\boldsymbol{\epsilon}^{B}}E_{\sigma^{A}}(\boldsymbol{n}^{A})$ can
only be a boundary term:

$\delta_{\boldsymbol{\epsilon}^{B}}E_{\sigma^{A}}(\boldsymbol{n}^{A})=\partial\delta_{\boldsymbol{\epsilon}^{B}}E_{\sigma^{A}}(\boldsymbol{n}^{A})$
(77)

Similarly,
$\delta_{\boldsymbol{\psi}_{\boldsymbol{n}^{B}}}E_{\sigma^{A}}(\boldsymbol{n}_{0})=\partial\delta_{\boldsymbol{\psi}_{\boldsymbol{n}^{B}}}E_{\sigma^{A}}(\boldsymbol{n}_{0})$.
Now because
$\partial\delta_{\boldsymbol{\epsilon}^{B}}E_{\sigma^{A}}(\boldsymbol{n}^{A})$
depends only on the field $\boldsymbol{n}^{A}$ evaluated along
$\partial\sigma^{A}$, which is far from $\vec{0}$, we can Taylor expand the
whole expression around $\boldsymbol{n}^{A}=\boldsymbol{n}_{0}$, using the
fact that
$\boldsymbol{\epsilon}^{B}\to\boldsymbol{\psi}_{\boldsymbol{n}^{B}}$:

$\partial\delta_{\boldsymbol{\epsilon}^{B}}E_{\sigma^{A}}(\boldsymbol{n}^{A})=\partial\delta_{\boldsymbol{\psi}_{\boldsymbol{n}^{B}}}E_{\sigma^{A}}(\boldsymbol{n}_{0})+\delta_{\boldsymbol{\psi}_{\boldsymbol{n}^{A}}}\partial\delta_{\boldsymbol{\psi}_{\boldsymbol{n}^{B}}}E_{\sigma^{A}}(\boldsymbol{n}_{0})+O(e^{-\frac{3}{2}mR})$
(78)

This means that

$V^{AB}=\delta_{\boldsymbol{\psi}_{\boldsymbol{n}^{A}}}\partial\delta_{\boldsymbol{\psi}_{\boldsymbol{n}^{B}}}E_{\sigma^{A}}(\boldsymbol{n}_{0})+\frac{1}{2}\delta^{2}_{\boldsymbol{\epsilon}^{B},\boldsymbol{\epsilon}^{B}}E_{\sigma^{A}}(\boldsymbol{n}^{A})-\frac{1}{2}\delta^{2}_{\boldsymbol{\psi}_{\boldsymbol{n}^{B}},\boldsymbol{\psi}_{\boldsymbol{n}^{B}}}E_{\sigma^{A}}(\boldsymbol{n}_{0})+(A\leftrightarrow
B)+O(e^{-\frac{3}{2}mR})$ (79)

Outside the soliton core, we can expand
$\frac{1}{2}\delta^{2}_{\boldsymbol{\epsilon}^{B},\boldsymbol{\epsilon}^{B}}E_{\sigma^{A}}(\boldsymbol{n}^{A})$
around $\boldsymbol{n}_{0}$, and the latter two terms will cancel to order
$\lvert\boldsymbol{\epsilon}^{B}\rvert^{2}\lvert\boldsymbol{\psi}_{\boldsymbol{n}^{A}}\rvert^{\prime}$,
where $\lvert\boldsymbol{\psi}_{\boldsymbol{n}^{A}}\rvert^{\prime}$ is the
maximum value of
$\lvert\boldsymbol{\psi}_{\boldsymbol{n}^{A}}(\vec{x})\rvert^{\prime}$ over
everywhere except the soliton core. This is still subleading. Inside the
soliton core, both terms are order
$\lvert\boldsymbol{\epsilon}^{B}\rvert^{\prime 2}$, where
$\lvert\boldsymbol{\epsilon}^{B}\rvert^{\prime}$ is the maximum value of
$\lvert\boldsymbol{\epsilon}^{B}(\vec{x})\rvert$ over just the soliton core.
Because of our assumption that $\boldsymbol{n}$ approaches
$\boldsymbol{n}^{A}$ in $\sigma^{A}$ as $R$ increases, the size of the soliton
core approaches a constant value, so it can be arbitrarily small in comparison
to $R$, and so in particular the distance between the centre of soliton $B$
and the closest edge of the core of soliton $A$ can be larger than
$\frac{3R}{4}$ for large enough $R$. Then
$\lvert\boldsymbol{\epsilon}^{B}\rvert^{\prime
2}\sim\lvert\boldsymbol{\epsilon}^{B}\rvert^{3}=O\left(e^{-\frac{3}{2}mR}\right)$
and we get our final expression for the general interaction energy at large
separation:

$V^{AB}=\delta_{\boldsymbol{\psi}_{\boldsymbol{n}^{A}}}\partial\delta_{\boldsymbol{\psi}_{\boldsymbol{n}^{B}}}E_{\sigma^{A}}(\boldsymbol{n}_{0})-\delta_{\boldsymbol{\psi}_{\boldsymbol{n}^{B}}}\partial\delta_{\boldsymbol{\psi}_{\boldsymbol{n}^{A}}}E_{\sigma^{A}}(\boldsymbol{n}_{0})+O\left(e^{-\frac{3}{2}mR}\right),$
(80)

where we use the fact that $\partial\sigma^{B}$ is equal to
$\partial\sigma^{A}$ with the opposite orientation.

We find in the main text that $V^{AB}\sim\frac{e^{-mR}}{\sqrt{mR}}$, so the
correction term is indeed subleading. By repeating the above calculations in a
more geometric language, we can view the term
$\partial\delta_{\boldsymbol{\epsilon}}E_{\sigma^{A}}$ as a one-form on the
configuration space of magnetisation fields, acting on the vector
$\boldsymbol{\epsilon}$,
$\partial\delta_{\boldsymbol{\epsilon}}E_{\sigma^{A}}(\boldsymbol{n}^{A})=\omega_{\boldsymbol{n}^{A}}(\boldsymbol{\epsilon})$.
Then our final expression can be seen to be the external derivative of this
one-form acting on the two tail vector fields,
$V^{AB}=d_{\boldsymbol{n}_{0}}\omega(\boldsymbol{\psi}_{\boldsymbol{n}^{A}},\boldsymbol{\psi}_{\boldsymbol{n}^{B}})$.

We can apply this formula to the chiral magnet energy functional (1). First we
find the boundary term from the first variation of the energy:

$\partial\delta_{\boldsymbol{\psi}}E_{\sigma^{A}}(\boldsymbol{n})=\int_{\partial\sigma^{A}}\boldsymbol{\psi}\cdot(\partial_{i}\boldsymbol{n}+\boldsymbol{D}_{i}\times\boldsymbol{n})dS^{i}$
(81)

then we vary this term with respect to a second field:

$\delta_{\boldsymbol{\psi}^{\prime}}\partial\delta_{\boldsymbol{\psi}}E_{\sigma^{A}}(\boldsymbol{n}_{0})=\int_{\partial\sigma^{A}}\boldsymbol{\psi}\cdot(\partial_{i}\boldsymbol{\psi}^{\prime}+\boldsymbol{D}_{i}\times\boldsymbol{\psi}^{\prime})dS^{i}$
(82)

Note that the addition of a boundary term to the energy like
$-\boldsymbol{D}_{i}\cdot(\boldsymbol{n}_{0}\times\partial_{i}\boldsymbol{n})$,
which can be motivated physically and in terms of analysis [30, 32], does not
affect this quantity and thus does not enter into the interaction potential.
Also, while for simplicity we chose $\partial\sigma^{A}$ to be the straight
line equidistant from the soliton centres, the derivation above only
fundamentally depends on the fact that $\partial\sigma^{A}$ divides the plane
into two halves with a soliton core in each half. We see this independence of
the exact boundary in the main section when we calculate the interaction
potential for the chiral magnet specifically. Finally we substitute this into
(80) to find (47).

## References

  * [1] A. Bogdanov and A. Hubert,  _Thermodynamically stable magnetic vortex states in magnetic crystals_ ,  Journal of Magnetism and Magnetic Materials 138(3), 255 (1994),  https://doi.org/10.1016/0304-8853(94)90046-9.
  * [2] N. Manton and P. Sutcliffe,  _Topological Solitons_ ,  Cambridge Monographs on Mathematical Physics. Cambridge University Press,  10.1017/CBO9780511617034 (2004).
  * [3] N. Manton,  _A remark on the scattering of bps monopoles_ ,  Physics Letters B 110(1), 54 (1982),  https://doi.org/10.1016/0370-2693(82)90950-9.
  * [4] A. A. Thiele,  _Steady-state motion of magnetic domains_ ,  Phys. Rev. Lett. 30, 230 (1973),  10.1103/PhysRevLett.30.230.
  * [5] K. Everschor, M. Garst, R. A. Duine and A. Rosch,  _Current-induced rotational torques in the skyrmion lattice phase of chiral magnets_ ,  Phys. Rev. B 84, 064401 (2011),  10.1103/PhysRevB.84.064401.
  * [6] O. A. Tretiakov, D. Clarke, G.-W. Chern, Y. B. Bazaliy and O. Tchernyshyov,  _Dynamics of domain walls in magnetic nanostrips_ ,  Phys. Rev. Lett. 100, 127204 (2008),  10.1103/PhysRevLett.100.127204.
  * [7] J. Müller and A. Rosch,  _Capturing of a magnetic skyrmion with a hole_ ,  Phys. Rev. B 91, 054410 (2015),  10.1103/PhysRevB.91.054410.
  * [8] D. M. A. Stuart,  _Analysis of the adiabatic limit for solitons in classical field theory_ ,  Proc. Roy. Soc. Lond. A 463, 2753 (2007),  10.1098/rspa.2007.0130.
  * [9] S. Demoulini and D. Stuart,  _Adiabatic limit and the slow motion of vortices in a chern-simons-schrodinger system_ ,  Commun. Math. Phys. 290, 597 (2009),  10.1007/s00220-009-0844-y,  0905.0998.
  * [10] A. Hubert and A. Bogdanov,  _The properties of isolated magnetic vortices_ ,  physica status solidi (b) 186 (1994),  10.1002/pssb.2221860223.
  * [11] L. LANDAU and E. LIFSHITZ,  _3 - on the theory of the dispersion of magnetic permeability in ferromagnetic bodies reprinted from physikalische zeitschrift der sowjetunion 8, part 2, 153, 1935._ ,  In L. PITAEVSKI, ed., _Perspectives in Theoretical Physics_ , pp. 51–65. Pergamon, Amsterdam,  ISBN 978-0-08-036364-6,  https://doi.org/10.1016/B978-0-08-036364-6.50008-9 (1992).
  * [12] T. L. Gilbert,  _A phenomenological theory of damping in ferromagnetic materials_ ,  IEEE Transactions on Magnetics 40, 3443 (2004).
  * [13] N. Manton,  _The force between ’t hooft-polyakov monopoles_ ,  Nuclear Physics B 126(3), 525 (1977),  https://doi.org/10.1016/0550-3213(77)90294-2.
  * [14] T. Skyrme,  _A unified field theory of mesons and baryons_ ,  Nuclear Physics 31, 556 (1962),  https://doi.org/10.1016/0029-5582(62)90775-7.
  * [15] B. Schroers,  _Dynamics of moving and spinning skyrmions_ ,  Zeitschrift für Physik C Particles and Fields 61(3), 479 (1994),  10.1007/BF01413188.
  * [16] N. Manton, B. Schroers and M. Singer,  _The interaction energy of well-separated skyrme solitons_ ,  Communications in Mathematical Physics 245, 123 (2004),  10.1007/s00220-003-1006-2.
  * [17] B. Piette, B. Schroers and W. Zakrzewski,  _Multisolitons in a two-dimensional skyrme model_ ,  Zeitschrift für Physik C Particles and Fields 65 (1994),  10.1007/BF01571317.
  * [18] J. M. Speight,  _Static intervortex forces_ ,  Phys. Rev. D 55, 3830 (1997),  10.1103/PhysRevD.55.3830,  hep-th/9603155.
  * [19] A. Bogdanov,  _New localized solutions of the nonlinear field equations_ ,  ZhETF Pisma Redaktsiiu 62, 231 (1995).
  * [20] S.-Z. Lin and S. Hayami,  _Ginzburg-landau theory for skyrmions in inversion-symmetric magnets with competing interactions_ ,  Phys. Rev. B 93, 064430 (2016),  10.1103/PhysRevB.93.064430.
  * [21] D. Foster, C. Kind, P. J. Ackerman, J.-S. B. Tai, M. R. Dennis and I. I. Smalyukh,  _Two-dimensional skyrmion bags in liquid crystals and ferromagnets_ ,  Nature Physics 15(7) (2019),  10.1038/s41567-019-0476-x.
  * [22] V. M. Kuchkin and N. S. Kiselev,  _Turning a chiral skyrmion inside out_ ,  Phys. Rev. B 101, 064408 (2020),  10.1103/PhysRevB.101.064408.
  * [23] M. Kameda, R. Koyama, T. Nakajima and Y. Kawaguchi,  _Controllable interskyrmion attractive interactions and resulting skyrmion-lattice structures in two-dimensional chiral magnets with in-plane anisotropy_ ,  Phys. Rev. B 104, 174446 (2021),  10.1103/PhysRevB.104.174446.
  * [24] D. Feist,  _Interactions of b = 4 skyrmions_ ,  Journal of High Energy Physics 2012 (2011),  10.1007/JHEP02(2012)100.
  * [25] F. N. Rybakov and N. S. Kiselev,  _Chiral magnetic skyrmions with arbitrary topological charge_ ,  Phys. Rev. B 99, 064437 (2019),  10.1103/PhysRevB.99.064437.
  * [26] V. M. Kuchkin, B. Barton-Singer, F. N. Rybakov, S. Blügel, B. J. Schroers and N. S. Kiselev,  _Magnetic skyrmions, chiral kinks, and holomorphic functions_ ,  Phys. Rev. B 102, 144422 (2020),  10.1103/PhysRevB.102.144422.
  * [27] I. Dzyaloshinsky,  _A thermodynamic theory of “weak” ferromagnetism of antiferromagnetics_ ,  Journal of Physics and Chemistry of Solids 4(4), 241 (1958),  10.1016/0022-3697(58)90076-3.
  * [28] T. Moriya,  _Anisotropic superexchange interaction and weak ferromagnetism_ ,  Phys. Rev. 120, 91 (1960),  10.1103/PhysRev.120.91.
  * [29] U. Güngördü, R. Nepal, O. A. Tretiakov, K. Belashchenko and A. Kovalev,  _Stability of skyrmion lattices and symmetries of quasi-two-dimensional chiral magnets_ ,  Physical Review B 93, 064428 (2016),  10.1103/PhysRevB.93.064428.
  * [30] C. Melcher,  _Chiral skyrmions in the plane_ ,  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 470(2172), 20140394 (2014),  10.1098/rspa.2014.0394,  https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.2014.0394.
  * [31] B. Barton-Singer, C. Ross and B. J. Schroers,  _Magnetic skyrmions at critical coupling_ ,  Communications in Mathematical Physics 375(3), 2259 (2020).
  * [32] B. Schroers,  _Gauged sigma models and magnetic skyrmions_ ,  SciPost Physics 7(3), 030 (2019).
  * [33] H. Georgi and S. L. Glashow,  _Unified weak and electromagnetic interactions without neutral currents_ ,  Phys. Rev. Lett. 28, 1494 (1972),  10.1103/PhysRevLett.28.1494.
  * [34] G. Hooft,  _Magnetic monopoles in unified gauge theories_ ,  Nuclear Physics B 79(2), 276 (1974),  https://doi.org/10.1016/0550-3213(74)90486-6.
  * [35] A. M. Polyakov,  _Particle spectrum in the quantum field theory_ ,  JETP Lett. 20, 194 (1974).
  * [36] B. Barton-Singer,  _Analytical predictions of skyrmion stability_ ,  Paper in preparation.
  * [37] N. M. Temme,  _Bessel Functions_ , chap. 9, pp. 219–255,  John Wiley & Sons, Ltd,  ISBN 9781118032572,  https://doi.org/10.1002/9781118032572.ch9 (1996), https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781118032572.ch9.
  * [38] C. Joslin and C. Gray,  _Multipole expansions in two dimensions_ ,  Molecular Physics 50(2), 329 (1983).
  * [39] A. Bogdanov and A. Hubert,  _The stability of vortex-like structures in uniaxial ferromagnets_ ,  Journal of Magnetism and Magnetic Materials 195(1), 182 (1999),  https://doi.org/10.1016/S0304-8853(98)01038-5.
  * [40] A. N. Bogdanov and D. Yablonskii,  _Thermodynamically stable “vortices” in magnetically ordered crystals. the mixed state of magnets_ ,  Zh. Eksp. Teor. Fiz 95(1), 178 (1989).
  * [41] R. Palais,  _The principle of symmetric criticality_ ,  Communications in Mathematical Physics 69 (1979),  10.1007/BF01941322.
  * [42] S.-Z. Lin, C. Reichhardt, C. D. Batista and A. Saxena,  _Particle model for skyrmions in metallic chiral magnets: Dynamics, pinning, and creep_ ,  Phys. Rev. B 87, 214419 (2013),  10.1103/PhysRevB.87.214419.
  * [43] P. Roget,  _V. explanation of an optical deception in the appearance of the spokes of a wheel seen through vertical apertures_ ,  Philosophical Transactions of the Royal Society of London pp. 131 – 140 (1824).
  * [44] B. Gregorcic and G. Planinsic,  _Why do photo finish images look weird?_ ,  Physics Education 47(5), 530 (2012),  10.1088/0031-9120/47/5/530.
  * [45] D. Walsh,  _Playing detective with rolling shutter photos_ ,  Accessed on 30.9.21 (2013).
  * [46] V. M. Kuchkin, K. Chichay, B. Barton-Singer, F. N. Rybakov, S. Blügel, B. J. Schroers and N. S. Kiselev,  _Geometry and symmetry in skyrmion dynamics_ ,  Phys. Rev. B 104(16), 165116 (2021).
  * [47] S. B. Gudnason and J. M. Speight,  _Realistic classical binding energies in the $\omega$-skyrme model_,  Journal of High Energy Physics 2020(7), 1 (2020).
  * [48] N. S. Manton,  _Unstable manifolds and soliton dynamics_ ,  Phys. Rev. Lett. 60, 1916 (1988),  10.1103/PhysRevLett.60.1916.
  * [49] C. Ross, N. Sakai and M. Nitta,  _Skyrmion interactions and lattices in chiral magnets: analytical results_ ,  Journal of High Energy Physics 2021 (2021),  10.1007/JHEP02(2021)095.
  * [50] J. Eells,  _On the geometry of function spaces_ ,  In _Symposium internacional de topologıa algebraica International symposium on algebraic topology_ , pp. 303–308. Citeseer (1958).