Near-BPS baby Skyrmions

We consider the baby-Skyrme model in the regime close to the so-called restricted baby-Skyrme model, which is a BPS model with area-preserving diffeomorphism invariance. The perturbation takes the form of the standard kinetic Dirichlet term with a small coefficient $\epsilon$. Classical solutions of this model, to leading order in $\epsilon$, are called restricted harmonic maps. In the BPS limit ($\epsilon\to 0$) of the model with the potential being the standard pion-mass term, the solution with unit topological charge is a compacton. Using analytical and numerical arguments we obtain solutions to the problem for topological sectors greater than one. We develop a perturbative scheme in $\epsilon$ with which we can calculate the corrections to the BPS mass. The leading order ($\mathcal{O}(\epsilon^1)$) corrections show that the baby Skyrmion with topological charge two is energetically preferred. The binding energy requires us to go to the third order in $\epsilon$ to capture the relevant terms in perturbation theory, however, the binding energy contributes to the total energy at order $\epsilon^2$. We find that the baby Skyrmions - in the near-BPS regime - are compactons of topological charge two, that touch each other on their periphery at a single point and with orientations in the attractive channel.


Introduction
Nuclear binding energies are relatively small compared to the total nuclear mass, i.e. roughly of the order 1/100 and slightly less than that for the light nuclei. A challenge for solitonic nuclear models, such as the Skyrme model, is to reproduce this result since they usually overestimate the binding energies. This sparked the interest, which has grown in the recent years, in looking for BPS solitonic models and small perturbations around them. In BPS models the energy is proportional to the topological charge, in this case the baryon number, so there is no binding energy between the nuclei. A small perturbation of the BPS model would give, presumably, a small binding energy. These "near-BPS" models have small binding energies already at the classical level. One model of this kind is the so-called BPS Skyrme model which consists of a Lagrangian with a potential plus a sextic term for the pion fields, which is also the squared topological charge density [1][2][3]. The ordinary Dirichlet (two-derivative kinetic) term and the Skyrme (four-derivative) term are considered as small perturbations. The BPS Skyrme model has the nice phenomenological feature of having an infinite-dimensional moduli space for static solutions consisting of all volume-preserving diffeomorphism maps. This resonates well with the liquid drop model of an incompressible fluid so successful in describing nuclear matter. Various theoretical and phenomenological studies have been done regarding the near-BPS Skyrme model [4][5][6][7][8] but a fundamental question still remains unanswered: When the non-BPS perturbations are gradually switched off, the solution in an arbitrary baryonic sector flows to a particular BPS solution, one of the many in the infinite dimensional moduli space, but which one? Clearly, it is the infinite number of zero modes -the very same feature which makes the model attractive -that makes the question hard to answer, both analytically and numerically. Attempts to answer this question by direct numerical study have not given a clear definite answer [9]. It was suggested, in modeling the nuclear stars, that the answer could be the axially symmetric multi-Skyrmion [6], but this, for reasons that we will explain, cannot be the true minimum of the energy. It is like searching for a needle in a haystack; clearly we first need a good and well-motivated analytical guess and then we can embark on the numerical study. In this paper we will provide, with the use of a toy model, an analytical guess, at least for a certain class of potentials.
The baby-Skyrme model, as the name suggests, is a toy model for the Skyrme model in one spacetime dimension less, i.e. 2+1 instead of 3+1, and with one target-space dimension less, S 2 instead of S 3 [10][11][12][13]. In this sense it has been used in the past as a test bed for various ideas and conjectures about the more difficult model. Among the various similarities between the two models, both of them possess a BPS restricted model. 1 For the baby-Skyrme model, the restricted version contains only a four-derivative term (the square of the topological charge) and a potential. This BPS model has area-preserving diffeomorphism invariance for the static solutions and it was actually discovered much earlier than the related model in 3 + 1 dimensions [14]. The soliton structure depends very much on the type of potential that is chosen, and in particular on the potential behavior close to the minimum. In this paper we will concentrate on the type of potentials that are quadratic at the minimum and for which the soliton solution of the BPS model is a "compacton" [15][16][17]. A compacton owe its name to the fact that all nontrivial behavior is contained in a compact region of space, outside of which the field is exactly at the minimum of the potential. The topic of baby Skyrmions is interesting in its own right; particularly, much attention has paid given recently to the applications in magnetic materials, see ref. [18] for a review and refs. [19][20][21] for some recent theoretical work. So what we describe in this work may hopefully be applied also in these cases, with opportune modifications.
Deformations of BPS models by small perturbations have been discussed in various cases, see for example refs. [22,23]. The general feature is that the moduli space of solutions is lifted by some effective potential and the solutions, when the perturbation is small, lie close to the minima of said potential. In the case of BPS models with volume preserving diffeomorphism symmetry, the moduli space is infinite dimensional and thus the problem becomes considerably more difficult. A rigorous mathematical definition of the problem has been put forward in ref. [24]. In particular, if the perturbation is the kinetic Dirichlet term, solutions have been coined "restricted harmonic". Harmonic because they minimize the Dirichlet term, restricted because the minimization is constrained to a sub-space of the functional space, namely the infinite-dimensional moduli space of solutions of the BPS equations. In this paper, we will show that the concept of restricted harmonic is still not enough to solve the problem in the case of compactons. In fact, for compactons, the restricted harmonic maps are still an infinite-dimensional subspace of the moduli space. Finding the correct solution thus requires to go to even higher order in the perturbation. This in particular affects the way the binding energy is going to zero as the perturbation is switched off.
In this paper, we focus on the baby-Skyrme model with a standard pion mass term and the perturbation is only the Dirichlet (kinetic) term with coefficient . When is switched off, the restricted baby-Skyrme model is recovered. We address the near BPS behavior of the solutions in the model with topological charges Q = 1, 2 and 4 by large brute-force numerical simulations and by developing a framework for perturbation theory around the BPS solution. We find that the correction to the energy (mass) of the axially symmetric compactons, to leading order in , comes from the kinetic (Dirichlet) term and remains inside the compacton domain. This is specific to the case of the compactons. This leadingorder correction thus cannot teach us anything about the binding energy of the solutions in the near-BPS regime. A main result of this paper, in the part of the perturbation theory, is that the binding energy does not appear at any particular order, if we strictly truncate the correction to a specific polynomial order. Nevertheless, we can calculate the binding energy using a linearized perturbation and it shows that the contribution to the binding energy is in some sense all-order, which we could choose to interpret as a "nonperturbative effect". This happens because the perturbation that gives rise to the binding energy takes the form exp(−r/ √ ) and thus does not have any nonvanishing Taylor coefficients upon expansion. For the choice of potential we work with in this paper, we find that the near-BPS solutions take the form of charge-2 baby Skyrmions placed very close to each other so that their periphery just touch.
The paper is organized as follows. In sec. 2 we give a brief review of the baby-Skyrme model and of the solutions in the BPS limit. In sec. 3 we develop a systematic framework for calculating the perturbations to the BPS soliton, which describes the solutions in the regime near the restricted BPS model with the Dirichlet term being a small perturbation. In order to calculate the binding energy, we find that it is necessary to go to the third order or the next-to-next-to-leading order in perturbation theory. We test the perturbative scheme on axially symmetric compactons, for which we have a comparison with the exact numerical solutions. In sec. 4 we present the results of the full numerical PDE computations for topological charges Q = 1, 2, 4 and calculate the total energy to be used for comparison with the results in sec. 5, where we calculate the binding energies in perturbation theory, by solving a linear PDE. In sec. 6 we contemplate what the solutions with large topological charges look like. Finally we conclude with a discussion in sec. 7. This paper studies primarily the standard pion mass term, but the cases of other potential are relegated to appendix A.

