Baby Skyrmion crystals

This paper describes a model for baby Skyrme crystal chunks with arbitrary potential by considering energy contributions from the bulk and surface of a crystal chunk. We focus on two potentials which yield distinct Skyrme lattices: the standard potential $V=m^2(1-\varphi^3)$ and the easy plane potential $V=\frac{1}{2}m^2 (\varphi^1)^2$. In both models, the static energy functional is minimized over all $2$-dimensional period lattices, yielding the minimal energy crystal structure(s). For the standard potential, the Skyrmions form a hexagonal crystal structure, whereas, for the easy plane potential, the minimal energy crystal structure is a square lattice of half-charge lumps. We find that square crystal chunks are the global minima in the easy plane model for charges $B>6$ with $2B$ a perfect square ($m^2=1$). In contrast, we observe that hexagonal crystal chunks in the standard model become the global minima for surprisingly large charges, $B>954$ ($m^2=0.1$).


I. INTRODUCTION
The Skyrme model [1] is a nonlinear field theory of pions which possesses topological solitons that describe baryons. It has been derived as a low-energy effective field theory of quantum chromodynamics (QCD) in the large colour limit [2,3] and, more recently, from holographic QCD models such as the Sakai-Sugimoto model [4]. One of the outstanding problems in the Skyrme model is the correct prediction of nuclear binding energies. One would like to be able to predict correct binding energies using the Bethe-Weizsäcker semi empirical mass formula, which is composed of five terms: the volume term, the surface term, the Coulomb term, the asymmetry term, and the pairing term. The coefficients in each term are normally determined empirically, and the problem at hand is: can the coefficients be estimated by using Skyrmions? In the Skyrme model, the classical mass of a Skyrmion roughly plays the same role as the volume and surface terms [5]. To be able to address these first two terms, we need to understand the phases of nuclear matter in the Skyrme model. An important question arises when studying phases of nuclear matter regarding the nature of high-density and low-density phases, and the transition between these phases. At high densities the Skyrmions form a crystal, whereas at low densities the Skyrmions are localised to their corresponding lattice points. As the ground state of nuclear matter has a crystalline structure in the classical approximation [6], understanding the infinite crystalline structure is key.
The baby Skyrme model [7] is a (2 + 1)-dimensional analogue of the (3+1)-dimensional Skyrme model, where interest in the baby Skyrme model has peaked again with the apparent prevalence of baby Skyrmions in condensed matter systems [8], quantum hall systems [9,10], chiral magnetic systems [11] and nematic liquid crystals [12]. In this paper we investigate the crystalline structure for * mmpnl@leeds.ac.uk baby Skyrmions and formulate our method in terms of an arbitrary potential. The choice of potential is crucial for baby Skyrmions as it determines the behaviour of the solitons and, thus, the underlying Skyrmion crystalline structure. For the first time, we propose a method to determine the surface energy contribution of a crystal chunk, once the minimal energy infinite crystal structure is determined. In order to predict the minimal energy of a charge-B crystal chunk with a fixed area, we then study isoperimetric problems for particular crystal symmetries.
For the standard baby Skyrme model [7] we find that the solitons form a hexagonal crystal structure with D 6 symmetry, which was first proposed by Hen and Karliner [13,14]. This hexagonal crystal structure is not unique to the baby Skyrme model; it also arises in quantum hall systems [10], chiral magnetic Skyrmion systems [15,16], Ginzburg-Landau vortices [17] and 3D Skyrmions in analogy with fullerene shells in carbon chemistry [18]. However, in the easy plane model [19][20][21] the optimal crystal structure is found to be a square lattice of half solitons, similar to that of the conjectured cubic crystal of half Skyrmions in the Skyrme model. This square crystal structure also arises in chiral magnets with easy plane anisotropy [22].
This paper is laid out as follows. We begin by discussing the general baby Skyrme model. From here, we introduce our numerical minimisation procedure and define the initial configurations that we use to initialise our algorithm. Then static multisoliton solutions are considered on the plane and a discussion of the possible global minima is presented. In Sec IV we investigate the lattice structure of baby Skyrmions and formulate a method to determine the optimal soliton crystal. Once these minimal energy infinite crystals are known, we construct a crystal slab model to numerically determine the surface energy of a crystal chunk. Finally, we study chunks of the infinite crystal in a bid to predict the classical energies of baby Skyrmion crystals.

