Skyrmion interactions and lattices in chiral magnets: analytical results

We study two-body interactions of magnetic skyrmions on the plane and apply them to a (mostly) analytic description of a skyrmion lattice. This is done in the context of the solvable line, a particular choice of a potential for magnetic anisotropy and Zeeman terms, where analytic expressions for skyrmions are available. The energy of these analytic single skyrmion solutions is found to become negative below a critical point, where the ferromagnetic state is no longer the lowest energy state. This critical value is determined exactly without the ambiguities of numerical simulations. Along the solvable line the interaction energy for a pair of skyrmions is repulsive with power law fall off in contrast to the exponential decay of a purely Zeeman potential term. Using the interaction energy expressions we construct an inhomogeneous skyrmion lattice state, which is a candidate ground states for the model in particular parameter regions. Finally we estimate the transition between the skyrmion lattice and an inhomogeneous spiral state.


Introduction
Magnetic skyrmions are topologically non-trivial solitons which occur in certain magnetic materials [1]. They were first predicted in [2] as minimum energy configurations of the magnetisation vector field and have been observed in a variety of magnetic materials [3][4][5]. In contrast to the case of baby skyrmions -field theoretical analogues of magnetic skyrmions stabilised by higher derivative terms in the energy [6,7] magnetic skyrmions are stabilised by a first order Dzyaloshinskii-Moriya (DM) interaction term (with coefficient κ) [8,9]. The fact that the DM term can give a negative contribution to the energy avoids a contradiction with the usual Derrick scaling argument [7]. A systematic mathematical study of magnetic skyrmions in the presence of an external magnetic field and without anisotropy is given in [10]. The various phases of a chiral magnet have been explored with easy plane anisotropy in addition to an external magnetic field [11][12][13][14]. The skyrmion lattice phase has been studied extensively in [15,16]. A model of a skyrmion lattice made up of non

JHEP02(2021)095
The interaction energy between a pair of skyrmions was discussed for a specific easyaxis potential in [32] and has recently been discussed [33][34][35], using the dipole approximation introduced in ref. [6]. Most previous work, including refs. [34,35], only considers the case of an external magnetic field without an anisotropy term. It has been found that the interaction energy decays exponentially at large separations in this case [36]. The case of easy plane anisotropy being more important than the Zeeman term was considered in ref. [37] where an attractive interaction was found between individual skyrmions. The skyrmion interaction energy for the solvable line does not appear to have been studied before. The boundary contributions to the energy caused by the slow, power law, decay of the skyrmion configurations complicate the asymptotic approach used for B > 2A. This necessitates a slightly different approach which is possible due to the presence of analytic skyrmion solutions.
In this paper we conduct a semi-analytic study of interactions between magnetic skyrmions along the solvable line, and apply it to construct a model of a skyrmion lattice. To achieve this, we make use of the dipole approximation from [6]. This approximation involves working with a superposition of two skyrmions, and is reliable as long as the skyrmions are well separated, becoming less reliable as the distinction between the constituent skyrmions becomes less pronounced. The accuracy of the approximation should be measured in terms of the separation between two skyrmions divided by the diameter of a single skyrmion. The diameter of the skyrmion is determined by three dimensionful parameters, the DM interaction strength κ, and the two parameters A (for magnetic anisotropy) and B (for a Zeeman field) in the potential. From the numerical study of the phase diagram of a chiral magnet with easy plane anisotropy, [11,12], it is known that large values of A and B correspond to homogeneous ferromagnetic phases. The solvable line, 2A = B, lies along the phase boundary between the polarised ferromagnetic phase and the spin canted ferromagnetic phase. At large values of the potential along the solvable line, magnetic skyrmions are obtained as positive energy excited states on top of the ferromagnetic ground state. The skyrmions have smaller radii as the potential parameter A becomes larger. Thus the dipole approximation, at a fixed separation between two isolated skyrmions, works best as the potential parameter A becomes large. For small values of A, on the other hand, the skyrmions along the solvable line have negative energy and the homogeneous ferromagnetic state becomes unstable. This leads to the possibility of the ground state being an inhomogeneous skyrmion lattice in agreement with numerical studies [15,16]. The other possible inhomogeneous state is a spiral or helical state [18,19] where the configuration is effectively one dimensional.
We apply the dipole approximation at smaller values of A to understand some of the qualitative behaviour of the triangular skyrmion lattice phase from [11] along the solvable line. We take two exact single skyrmion solutions and consider their superposition as a function of their separation. Our results show that the dipole approximation can be used reliably to find the transition between the ferromagnetic phase and the lattice phase. However, the dipole approximation breaks down before we can study the phase transition between the lattice phase and the spiral phase. We find evidence that a model of a skyrmion lattice as, effectively, non interacting closely packed skyrmions [17] becomes JHEP02(2021)095 a better approximation. Using this we estimate the critical value of the external magnetic field where the transition to the spiral phase occurs.
The key feature of the present work is that by working on the solvable line there are exact expressions for the skyrmion configurations enabling us to make exact statements about the energy of a single skyrmion. This makes it possible for us to identify a specific value of the anisotropy below which the single skyrmion solution has negative energy and we have a state lower in energy than the ferromagnetic state. Combining the exact single skyrmion configurations with the dipole approximation we arrive at an analytic expression for the interaction energy density of the superposition without the need to work asymptotically. After numerical integration we find that the interaction energy is positive and has a power law decay in contrast to the exponential decay of the pure Zeeman potential [34][35][36]. An exponentially decaying interaction was also found in [32] for the case of a specific easyaxis potential. Comparing these previous studies with our appendix A, we see that the interaction energy will exponentially decay when the skyrmion field decays exponentially. Two-body and three-body skyrmion interactions in the absence of an anisotropy term in the potential have been studied in [38]. We show that the power law decay of the skyrmion fields necessitates the inclusion of a boundary term in the energy functional in order to have a well defined variational problem when deriving the equations of motion. Along the solvable line, we find a repulsive interaction. This paper is organized as follows. In section 2, we review solvable models for magnetic skyrmions and clarify the necessity of a boundary term. In section 3, we calculate the interaction energy and force between two isolated skyrmions. In section 4, the interaction energy is applied to describe a skyrmion lattice, and the phase transition to the spiral phase is estimated. Section 5 is devoted to a summary and discussion. Appendix A gives the details of the boundary term in the energy functional and its relation to the asymptotic behaviour of the magnetisation vector field. In appendix B we express the interaction energy as a function of dimensionless variables. Spiral states are described in appendix C.