The model
The model is based on the BPS baby-Skyrme model with a non-BPS deformation which is manifested as the kinetic term with coefficient where the kinetic term and the Skyrme terms, respectively, are given by and φ = (φ 1 , φ 2 , φ 3 ) is a real 3-vector on which the nonlinear sigma model constraint φ · φ = 1 is imposed by the Lagrange multiplier, λ, and the metric signature is here taken to be the mostly positive one. The symmetry of the Lagrangian for V = 0 isG = O(3). In the vacuum, this symmetry is spontaneously broken toH O(2), which gives rise to the If we choose the potential of the form with 1 2 ≤ p < 2, then the soliton solution, called a baby Skyrmion, has support on a compact domain in the limit of → 0 and is then dubbed a compacton [16]. In this limit the model is then called the BPS baby-Skyrme model or restricted baby-Skyrme model.
In this paper, we will use the pion mass term which turns the baby Skyrmion into a compacton in the limit → 0.
The topological degree of the baby Skyrmion or compacton is given by which also counts the number of baby Skyrmions in R 2 .
For the analytic calculations, it will prove convenient to use the following parametrization in terms of which the nonlinear sigma model constraint φ · φ = 1 is manifest and we can thus write the Lagrangian components as 9) and the potential as In this parametrization, the topological charge reads We will keep both formalisms in this paper, as ω is useful for analytical calculations and φ is useful for the numerical analysis.

BPS solutions
The model contains a BPS submodel, which is given by setting = 0 in the Lagrangian (2.1): In order to find the BPS equation and corresponding solutions for each topological sector, Q, it proves convenient to rewrite the static energy density (ignoring the Lagrange multiplier term) where we have performed a Bogomol'nyi completion in the last step. The BPS equation is thus 14) and when satisfied, the total energy is proportional to the topological degree of the baby Skyrmion, Q.
Changing parametrization of φ to stereographic coordinates (2.7), the BPS equation now reads Inserting the Ansatz ω = e iN θ ζ(r), we get Choosing the lower sign, we get which with κ = −ξ 2 0 can be written as .

(2.19)
If we set the integration constant ξ 0 = 1, we can move the coordinate singularity to ξ = 0, obtaining the solution where ξ ≡ r R and the compacton radius is It will prove useful to calculate the BPS mass In the fourth line we have chosen the lower sign, corresponding to the boundary conditions ζ(0) → ∞ and ζ(R) = 0.
Note that the topological charge of this axially symmetric configuration is (2.23) For axially symmetric baby Skyrmions, we will use N (which is equal to Q) to denote the topological charge, whereas for more complicated configurations, Q is the total topological charge. For instance, later it will be useful to consider a Q = 4 baby Skyrmion that is composed by two N = 2 axially symmetric solutions.

Energy bound
In the previous section, we have shown that the energy is bounded from below by the Bogomol'nyi type bound for the BPS sector = 0, In addition, it can also be shown that there is a bound in the energy for the first term in the Lagrangian (2.1), The total energy in the model (2.1) is thus bounded from below by This bound is, however, only satisfied in the two limits: → 0 (BPS baby-Skyrmion limit) and → ∞ (BPS lump limit) [23].

Perturbation in
We will now consider making a perturbation in around a background solution. That is, the background, ϕ, is a BPS compacton, which minimizes the BPS submodel (2.12): For reference, it will be useful to write the explicit form of the charge-N axially symmetric compacton solution (2.20) in the vector coordinate ϕ: with R the compacton radius (2.21) and x + iy = re iθ the standard polar coordinates in We will consider the corrections to the energy order by order in the following sections.

Leading-order correction
The leading order correction, proportional to , comes from inserting the background solution into the kinetic term. However, since ϕ depends on the moduli parameters of the BPS sector, not all the possible background solutions are equivalent choices at this order.
Recently, it has been conjectured that the right choice among all the possible maps is given by the minimizer of E 2 on the moduli space of the BPS solutions [23,24]. A map respecting this request is said to be restricted harmonic and, in the case of a single compacton, it is identified by the axial symmetric solution (2.20). The complete proof of this statement is given in the next section.
The first order correction in terms of is therefore The parenthesis on the last line takes the value 1.678 for N = 1, which is about 5/3 of the energy bound for the kinetic term.
To this order, the energy reads We can now consider the energy per N as a function of N . In particular, this function has a minimum d dN 5) which is N = 4 3 log 2 7 2.180. (3.6) This means that for an axially symmetric configuration with topological charge N = 1, 2, the leading-order energy (i.e. to O( )) per charge, Q = N , decreases as a function of N . This implies that an axially symmetric 2-Skyrmion has lower energy than 2 well separated 1-Skyrmions. 2 This, in particular, means that for parametrically small , the lowest energy configuration with topological charge Q = 2M will consist of M 2-Skyrmions at a separation distance that is not determined at this order in the calculation.