II. BABY SKYRME MODEL
The general static baby Skyrme model consists of a single scalar field ϕ : Σ → S 2 where (Σ, g) is a 2-dimensional Riemannian manifold, and (S 2 , h, ω) is the 2-sphere embedded in R 3 with the induced flat Euclidean metric h and area 2-form ω. We will often write the baby Skyrme field as the 3-vector ϕ = (ϕ 1 , ϕ 2 , ϕ 3 ). The static energy functional of this model on Σ is given by where V : S 2 → R is the potential of the baby Skyrme model, | · | denotes the Hilbert-Schmidt norm and vol g is the volume form on Σ associated with its metric g. The parameter κ is a standard coupling constant for which we will set κ = 1 for our numerical analysis. Note that the differential form dϕ ∈ Ω 1 (Σ) is a linear map dϕ x : T x Σ → T ϕ(x) S 2 and so the Hilbert-Schmidt norm |dϕ x | depends on both the domain metric g and target metric h. However, the pullback of the area form ω ∈ Ω 2 (S 2 ) is ϕ * ω ∈ Ω 2 (Σ) and so its norm only depends on the metric g.
We will follow the terminology of harmonic map theory and refer to the first term in (1) as the Dirichlet energy. The Dirichlet term is also commonly referred to as the σ-model term. The second term is known as the Skyrme energy. It is conventional to label the three terms in (1) as E 2 , E 4 , and E 0 respectively, where each term is thought of as a polynomial in spatial derivatives with the subscript denoting the degree.
Let us introduce oriented local coordinates (x 1 , x 2 ) on the domain Σ and let {∂ 1 , ∂ 2 } be a local orthonormal basis for the tangent space T x Σ at x ∈ Σ. Then the Dirichlet energy in local coordinates is given by [23] where ∂ i := ∂ ∂x i . The Skyrme energy in local coordinates is It is easy to show that the pullback of the area 2-form ω on the 2-sphere to Σ is given by If the domain Σ is compact, then the baby Skyrme map ϕ : Σ → S 2 has an associated topological degree given by the pullback of the normalised area 2-form of the target space S 2 , In terms of the local coordinates (x 1 , x 2 ) on Σ, the topological degree is explicitly We refer to minimisers of the static energy functional E for fixed degree B as baby Skyrmions. The topological degree B is also referred to as the topological charge, or just charge, which we adopt throughout. Finding baby Skyrmions involves numerically solving partial differential equations. We do this using an accelerated gradient descent algorithm for second order dynamics, detailed in Sec II B.
In order for static (multi)soliton solutions to exist in the baby Skyrme system, we must evade Derrick's nonexistence Theorem. Consider a variation ϕ λ : Σ×R → S 2 of the baby Skyrme field ϕ such that ϕ λ=0 = ϕ. This has infinitesimal generator ∂ λ ϕ λ | λ=0 ∈ Γ(ϕ −1 T S 2 ), where ϕ −1 T S 2 is the vector bundle over Σ with fibre T ϕ(x) S 2 over x ∈ Σ. Explicitly, if we consider the spatial rescaling x → e λ x, then we have a one-parameter family of maps ϕ λ = ϕ(e λ x) such that ϕ λ=0 = ϕ. The rescaled static energy functional is then If the baby Skyrme field configuration ϕ is a minimiser of the energy E, then we require which yields the familiar virial constraint E 4 = E 0 . Unlike the (3 + 1)-dimensional Skyrme model, the potential E 0 = Σ V (ϕ)vol g is necessary in the baby Skyrme model otherwise the energy E[ϕ] can be lowered by spatial rescaling and thus cannot have minima. So the baby Skyrmions have a preferred size determined by the ratio κ/m. There also exists a lower topological Bogomol'nyi bound on the (static) energy given by [24] In comparison to the σ-model, the addition of the Skyrme term stabilizes the σ-model to spatial rescalings. The addition of any term that is cubic, or more, in spatial derivatives would stabilise the model (for example, the O(3) σ-model coupled to a massive vector meson [25]), however the Skyrme term is the lowest order expression that retains the second order nature of the equations of motion in terms of time derivatives. Throughout this paper, there are three choices of the physical space Σ that we will consider. The first physical space we will consider is the plane Σ = R 2 . For the solitons to have finite energy, it is necessary to impose the boundary conditions and select ϕ ∞ from the vacuum manifold of the model, i.e. such that V [ϕ ∞ ] = 0. Without loss of generality, we choose the vacuum ϕ ∞ = (0, 0, 1) throughout. This gives us the one-point compactification of space R 2 ∪{∞} ∼ = S 2 . The baby Skyrme field can then be viewed as the map ϕ : S 2 → S 2 , which has a conserved topological charge B ∈ π 2 (S 2 ) = Z, characterized as the winding number of the map and given explicitly by (6).
The second physical space we will consider is that of the 2-torus Σ = R 2 /Λ, in which our field satisfies the doubly-periodic boundary conditions ϕ(x) = ϕ(x + n 1 X 1 + n 2 X 2 ). Here n 1 , n 2 ∈ Z and X 1 , X 2 ∈ R 2 are a fundamental pair of periods that generate the lattice Λ. The maps ϕ : R 2 /Λ → S 2 have an associated integer degree, and so admit topological solitons.
Finally, the third physical space we consider is the infinite cylinder Σ = S 1 × R. This corresponds to a Dirichlet boundary condition in the x 2 -direction, lim |x 2 |→∞ = ϕ ∞ , and a periodic boundary condition in the x 1 -direction, ϕ(x) = ϕ(x + n 1 X 1 ), where n 1 ∈ Z and X 1 ∈ R 2 is a vector in the x 1 -direction. The maps ϕ : S 1 × R → S 2 also have a conserved integer topological degree and admit topological solitons.