Solvable models for magnetic skyrmions
In this section, we review solvable models for magnetic skyrmions. In subsection 2.1, we introduce the model and classify homogeneous ground states, when the DM term can be neglected. Then in subsection 2.2, we present exact skyrmion solutions. Finally in subsection 2.3, we discuss the instability of the ferromagnetic phase.

The model
This subsection gives a concise summary of previously known results that are directly relevant for our study and clarifies our conventions at the same time. We study chiral magnets in two spatial dimensions described in terms of a real three component magnetisation vector field n : R 2 → S 2 , constrained as n 2 = 1. In works such as [2,11,15] the following energy functional is considered where the potential U consists of a Zeeman term with the magnetic field B as coefficient, and an anisotropy term with the coefficient A The term corresponding to the parameter κ is the Dzyaloshinskii-Moriya (DM) interaction energy [8,9], and the rotated derivative is defined as where the parameter α ∈ [0, 2π) denotes a one parameter family of models corresponding to different types of DM term, for example α = 0 corresponds to the Bloch type DM and α = π 2 corresponds to the Néel type DM interaction. The subscript "bulk" indicates that we need to add a boundary term to properly define the DM interaction term (to apply the variational principle), as we describe subsequently.
For large values of the potential parameters |A|, |B| κ, for which the DM term can be neglected, the ground state is determined by the minimum of the potential. There are three types of minima of the general potential U (n 3 ) depending on values of A, B:

JHEP02(2021)095
Let us note that the model is invariant under the transformation We consider the parameter region B ≥ 0 in the following, since we can obtain identical result for B < 0 by just changing the sign of n 3 . If we tune the parameters A, B in the potential as (and add a constant A), we obtain a potential V ( We call this choice of parameters, B = 2A, the solvable line, since explicit skyrmion configurations can be constructed as shown in the next subsection. Instead of the constrained vector field n with n 2 = 1, it is often more convenient to use the unconstrained complex field w defined by the stereographic projection The energy functional in eq. (2.1) with the potential replaced by V in eq. (2.6) is rewritten as (2.8)

Exact single skyrmion solutions
Defining isolated skyrmion states requires knowledge of the minima of the potential. This is because the magnetisation vector of the skyrmion configuration is required to go to the ground state asymptotically.
In the case of derivative interactions like the DM term, particular care is needed to properly define the energy functional. When varying fields to obtain field equations, we need to integrate by parts, which produces a surface term. If the fields do not decay sufficiently fast asymptotically, the surface term obstructs the variational principle used to obtain the field equation. Hence it is mandatory to subtract a boundary term to remove this obstruction and allow us to apply the variational principle. In the case of the DM interaction, this is accomplished by adding the following boundary term [30,31] to the energy functional (2.1) Here n ∞ is the asymptotic value of n determined by the homogeneous state. Explicitly, for the ferromagnetic phases, n ∞ = ±e 3 with the top sign when B is positive and the lower sign when B is negative.

JHEP02(2021)095
We wish to stress that only by adding the boundary term, does the DM term become well-defined in terms of the variational problem for the energy functional. More details about the necessity of the boundary term are discussed in appendix A. Some properties of the boundary term in eq. (2.9), are discussed in refs. [10,20].
Taking into account the boundary term in eq. (2.9) and the potential in eq. (2.6), the model of a chiral magnet that we consider is The tuning of the potential allows one to obtain exact solutions of a hedgehog type [20]. In terms of the polar coordinates x 1 = r cos ϕ, x 2 = r sin ϕ, one finds the exact solution by imposing a boundary condition n → (0, 0, 1) at infinity r → ∞, with phase γ = π 2 − α and profile function The hedgehog solution is called a skyrmion. It gives a mapping from the base space R 2 to the target space S 2 . As n tends to a constant as r → ∞ we can regard R 2 as S 2 by adding the point at infinity, thus we can define the degree Q corresponding to π 2 (S 2 ) This topological charge is Q = −1 for a skyrmion configuration. It is worth repeating that in this context a skyrmion has topological charge Q = −1 while an anti-skyrmion has topological charge Q = +1. This is in contrast to the case of baby skyrmions as discussed in [7]. The skyrmion solution has energy (2.14) Here, when calculating the energy of the skyrmion, taking into account the boundary term (2.9) is crucial. This is an exact expression for the skyrmion energy as a function of A κ 2 along the solvable line. The analytic skyrmion solutions, eqs. (2.11) and (2.12), along the solvable line enable us to compute this exact energy expression. As highlighted in table 1, subtracting the boundary term leads to the skyrmion having negative energy which is very important for our approach to constructing the skyrmion lattice. The profile function in eq. (2.12) shows that n 3 = 0 (the magnetisation vector lies in the x − y plane) at which is defined as the radius of a skyrmion.

JHEP02(2021)095
Skyrmion BPS skyrmion BPS anti-skyrmion Table 1. A table comparing the energy of the exact skyrmion solution and BPS anti-skyrmion with and without the boundary term, eq. (2.9). The BPS case is B = 2A = κ 2 and the anti-skyrmion is constructed from eq. (2.21) with f (z) = z. Subtracting the boundary term is necessary for the skyrmions to have negative energy.
In terms of the complex field w of the stereographic projection, the energy functional and topological charge can be expressed as and respectively. The single skyrmion exact solution can be rewritten in terms of the stereographic complex field w in eq. (2.7) as where the centre of the skyrmion is denoted as z = b, here n = (0, 0, −1) T . In ref. [20] it was shown that for the BPS case (B = 2A = κ 2 ) we can complete the square in the energy functional (2.10) and find first order BPS equations In terms of the stereographic field w in eq. (2.7), the BPS equations become and there are infinitely many BPS solutions, with topological charge Q ≥ −1, of the form for f (z) an arbitrary holomorphic function. Solutions of these BPS equations with the boundary condition n 3 → 1 at r → ∞ have energy given by their topological charge x. This is in contrast to the usual story for soliton solutions in the familiar O(3) sigma model [7] where the energy is given by 4π|Q[ n]|. The exact Hedgehog solution in eq. (2.18) at 2A = κ 2 coincides with the BPS solution with Q = −1. Full details of these BPS solutions are given in refs. [20,30].

Instability of ferromagnetic phase
We are now in a position to consider the ground state of the chiral magnet described by the energy functional in (2.10) along the solvable line B = 2A > 0. The solvable line coincides precisely with the boundary between the ferromagnetic (n 3 = +1) and canted ferromagnetic phases. Along the solvable line, we obtain the ferromagnetic ground state n 3 = +1 in the limit A → ∞. If we decrease the potential strength parameter A, the DM interaction becomes more important and can result in a negative energy. The exact single skyrmion energy eq. (2.14) shows that a single skyrmion is a positive energy excited state above the ferromagnetic ground state for A > κ 2 . On the other hand, the negative energy of the exact single skyrmion solution for A < κ 2 reveals that the homogeneous ferromagnetic state is no longer the lowest energy state. As A decreases down to κ 2 , we expect that skyrmions will start to appear as defects in the homogeneous ferromagnetic state. Therefore we expect that there is a phase boundary at A = κ 2 , below which an inhomogeneous state with many skyrmions is favored over the homogeneous ferromagnetic state. Repeating this for emphasis, the analytic skyrmion solutions become negative energy states for A < κ 2 and thus the ferromagnetic state is no longer the lowest energy solution of the equations of motion. This suggests that there will be a phase transition to a new inhomogeneous ground state constructed from the negative energy skyrmions. This is our first observation of the precise value of the phase transition point in a chiral magnet. This transition point cannot be established precisely by numerical methods and shows the power of the analytical single skyrmion solutions.
Detailed numerical studies [11,13,14] showed that the inhomogeneous phases are the skyrmion lattice phase, at intermediate values of A below κ 2 , and the spiral phase, for the region much below κ 2 . In subsequent sections, we wish to shed more light on the nature of these inhomogeneous ground states and their phase boundaries by analytic means, rather than numerically.
As we have an exact expression for the energy (2.14) and the radius (2.15) of a skyrmion, we can compute the energy density inside a skyrmion. We assume that the skyrmion energy (2.14) is approximately contained inside the circle of the skyrmion radius r = κ/A with the area π (κ/A) 2 . This leads to the skyrmion energy density being We observe that the energy density is symmetric around A = κ 2 2 . For values of A higher than κ 2 , the energy density is positive. As A → ∞, the single skyrmion energy increases towards 4π and the skyrmion radius goes to zero. Conversely, the skyrmion radius diverges in the limit A → 0, indicating that the n 3 = 0 circle becomes very large and that the region where n 3 = −1 covers most of the interior. In the same limit the energy of a single skyrmion diverges. However, the average energy density (2.22) goes to zero.

JHEP02(2021)095 3 Interactions of magnetic skyrmions
Although it is energetically more favorable to produce skyrmions below the phase transition point, one cannot create infinitely many skyrmions because of their interaction. Therefore, it is important to study the interactions between skyrmions in order to understand the inhomogeneous ground state. We can use the method from refs. [6,36] and consider the case of two Q = −1 skyrmions. The centres of the skyrmions are defined to be the points where n = 0, 0, −1 T . In this section we explain the approximation and use it to fit the interaction energy as a function of the ratio of separation to the skyrmion diameter for various values of A κ 2 .

Superpositions of skyrmions
We study the interaction between two magnetic skyrmions of degree Q = −1 by constructing a degree Q = −2 skyrmion configuration, which approaches a solution of the field equations at infinite separations. In refs. [6,36], the magnetisation field n is used to superpose two skyrmion solutions, n A and n B . The degree of the superposition is the sum of the degrees of n A and n B , in this case Q = −2. Here, we work at the level of the complex field in eq. (2.7) where the superposition is The superposition in terms of the complex field has the advantage of automatically satisfying the constraint n 2 = 1, in contrast to the superposition in terms of n. The two degree Q = −1 fields, w A and w B , both solve the equations of motion. However, the superposition w does not as the equations of motion are non-linear. As stated above the superposition becomes a solution only in the limit of infinite separation. Hence the accuracy of the approximation improves as the centres of the skyrmions (the poles of the w's) get further apart. Some examples of these superposition configurations at different values of A are included in figure 2. These figures show that as A decreases the approximation becomes less reliable. In fact the key quantity to consider is the ratio of the separation to the skyrmion diameter, this will be made explicit below. For the case B > 2A where the skyrmion fields are exponentially decaying, the decay exponent gives a characteristic length scale which controls the validity of the dipole approximation. Along the solvable line the power law decay complicates matters and it is not as obvious how to extract a characteristic length scale. However, it is still known that the superposition approaches a true solution for infinite separation. This means that the accuracy of the approximation will improve with increasing separation. In fact we show below that it is best to work in terms of a dimensionless separation which demonstrates that the accuracy of the approximation also depends on the value of A.
We take a skyrmion w A centred at the origin and another w B centred at The magnetisation vector field for the superposition of two Q = −1 hedgehog skyrmions with one skyrmion centred at zero and the other skyrmion centred at (30κ −1 , 0). On the left the plot is for A = κ 2 2 , in the middle is A = 2 10 κ 2 , and on the right is A = 1 10 κ 2 . This shows that as A decreases, at fixed separation, the interpretation of the superposition as two distinct skyrmions disappears and the approximation becomes worse.
The superposition is then given by This has degree Q = −2 for all values of b, except b = 0. At b = 0 there is a cancellation between the numerator and the denominator and we end up with w = 2iκe −iα Az , a Q = −1 configuration that is not a solution of the equations of motion.
We define the interaction energy of the pair of skyrmions w A , w B to be and regard it as the potential energy of the two-body force between the skyrmion pair. Our energy functional E is dimensionless (a dimensionful constant is divided out from the physical energy). With this convention, the interaction energy E int [w] is also dimensionless. For dimensional reasons, E int [w] is a function of only two dimensionless ratios of the three dimensionful parameters, κ, A, b. We define dimensionless coordinates ζ, β as The physical (dimensionful) separation |b| between two skyrmions is given in terms of the dimensionless separation |β| as |b| = |β|r/2. In appendix B, we show that the interaction energy has the following simple structure as a function of two dimensionless parameters β and κ 2 /(2A)

JHEP02(2021)095
We observe that the interaction energy is independent of the angle parameter α of the DM interaction, and depends on κ 2 /(2A) only linearly. The coefficient I sup 2 depends only on the magnitude of the dimensionless separation, |β|.

Qualitative analysis of interactions
At the BPS point B = 2A = κ 2 , the interaction energy of the superposition w is positive. This is demonstrated in the following way. The completing the square argument from ref. [20] can be applied to eq. (2.16) with A = κ 2 /2. Using the superposition w with the degree Q = −2, we obtain The second term is zero when the BPS equations are satisfied and is positive for every other configuration. This means that the energy in the degree Q = −2 sector is bounded below by −8π. Now from the definition of the interaction energy in (3.4) we have where we have used that both w A and w B are BPS skyrmions with energy E = −4π. This means that the interaction energy, E int must be positive to respect the completing the square argument. We also know that the interaction energy must vanish in the limit of infinite separation. Combining these facts, we find that the interaction is repulsive at least for asymptotically large separations. It was not known if the interaction energy is positive for other values of A κ 2 . For the range of A which we consider, we find that the interaction energy between the exact skyrmion solutions is always positive.
For general external magnetic fields and anisotropies there are no known exact solutions so a different approach is needed to study the interaction energy. One way to approach this is to consider the asymptotic form of the equation for the hedgehog profile, Θ(r), as is done for A = 0 in [35,36]. Another approach used in [34], again for A = 0, is to proceed numerically. Both of these approaches find that the interaction energy is repulsive and decreases exponentially for increasing separation.

Computation of the interaction energy
We are now ready to compute the interaction energy between two exact single skyrmion solutions w A and w B in eq. (3.2). Our approach proceeds by evaluating eq. (3.6) for the superposition w = w A + w B . We can assume without loss of generality that β is a real parameter. This is because the invariance of the energy under rotations, w(ξ) → e iθ w e −iθ ξ , means that the interaction energy just depends on the modulus of the dimensionless separation, |β|. We compute the integrand of eq. (3.7) and numerically integrate it using Mathematica, to find an expression for E int [w] using eq. (3.6). The integration is carried out in a 2 × 10 6 by 2 × 10 6 region so that we can approximate the boundary as being at infinity compared to the skyrmion separation. After repeating the integration at many separations, between |β| = 30 and |β| = 500 we fit the interaction energy as a function of the dimensionless separation |β| using a log-log plot and find: The interaction energy is then (3.11) The power law fit for I sup 2 (|β|) + 16π is plotted in figure 3. There is good agreement between the fit and the numerical computation for larger separations. At small separations there is a slight deviation with eq. (3.11) marginally over estimating the strength of the interaction.
The general trends that we observe are that skyrmions are small, weakly interacting, and have positive energy for large values of A κ 2 . As A decreases towards κ 2 , the interskyrmion interaction gets stronger and the skyrmions increase in size and decrease in energy, with the energy of an individual skyrmion vanishing at A = κ 2 , the phase transition point between the ferromagnetic and skyrmion lattice phases.
The interaction gets stronger and the skyrmion energy becomes more negative as A decreases. On the other hand, the relative separation |β| becomes smaller as A κ 2 decreases. Therefore our results get less reliable as A decreases, and it is likely that we are overestimating the repulsive interaction energy. This can be observed directly from figure 2 where for smaller values of A the distinction between the two single skyrmions becomes less pronounced. At short distances, the skyrmions will have a significant effect on each others shape and a more sophisticated method is needed to account for this.

Average energy density of skyrmion lattice state
For A < κ 2 , the homogeneous ferromagnetic state becomes unstable due to the emergence of negative energy single skyrmion solutions. The negative energy of individual skyrmions results in a preference to produce as many skyrmions as possible. However, the repulsive interaction prevents the skyrmions from becoming too densely packed. This competition is expected to lead to a lattice of skyrmions. Depending on symmetries, either a triangular or square lattice configuration will have the lowest energy. An extensive numerical study in ref. [11] indicated that a triangular lattice configuration is the ground state for the parameter range where the exact skyrmion solutions occur, as in the case of Abrikosov lattices of vortices in superconductors and superfluids. A cubic lattice of merons, half skyrmions, was also found numerically in another parameter region with larger anisotropy, B < 2A. In light of this numerical work we construct a simple model for a triangular lattice of skyrmions. Using this model, we compute the average energy per unit cell and find the lattice spacing which minimises this.
Assuming that the two-body force between skyrmions is the dominant interaction, we estimate the average energy density of the triangular lattice state. Let us identify the lattice spacing to be the separation |b| in the calculation of the two-body force between skyrmions in eq. (3.2). We find it convenient to use a dimensionless separation defined as the ratio of |b| to the diameter of a single skyrmion 2r = 2κ/A Using the interaction energies computed in the previous section, we now compute the approximate energy density for a triangular skyrmion lattice. Our approach is to compute the energy per unit cell. The unit cell contains one skyrmion, with 6 nearest neighbours, and has area √ 3 figure 4. The average energy density is then given by where E Sky is the skyrmion energy given in eq. (2.14), and C 1 . On the left is the critically coupled case A = κ 2 2 with its minimum at 11 and on the right is A = κ 2 whose minimum is for infinitely separated skyrmions.
Depending on the value of A κ 2 , eq. (4.2) will have different behaviour as a function of . In figure 5 we plot e Lattice as a function of for the BPS case, A κ 2 = 1 2 , and at the phase transition to the ferromagnetic phase, A κ 2 = 1. In the BPS case e Lattice is negative at its minimum and there the triangular skyrmion lattice is the ground state. While for A κ 2 ≥ 1 the minimum of e Lattice is positive and reached for asymptotically separated skyrmions. This is because the ferromagnetic state is the ground state and skyrmions arise as positive energy excitations above it.
When A κ 2 < 1 minimising eq. (4.2) with respect to leads to = min with where δ = 1.735. The average energy density e Lattice ( min ) and the dimensionless lattice spacing at the minimum, min , are plotted against A κ 2 in figure 6 and figure 7, respectively. These figures show that as A κ 2 decreases the dimensionless lattice separation at the minimum decreases. While at first the lattice energy density decreases with decreasing A κ 2 it reaches a minimum near A κ 2 = 0.3 before increasing towards zero. Let us make a few observations. Firstly we observe that a finite lattice spacing is achieved by the balance between the negative energy of constituent skyrmions and the positive energy due to the repulsive interaction between skyrmions. Secondly, this semianalytic understanding is most reliable quantitatively for regions near the phase transition point A κ 2 = 1. Thirdly for small A κ 2 , |β| becomes small. Hence our calculation of the interaction energy becomes less reliable. Since we are using a superposition of single skyrmion solutions as a variational Ansatz for a Q = −2 configuration, we naturally expect that the true total energy at a fixed separation, for a Q = −2 configuration, may be lower. Namely it is quite likely that our approximate result for the interaction energy is overestimating the repulsive interaction. Therefore we expect that the skyrmion separation in the lattice will be smaller than our estimate for small values of A κ 2 .

Transition to the spiral state
For smaller values of the parameter A, compared to the DM interaction parameter κ 2 , we expect from [11] that spiral states have lower energy than the lattice. Thus we expect that a phase transition between the skyrmion lattice and spiral states occurs. As we have seen in the preceding section, our estimate of interaction energy is likely to be an overestimate for smaller values of A κ 2 . On the other hand, the exact single skyrmion solution has energy E Sky = 4π(1 − κ 2 A ) → −∞ and radius r = κ A → ∞ as A → 0, with the average energy density e Sk inside the radius tending to a finite value.
From figure 7 we see that for small A κ 2 the lattice spacing and the skyrmion diameter become the same order of magnitude. This suggests that near the transition between the lattice and spiral phases we need a better model of the lattice. As a candidate model of this kind, we consider the case of skyrmions separated by a multiple of the skyrmion diameter, with c a positive constant. The case of c = 1 would have skyrmions separated by 2 κ A , the skyrmion diameter, which means that the magnetisation vector n will not attain the value JHEP02(2021)095 On the left for larger A the lattice spacing, |b| is much larger than the diameter of the skyrmion, 2r. On the right for lower A κ 2 both the skyrmion diameter and lattice separation increase but their ratio, = |b| 2r , decreases. The scale in the figure on the right has been exaggerated for effect. e 3 between the skyrmions. A situation like this is called a hexagonal meron lattice in [39]. As such we are interested in c > 1 and will in fact focus on c = 2. We are also assuming that the repulsive skyrmion interaction can be neglected, except to keep the skyrmions from merging, such that they form a lattice. This leads to the situation depicted on the right in figure 8. This type of triangular lattice made from nearly touching skyrmions was considered in [17]. This scenario seems to be more realistic for small values of A κ 2 , since we are overestimating the interaction energy in that region.
Using this alternative model of a skyrmion lattice we can estimate the value of A that the transition between the lattice phase to the spiral phase occurs at. To do this we note that for V (n The most general spiral solution at A = B = 0 is obtained by rotating and/or translating this solution in the x-y plane with the simultaneous rotation of the magnetization vector n. This has an identical energy density so we choose to work with the simplest case, eq. (4.5). When the potential is non-zero the spiral configuration given by eq. (4.5) has average energy which is shown in appendix C. While it is not a solution of the equations of motion it provides an approximation for a spiral state when the potential is small. In [10] for A = 0, B = 0 a spiral state like eq. (4.5) was motivated by comparing the energy functional eq. (2.1) to an energy functional whose ground states are Beltrami fields, solutions to ∇ × n + κ n = 0. A rough estimate of where the transition between the skyrmion lattice and the spiral takes place can be found by comparing this to the average energy density of a single skyrmion. The skyrmion average energy is given by eq. (C.5). The transition will JHEP02(2021)095 occur when these two energy densities are equal. When the separation is |b| = 4 κ A (taking c = 2 in eq. (4.4)) this gives the phase transition at this value of A matches the location of the phase transition found numerically in [11]. We refer the reader to appendix C for the explicit details of this computation. At this value of the coupling the spiral state has energy density This energy density is much lower than the energy density that we found for the skyrmion lattice in figure 6 using the dipole approximation. This is further evidence that our dipole approximation underestimates the energy density because it is over estimating the interaction energy.

Conclusion
In this paper, we have studied the interaction between two well-separated skyrmions and the phase structure of a chiral magnet. We have focused on the solvable line B = 2A, where A is the anisotorpy parameter and B is the magnetic field parameter in the potential, so that we can make use of the exact skyrmion configurations from ref. [20]. For A κ 2 1, where the DM interaction can be neglected compared to the potential, we have found the homogeneous ferromagnetic state n 3 = 1 as the ground state. This ferromagnetic phase continues for smaller values of A κ 2 until we reach the critical point at A κ 2 = 1, where the exact skyrmion solution has zero energy. Below this point A κ 2 < 1, the exact skyrmion solution has negative energy. The presence of a negative energy solution to the equations of motion means that the homogeneous ferromagnetic state is no longer the lowest energy state. This suggests that the true ground state becomes an inhomogeneous skyrmion lattice state as is found numerically in refs. [11,13,14]. As the skyrmion energy is known analytically the transition point, A = κ 2 , is found exactly without the ambiguity of numerical simulations. This is only possible on the solvable line where exact skyrmion solutions are known.
We have studied the properties of this skyrmion lattice state by working out the interaction between exact skyrmion configurations along the solvable line B = 2A. To achieve this we used the dipole approximation for the superposition of two, exact, single skyrmion solutions centred with a finite separation. and have found that the interaction energy is independent of the angle parameter of the DM interaction which distinguishes Bloch, Néel or other types of skyrmion. Exact expressions for the interaction energy density were obtained and after numerically integrating these a power law fit is found for the interaction energy of well separated skyrmions. Using the computed interaction energies, we have been able to predict that the phase transition from a skyrmion lattice state to a ferromagnetic state will happen for A κ 2 = 1. Our results have also demonstrated that in the skyrmion lattice phase as A κ 2 decreases the ratio of lattice separation to the skyrmion diameter decreases. This suggests that near the boundary between the lattice and spiral phases the lattice can be approximately computed by taking the lattice density to be the energy density of a single skyrmion.

JHEP02(2021)095
This allows us to gain a semi-analytic understanding of the phase diagram of a chiral magnet along the solvable line. In particular we determine the location of the phase transition between the homogeneous ferromagnetic phase and the skyrmion lattice phase precisely, and estimated the transition between the skyrmion lattice phase and the spiral phase.
We have emphasized the need for the boundary term in the energy functional to make the DM interaction term well defined, and have clarified its close relation to the asymptotic behavior of the magnetisation vector.
We need to emphasise that our work deals only with static configurations, since the exact skyrmion configurations from ref. [20] that we have used are minimisers of the static energy functional. Understanding the effect of time evolution of excited states above the ground state is an open problem, which requires an understanding of the dynamics of skyrmion configurations and their superposition. Another interesting open problem is to find the values of A and B at which the single skyrmion energy becomes negative away from the solvable line.
We have studied skyrmions in two dimensions, applicable to a thin film. In three spatial dimensions, skyrmions are lines for which one can study collective modes propagating along a skyrmion line [40]. Extension of the present work to three dimensions remains an open problem.

A DM interaction and boundary terms
In this appendix, we explain why the inclusion of the boundary term in the energy functional is mandatory. Let us consider the energy functional E[ n] bulk in the bulk with the DM interaction in eq. (2.1). To obtain the equation of motion by the variational principle, we consider a variation of the bulk energy functional δE[ n] bulk through the variation of the magnetization vector δn a

JHEP02(2021)095
where we have performed a partial integration resulting in the total derivative term (the last term). We can obtain the equation of motion if and only if the total derivative term (surface term in the variation of the total energy) is absent. We can achieve this goal by adding a boundary term to the energy functional whose variation should satisfy where the booundry value of the magnetization vector is denoted as n α,bound a . We can take the following boundary energy functional as a solution of the above requirement where we denote the boundary value of the variable n α as n α,bound . When we evaluate the contribution from the first term, we note that n α a (∇ i n α,bound a ) → n α,bound a (∇ i n α,bound a ) = 0 because of (n α,bound a ) 2 = 1. Therefore we need to consider only the second term. Thus the boundary energy functional is finally given by For definiteness, let us consider the case of single skyrmion (Q = −1) solution in the parameter region 2A ≤ B. Since the minimum of the potential is at n α = (0, 0, 1) ≡ n ∞ , we should impose a boundary condition fixing the boundary value as n α,bound a = n a,∞ . It is convenient to rewrite the boundary energy functional in eq. (A.4) using an identity in ref. [20] where ω stands for the vorticity density in the n 1 , n 2 components. The vorticity Ω is defined as in order to make variational principle well-defined. Thus we find that the correct energy of the skyrmion solution for 2A ≤ B is given by By choosing the potential V (n 3 ) for the solvable line (B = 2A) instead of the generic U (n 3 ), we obtain our energy functional in eq. (2.10).
In order to study the contribution from the boundary energy functional in the case of boundary at infinity precisely, we need to regularize the boundary by taking a large circle with the radius R and take the R → ∞ limit. The volume aib n α a n α,bound b of the parallerogram formed by the three unit vectors n α a , n α,bound b and the vector e i tangent to the boundary vanishes as R → ∞. However, the circumference of the boundary increases JHEP02(2021)095 as R → ∞. As a total derivative, E boundary can be rewritten into an integral along a circle C with a large radius R centered at the location of the skyrmion (defined by n 3 = −1) where n α ϕ is the tangential component of the magnetization vector at the large circle C in the x 1 , x 2 plane. Therefore the boundary energy functional contribute a finite result if the magnetization vector decays only as a power of 1/R: n α ϕ ∼ 1/R. It has been noted [10] that the value of the boundary term ∇ · ( n × δ n) depends on the fall off conditions on n and δ n. This is discussed from the more general gauged sigma model point of view in ref. [31]. In the familiar case of a skyrmion energy including just a Zeeman term studied in ref. [10], the A = 0 case of eq. (2.10), the fields and their fluctuations fall off exponentially fast and the boundary term is zero. However, along the solvable line, B = 2A, the fields and fluctuations fall off more slowly and the boundary term is not zero and may not even be well defined [31].
To showcase the explicit decay properties of the hedgehog configurations consider the radial equation of motion We focus on the case of a magnetic field aligned with the positive z-axis, B > 0. In particular, we consider the case of B ≥ 2A > 0, where the boundary condition is given by n → n ∞ = (0, 0, 1) for large values of A, B. In the limit r → ∞ the boundary condition on the profile function is lim r→∞ Θ(r) = 0. Considering the limit of (A.10) and keeping only terms linear in Θ we have decaying exponentially fast as r → ∞.

JHEP02(2021)095
This means that for B = 2A, the magnetisation vector field decays slowly and the boundary term is needed to make the variational problem well defined. For B > 2A, however, the magnetisation vector field decays exponentially fast and the boundary term in the energy is zero. Therefore the boundary term is mandatory for our solvable line B = 2A (the phase boundary between the canted and polarised ferromagnetic phases), from the viewpoint of the variational principle. For 0 < B < 2A (the canted ferromagnetic phase), we have to analyze the asymptotic behaviour of the magnetization vector using different differential equations from eqs. (A.10) and (A.11), since the boundary condition is now different, n ∞ = (0, 0, 1).

B Expression for the energy density
In this appendix, we calculate an expression for the skymion energy in terms of dimensionless parameters. The energy density with the boundary term removed is given in eq. (2.16). This can be written in terms of two dimensionless parameters in the following way. We change coordinates to the dimensionless ones ζ = 2A κ ze −iα , and β = 2A κ ze −ib in eq. (3.5 One should note that the functional dependence on the dimensionless parameter κ 2 /(2A) is found to be linear. This structure is valid for any field configurations including our superposition of two skyrmions. The two integrals I 1 , I 2 for the superposition are only functions of the magnitude of the dimensionless separation, |β|.
Applying the formula (B.2) to the case of a single skyrmion, such as in eq. (2.18), these integrals can be evaluated explicitly to give I 1 (β) = 8π, which is consistent with the expression for the energy in eq. (2.14).

JHEP02(2021)095
Applying the formula (B.2) to the superposition w = w A + w B and using the definition of the interaction energy (3.4), we can write the interaction energy as where I sup 1 (β) and I sup 2 (β) are the dimensionless integrals for the superposition field w = w A + w B .
For purely anti-holomorphic functions we can take this a step further and explicitly compute I sup 1 (β). The key observation for this is that in eq. (2.17) the degree Q[w] is written as When w is anti-holomorphic, such as for a superposition of skyrmions, it becomes The superposition w = w A + w B of two single skyrmions is an anti-holomorphic function with Q[w] = −2, as long as β = 0, and thus The period is defined as L = 2π κ so that n 3 = 0. If A is not zero then for eq. (4.5) the boundary conditions imply that The energy density of a single skyrmion in a disc of radius r = c κ A is Setting this equal to the spiral energy density we find which leads to the following quadratic equation for A, solving this for A results in A = κ 2 16 8 + 3c 2 ± 64 + 16c 2 + 9c 4 . (C.8) As c is positive and we know that A ≤ κ 2 , as above A = κ 2 we are in the ferromagnetic phase, we take the negative sign and get A = κ 2 16 8 + 3c 2 − 64 + 16c 2 + 9c 4 . (C.9) An interesting choice is to take c = 2, that is our lattice consists of a single skyrmion in a unit cell which is the disc with radius r = 2 κ A , for which we find A = κ 2 16 20 − 4 √ 17 0.22κ 2 (C. 10) as the value of A at which we expect a transition between the ferromagnetic phase and the skyrmion lattice. This value of A matches the value that can be read off the phase diagram that was numerically found in [11].
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.