Restricted harmonic maps
In the previous section, the notion of restricted harmonic maps [24] was briefly discussed. In what follows, we will formally define such a map that is used as the zeroth-order background of a field expansion in a near-BPS model. After the general definition, we focus on the near-BPS baby-Skyrme model, giving a review of the criterion obtained in ref. [24] to identify the correct restricted harmonic map. The axially symmetric solution used in eq. (3.3) is verified to respect this criterium. In the last part of this section, we further verify that even two axially symmetric compactons sitting side-by-side with a random orientation are restricted harmonic. The last result will be useful in sec. 5 in which we discuss the interaction between two near-BPS baby Skyrmions and calculate the binding energy.
Given a set of scalar fields φ a mapping d + 1 dimensional Minkowski spacetime (M, η) with flat metric η to the target space (N , h) with N = S 2 being the 2-sphere with metric h, a near-BPS Lagrangian L can be written as where L BPS describes the pure BPS sector and is a small parameter multiplying the Dirichlet term with h ab representing the metric of the target manifold N .
Taking into account only the BPS sector, we denote by ϕ a (x, λ) the generic static solution of the model that depends on the moduli parameters λ. Explicitly, we consider a BPS system for which λ consists of all the possible orbits in the group of volume-preserving diffeomorphisms.
In the limit of → 0, the total static energy E of the system (3.7), at the lowest order in , takes the form where λ * ⊂ λ is the subset of λ that minimizes the Dirichlet energy E 2 The map ϕ a (x, λ * ) that locally minimizes E 2 within the domain of the group of volumepreserving diffeomorphisms is said to be restricted harmonic.
A possible proof for the expression (3.9) can be performed assuming a Taylor expansion for the field φ in terms of , All orders of φ can be solved by minimizing the total action S order by order. It is important to note that a perturbative expansion in terms of is not always possible for the exact static solution φ, since the dependence on such a parameter could be nonanalytic. We therefore assume this hypothesis here only to give a simple justification of the expression (3.9); for a more complete treatment of the topic, see refs. [23,24].
Since we are dealing only with static configurations, we consider the static energy E to be minimized instead of the action S and we analogously write a Taylor series for the energy Hence, the zeroth-order of the static energy is given by whose static solution is by definition 14) The first-order in the expansion reads where the first term in the parenthesis on the first line vanishes due to the equations of motion. Since the functional E (1) contains only the background solution, the minimum of E (1) must be sought not among all the field configurations but only within the moduli space of the BPS solution (3.14). Therefore, following the principle of least energy (action) whose solution ϕ a (x, λ * ) is by definition the restricted harmonic map discussed above. We clarify that for an infinite moduli space, such as the group of volume-preserving diffeomorphism, the derivative with respect to λ does not have a mathematically well-defined meaning and we use it here only to simplify the notation. Summing the zeroth-order (3.13) and the first order (3.15) of the energy calculated with the solution ϕ a (x, λ * ) we obtain the expression (3.9), as we wanted to prove.
After the formal definition, we need some practical tools to calculate the restricted harmonic maps of a given system. A mathematical criterion that is able to identify whether a map is restricted harmonic (or not) has been proposed in ref. [24]. In the following, we review the theorem given there, which we use for the near-BPS baby-Skyrme model in the case of a single and two compactons.
Given a smooth map φ from the manifold M with metric g = g ij dx i ⊗ dx j to the manifold N with metric h = h ab dφ a ⊗ dφ b , the Dirichlet energy is defined as where d is the number of dimensions of the space M. Using the map φ we can construct the pull-back φ * h of the metric h to M which is necessary for what follows. Among all the maps φ with finite Dirichlet energy connected by a volume-preserving diffeomorphisms, a mapφ is restricted harmonic if and only if the one-form divφ * h on M is exact [24]. It is useful to recall that the divergence of a symmetric (0, 2) tensor ω = ω ij dx i ⊗ dx j on M it is nothing but with the connection Γ (Christoffel symbols) defined as Using this theorem we verify that the axially symmetric solution used in eq. (3.3) is restricted harmonic. Using polar coordinates on M = R 2 and the vector notation φ a for the field, with the constraint φ a φ a = 1, we write a generic Q = N axially symmetric compacton in the form where the function f depends only on the radial coordinate and α is the orientation phase [14]. With this choice, the metric h reduces to the standard euclidean metric and the pull-back φ * h of h can be written as with dx i = (dr, dθ). Note that there is no longer any dependence on the orientation phase α. Taking the divergence of this tensor yields resulting in a one-form whose exactness we must prove. Here, it is useful to recall that according to Poincaré's lemma all closed forms on a contractible manifold are exact. Hence, to complete the proof we must verify that d div φ * h = 0, with d being the exterior derivative. Explicitly, which gives the necessary and sufficient condition for eq. (3.21) to be restricted harmonic.
The same proof can be extended to the case of two axially symmetric compactons sitting side-by-side without overlap. A solution for two separated compactons with Q = N + N and random orientations can be written as where φ a 1 , φ a 2 are axially symmetric with respect to the points (−x 0 , 0) and (x 0 , 0) respectively, i.e., with and α, β are two independent phases. Note that in order to have two separated compactons we require |x 0 | > R, where R is the compacton radius (2.21).
Following the same procedure as adopted above, we calculate the pull-back of the metric h induced by φ: where dx i = (dx, dy) are the Cartesian coordinates in M = R 2 . In the last equality of eq. (3.28), we used the propriety of the compacton that the fields are constant outside of its radius so that At this point it is useful to manipulate the expression (3.28) by a series of change of variables in order to write φ * h as a sum of symmetric (0, 2) tensors of the form (3.22). Hence, where x i ± = (x ± , y ± ) are defined by whilex i ± = (r ± , θ ± ) are the polar coordinates Defining for simplicity the final result for the two tensors is which takes the same form as that of eq. (3.22), but with translated coordinates.
Applying the divergence to eq. (3.33) and using the linearity we obtain a one-form whose exactness we must verify. Using Poincaré's lemma again as well as the linearity of the exterior derivative, we get which finally proves that the map (3.25) is restricted harmonic.
The results of this section show that both an axially symmetric compacton with random orientation and two nonoverlapping axially-symmetric compactons, with random orientations are restricted harmonic maps. In this paper, we therefore use the first solution as the zeroth-order background for a single near-BPS baby Skyrmion, as already anticipated in sec. 3.1, and the second solution for two near-BPS baby Skyrmions.
With this last result, we are now able to calculate the energy of two near-BPS baby Skyrmions with topological charge Q = N + N at the leading order in , as already performed for the single case Q = N in eq. (3.4). Using the restricted harmonic map (3.25) and the properties (3.29), we get (3.38) where N ⊕ N means an axially symmetric charge-N compacton side-by-side of another axially symmetric charge-N compacton. At this order, the energy of the Q = N + N configuration is therefore simply the sum of the energies of the two components with no information about the interaction between two near-BPS baby Skyrmions. An investigation of the next-to-leading order (NLO) is required to calculate the binding energy of such a configuration.