A. Initial configurations
To initialise the numerical minimisation procedure, the gradient descent algorithm requires an initial configuration, or approximation to the static soliton. Consider the axially symmetric field configuration ϕ = (sin f (r) cos Bθ, sin f (r) sin Bθ, cos f (r)) , with the monotonically decreasing radial profile function f (r) satisfying f (0) = π and f (∞) = 0. Equivalently, the profile function f vanishes at the boundary of the grid.
Here, r and θ are polar coordinates in the plane and there exists an internal phase that has been set to zero by applying the global symmetry that rotates the ϕ 1 , ϕ 2 field components [26]. For our numerical minimisation procedure, the profile function is taken to be [27] f (r) = π exp(−r).
The initial field configuration is a linear superposition of static solutions; typically we use a set-up of N charge-1 skyrmions with a favourable relative phase-shift between each other (for maximal attraction). This is known as the attractive channel, and is dependent upon the choice of potential. The superposition is justified because the profile function decays exponentially. The superposition is done in the complex field formalism, i.e. where W : S 2 → CP 1 is the stereographic projection of the ϕ field of S 2 . We use the profile function of a static solution (typically of topological charge one) to obtain Using the radial ansatz (11), this is We can then assume that if the solitons are well separated in relation to their size, then we can approximate the resulting solution as W = N i W i , where N is the total number of solitons in the system, and B = N i B i is the total baryon number of the system. In terms of the stereographic coordinate W , the baby Skyrme field is B. Numerical minimisation procedure In order to find local minima of the static energy, we must numerically relax the baby Skyrme field. The numerical methods are carried out on a N 1 × N 2 grid with lattice spacings ∆x 1 , ∆x 2 . The baby Skyrme energy is then discretised using a 4th order central finite-difference scheme. This yields a discrete approximation E dis [ϕ] to the static energy functional E[ϕ], which we can regard as a function E dis : C → R, where the discretised configuration space is the manifold C = (S 2 ) N1 N2 ⊂ R 3 N1 N2 [28,29].
To compute the minima of the discretised static energy, initially a gradient descent method was chosen. However, gradient descent can be particularly slow method when the Hessian is of poor condition. A more efficient way to is to simulate the time development using an accelerated gradient descent algorithm known as arrested Newton flow [30]. The essence of the algorithm is as follows: we solve Newton's equations of motion for a particle on the discretised configuration space C with potential energy E dis . Explicitly, we are solving the system of 2nd order ODEsφ with initial velocityφ(0) = 0. Setting ψ :=φ as the velocity with ψ(0) =φ(0) = 0 reduces the problem to a coupled system of 1st order ODEs. We implement a 4th order Runge-Kutta method to solve this coupled system. The main advantage in implementing the arrested Newton flow algorithm is that the field will naturally relax to a local minimum. After each time step t → t + δt, we check to see if the energy is increasing. If E dis (t + δt) > E dis (t), we take out all the kinetic energy in the system by setting ψ(t + δt) = 0 and restart the flow. The flow then terminates when every component of the energy gradient δE dis δϕ is zero to within a given tolerance (we have used 10 −4 ). Unless stated otherwise, the plots shown throughout were simulated on a grid with 0.05 lattice spacings and grid sizes 1000 × 1000.
It is essential that we enforce the constraint ϕ · ϕ = 1. This is normally done by including a Lagrange multiplier term into the Lagrangian, and the form for the Lagrange multiplier can be found taking the dot product of the field with the resulting Euler-Lagrange equations. However, to do this numerically we have to pull our target space back onto S 2 . This is done by normalizing the Skyrme field ϕ each loop, We also need to project out the component of the energy gradient, and velocity, in the direction of Skyrme field, that is δE and

III. BABY SKYRMIONS ON R 2
Consider the plane R 2 with the usual flat Euclidean metric g ij = δ ij . The static energy functional of the baby Skyrme model on R 2 takes the familiar form (20) For numerical analysis, it proves convenient to express the static energy functional using Einstein's summation notation, that is where i, j ∈ {1, 2} and a, b ∈ {1, 2, 3}. The energy functional has a continuous O(3) symmetry before the symmetry is broken by the choice of potential term V (ϕ). To carry out arrested Newton flow, or a numerical relaxation method using the gradient of the energy, we need to calculate the energy gradient explicitly. We will do this in index notation for numerical convenience. The variation of the energy density with respect to field ϕ a is δE (a) Standard charge-1 energy density.
FIG. 1: Plots of the energy density of (a) the axially symmetric charge-1 baby Skyrmion, for the standard potential V (ϕ) = m 2 (1 − ϕ 3 ), and (b) the charge-1 baby Skyrmion for the easy plane potential