NLO and N 2 LO corrections
We will now consider the next-to-leading order (NLO) and next-to-next-to-leading order (N 2 LO) corrections to the energy, which corresponds to taking into account the corrections of order O( 2 ) and O( 3 ), respectively. To this end, we will perform a linear perturbation of the model For the NLO and N 2 LO corrections, we need to calculate the variation up to third order (in the fields) of the Lagrangian (2.1): where the | denotes that the expression (to the left of the bar) is evaluated on the background by setting φ = ϕ and we have defined (3.41) We have furthermore replaced the Lagrange multiplier constant λ with an expansion λ → λ 0 + δλ, where it is understood that λ 0 is the Lagrange multiplier that solves the sigmamodel constraint for the background BPS solution. δλ is then a Lagrange multiplier that ensures that the perturbation does not bring the total field φ = ϕ + δφ out of the O(4) group, i.e. it should still preserve the unit length constraint φ · φ = 1 + O( 4 ). More precisely, if we vary the above Lagrangian density with respect to δλ we obtain which is solved by [12] which satisfies the constraint (3.42) up to O(∆ 4 ), which we shall verify is of order O( 4 ) and hence will give a contribution to the energy of order O( 5 ). Notice that the above form automatically restricts ∆ to be orthogonal to the background solution ϕ to leading order.
The Hessian of the Lagrangian density with respect to the derivatives of the fields, V µν ab , is given by the derivative of the Hessian is and the Lagrange multiplier of the background, λ 0 , reads Substituting the form of the variation (3.43) into the perturbation Lagrangian density (3.40) yields which is the complete Lagrangian up to third order in , (i.e. discarding terms of O( 4 )).
Notice that the term with δλ is first needed at the fifth order in and can thus safely be ignored -its job was to produce the form (3.43).
The reason for keeping terms up to third order in , as we shall see, is to retain the last term in eq. (3.49), which will enable the fluctuation to propagate outside of the compacton and hence give rise to binding energy. In order to calculate the energy consistently, we have thus kept all other terms up to third order in . However, for calculating ∆ we will linearize its equation of motion where we have used that V µν cb is symmetric under the simultaneous exchange of µ ↔ ν and c ↔ b.
Let us now consider what happens outside of the compacton. Since ϕ a → δ a3 is at its vacuum and all derivatives of the background field vanish identically, most terms including the source term (right-hand side of eq. (3.50)) switches off. In fact, the only remaining terms outside of the compacton background are V µν 1ab ⊂ V µν ab and the λ 0 terms in eq. (3.50). Thus the linearized equation of motion of the perturbation, ∆, reduces outside of the compacton background to It is interesting to see that the kinetic term here (outside the compacton) is given by V µν 1ab which is third order in (i.e. O( 3 )). Thus, had we only kept terms up to second order in in the Lagrangian (3.40), the equation of motion would have been a constraint setting ∆ = 0. If ∆ = 0 outside the compacton, then there is no information propagating between two compactons and hence there is no binding energy.
A natural question then arises: At which order in is the binding energy of baby Skyrmions captured? If we do not restrict to the linearized equation of motion for the fluctuation, it is clear that the above conclusion about the equation outside of the compacton does not change (although we do not prove this). It is instructive to look at the solution for the fluctuation outside of the compacton as the solution to eq. (3.51), which reads where c 1,2 are constants. This solution illustrates perfectly the problem of describing the binding energy at a specific order. In fact, although the tail does not vanish for a nonvanishing , any order in vanishes due to the exponential. In some sense, this solution is all-order in or "nonperturbative".
A word of caution is that although we include this "nonperturbative" effect of the allorder tail for the linear perturbation, we solve a linearized equation for the perturbation which only captures every effect to second order in and some effects -like the crucial tail (3.52) -at third order in . Nevertheless, we calculate the full energy to third order in using this solution to the linearized equation of motion.
Writing out eq. (3.50) explicitly in the static case, we get where we have defined The energy density of the perturbation can be written as with the second order (NLO) in : where the former expression contains only quadratic terms in ∆ which contribute to the linear equation of motion and the latter gathers the cubic terms.
Notice that the entire perturbation energy vanishes manifestly for ∆ ∝ ϕ. Due to the nonlinearity in the form of δφ in eq. (3.43), a mixing can occur if both a transverse ∆ ⊥ and a longitudinal perturbation δc(x)ϕ is turned on at the same time. Inserting ∆ = δcϕ + ∆ ⊥ into quadratic part of the perturbation Lagrangian yields The solution for δc -which is sourceless -with vanishing boundary condition at r = 0 and r → ∞ (as measured from the compacton origin) yields δc = 0. Because of the absence of a quadratic term for δc, variation neither gives an auxiliary equation (algebraic) nor a dynamic equation for δc. We checked that the quadratic terms in δc do not appear in the cubic part of the Lagrangian either and hence can first appear at order O( 4 ), which we shall not consider in this paper 3 . We shall thus only consider ∆ = ∆ ⊥ which is transverse to the background solution ϕ, i.e. for which it holds ϕ · ∆ ⊥ = 0.
At this point, it will prove useful to specialize to the case of the background BPS solution for ϕ, the transverse perturbations for ∆ = ∆ ⊥ and switch to polar coordinates in R 2 , for which the perturbation energy reads for the NLO terms, for the NNLO terms quadratic in ∆ and for the NNLO terms cubic in ∆, and we have defined the transverse ∆ perturbations and the background BPS compacton solution is described by the radial profile with R being the compacton radius of eq. (2.21).
The corresponding equations of motion can be written as where we have defined the matrices First notice that there is only a source term (right-hand side of eq. (3.72)) for the upper equation, viz. for δf . Notice then that the off-diagonal components reside only in the matrices X θ , X rθ and hence couple the second equation to δf θ and δf rθ . The source for the δθ perturbation is thus the coupling to the nontrivial θ behavior in δf . It is now clear that, if we restrict to axially symmetric perturbations such that δf = δf (r), the second equation decouples and is sourceless; hence it is satisfied by the trivial solution δθ = 0, everywhere.
The perturbation governing axially symmetric baby Skyrmions is thus simply described by the equation and the corresponding energy for the axially symmetric perturbation is for the NLO terms, for the NNLO terms quadratic in δf and for the NNLO terms cubic in δf .

Axially symmetric solutions
In order to verify our perturbative scheme, we start with axially symmetric baby Skyrmions. It will prove useful to interpret the perturbation δf as follows where we have used eq. (3.70) and set δθ = 0. It is now clear that δf indeed is an additive correction to the BPS background profile function f .
The naive attempt is to solve eq. (3.83) with boundary conditions Had the soliton background solution been smooth, it would probably have worked out. However, the first derivative of the background solution is discontinuous at r = R with the compacton radius R given by eq. (2.21). Since the full soliton solution is smooth, the perturbation should counteract this cusp (jump in the first derivative), which we can write as a new condition at r = R: which with the perturbation δf of eq. (3.70) yields (3.90) We will refer to this extra condition in addition to the conditions (3.88) as the "new boundary conditions." In fig. 1 is shown the profile function f = arccos(φ 3 ) for the N = 1 axially symmetric baby Skyrmion with = 0.01 obtained by different methods. Since the BPS solution is a good starting point, we subtract off the BPS profile to better see the differences in the panel on the right-hand side. The naive implementation with a smooth perturbation is shown as the dark-green curve and does not resemble an improvement with respect to the BPS background solution compared to the exact solution (numerical). The perturbation subject to the new boundary condition (3.90) (orange curve) on the other hand gives an incredible improvement over the BPS solution and the discrepancy with respect to the exact solution (numerical) is really tiny and is expected to be due to nonlinearities (since we have linearized its equation of motion (3.83)).  Figure 2. The mass of the baby Skyrmions in the perturbative scheme as a function of (on a logarithmic scale): The red line is the BPS bound, the orange curve is the leading order (LO) correction, the black pluses are the corrections calculated at next-to-leading order (NLO) using the linear perturbation, the green-dashed curve is a fit to the latter points, the red crosses are the corrections calculated at next-to-next-to-leading order (N 2 LO) using the same linear perturbation, the dark-blue curve is a fit to the latter points and the black curve is the exact energy calculated using the full (nonlinear) equations of motion. The left columns show the total energy while the right columns show the energy relative to the LO correction.
We are now ready to calculate the perturbative corrections to the energy of the N = 1, 2, 4 axially symmetric baby Skyrmions and the results are shown in fig. 2. The corrections are calculated using the linear perturbation, which is a solution to equation (3.83) subject to the new boundary condition (3.90). The NLO energy is calculated using the energy to order 2 using eq. (3.84), whereas the N 2 LO energy is calculated as the sum of eqs. (3.85) and (3.86). We can see that the perturbative LO corrections overshoot the energy for every N and the NLO corrections correct this overshooting too strongly, yielding an undershoot. Including finally the N 2 LO corrections, the result is extremely close to the exact numerical one. Notice by careful inspection of the figure, that the red crosses actually match the exact result (black curves) better than the fit (dark-blue curves). Although we write the expansion to fourth order in , it should not be trusted beyond the third order. The higher orders simply represent the all-order contributions from the tail that are important for capturing the binding energies. N is shown in fig. 3. As can be seen

Numerical calculations
In this section, we present numerical solutions to the full equations of motion. This is a computationally extremely expensive work and can be carried out only for the 2-dimensional case unless very clever adaptive methods are being used. The numerical calculations performed here are carried out with a fourth-order 5-point stencil finite-difference method using the arrested Newton flow [26] on a square grid of sizes up to ∼ 3226 2 ∼ 10 7 with lattice spacing down to ∼ 0.0037. Our numerical accuracy is about 10 −6 or better. Such expensive grids would take a very long time on conventional CPU clusters, so our code is written in CUDA C and is executed on a GPU cluster. In order to avoid warp divergence, the code is made such that the bulk of the lattice is launched as one kernel without any conditionals and the edges are launched as two different additional kernels corresponding to the vertical and horizontal edges of the lattice (also without any conditionals).
In this section we explore the numerical solutions of baby Skyrmions that do not possess axial symmetry, however, we also check the axially symmetric configurations both for ensuring that they exist in the given range of parameter space and for comparing our numerical accuracy to that obtained using ODEs for the axially symmetric system.
Our aim here is to confirm whether the energetically most favorable configuration -in the case of small but finite 0 < 1 -consists of N = 2 baby Skyrmions sitting next to each other at a finite or vanishing distance between them. Such a scenario is supported by the notion of restricted harmonic maps.
To this end, we start with two N = 2 baby Skyrmions (hence with a total topological charge Q = 2 + 2 = 4), situated next to each other. Some of the numerical results are shown in fig. 4. The figure is organized into four columns showing the topological charge density, the field configuration using a color scheme, the total energy density and finally the energy density of the kinetic term (− L 2 ). The colors on each graph is relatively rescaled to the values of the graph. The color scheme for the second column of fig. 4 illustrates the orientation of the fields at each point of the configuration space (R 2 ). The scheme is defined such that the hue is correlated with the complex phase arg(φ 1 + iφ 2 ) and the lightness is determined by the value of φ 3 . Explicitly, φ 3 = 1 is white, φ 3 = −1 is black, arg(φ 1 + iφ 2 ) = 0 is red, arg(φ 1 + iφ 2 ) = 2π/3 is green and arg(φ 1 + iφ 2 ) = 4π/3 is blue.
In particular, we find that the two N = 2 baby Skyrmions next to each other is always stable -even for = 1. What happens when is decreased, is that the hole in each N = 2 soliton is shrunk to a point and asymptotically, the two Skyrmions become axially symmetric -almost unaware of each others presence. Notice that the orientation of the two baby Skyrmions is always in the attractive channel, i.e. such that e.g. two patches of red touch each other, see the second column of fig. 4. At larger values of , the two N = 2 solitons are elongated along the axis joining their centers. This can be interpreted as an effect due to the strong binding force between them. For small , this effect goes away and the binding energy is also drastically reduced. We will discuss the binding energies separately in the next section.
Baby-Skyrmion solitons are free to deform their shape and are able to split up to the  . The Q = 2 + 2 baby Skyrmion solution as a function of , which is the stable Q = 4 solution. The columns display the topological charge density, the field orientation using the color scheme described in the text, the total energy and finally the kinetic term − L 2 . We note that the two stretched N = 2 baby Skyrmions become almost perfectly round in the last row (for = 0.0207). lowest energy state. Nevertheless, many local minima of the energy functional exists and they correspond to metastable states. In order to know which configurations are the stable ones, we start with many different initial guesses and see what they converge to by the numerical flow.  The columns display the topological charge density, the field orientation using the color scheme described in the text, the total energy and finally the kinetic term − L 2 . The last row (i.e. for = 0.144) shows that the solutions has decayed into a lower-energy state with tetrahedral symmetry instead of axial symmetry.
We expect from the fact that N ≈ 2, that the N = 4 axially symmetric soliton will be only metastable. As an explicit check, we perform the numerical calculations which are shown in fig. 5. Perhaps surprisingly, it turns from metastable at ≈ 1 into unstable around ≈ 0.15, see the last row of fig. 5. The numerical algorithm finds a lower-energy configuration, which is composed by four N = 1 Skyrmions attached to each other in a tetrahedral arrangement. This is somewhat surprising, because as we shall see shortly, the N = 1 soliton is actually not stable for = 0.15 unless it is alone. Once it is near another N = 1 soliton, they will merge into a deformed N = 2 soliton at that . The high degree  . The Q = 1+1+1+1 tetrahedral baby-Skyrmion solution as a function of . The columns display the topological charge density, the field orientation using the color scheme described in the text, the total energy and finally the kinetic term − L 2 . We note that the tetrahedrally symmetric baby Skyrmion exists also for large (i.e. = 1). of discrete symmetry -tetrahedral symmetry -somehow prevents that from happening.
There are two possibilities: Either there is a crossover, so that the tetrahedral config-  . The Q = 1 + 1 baby-Skyrmion solution as a function of . The columns display the topological charge density, the field orientation using the color scheme described in the text, the total energy and finally the kinetic term − L 2 . uration becomes stable and the N = 4 axially symmetric soliton becomes unstable around = 0.15 or both the tetrahedral and axially symmetric solitons exist for 0.15 and the axially symmetric solutions simply becomes unstable around = 0.15 such that the decay to the nearest (in energy) metastable state is the tetrahedral baby Skyrmion. In order to determine the phase diagram of the Q = 4 sector, we need to know the energies of the tetrahedral soliton also for larger values of . Thus we use the tetrahedral configuration as a seed (initial guess) in the calculation and calculate its energy for larger all the way to = 1, see fig. 6. This reveals that the tetrahedral Q = 4 baby Skyrmion exists also for large values of (even though two N = 1 baby-Skyrmions next two each other are unstable, see  The columns display the topological charge density, the field orientation using the color scheme described in the text, the total energy and finally the kinetic term − L 2 . below).
Then to determine the phase diagram of the Q = 4, we plot the total energies of the different solutions in fig. 9. From the figure -especially the right-hand side panel -we can see that the N = 4 axially symmetric configuration is metastable for 0.15 with the highest energy. The tetrahedral soliton is metastable with an intermediate energy, but not too far above the lowest-energy state. Finally, the two N = 2 Skyrmions side-by-side is the most energetically favorable solution in the entire range of the considered here. For = 0.0207 it is rather convincing that the → 0 limit turns the solitons into a lattice of axially symmetric solutions sitting next to each other with the most stable axially symmetric component being the N = 2 baby Skyrmion. The contact to neighboring solitons becomes point-like in the → 0 limit. This is consistent with the restricted harmonic property that the near-BPS (baby) Skyrmions should possess.
We will now turn to the Q = 2 sector. Again, as N ≈ 2, we know that the N = 1 axially symmetric baby Skyrmion is stable only in isolation. When exposed to more baby Skyrmions, it has to be either metastable or unstable; that is, it should be energetically favorable to combine into N = 2 Skyrmions. It turns out that for = 1, two N = 1 axially symmetric baby Skyrmions next to each other are unstable and they immediately combine into a single N = 2 soliton with axial symmetry [25]. However, for small enough , the two N = 1's side-by-side become metastable, see fig. 7. The critical 2 where the baby Skyrmion departs from axial symmetry is around ≈ 0.2, whereas the critical 1 where the baby Skyrmion is clearly composed of two individual solitons is around 1 ≈ 0.1. For < 0.06 the two individual solitons are almost only connected by a "single point." As a check on our numerical calculations, we compute the N = 2 axially symmetric baby Skyrmions as well, see fig. 8. The figure shows that the hole in the middle of the soliton is shrunk to a point in the limit of → 0.
To check how these solutions fit into our picture, we show their total energies in fig. 10. We can thus confirm that the phase diagram in the Q = 2 sector is determined by the LO correction to the energy, which makes the two N = 1 baby Skyrmions sitting side-by-side metastable for 0.15 and unstable above. We can also confirm by comparing the 2dimensional PDE calculations with the 1-dimensional ODE ones (for the N = 2 case), that our numerical accuracy is incredibly good for the entire range of s considered here.
In order to complete our perturbative scheme, we need to calculate the binding energies between two Skyrmions, which will be the topic of the next section.