A. Standard baby Skyrmions
Numerous potentials have been proposed [7,19,26,[31][32][33][34][35] and studied extensively in the literature. However, there are two choices of potential that we are particularly interested in: the standard potential and the easy plane potential. These two theories are quite distinct and, as we will describe below, we should expect different phenomena. In the standard baby Skyrme model [7], the standard potential is an analogue of the pion mass term in the Skyrme model, and takes the form If we consider excitations around our unique choice of vacuum ϕ ∞ = (0, 0, 1), then the fields ϕ 1 and ϕ 2 acquire a mass m. The standard potential (22) spontaneously breaks the O(3) symmetry to an O(2) symmetry that acts on the field components ϕ 1 , ϕ 2 . For this potential, the charge-1 baby Skyrmion is axially symmetric and exponentially localised, see Fig. 1a. Piette et al. [7] studied the asymptotic interactions of standard baby Skyrmions and found that two well-separated charge-1 solitons have an interaction energy that can be calculated using a dipole approximation, such that where χ 1 − χ 2 is the relative phase. The attraction between these two well-separated baby Skyrmions is greatest when χ 1 − χ 2 = π. This is known as the attractive channel.
There is a rather nice way to graphically represent the phase of a baby Skyrmion [36] using a HSV color model, which is almost analogous to the Runge color sphere coloring in the Skyrme model. We begin by plotting the energy density and color it using the stereographic coordinate W , given in (13). The phase of W , arg(φ 1 + iφ 2 ), gives the hue of the color and is defined such that arg(φ 1 + iφ 2 ) = 0 is red, arg(φ 1 + iφ 2 ) = 2π/3 is green and arg(φ 1 + iφ 2 ) = 4π/3 is blue. We use the value of φ 3 to determine the lightness, such that φ 3 = +1 is white and φ 3 = −1 is black [37]. The coloring scheme detailed above is shown in Fig. 2.
For both potentials, multi-charged baby Skyrmion solutions have an underlying modular structure. One such structure of interest for the standard potential (22) is that of chains of solitons. This was first investigated by Harland [38], in the context of baby Skyrmions, and then later by Foster [39] and also Shnir [40]. Each chain has its ends capped by charge-2 solutions and the chain links are built from either charge-1 baby Skyrmions, with a relative phase of π with each neighbour, or charge-2 solitons. Shnir [40] showed that a chain with charge-1 links has a lower energy than a chain with charge-2 links within each homotopy class. For low charge solutions, chains appear to be good candidates for the global minima. A typical chain configuration for the standard potential (22) is displayed in Fig. 3a.
Foster [39] also investigated baby Skyrmions on a cylinder R × S 1 , and calculated the minimum energy per charge of an infinite chain to be E chain = 1.4549. We have carried out the same calculation using the lattice variation method detailed in Sec IV. This is done by imposing a Dirichlet boundary condition in the x 2direction, lim |x 2 |→∞ = ϕ ∞ , and a periodic boundary condition in the The periodic cell length L is then varied to minimise the energy, and a minimum energy of E chain = 1.4548 was found for a periodic cell length of L = 8.53. Thus, our results provide excellent fidelity and are displayed in Fig. 3e.
Later, it was realised by Winyard [41] that the energy density peaks at the ends of the chains could be reduced by joining the two ends into a ring-like solution with an added energy correction for the curvature of the ring. Using a least squares fitting, they were able to obtain values for the energy contributions from the chain caps and the ring curvature. They showed that ring solutions are a better candidate for the global minima for B > B ring ∈ Z. For the mass m 2 = 0.1, this transition from chains to rings is numerically found to occur at B ring = 15. A typical ring configuration is displayed in Fig. 3c. Soliton crystals in the standard baby Skyrme model, with the standard potential (22), were studied by Hen and Karliner [13,14]. Through their work they observed that the minimal energy soliton crystal was almost hexagonal by use of simulated annealing. So one would expect chunks of the infinite hexagonal crystal to be global minima for some B > B crystal ∈ Z. This prompts the basis of this paper: at what charge do chunks of the infinite soliton crystal become the global minima?

B. Easy plane baby Skyrmions
The second potential of particular interest is the easy plane potential, proposed by Jäykkä and Speight [19]. As with the standard potential, the easy plane potential leaves an unbroken O(2) symmetry. However, the canonical choice of vacuum ϕ ∞ = (0, 0, 1) distinguishes a point on the O(2) orbit and breaks the symmetry further to a discrete D 2 symmetry. Unlike the standard model, the charge-1 baby Skyrmion is not axially symmetric but rather is composed of two charge-1/2 baby Skyrmions. This is shown in Fig. 1b. As we did before, let us consider elementary excitations around our canonical choice of vacuum ϕ ∞ = (0, 0, 1), then the field ϕ 1 acquires a mass m and the ϕ 2 field is massless. Adapting the dipole approximation proposed by Piette et al. [7], and assuming that ϕ 2 mediates the dominant interaction asymptotically [19], gives us an interaction energy This shows that the interaction energy depends only on the average phase of the dipoles, which is exactly opposite of the situation in the standard model. While the phase coloring is particularly useful for the standard potential, it is actually more instructive to use the field structure of the field component ϕ 1 for the easy plane model. The red peaks and blues peaks are attracted to one another, but peaks of the same color are repelled by each other. As each individual peak in the energy density resembles a CP 1 model lump, then we refer to each lump as a half charge lump. Each of these half lumps are located at the red and blue peaks of ϕ 1 and come in pairs. So more information can be gained by studying plots of the ϕ 1 density than the energy density itself.
In contrast to the standard model, chains do not appear to be the global minima for low charges in the easy plane model. For charges B ≤ 6 with mass m 2 = 1, the global minima are 2B-gons or ring-like solutions. Chunks of an infinite crystal with a square/rectangular crystalline structure seem to be the global minima for almost all charges B > 6. An example of such a global minimum for B = 8 can be seen in Fig. 4c. The easy plane model also exhibits a modular structure with some more exotic local minima consisting of square and polygonal building blocks. One such solution is the B = 10 easy plane baby Skyrmion built from square and hexagonal units in Fig. 4e.
Although chunks of the assumed infinite crystal are prevalent, ring-like solutions and chain solutions do exist as other local minima. Jäykkä and Speight [19] showed that 2B-gon rings are the global minima for low charges, that is a single ring of 2B half lumps. For higher charges, it is energetically favourable for the ring solutions to form a double ring structure with some discrete symmetry. Example solutions for the easy plane model are shown in Fig. 4. The charge-5 chain in Fig. 4a is coincidentally a chunk of the assumed infinite soliton crystal and is only a local minimiser for B = 5, rather a 10-gon of half lumps is the global minimiser. The particularly interesting aspect of the charge-5 chain is that its shape is closer to that of a square soliton crystal chain than that of a double hexagon. This suggests that the square soliton crystal is a lower energy crystalline structure than the hexagonal soliton crystal.

IV. LATTICE STRUCTURE OF BABY SKYRMIONS
In a series of papers by Hen and Karliner [13,14], they determine the minimal energy soliton crystal for the standard model to be hexagonal. They scanned the parallelogram parameter space at a constant Skyrmion density to find the parallelogram that minimises the static energy. Once they found the optimal parallelogram, they then varied the Skyrmion density to find the minimal energy Skyrmion structure. In what follows, we refer to the shape of the lattice Λ as the lattice structure and the energy minimising Skyrmion as the soliton/Skyrmion crystal. Note that for an energy minimiser ϕ to be a soliton crystal it has to satisfy the extended virial constraints detailed below. We use a more robust method based on the work done by Speight [42] and propose an analytic method to determine the optimal lattice structure for an arbitrary potential. We then apply this method to study Skyrmion crystals in the standard model and in the easy plane model.
The physical space of interest is the 2-torus R 2 /Λ, where Λ is the set of all 2-dimensional period lattices α is a scaling parameter and {X 1 , X 2 } is a basis for R 2 .
We have written the fundamental pair of periods in the form Y i = αX i ∈ R 2 for later convenience, where we will introduce a constraint such that the area of the period lattice is α 2 . The crystallographic restriction theorem states that there are 5 Bravais lattice types in 2dimensions [43]. In each of these lattice types the fundamental unit cell is a certain type of a parallelogram.
To find the Skyrmion crystal we minimise the static energy functional over all period lattices. Equivalently, we fix our domain of ϕ to be the unit 2-torus R 2 /Z 2 and identify every other torus R 2 /Λ with R 2 /Z 2 , but with a non-standard Riemannian metric g. This metric g on R 2 /Z 2 is the pullback of the flat Euclidean metricḡ on R 2 /Λ via the diffeomorphism R 2 /Z 2 → R 2 /Λ. As we vary the period lattice Λ then the metric g varies [42]. Now, let F : T 2 → R 2 /Λ be a diffeomorphism with F ∈ GL + (2, R) and T 2 = R 2 /Z 2 , as shown in Fig. 5. Using the identification GL + (2, R) = SL(2, R) × R * /Z 2 , let A = [X 1 X 2 ] ∈ SL(2, R) and α ∈ R * such that F = αA ∈ GL + (2, R). We will now identify the domain of ϕ as Σ = T 2 , so that the Skyrme field is a map ϕ : T 2 → S 2 . The metric on T 2 is the pullback g = F * ḡ of the flat Euclidean metricḡ on R 2 /Λ. Explicitly, this is with inverse The Riemannian volume form is simply vol g = √ det g dx 1 ∧ dx 2 = α 2 dx 1 ∧ dx 2 . Then, using the local form for the Dirichlet term (2) and the inverse metric (28), we can compute the Dirichlet energy on T 2 to be given by Since the Dirichlet energy is conformally invariant, it does not have a dependence on the scaling parameter α. Likewise, using the local form for the Skyrme term (3), the inverse metric (28) and the pullback of the area 2form (4), the Skyrme energy on T 2 is and the potential energy is simply Putting this together, we see that the static energy functional for baby Skyrmions on the unit area 2-torus T 2 with the non-standard Riemannian metric g is As before, we need an explicit description of the energy gradient for our numerical analysis. The variation of the energy density with respect to field ϕ a can be obtained from the Euler-Lagrange field equations, that is where i, j ∈ {1, 2} and a, b ∈ {1, 2, 3}. To find the optimal lattice structure, we must vary the static energy functional (32) with respect to the period lattice parameters X 1 , X 2 and α. Firstly, taking the variation of the static energy functional (32) with respect to the scaling parameter α, yields the following relation for the scaling parameter: Thus, the area of the period lattice is determined by the ratio of the flat Skyrme term to the flat potential term. Determining the fundamental pair of periods X 1 , X 2 which minimise the Dirichlet energy E 2 is a constrained quadratic optimization problem with the nonlinear constraint det([X 1 X 2 ]) = 1. For notational convenience, let us write Then the Dirichlet energy (29) can be expressed in the quadratic form where x = X 1 X 2 is a 4-vector and Q is a 4 × 4symmetric matrix. This constrained quadratic optimization problem can be solved by including the Lagrange term γ(det([X 1 X 2 ]) − 1) into (35), where γ ∈ R * is a Lagrange multiplier. This reduces the problem to an eigenvalue problem By definition, for an energy minimiser ϕ : R 2 /Λ → S 2 to be a soliton lattice, its stress tensor S[ϕ] must be L 2 orthogonal to the space of parallel symmetric bilinear forms E (a 3-dimensional subspace of the space of sections of the rank 3 vector bundle T * R 2 /Λ T * R 2 /Λ). Furthermore, if the Hessian of the soliton lattice is positive definite then it is a soliton crystal. In fact, Speight [42] showed that every baby Skyrme lattice is a soliton crystal. The stress tensor of ϕ is given by [32] Let E 0 be the 2-dimensional space of traceless parallel symmetric bilinear forms. Then the Skyrme field ϕ is a soliton lattice if and only if ϕ is L 2 orthogonal toḡ and E 0 , where we recall thatḡ is the metric on R 2 /Λ and not T 2 . This gives us the familiar virial constraint Let (y 1 , y 2 ) be local orthonormal coordinates on R 2 /Λ and ε ∈ E 0 . Then for S to be L 2 orthogonal to E 0 , we require As E 0 is spanned by ε 1 = (dy 1 ) 2 − (dy 2 ) 2 and ε 2 = 2dy 1 dy 2 , we get the following additional virial constraints and R 2 /Λ ∂ϕ ∂y 1 · ∂ϕ ∂y 2 dy 1 dy 2 = 0.
These additional virial constraints state that the Skyrme map ϕ must be conformal on average. We have shown above that for the soliton lattice ϕ to be critical with respect to variations of the period lattice Λ, it must satisfy the extended virial constraints in each lattice cell. The extended Derrick virial constraints (38)-(40) can be imposed as a consistency check when implementing the lattice optimisation method detailed above. For numerics, the domain manifold of interest is the 2-torus T 2 = R 2 /Z 2 with metric g = F * ḡ , where we previously introduced the diffeomorphism F : T 2 → R 2 /Λ. Recall that we used the identification F = αA ∈ GL + (2, R) for α ∈ R * and A = [X 1 X 2 ] ∈ SL(2, R). Thus, using the compact notation (34), the generalised virial constraints (39) and (40) on T 2 are given, respectively, by and In this section, we have shown that the problem of determining the optimal lattice structure that minimises the baby Skyrme energy (32) amounts to solving an eigenvalue problem (36). During each iteration of our numerical minimisation algorithm, we perform an accelerated gradient descent then we compute the scaling parameter α via (33) and solve the eigenvalue problem (36) to give us the pair of periods X 1 , X 2 . We also check that the generalised virial constraints (38), (41) and (42) are satisfied each iteration, showing that the energy minimiser ϕ : R 2 /Λ → S 2 is indeed a soliton lattice and thus a soliton crystal. This determines the lattice Λ and the algorithm in turn determines the Skyrmion crystal. The numerics detailed throughout this section were carried out initially on a 200 × 200-grid with lattice spacing ∆x = 0.005, with a final higher accuracy simulation carried out on a 500 × 500-grid with lattice spacing ∆x = 0.002. Finer meshes were tried but there were no considerable changes in the final energy. Note that the coarser 200 × 200-grid would be sufficient as this gives approximately the same accuracy as the numerics for the baby Skyrmions on R 2 . It is also worth noting that the lattice spacings are fixed sizes on the discretised unit area 2-torus T 2 , whereas the equivalent lattice spacings on the discretised 2-torus R 2 /Λ vary as the lattice Λ varies.

A. Standard baby Skyrmion crystals
Employing the lattice optimisation method detailed in Sec IV for the standard potential (22) with m 2 = 0.1 (and κ 2 =1), the optimal lattice is found to be an equianharmonic lattice with the baby Skyrmions forming a hexagonal Skyrmion crystal with D 6 symmetry. This is found for almost all B = 2 initial configurations on random initial lattice geometry (with the exception of relaxing to the infinite chain solution occasionally). Each unit  cell contains a charge of B = 2 and has sides of equal length L crystal = 9.60 with the angle between the two periods X 1 , X 2 being θ = 2π 3 , giving a unit cell area of A = 79.84. We find that the Skyrmion crystal has energy E crystal = 1.4543, which is lower than the infinite chain energy E chain = 1.4548. Note that when we refer to energy values, we have normalised them by the Bogomolny bound, i.e. E := E/(4πB). The hexagonal Skyrmion crystal can be seen in Fig. 6. As a fidelity check with Hen and Karliner's work, we also determined the optimal lattice to be equianharmonic with a hexagonal Skyrmion crystal for m 2 = 0.1 and κ 2 = 0.03. The energy for this crystal is found to be E crystal = 1.0799, which is in excellent agreement with their numerically determined value of E crystal = 1.08.
Other soliton crystals were searched for at numerous topological charges for various initial configurations and initial lattice geometry. However, they all had a tendency to relax into a chain or rows of separate chains with the infinite chain energy E chain = 1.4548. A slightly lower energy configuration was found for rows of adjacent chains with all the charge-1 links rotated by π in one chain relative to the other. This attractive chains configuration has an energy of E 2-chains = 1.4545 and is shown in Fig. 7.

B. Easy plane baby Skyrmion crystals
As previously proposed in Sec III B, it seems likely that there may possibly be a few soliton crystals for the easy plane model. This prompts the search for Skyrmion crystals for a range of charges with various initial configurations. The lowest energy Skyrmion crystal is a square of half lumps with D 4 symmetry in a square lattice for

V. BABY SKYRMION CRYSTAL CHUNKS
The soliton crystal is the lowest energy solution, so one would expect chunks of the soliton crystal to be the global minima for charges past a critical charge B crystal . A starting point would be to split the crystal chunk energy into a bulk volume, or area, term and a surface term. For a given charge B, we know the minimal energy soliton crystal, the corresponding lattice Λ and the lattice area. So, the bulk area term is easy to calculate. However, the problem lies in minimising the surface energy contribution for a fixed area, which corresponds to minimising the crystal perimeter for a fixed area. This is known as an isoperimetric problem. Even once the minimal energy crystal chunk shape has been found, we still require an estimate of the surface energy (per unit length) to determine the surface energy of the crystal chunk.

A. Surface energy of a baby Skyrmion crystal chunk
The surface energy per unit length of a Skyrmion crystal chunk can be predicted by using a crystal slab model. Skyrmion crystals are layered on an infinite cylinder Σ = R × S 1 of width L = L crystal , and the number of layers n ∈ N are increased to estimate the surface energy contribution. As stated in Sec II, this corresponds to a Dirichlet boundary condition in the x 2 -direction, lim |x 2 |→∞ = ϕ ∞ , and a periodic boundary condition in the x 1 -direction, ϕ(x) = ϕ(x + n 1 X 1 ). Each layer contributes a charge of 2 in this model, giving an n-layer crystal slab a total charge of B = 2n. The crystal slab layering can be seen in Fig. 9.
For the standard potential, each hexagonal baby Skyrmion has 6 bonding sides (or nearest neighbours) and each crystal slab edge has 2 unbonded sides. Similarly, for easy plane crystal chunks, each baby Skyrmion has 4 bonding sides and each crystal slab edge has 2 unbonded sides. We can express the crystal slab energy as where N free is the number of free bonds (or unbonded sides) and E bond is the energy of each free bond. For the standard and easy plane crystal slabs in Fig. 9 we have N free = 4. We can approximate the surface energy contribution by applying a least squares fitting to the crystal slab energy normalised by the Bogomolny bound, where E bond is the normalised free bond energy such that E bond = 4πE bond . For approximating the surface energy, we computed the energies of various n-layer crystal slabs with n ∈ {3, . . . , 11}. Using a trust region reflective algorithm, and the crystal slab approximation (44), we find that for the standard potential E bond = 0.0031 (with m 2 = 0.1) and for the easy plane potential E bond = 0.0103 (with m 2 = 1).

B. Standard crystal chunks
To model a Skyrmion crystal chunk, we split the crystal chunk energy into a bulk volume term and a surface term, E chunk = E bulk + E surface . The surface energy contribution of a baby Skyrme crystal is determined by the number of unbonded sides. As stated before, for the standard potential, each hexagonal baby Skyrmion has 6 bonding sides (or nearest neighbours) which means there are many possible arrangements for crystal chunks for a given charge B. For easy plane crystal chunks, the square lattice is the minimal energy crystal configuration so we only consider each easy plane baby Skyrmion to have 4 bonding sides. Our aim is to determine the shape of an equilibrium crystal by minimising the total surface energy associated to the crystal-vacuum interface. In crystallography one normally employs the Wulff construction method to determine the equilibrium shape of a crystal chunk. However, we take a simpler approach and only consider the perimeter of the crystal chunk boundary, not its shape. Equivalently, we are considering the number of free bonds in a given crystal chunk. This enables us to write the energy in the form Therefore, we want to find the crystal chunk that minimises the number of free bonds, and hence its surface energy contribution, for a fixed charge B and crystal area A.
The infinite standard crystal has a discrete D 6 symmetry, so we propose that minimal energy chunks of the infinite crystal take the form of layered hexagonal solitons, as can be seen in Fig. 11. The number of charge-2 units in each layer is precisely 6n. As we consider each charge-2 baby Skyrme unit to have 6 bonding sides, we can determine the number of free bonds in an n-layer crystal chunk to be N free = (12n + 6). This accounts for the 2 free bonds on each outer charge-2 soliton plus the additional free bond at each vertex of the crystal chunk. The total charge of a crystal chunk is B = 2(3(n+1)n+1), such that Thus, we can approximate the normalised energy of a hexagonal standard baby Skyrmion crystal chunk to be given by To determine the transition charge B crystal where chunks of the infinite soliton crystal become the global minima, we need to compare the crystal chunk model (47) to chain and ring models. Using the models proposed by Winyard [41], and our numerically determined value for the infinite chain, we are able to approximate ring and chain solutions. The results are plotted in Fig. 10, which includes data points from crystal chunk solutions for numerous charges B up to B = 2054. It can be observed that the crystal chunk approximation (47) fits the data very well. We find that crystal chunk solutions become global minima for charges B > B crystal = 954. Energy density plots of the Skyrmion crystal chunk solutions are shown in Fig. 11. All of these crystal chunks were found numerically on grids with lattice spacing 0.05, with grid sizes ranging from 800×800 to 2400×2400. Crystal chunk solutions for B = 1262 and B = 2054 are not shown but were obtained on grids with lattice spacing 0.1 and grid sizes 3000 × 3000 and 3800 × 3800, respectively.

C. Easy plane crystal chunks
To predict the energy of an easy plane crystal chunk is somewhat more challenging. There exists three soliton crystals, all relatively close in energy with one other. For charges B such that √ 2B ∈ Z, the minimal energy crystal chunk is the minimal perimeter √ 2B × √ 2B-square of half lumps. For non-square 2B, it becomes exceedingly difficult to predict the global minima. This is because there exists a smorgasbord of local minima for a given charge, which increases with the charge number. Nevertheless, some progress can be made if we consider rectangular crystal chunks built from the square Skyrmion crystal. In the first instance, if 2B has factors other than 2 and B, say a ∈ Z and 2B/a ∈ Z, then the (local) minimal energy crystal chunk will be an a × 2B/a-rectangle of half lumps such that the sum a + 2B/a is minimal with respect to the other pairs of possible factors. We ignore the trivial factors 2B and 1 as the 2B × 1-chain is most likely not a local minimiser for the easy plane model. Random initial configurations do not relax to a single linear chain, unlike the standard model. Even starting with a single chain initial configuration in attractive channel orientations does not result in a relaxed final state of a single chain, it normally relaxes to the double chain. To find the pair of factors with minimum sum, one would find the factor a ∈ Z of 2B that minimises the perimeter function f (a) = 2(a + 2B/a).
Using the above information, we are able estimate the energy of a crystal chunk for a given charge B. Similar to the standard model, we split the crystal chunk energy into a bulk term and a surface term. For charges B, we need to determine the pair (a, 2B/a) of minimal sum factors of 2B. Then we can calculate the surface energy and determining the bulk energy is straightforward. Explicitly, the normalised energy for a charge-B crystal chunk is given by This approximation for square crystal chunks (2B = a 2 ) is plotted in Fig. 12, alongside the corresponding true numerically determined (normalised) energies. Clearly, the further a deviates from √ 2B the higher the surface energy contribution. So one would expect there to be normalised energy bands at high charges for this rectangular approximation. These bands can be determined in the limit B → ∞, and for such highly rectangular pairs (a, 2B/a) the bands in the limit B → ∞ are given by Since there are three soliton crystals, there are obviously better crystal chunk solutions for highly rectangular factors (a, b). As an example, lets consider local minima for the charge-13 easy plane soliton. There are three solutions that one might expect to be contenders for the crystal chunk for this charge. Firstly, the natural choice is the minimal perimeter rectangle, which will be a double chain, or simply a 2 × 13-rectangle of half lumps. This is depicted in Fig. 13a. The next idea is to consider the minimal perimeter rectangle of half lumps for a B − 1 = 12 crystal chunk, then add a half lump pair to one of its corners to create a defect. This results in a 6×4-rectangular crystal chunk with one distorted hexagonal corner, as shown in Fig. 13b. The third contender is akin to the corner cutting method used in the Skyrme model [44,45], in which we try to remove half lumps (one blue and one red) from two corners of the minimal perimeter 7 × 4-rectangular B = 14 crystal chunk. This does not have the intended effect of missing half lumps on corners, rather it pulls the rest of the row, and the adjacent row, away from the chunk to form an arc with two half lumps more than the half lump height of the rectangular chunk. This can be seen in Fig. 13c. Out of these three most likely crystal chunk candidates for a charge-13 soliton, the rectangular crystal chunk with one distorted hexagonal corner is the minimal energy solution. So for a charge-13 baby Skyrmion the minimal perimeter rectangle model fails as a candidate for the global minimal energy crystal chunk in the easy plane model. So even at charge-13 we have already found a lower energy crystal arrangement than the rectangular crystal chunk. One would expect that adding hexagonal/octagonal defects to nearly square crystal chunks would result in lower energy solutions than rectangular crystal chunks. However, for square crystal chunks, such that B ∈ √ 2B, the rectangular crystal chunk model (48) is an excellent approximation.
In brief, we conjecture that the prevalent minimal energy crystal chunks for the easy plane model are squares of half lumps if √ 2B ∈ Z, otherwise they are minimal perimeter rectangles of half lumps, with some crystal chunks having distorted hexagonal defects. We have also proposed an empirical model to determine the energy for a given square/rectangular crystal chunk.

VI. CONCLUDING REMARKS
In this paper, we have presented a method to determine soliton crystals on an optimised lattice for arbitrary potentials in the baby Skyrme model. Once the minimal energy soliton crystal is known, the solitons can be layered by use of a crystal slab model and the surface energy per unit length obtained numerically. Using insight obtained from the soliton crystal and the surface energy, chunks of the Skyrmion crystal can be constructed and their corresponding energies determined. For the standard potential, we have demonstrated that the global minimum energy Skyrmion crystal exhibits a clear hexagonal D 6 symmetry. This hexagonal soliton crystal has a lower normalised energy than the infinite chain solution proposed by Foster [39]. We propose that the global minima are layered hexagonal crystals for B > 954 with m 2 = 0.1.
We determined that rows of adjacent infinite chains in attractive orientations have normalised energies close to that of the the infinite hexagonal crystal. So it is quite possible that concentric ring solutions in attractive orientations could be the global minima for charges B with B ring < B < B crystal , this would need to be investigated. Likewise, Winyward [41] showed that chain solutions could also intersect to form junctions, and proposed that networks of standard baby Skyrmions could be the global minima between rings and crystal chunks. This too would need to be studied.
Solitons in the easy plane model take the form of configurations of half lumps. This model has three soliton crystals all relatively close in energy: square, hexagonal and octagonal. Of these, the square Skyrmion crystal of half lumps is the global minimum. This is more reminiscent of the three dimensional Skyrme system and, in a manner of respect, a better analogue. The easy plane model exhibits a plethora of local minima with various different types of symmetries. We conjecture that, when 2B is a perfect square, square crystal chunks are the global minima. For rectangular 2B, the minimal energy crystal chunks are, as close to square as possible, rectangular crystal chunks of half lumps with some chunks having hexagonal surface defects. The study of internal anomalies has not been carried out, so it is possible that FIG. 14: Baryon density plot of the cubic lattice of half Skyrmions in the 3D Skyrme model. The 3D crystal is colored using the Runge coloring scheme detailed in [45]. the inclusion of a defect into the bulk is more energetically favourable over a surface defect.
The obvious extension of the work detailed in this pa-per is to the three dimensional Skyrme model. The crystal structure has been studied extensively in the literature [6,18,[45][46][47][48][49][50]. All the lattice variations that have been studied were for a cubic lattice of side length L, in which only L is varied. In this case, the "optimal" soliton crystal is the cubic arrangement of half Skyrmions [6] as shown in Fig. 14. It is believed that the simple cubic arrangement of half Skyrmions is the lowest energy Skyrmion crystal. However, if the Skyrme model and the baby Skyrme model truly are analogous, it is possible that a hexagonal Skyrmion configuration could prevail (or even an entirely different crystalline structure).