Binding energies
In this section we will discuss the binding energies of two N = 2 baby Skyrmions sitting side-by-side and, in the last part, also the case of two N = 1 side-by-side. We have performed very high-resolution numerical calculations that gives us the answers we seek for the baby Skyrmion case. However, for the 3-dimensional Skyrmions, the numerical calculations become much more difficult -even with the utilization of adaptive methods. Therefore we want to push the perturbative approach to capture the physics and respective energies of the composite near-BPS solitons.
Since we have the compacton solutions, the background is analytically known. The difficulty is to impose the boundary or cusp condition (3.90) that the derivative of the perturbation fields obey at the boundary of the compacton, but this can easily be done in polar coordinates. Furthermore, the composite configuration of two baby Skyrmions side-by-side will need a gluing condition, which we shall discuss shortly. Finally, the perturbations must go to zero at spatial infinity, which is the easiest condition to impose.
In section 3.2 we showed that the restricted harmonic map for the Q = 2 + 2 configuration consists of two separated axially symmetric compactons with random relative orientation. However, if for a single baby Skyrmion the problem of the background solution is solved by the restricted harmonic condition, for the multi-soliton case further considerations are needed. To understand the reason for this difference, we must consider the physical meaning of the zeroth-order background field. This field represents the limit of the exact solution φ of the Lagrangian (2.1) when tends to zero, i.e. the final configuration that we obtain if we adiabatically switch off the interaction among the solitons. The static solution φ exists only if the baby Skyrmions interact with an attractive force that depends on the relative orientation among them [12]. Therefore, we expect that even in the limit of vanishing interaction ( → 0) the final configuration will conserve a particular relative orientation. Then, we conclude that for the near-BPS multi-baby Skyrmion case, the restricted harmonic condition is not sufficient to identify the correct zeroth-order solution since it does not yield any restriction on the relative orientation (or distance) among axially symmetric compactons.
In order to fix the background configuration among all the restricted harmonic maps, we propose to identify the right map by looking at the form of the interaction between two well-separated baby Skyrmions. The potential V among these solitons has been calculated in ref. [12] for the baby-Skyrme model (the same Lagrangian (2.1) with = 1) and for the Q = 2 + 2 case it takes the form where R is the distance between the soliton centers, α and β are the two respective phases and m is the (perturbative) pion mass. From this expression we recognize that, for fixed distance R, the minimum of the potential is obtained for α − β = 0, i.e. when the two baby Skyrmions have the same orientation. Therefore, even if the form of the potential (5.1) holds only at large distances, we guess that in the limit of vanishing interaction (when → 0) the two solitons keep their relative orientation unchanged.
Instead, for what concerns the relative distance, we expect that in the limit of vanishing interaction, the two baby Skyrmions must flow to a configuration in which they overlap their tails less and less. However, since the tail of a near-BPS baby Skyrmion vanishes rapidly with as is clear from eq. (3.52), the two solitons flow to a distance that is very close to the sum of their radii.
Summarizing, we guess that in the limit → 0, the background solution for the Q = 2+2 configuration is given by two axially symmetric compactons with the same orientation touching each other at a single point. This guess seems to be confirmed by numerical calculations in fig. 4 and we will thus use it in the following.
Since the zeroth-order configuration for the two compactons is such that the other baby Skyrmion is a copy of itself, spatially translated by a distance of its own diameter (roughly), it is not necessary to calculate the perturbation over the entire plane, but we can reduce the problem by identifying appropriate boundary conditions for a field living in an appropriate half of R 2 .
The boundary conditions are sketched in fig. 11. In addition to setting the perturbation to zero at r = 0 and imposing the condition (3.90), we need to impose the following Figure 11. Coordinate system (r, θ) for perturbation of the two N = 2 baby Skyrmions sitting side-by-side. Boundary conditions need to be imposed at O and at the compacton boundary r = R as well as at x = 0 (at re iθ = R + a + iy) for attaining the ability of gluing the two solitons together.

boundary (gluing) conditions
which ensure that the soliton can be glued together with its partner. The second boundary condition is "odd" because this field is odd under parity transformations and the only way to connect the perturbation to a field that is not flipped in x → −x (because that would make it an anti-Skyrmion) is to impose the boundary condition as given in the above equations. Because the condition (3.90) is easier in polar coordinates, we prefer to impose the x = 0 gluing conditions on the equation in polar coordinates.
In terms of the perturbation fields given in eq.  Fig. 12 shows the boundary conditions for the lattice problem (the discretized PDEs) at the x = 0 boundary. The lattice point (I + 1, j) is determined by the condition (5.9) which is written in discretized form as with β the offset of the x = 0 boundary from the second-last lattice point in the calculation.
It will be instructive to write the perturbations δf, δθ as follows  where we have used eq. (3.70).
We now turn to solving the coupled PDEs (3.72) with the boundary conditions (5.5) and (5.9). As a good starting point, we take the perturbation of the axially symmetric compacton for given and N as the initial condition. What the PDE problem then boils down to, is to implement the gluing conditions (5.9) at x = 0 (i.e. midway between the two compactons). The gluing conditions break the axial symmetry which provides nontrivial θ dependence for δf which in turn acts as a source for δθ.  The only θ-dependence is induced by the gluing condition (5.9) near the right-hand side boundary of each panel; the condition on δf is slightly cumbersome as δf 2 / cos 2 (N θ) must obey the Neumann boundary condition at x = 0. The effect is most visible for small values of a and as soon as a is of order of the side of the "yellow" ring, which is the perturbation deformation due to the cusp condition, the θ-dependence is almost negligible. This can also be seen in the second column of the figures, where a nontrivial solution for δθ is most pronounced for a = 0. The θ-dependence in δθ at the x = 0 boundary is visible as positive (negative) values for y negative (positive). This effect persists when a is increased, but the magnitude of the values of δθ at the x = 0 boundary become exponentially suppressed in line with the tail of the perturbation around the compacton does. The fact that the solutions are N = 2 baby Skyrmion is clearly visible in columns 3 and 4, where the cusp in δf switches sign four times as one goes around the compacton border. The induced θdependence is also very clear for a = 0 and a = 0.0221 in column 5 in the third component of the vector perturbation, δφ 3 . The last column in the two figures shows the energy density of the perturbation and it is not easy to the result of the integration with the naked eye. In fact, the largest effect of the gluing is geometric, meaning that cutting part of the tail of the perturbation at x = 0, diminishes the contribution to the energy. One could anticipate that the gluing of the two compactons would give a positive contribution counteracting the geometric decrease in the perturbation energy; although that happens for the field δf , the induction of θ-dependence turns on nontrivial behavior for δθ which lowers the energy slightly. At the distance a = 0.243, there is almost no visible effect of the gluing conditions and the energy is very close to the sum of two compactons. Fig.13 shows the case for very small = 0.01 and in fig. 14 = 0.0428 is slightly larger. The effect is simply an amplification of perturbations and in particular, the thickness of the perturbation in δf due to the cusp condition is far larger for = 0.0428 than for = 0.01. Fig. 15 shows the perturbative contribution to the energy at next-to-next-to-leading order (N 2 LO) for two N = 2 baby Skyrmions as function of the separation distance 2a, see fig. 12. The four panels show = 0.01, = 0.0207, = 0.0428 and = 0.0886, respectively and for all panels, the smallest energy (meaning the largest negative contribution to the energy) is at a = 0, which means when the two compactons exactly touch each other at a point. For the three smallest values of (i.e. = 0.01, 0.0207, 0.0428), the perturbation energy turns out to be lower than the geometric energy for small a (and equal at large a). The geometric energy is simply the perturbation energy of the axially symmetric compacton cut off at x = 0. However, for = 0.0886 the perturbation energy crosses over the geometric line at a ∼ 0.12 and is slightly above the geometric one at a = 0, indicating that the gluing condition builds up some tension or excess energy in the perturbation for this "large" value of . Of course, our perturbation scheme is best trusted at small , so this may well be an indication of the approximation starting to deteriorate.
We finally, compare the N 2 LO result of the perturbation theory of the energy of two N = 2 baby Skyrmions sitting side-by-side -taking into account the binding energywith the full numerical (brute-force) PDE calculations, see fig. 16. The two results are incredibly close to each other, in particular for small values of . Since the binding energy requires a PDE calculation as well, it is not clear which result is more accurate. We expect the perturbative result to be most accurate of the two for small ∼ 0.01, whereas the full numerical PDE calculations are most accurate for large values of ∼ 0.1.
Before summarizing the results for the binding energies, we will repeat the perturbative N 2 LO calculation of the two baby Skyrmions side-by-side, but this time for N = 1 compactons. This solution is only metastable for small while for 0.15 it is unstable. Nevertheless, since we have performed the full PDE computations of these composite baby Similarly to the case of Q = 2 + 2, for the case of two N = 1 baby Skyrmions side-byside, we must specify the background solution around which we calculate the perturbation (5.11). Again the restricted harmonic condition does not identify uniquely the right compacton configuration that is needed since it does not specify any condition on the relative orientation (or relative distance). Following the discussion of the Q = 2 + 2 case, we look at the long-range potential of two well-separated N = 1 baby Skyrmions [12]: where again R is the relative distance, α and β are the two respective phases and m is the pion mass. The potential (5.13) generates the maximum attractive force in case of two baby Skyrmions with opposite orientation, i.e. α = β + π. Therefore, in contrast to the Q = 2 + 2 case, among all the axially symmetric compactons side-by-side, we choose two N = 1 solitons with opposite orientation as the background solution. Our guess is confirmed by the numerical simulation in fig. 7.
Similarly to the case of two N = 2 compactons side-by-side, figs. 17 and 18 which are for two N = 1 compactons side-by-side, show the left-hand side of the composite soliton configuration with the gluing condition (5.9) imposed on the right-hand boundary of each panel for = 0.01 and = 0.0428, respectively. The other edges of the panels are not the limit of the calculations, which have merely been cropped so as to render the content as clear as possible. The two figs. 17 and 18 again display the 6 columns: the perturbations δf and δθ as well as the three components of the vector perturbations δφ, and finally the perturbation energy density at N 2 LO. Similarly to the N = 2 case, the largest θdependence is induced for very small separation distances 2a and especially for a = 0. It is clearly visible from the third and fourth columns of the figures that the compacton background solutions are N = 1 baby Skyrmions, where cusp condition only changes sign twice around the perimeter of the compacton. As before, the width of the perturbation due to the cusp condition grows with and hence is more pronounced in fig. 18   symmetric compacton at x = 0. In particular, there is no crossing of the N 2 LO result and the geometric line, even for = 0.0886. This is probably because the gluing condition is more relaxed for the N = 1 compactons, since the winding of the background solution is less (i.e. minimal). In particular, this eases the gluing condition (5.7) which is more strict for the N = 2 compacton's δφ 2 component.
We will now compare the N 2 LO result of the perturbative contribution to the energy of the two N = 1 baby Skyrmions sitting side-by-side (taking into account the binding energy between them) with the full numerical (brute-force) PDE calculation in fig. 20. to calculate any binding energy. Hence this calculation is simply a test of the perturbative scheme on the two N = 1 baby Skyrmions sitting side-by-side, even though they are only metastable. For large 0.15 we know from the full numerical PDE computations of sec. 4 that the two N = 1 baby Skyrmions side-by-side become unstable and merge into an axially symmetric N = 2 baby Skyrmion. This fact is known already from the leading order (LO) result of the perturbation theory, see eq. (3.6). We can thus anticipate that the N 2 LO perturbative scheme will become inaccurate for large , since the instability is due to the LO result and, in principle, not known to the NLO and N 2 LO computations. For this reason, for reference, we show also the single N = 2 baby Skyrmion (with axial symmetry) in fig. 20 as magenta pluses (full PDE results) and a black solid line (ODE results). From fig. 20(b) we can see that the perturbative result for the N 2 LO contribution to the energy works impressively well for 0.033 but then deviates and smoothly tends to the curve for the single axially symmetric N = 2 solution (black solid line).
We have now confirmed by comparison with full numerical PDE computations, that our perturbative scheme at N 2 LO works very well for the two N = 2 baby Skyrmions sitting side-by-side at least for 0.1 and also for the two N = 1 baby Skyrmions side-by-side for 0.033 (due to the instability of this solution). Since we already have a quite precise perturbative formula for the baby Skyrmion energy for small N and small at N 2 LO, see eq. (3.91), it will henceforth be worthwhile to separate out the binding energy of a "single bond" between two baby Skyrmions at N 2 LO. This is straightforwardly extracted from the data used in figs. 16 and 20. The result for two baby Skyrmions with total charge 2N sitting side-by-side is hence given by where we have not extracted an N -behavior, since for N = 1 and N = 2 data only we cannot determine whether the is a quadratic or a linear behavior (in N ) in the coefficients. Obviously, the constant term in both of the above expressions is expected to vanish. We have, nevertheless, included it because forcing it to be zero significantly worsens the fit. At this stage it is not clear whether the nonvanishing constant term is numerical error or there is a contact term in the binding energy, although we did not predict such a term.

Composite solutions
One could now in principle construct "nuclei" of any Q by rotating the constituent compactons such that the colors are matching at the "bond" and to a first approximation, the total energy can be calculated from eq. (5.15) with the addition of E binding ( , N ) for each "bond" binding the nucleus. The near-BPS property allows one to tune until realistic binding energies are obtained. Of course, this only provides one with the ground state energies and further development is needed for calculating the excitational spectrum. And of course, the baby Skyrme model is just a toy model. Figure 22. A sketch of a possible lattice for nuclear matter in the baby-Skyrme model. Each constituent baby Skyrmion is placed such that the colors match between the latter and the neighboring baby Skyrmion match.
We could further contemplate how the limit of large nuclei or nuclear matter would look like. In fig. 22 we show a possibility for an infinite "crystal" lattice of which large nuclei could be cut out from or nuclear matter could be made of. This particular lattice could be continued indefinitely. Nevertheless the shown lattice is made of ten baby Skyrmions and the energy to a first approximation could be calculated for any small value of by taking into account 19 bonds in the formula (5.15).
Or course, the lattice phase would be determined by minimizing the energy and it could be for some that the square lattice is energetically preferred to this triangle lattice, displayed in fig. 22. Nevertheless, it is expected that for sufficiently large , the large-Q solution becomes a chain [25] and the lattice would then only exist at finite density.

Conclusion and Discussion
In this paper, we have studied the near-BPS regime of the baby Skyrme model with the standard pion mass term as the potential and the kinetic term as the BPS-breaking perturbation with coefficient . The BPS solutions attain the BPS bound and hence the largest contribution to the energy is the BPS mass. The leading order contribution comes from the kinetic term and is of order . To this order, there is no contribution to the binding energy and no forces between the compactons. Nevertheless, we find from the leading-order result that the N = 2 baby Skyrmion is the stable solution for small . In order to get a nonvanishing tail outside the compactons, we must go to third order in and not truncate the solution to a finite order in , that is, the solution is of the form exp(−mr/ √ ), which in some sense is an all-order solution in . We test our perturbative scheme on the axially symmetric compacton solutions and calculate the N 2 LO corrections to the energy obtaining impressively good results. A key to making the perturbation around the compactons capture the physics of the baby Skyrmion is a delicate cusp condition that must be imposed on the compacton boundary. We then turn to a large-scale brute-force computation of the full PDEs and obtain numerical solutions for the baby Skyrmions for values of in the range [0.01, 1] and clearly observe that the baby Skyrmions tend to almost unperturbed compactons that touch each other at a single point in the attractive channel. We also find a new solution, by studying the N = 4 axially symmetric baby Skyrmion which for 0.15 collapses to a new metastable solution composed of four N = 1 solutions in a tetrahedral arrangement -this solution is to the best of our knowledge new. The surprise about this solution is that the neighboring N = 1 solutions should merge into an N = 2 solution, but this does not happen for the tetrahedral solution, probably because of the discrete symmetry yielding a delicate balance of the solution. We also study the two N = 1 baby Skyrmions side-by-side and find they are metastable for small values of 0.15. Finally, we turn to calculating the binding energies of two N = 1, 2 baby Skyrmions sitting side-by-side using our perturbative scheme. This is possible by imposing a gluing condition on the perturbation field and the results are in very good agreement with the brute-force numerical computations for small values of 0.1. Finally, we fit the results we have obtained, yielding an energy formula as a function of and N and the binding energy for each bond is fitted for the cases of N = 1 and N = 2.
In this paper, we considered the case of the pion mass term as the potential (2.5). For this potential, in the near-BPS limit, the mass of the "pions" (the perturbative particles) goes like 1 √ due to the in front of the Dirichlet (kinetic) term; thus at a sufficiently small they become even more massive than the baby Skyrmion itself. Clearly, in QCD terms, this would be phenomenologically unacceptable. There are other near-BPS limits that resolve this problem. For example, by choosing a different type of potential for the main BPS part of the model, for which the contribution to the pion mass vanishes, e.g. that of eq. (2.4) with p bigger than one. One could then add the massive term in the perturbation together with the kinetic term (so that it also is of order in the Lagrangian). These types of near-BPS limits have been discussed in the Skyrme literature, see for example refs. [9,23,27,28]. The main new observation we wish to make here is the following. We showed that the binding energy for the compacton case goes like 2 and not as as would be expected from the restricted harmonic argument. Thus a considerable small binding energy could possibly be obtained without pushing the pion mass too high.
Both the analytical guess and the numerical type of work and in particular the semianalytic numerics we performed in the present paper are tailored to the compacton case, so the non-compacton case will require a different approach or at least crucial modifications.
So far we know of several types of near-BPS solitonic models which can produce small classical binding energies. There are big differences among them, and it is not clear which one is most suitable for phenomenological applications to nuclear physics. For example there are cases in which the small binding energy is achieved by nuclei with size much smaller than the separation, and a relative position fixed by a potential which can be computed in the linear approximation [9,29], i.e. quite different from the type of bound state studied in the present paper. Sometimes it is possible to interpolate between these two regimes by dialing a parameter, for example in holographic QCD with the 't Hooft coupling [30]. It may be that by studying different types of potentials of the baby Skyrme model even more diverse near-BPS behaviors will be discovered.
A crucial point of the perturbative scheme deployed in this paper is that the leading order (LO) energy is finite. In app. A, we have extended the BPS solutions and the leadingorder mass correction to the case of various generalized potentials. For the straightforward generalization of the pion mass term, which is simply the latter to the power s (i.e. V ∼ (1−φ 3 ) s ), the BPS solution yields a finite integral for the kinetic energy for all the values of s studied in the appendix. Changing the potential to the modified pion mass, which possesses domain walls, (i.e. V ∼ (1 − (φ 3 ) 2 ) s , the BPS solutions can be found only indirectly, except for s = 1. But yet worse, for the solution with s = 1, the leading order mass correction diverges. This demonstrates that the problem of perturbation around the compacton is not due to the cusp, that one might naively expect, but is due to the behavior of the BPS solution at the origin (i.e. at r = 0). For the BPS Skyrme model in 3+1 dimensions, a similar issue with the leading-order mass correction is known, for example for the solutions in ref. [1]. As illustrated in the app. A, it may well be that one must find an appropriate class of potential for the BPS Skyrme model, in order to be able to port the perturbative scheme to the 3+1 dimensional model. potential. The BPS equation (2.16) for a generic potential takes the form We will consider a few cases in turn in the next subsections.
Starting with the simplest generalization of the standard pion mass term, we consider the latter to the power s: Finally, we move the singularity to ξ = 0 by setting ξ 0 = 1: The BPS bound is still attained, but now we can calculate the LO contribution to the energy using the above solution with 0 < s < 2. Unfortunately, we do not know how to perform this integral for any value of s in the given interval, but it can be carried out for certain specific values    When N is not close to an integer, it is less clear which N actually provides the minimal energy at LO. Thus we provide the data for F (N )/N in tab. where we have set ξ 0 = 1. This solution, however, does not yield a finite result for the leading-order correction to the mass.