Near-BPS baby Skyrmions with Gaussian tails

We consider the baby Skyrme model in a physically motivated limit of reaching the restricted or BPS baby Skyrme model, which is a model that enjoys area-preserving diffeomorphism invariance. The perturbation consists of the kinetic Dirichlet term with a small coefficient $\epsilon$ as well as the standard pion mass term, with coefficient $\epsilon m_1^2$. The pions remain lighter than the soliton for any $\epsilon$ and therefore the model is physically acceptable, even in the $\epsilon \to 0$ limit. The version of the BPS baby Skyrme model we use has BPS solutions with Gaussian tails. We perform full numerical computations in the $\epsilon\to 0$ limit and even reach the strict $\epsilon=0$ case, finding new nontrivial BPS solutions, for which we do not yet know the analytic form.


Introduction
The Skyrme model [1,2] is an attractive model for a field-theoretic approach to nuclei and it is the low-energy effective theory of the Witten-Sakai-Sugimoto model [3,4]. Its simplest version overestimates the nuclear binding energies at the classical level by about an order of magnitude compared with those observed in experiments. BPS solitons, on the other hand, saturate a Bogomol'nyi bound which implies that such solitons have vanishing binding energies. A BPS Skyrme model was proposed in ref. [5] by Adam-Sanchez-Guillen-Wereszczynski (see also refs. [6][7][8][9]) which is a BPS theory that has BPS solutions for any value of the topological charge or baryon number, B. Depending on the potential used in this BPS theory, the solutions are often of a compacton type, that is, the solutions have support on a compact region of space and has a discontinuous derivative at the border of the compacton region. This property can in part be traced back to the fact that the BPS Skyrme model does not have a kinetic term, but its only derivative term is of sixth order and is the topological charge-density squared.
Real-world nuclei have binding energies at the one-percent level of their rest mass, and it is expected that a small perturbation of a BPS model will yield multi-Skyrmion solutions with low binding energies -at least at the classical level. Deforming the BPS Skyrme model is easy enough: a natural deformation is to try and revert back towards the chiral Lagrangian and include the kinetic (Dirichlet) term with a small coefficient [10], see also ref. [11] for a deformation in holography. Unfortunately, the numerical solutions become unprecedentedly difficult in the tiny limit, whereas for of order one, there are no severe difficulties with the model. This fact is due to several factors. The main reason is that the BPS model has an infinite moduli space, i.e. of volume-preserving diffeomorphisms; the theory describes the Skyrmions as an incompressible fluid [12], akin to the liquiddrop model description of nuclei. Indeed this fact about the model is very welcome for phenomenology, and especially for finite density applications [13][14][15]. However, this means that the = 0 theory has infinitely many solutions, whereas the 0 < 1 theory picks out a (possibly finite) subset of these infinitely many solutions. The mathematical problem of picking out this subset of solutions has been coined restricted harmonicity by Speight [16]. The idea is that the solution to a near-BPS theory deformed by the Dirichlet term, whose BPS solutions can have all shapes with a fixed volume, must be the subset of the solutions that minimize the Dirichlet energy. It is thus a "constrained" harmonicity problem.
Because the near-BPS limit of the BPS Skyrme model is very challenging, we have in ref. [17] considered the analogous problem in the toy model, known as the baby Skyrme model [18][19][20][21]. The baby Skyrme model is "smaller" than the Skyrme model in two senses: it is formulated in 2+1 dimensional spacetime instead of 3+1 and it has the target space S 2 instead of S 3 . A further similarity with the full Skyrme model, is that it also possesses a BPS version, which we shall call the BPS baby Skyrme model or restricted baby Skyrme model. This BPS version of the model consists, again analogously to the BPS Skyrme model in 3+1 dimensions, of the topological charge-density squared as well as a potential. Furthermore, the BPS baby Skyrme model also enjoys area-preserving diffeomorphism invariance for static solutions [22], and it was in fact discovered earlier in this context, than in the full 3+1 dimensional Skyrme model. In our previous paper [17], we chose to stick with the most simplistic potential, namely the pion mass term. This had the consequence that the BPS solutions in that model are of the compacton type [23][24][25]. Other potentials, however, will provide BPS solutions with either exponential (Gaussian) tails or power-law tails [24]. Baby Skyrmions are also interesting in their own right, and in fact they are quite similar to the Skyrmions being studied heavily for the moment in magnetic materials, see ref. [26] for a review, and refs. [27][28][29][30][31] for some recent theoretical work.
The BPS property originally appeared in theories where BPS solitons were solutions preserving a fraction of supersymmetry. For the Skyrme model in 3+1 dimensions, this is not the case simply because its target space is not Kähler [32]. The first attempts at supersymmetrizing the Skyrme model ended up with a model with CP 1 target space [33,34], which is the right target space for the baby Skyrme model (but not the Skyrme model), although the bosonic Lagrangian was slightly different than that of the baby Skyrme model. The exact supersymmetric version of the baby Skyrme model was later constructed with N = 1 supersymmetry [35,36], but it does not contain the bosonic Lagrangian of the restricted baby Skyrme model either. It turns out that the restricted baby Skyrme model possesses N = 2 supersymmetry [36,37] and interestingly, supersymmetry forbids the presence of the kinetic Dirichlet term, which exactly corresponds to the → 0 limit mentioned above. The baby Skyrmion solitons preserve only a quarter of supersymmetry (i.e. one supercharge) [38,39]. The Skyrme model was successfully supersymmetrized in ref. [40] by complexifying the target space from SU (2) to SL (2, C), and like the 2-dimensional case, supersymmetry also forbids the presence of the kinetic Dirichlet term. This supersymmetric Skyrme model does possess soliton solutions, but they are not BPS states [41].
In this paper, we construct a BPS sector which is composed by the Skyrme term (which is also the topological charge-density squared) as well as a potential that does not give the pion a mass. Specifically, we will choose the pion mass term squared, which gives the BPS solitons a Gaussian tail. The deformation of the BPS sector is given by the kinetic (Dirichlet) term, with coefficient , and the normal pion mass term, with coefficient m 2 1 . This has the nice scaling with that keeps the mass of the perturbative pions constant and equal to m 1 . Because of the deformation being in the form of both the kinetic (Dirichlet) term and a potential term (the standard pion mass term), the solution of the deformation is not restricted harmonic, but some generalization which we will call generalized restricted harmonic (GRH). 1 We perform large-scale high-definition brute-force numerical computations for the topological sectors Q = 2 and Q = 4 and are able to dial all the way to = 0 -finding new nontrivial BPS solutions, for which we do not know an analytic expression. We know that for Q = 1, the restricted harmonic solution, i.e. the solution for > 0 is axially symmetric. We study the axially symmetric baby Skyrmions in the perturbative scheme, put forward in our previous paper [17] and find very good agreement between the perturbative results and the full numerical ODE results. The main difference here with respect to ref. [17] is that the BPS solution here is not a compacton and therefore there is no need for imposing special cusp conditions at the boundary of the would-be compacton. The most stable axially symmetric baby Skyrmion has charge Q = N = 2 and therefore it is expected that for large topological charge Q, the most stable solution is made of Q 2 almost axially symmetric N = 2 components -stitched together in some fashion. Throughout the paper, we are using the notation that Q denotes the total charge of the baby Skyrmion, regardless of its shape, and N denotes the charge of an (almost) axially symmetric baby Skyrmion. We then turn to calculating the bound states in the perturbative scheme, using the expansion up to O( 3 ) or next-to-next-toleading order (N 2 LO). To leading order (LO), the correction to the energy is given by inserting the BPS solution into the deformation terms -this cannot explain a bound state, because in order to avoid overcounting, we simply cut off the tails of the BPS and LO energy densities, that otherwise would overlap. This does explain an attractive force at LO. In order to understand the separation distance of two N = 2 (or even two N = 1) constituents, we need also a repulsive force at shorter distances. In our perturbative scheme, we do find such a force and it is due to energy accumulation near the middle of the bound state, caused by our imposed gluing conditions. Unfortunately, this solution based on the perturbative scheme does not match with the precise full numerical calculations, neither for the separation distances nor for the binding energies and we conclude in this paper that the perturbative scheme fails for solitons with tails -as opposed to compactons [17].
The paper is organized as follows. In sec. 2, we introduce our version of the deformed restricted baby Skyrme model, its BPS solutions, topological energy bounds and define a physical length scale for the BPS solutions. In sec. 3, we review and modify the perturbative -expansion scheme introduced in ref. [17]. We then test the scheme on axially symmetric baby Skyrmion solutions. In sec. 4, we perform full high-resolution numerical PDE computations and find nontrivial BPS solutions as well as near-BPS solutions. In sec. 5, we review the long-range interactions of ref. [20] but with the notation of our model. In sec. 6, we attempt at calculating the binding energies perturbatively, but discover that the method fails for solitons with tails. We conclude the paper with a discussion in sec. 7. We have relegated some details of a(n -dependent) BPS solution that contains the potential of the deformation Lagrangian to app. A.

The model
The model is based on the restricted (BPS) baby Skyrme model with a non-BPS deformation in the form of the kinetic term as well as the standard pion mass term, both with coefficient , is a real 3-vector, the kinetic term and the Skyrme terms are given by respectively, and the two potentials are with p = 1, 2. The nonlinear sigma model constraint φ · φ = 1 is enforced by means of the Lagrange multiplier, λ. In this paper we will use the mostly positive metric signature.
We have kept the pion mass, m 1 as a free parameter, but in order to avoid clutter, we have scaled away the coefficients in front of the BPS potential V 2 (φ) and the Skyrme term by an appropriate choice of energy and length units (without loss of generality).
The requirement of finite total energy effectively point compactifies 2-space from R 2 to R 2 ∪{∞} S 2 , i.e. a 2-sphere. It would also induce a spontaneous symmetry breaking from O(3) to O(2), but we break this symmetry explicitly with the two potentials. Therefore the target space is given by O(3)/O(2) S 2 which is also a 2-sphere, hence allowing for topologically nontrivial solutions with degree which is also equal to the number of baby Skyrmions in R 2 . The solutions that minimize the static energy are called stable baby Skyrmions, whereas local minima exists which we shall call metastable baby Skyrmions.
In ref. [17] we have studied the case with only V 1 as the potential, which in the limit → 0 has solutions of the compacton type [24], meaning that the soliton has only support on a compact domain D ⊂ R 2 , which could be of the shape of a disc, but not necessarily. This was the case because the potential V 2 was absent. In fact, any potential with 0 < p < 2 would give rise to a compacton [24]. The potential V 2 gives rise to solutions with a Gaussian tail, even in the limit of → 0. We will see this explicitly in the next subsection.
For the analytic calculations, it will be useful to use the following parametrization for which φ · φ = 1 is manifest. The Lagrangian components in terms of ω can be written as the potentials as with p = 1, 2, and the topological charge as The ω parametrization will prove useful for finding the BPS solutions, which we will turn to next, but we will turn back to the φ parametrization for the numerical analysis.

BPS solutions
The model contains a BPS submodel, which is found by sending to zero ( → 0) in the Lagrangian (2.1): It contains also a different BPS submodel which includes the potential m 2 1 V 1 , but the solutions are obviously dependent on , see app. A. The BPS equation can be found by performing a so-called Bogomol'nyi trick (ignoring the Lagrange multiplier term) The BPS equation is thus 14) and when satisfied, the total energy is proportional to the topological charge or degree of the baby Skyrmion, Q, which is the minimum of the energy in the given topological sector (for given Q).
We will now change parametrization of φ to stereographic coordinates (2.7) and then the BPS equation reads (2.15) Inserting the Ansatz ω = e iN θ ζ(r), we obtain the equation where we have chosen the lower sign. It will now prove convenient to change variables as [24] 1 for which the differential equation becomes which has the solution [24] γ = e −ξ 2 −κ , (2.20) where ξ ≡ r R , the characteristic radius is and κ is an integration constant. Changing variables back to ζ, we finally obtain the solution where in order to ensure that the target space is fully covered, we must make sure that the function is singular at ξ = 0 and hence we have set κ = 0. For large values of ξ, we have It will prove useful to calculate the BPS mass (2.24) In the fourth line we have chosen the lower sign, corresponding to the above-found BPS solution with the boundary conditions ζ(0) → ∞ and ζ(∞) = 0.
The topological charge of such an axially symmetric configuration is thus given by (2.25) We will call the topological charge N only for axially symmetric baby Skyrmions, whereas for generically shaped multi-baby Skyrmions we will denote it by Q.

Energy bound
In the previous section, the energy was shown to be bounded from below by the Bogomol'nyi bound in the BPS limit, = 0, There is another BPS limit, which is the double limit → ∞ and m 1 → 0. This limit is also bounded from below due to the model becoming the pure O(3) sigma model (consisting only of the kinetic term) [36], where the topological solitons are instead called lumps.
The total energy in the model (2.1) in the limit m 1 → 0 therefore obeys the energy bound This bound is a composite of two limits and it is indeed only saturated in those two limits: i.e. in the limit of → 0, the model contains BPS baby Skyrmions and in the limit of → ∞, m 1 → 0 the model contains BPS lumps [36].
Now we consider turning on m 1 , for which the composite energy bound (2.28) no longer can be saturated. It is however possible to calculate a bound similar to that of eq. (2.26), but instead of including only V 2 as the potential we include both the BPS potential and the deformation potential, that ism 2 Writing down the would-be BPS mass which is the bound for the model (2.1) in the limit of → 0 and m 2 1 =m 2 1 fixed. The expression in the parenthesis on the last line of eq. (2.29) has the limiting value of unity whenm 1 → 0 tends to zero; it thus coincides with the bound of eq. (2.26) as it must. The BPS solution corresponding to the above energy bound is given in app. A.
We can now combine the energy bound (2.27) for the sigma model and for the remaining terms in the Lagrangian (2.1) to obtain the energy bound for the full model as The BPS bound in the limit → 0 is clear. The first order correction to the energy bound, in , is due to the mass term and the lump mass. The higher-order corrections to the energy bound come from the would-be BPS bound on the Skyrme term together with both potentials, i.e.m 2 1 V 1 + V 2 .

Length scale
It will prove useful to know the length scale of the baby Skyrmion in the BPS limit, a.k.a. its radius. The short answer is R = √ 2N , which depends only on N (since we have fixed the coefficient of the BPS potential to be unity, which can always be done by a rescaling of lengths).
For a more precise estimate of the radius, we can first calculate at which value of ζ the fraction β of the energy is contained. This in turn translates to a radius and hence yields the coefficient multiplying R for the radius. The entire BPS energy is contained if ζ tends to zero in the integral (2.24), Therefore a fraction β of the BPS energy is contained by Using that ξ ≡ r R , we find the radius that contains the fraction β of the total BPS energy reads The radii corresponding to β = 0.9, 0.95, 0.99 are r 0.9 1.073R, r 0.95 1.224R and r 0.99 1.517R, respectively.

Perturbation in
We will now use the framework of perturbation theory around a background soliton solution, ϕ, developed in ref. [17] where The BPS solution in this model is given by eq. (2.22) in the ω parametrization. Since the numerical calculations will be performed in the φ parametrization, it will be useful to write the charge-N axially symmetric baby Skyrmion solution in that parametrization: with ξ ≡ r/R and R is the characteristic radius (2.21), x + iy = re iθ are the standard polar coordinates in R 2 and α is a phase modulus. The phase modulus is an internal parameter of the solution, but for a single solution is equivalent to a rotation in the plane by −α/N .
Next we will consider the corrections to the energy order-by-order in , in the following sections, starting with the leading order.

Leading-order correction
The leading order (LO) correction is the first and linear order in , and in contradistinction to the compacton case of ref. [17] comes from both the kinetic term as well as the (perturbative) pion mass term. Inserting the BPS solution into these terms thus gives where in the following we will only consider N > 0 positive. Of course, by including the pion mass in the perturbation, the leading-order contribution to the energy depends on the mass parameter m 1 . The parenthesis on the last line takes the value 1.169 for N = 1 and m 1 = 0, which is about 17% above the energy bound for the kinetic term. Turning on the pion mass, m 1 = 0.5, increases the value of the parenthesis to 1.419.
To this order, the energy reads Note that this is strictly above the BPS bound (2.31) considered at the linear level in . Considering now the energy per N as a function of N , we can determine which solution has the lowest energy per baby Skyrmion (nucleon). We have d dN which is the N with the minimum energy per N . We notice that the formal minimum of the energy per N , to leading order in , does not depend on m 1 .
In fig. 1 is shown the leading-order energy correction, divided by N . This is the energy correction per nucleon, to be multiplied by . We can see from the figure that the stable axially symmetric baby Skyrmion, to leading order in , will have N = 2. For example, two separated 1-Skyrmions will have a higher energy than an axially symmetric 2-Skyrmion and an axially symmetric 4-Skyrmion will have higher energy than two separated 2-Skyrmions. Also an axially symmetric 3-Skyrmion will have higher energy than three 1-Skyrmions or a 1-Skyrmion and a 2-Skyrmion.

NLO and N 2 LO corrections
The next corrections are the next-to-leading order (NLO) and next-to-next-to-leading order (N 2 LO) corrections to the energy which are of order 2 and 3 , respectively, and they will be calculated by introducing a linear perturbation around the background field where ϕ = (ϕ 1 , ϕ 2 , ϕ 3 ) is the background solution around which δφ = (δφ 1 , δφ 2 , δφ 3 ) is a small perturbation. We assume that δφ is of order , although this is not clear a priori ; we have to check it a posteriori.
In order to capture the NLO and N 2 LO corrections, we have to calculate the variation of the model (2.1) up to third order in the fields (assuming that δφ = O( )) L perturb [ϕ, δφ] = ∂L ∂φ a δφ a + 1 2 where the vertical bar " | " denotes evaluation on the background by setting φ = ϕ (to the left of the bar) and we have defined the following quantities 11) and the Lagrange multiplier for the background fields, λ 0 , is given by The Lagrange multiplier λ has been replaced by the expansion λ 0 + δλ, where λ 0 enforces the sigma model constraint for the background solution ϕ and δλ makes sure the total field φ = ϕ + δφ remains inside the O(4) group, i.e. it preserves the constraint φ · φ = 1 + O( 4 ). Varying δλ yields which has the solution [20,17] This form of the perturbation satisfies the constraint (3.13) up to O(∆ 4 ) and since ∆ will turn out to be of order , that makes it of order O( 4 ). Notice that the perturbative contribution to the energy from the terms multiplying δλ is exactly the bracket in eq. (3.13) and hence is of order O( 4 ), which we can safely disregard. Therefore, we do not need to find the explicit expression for δλ to this order in perturbation theory; its job was to produce eq. (3.14).
We will now substitute the form of the variation (3.14) into the perturbation Lagrangian (3.8): where the NLO Lagrangian is given by which is the complete perturbation Lagrangian to third order in .
In principle, we should solve this nonlinear problem for ∆, which however is almost as difficult as the original problem, without introducing the perturbation theory on top of the soliton background. Therefore, we will linearize the above Lagrangian, i.e. we will only use the linear and quadratic parts in ∆ (i.e. eqs. (3.17) and (3.18)) to determine its equation of motion.
One may then wonder why go to the third order in if it results in cubic terms in ∆, that we anyway will discard once we turn to solving the equation of motion for the perturbation. The answer is that we need the last term of eq. (3.18), as it gives a term ∂ 2 in the equation of motion, whereas the last term in eq. (3.17) will give a term f (ϕ, ∂ µ ϕ)∂ 2 in the equation of motion for ∆. Although we are preparing for the study of the limit of very small (say about 1/100), it is still much larger than the exponentially suppressed operator coming from eq. (3.17). Therefore, neglecting it will yield an incorrect behavior at distances of the order of the Skyrmion size away from the Skyrmion (dependent on of course).
The linearized static equation of motion for ∆ can thus be written down where we have defined the quantities It will now prove instructive to consider the equation of motion for the fluctuation ∆ at an asymptotic distance from the BPS baby Skyrmion background. In that case, the background solution (3.2) can be approximated as Inserting this approximation into eq. (3.20), we can write the equation of motion for the fluctuation, to second order in e − ξ 2 2 as: where we have defined Clearly, only the first line of eq. (3.28) is not exponentially suppressed and therefore guarantees the propagation of the fluctuation ∆ in the asymptotic regime (away from the background baby Skyrmion). Now, importantly, all the terms with a factor of on the left-hand side of eq. (3.28) come from the O( 3 ) Lagrangian (3.18) and therefore had we only gone to second order in the expansion, the linearized equation for the fluctuation would have looked like this: where we have kept only the terms with the highest powers of ξ on both sides. To this order, the equation of motion for the fluctuation is badly behaved, because it schematically takes the form whose right-hand side diverges exponentially. We expect of course that the fluctuations go to zero at asymptotic distances, and this will also be the case for the third-order in equation of motion (3.28), because of the terms (∂ 2 − m 2 1 )∆ a , which have the expected exponential falloff compatible with the boundary conditions.
We have assumed that the perturbation field is proportional to . In order to attempt to address this point, let us consider (for simplicity) a few of the terms of the third-order in equation of motion (3.28) for the fluctuation: Since w and w are of order one, we can estimate at which distance from the background baby Skyrmion the third order terms become dominant: The smaller values of , the larger distances are needed before the third-order terms take over. Since the exponential of −r 2 quickly becomes infinitesimally small, the tail of the perturbation at asymptotic distances will take the form where c 1,2 are constants. Notice that this is seemingly independent of . All thedependence is contained at distances for which the -dependence becomes quite complicated. Indirectly, the dependence on in eq. (3.34) is possessed by c 1,2 by gluing it together with the solution at distances given by eq. (3.35). Importantly, in the limit of = 0, c 1,2 = 0 because there is no tail correction in the BPS limit. To leading order, the coefficients must behave like c 1,2 ∝ p , with p ∈ Z >0 a positive integer. We have, however, not been able to prove rigorously that p = 1 (i.e. it could be larger than one).
It is worthwhile to compare the situation of the setting at hand -a baby Skyrmion having a Gaussian tail in the BPS limit and where the perturbation of the model entails both the Dirichlet kinetic term as well as the standard pion mass term, see eq. (2.1) -with the compacton case studied in ref. [17].
In the compacton case, as the name suggests, the BPS limit of the baby Skyrmion is a compacton (hence no tail in the BPS limit) and the perturbation in that case was done solely by adding the Dirichlet kinetic term (with coefficient ). Let us recapitulate the situation in the compacton case of ref. [17]. The tail of the perturbation in that case became non-dynamic unless we went to the third order in (which would prevent the compactons of knowing of one another and hence prevent the calculation of binding energies). Going to the third order in the expansion gave a dynamic tail to the perturbation, which however was nonanalytic in . The nonanalycity was due to the fact that the (inverse) propagator had the form ∂ 2 − m 2 and the tail thus took the form e − mr √ , which when Taylor expanded in vanishes at any finite order.
In this case of a baby Skyrmion with a Gaussian tail in the BPS limit, the situation draws some similarities to the compacton case, but nevertheless is dramatically different. The similar property of the expansion, is that we still have to go to the third order in , because otherwise -as demonstrated above -the equations of motion for the perturbation become ill-defined at asymptotic distances or at least incompatible with suitable boundary conditions. A huge difference is that the baby Skyrmions, in the BPS limit, themselves have Gaussian tails, and therefore already "feel" each other when two or more of them are placed at a finite distance from each other. The BPS solution is thus nontrivial for any separation distance. If we now turn on a small but finite according to the Lagrangian at hand (2.1), the perturbation field ∆ adds an exponential tail to the background solution. The field configuration now flows to the nearest GRH solution in field space. The solution to the generalized restricted harmonicity is, however, extremely difficult compared to the compacton case, which is nearly trivial. That is, axially symmetric baby Skyrmions placed at distances such that their compacton regions do not overlap -for details, see ref. [17]. Now, since the background fields (meaning the fields of the background BPS baby Skyrmion) tend to the vacuum exponentially (or rather like the Gaussian), the governing equations of motion for the perturbation are simply and thus the tail of the perturbation (not the BPS solution) is exponential too, see eq. (3.34), but in stark contradistinction to the compacton case, it is independent of (it is nevertheless "glued together" with a solution closer to the baby Skyrmion that is dependent on and hence a suppression of the tail is still in effect, but we shall check this a posteriori ). Of course, this is not an accident, but a consequence of the construction of this more elaborate (and aimed at being more physical) model, compared with the compacton case of ref. [17]. We include the Dirichlet kinetic term with a coefficient , because we want small binding energies (that is a small perturbation added to the BPS model). Then we include the mass term, also with coefficient , in order to prevent that the pion mass becomes unrealistically large in the small-limit. Another feature -also by constructionis that the potential in the BPS sector (V 2 ) does not give rise to a pion mass and therefore the coefficient of the potential is not restricted by the value of the pion mass.
The picture that forms, which we will elaborate on in a later section, is that at asymptotically large distances, the governing equation for the perturbation is eq. (3.36) and the tails of two or more baby Skyrmions attract each other (in the attractive channels) with a force that is at least cubic in , but the mass of the fluctuations is independent of . The pion mass can thus be set to any realistic value 2 .
A word of caution is that we calculate the perturbation using the linearized equation of motion that contains some (crucial) terms (3.18) at O( 3 ), but not the cubic terms (3.19) at the same order. As for the calculation of the energy, we will compare the energy calculated at NLO and N 2 LO to the exact numerical solutions in the axially symmetric case in the next section.
The static energy of the perturbation is simply given by (3.37) A comment in store about the perturbation, ∆, is that the entire static energy vanishes for ∆ ∝ ϕ. However, since the expression (3.14) for δφ is nonlinear in ∆, a few cross terms survive if we take ∆ = δcϕ + ∆ ⊥ , with ∆ ⊥ · ϕ = 0 a perpendicular perturbation. The cross terms only give rise to a linear first-order PDE for δc, which must vanish once subject to the boundary conditions lim |x|→∞ δc = 0 and δc(0) = 0, where the origin is at each background baby Skyrmion center. This justifies setting ∆ = ∆ ⊥ . For more details, see ref. [17].
It will now prove useful to specialize to the case of the background BPS baby Skyrmion solution for ϕ, using transverse perturbations for ∆ = ∆ ⊥ and switching to polar coordinates in the plane, for which the static perturbation energy reads for the NLO terms, 40) for the N 2 LO terms quadratic in ∆ ⊥ and for the N 2 LO terms cubic in ∆ ⊥ , and we have defined the transverse perturbations ∆ ⊥ as 42) and the background BPS baby Skyrmion solution is given by The corresponding static equations of motion read where we have defined the matrices (3.47) as well as the functions Firstly, the source term, i.e. the right-hand side of eq. (3.44) only exists for the upper equation, that is for the equation of motion for δf . Secondly, the mixing between the upper and the lower equations only appears in terms involving a θ derivative and a mixed r and θ derivative. Thus, δθ is a homogeneous source-free equation of motion, unless δf has nontrivial θ dependence. Therefore, if we restrict to axially symmetric background BPS solutions and turn on only axially symmetric perturbations, δf = δf (r), then δθ decouples and is trivially satisfied (δθ = 0 everywhere).
The equation of motion for the perturbation of axially symmetric baby Skyrmions thus reduces to and the corresponding static energy for the perturbation is for the NLO terms, for the N 2 LO terms quadratic in δf and for the N 2 LO terms cubic in δf .

Axially symmetric solutions
In order to verify the accuracy of our perturbative scheme, we start with axially symmetric baby Skyrmions. As was shown in ref. [17], the perturbation can in this case be written as  Notice that the perturbation qualitatively captures the correct behavior of the exact solution everywhere. For r 1.5R the perturbative correction to the BPS solution matches the exact solution extremely well, see the right-hand side panels of fig. 2 (the characteristic radius R = √ 2N is marked with a vertical dashed black line on each graph) and this is the radius within which 99% of the energy of the BPS baby Skyrmion is contained. However, for r 2R the perturbative correction overshoots (undershoots) the exact solution for N = 1, 2 (N = 4), which indicates the level of precision of the perturbative method at the linear order in the perturbation field (corresponding to second and some third order terms in ).
One could think that it would be necessary to include all terms at the third order in the expansion, hence making the problem a nonlinear one. However, as the discrepancy between the perturbative solution and the exact solution appears only at distances r 2R, the exponential suppression of the tail of the BPS background solution makes the perturbative contribution to the energy subdominant in that (asymptotic) region; in fact, less than 1% of the energy density of the BPS background baby Skyrmion resides beyond the distance of ∼ 1.5R from the center of the soliton. In order to demonstrate that this statement is true, we now turn to calculating the perturbative corrections to the energy of the axially symmetric baby Skyrmions. The red line shows the BPS bound, the orange line is the leading-order (LO) correction, the black pluses are the next-to-leading order (NLO) corrections in calculated using the linear perturbation, the green dashed line is a fit to the NLO points, the red crosses are the next-to-next-to-leading order (N 2 LO) corrections in calculated using the same linear perturbation, the dark-blue dashed line is a fit to the N 2 LO points and the black line is the exact energy calculated using the full (nonlinear) equations of motion. The left columns show the total energy, E, the middle and right columns show the energy with the BPS and the LO correction subtracted off for m 1 = 0.5 and m 1 = 1, respectively, to better see the differences between the exact, NLO and N 2 LO corrections.
The perturbative corrections to the energy of the N = 1, 2, 4 axially symmetric baby Skyrmions for m 1 = 0.5, 1 are shown in fig. 3 with the left-hand side panels showing the exact energies (solid black lines), the BPS energies (solid red lines), the NLO corrections (black pluses), the N 2 LO corrections (red crosses). Both the NLO (i.e. O( 2 )) corrections (3.54) and the N 2 LO (i.e. O( 3 )) corrections (3.55)+(3.56) are calculated using the perturbation, which is a solution to the linearized equation of motion (3.53), which contains all terms of the second-order-in-Lagrangian (3.17) and the quadratic terms of the third-order-in-Lagrangian (3.18). The green dashed line shows a fit to the NLO corrections and the dark-blue dashed line for the N 2 LO corrections. The middle and right panels of the figure show the same information, for m 1 = 0.5 and m 1 = 1, respectively, but with the BPS and LO corrections subtracted off, so as to better see the differences between the NLO corrections, the N 2 LO corrections and the exact energies.
First of all, we can see from the total energies, shown in the left-hand side panels of fig. 3, that the precision of the perturbative scheme is extremely good. More precisely, the LO correction overshoots the exact result slightly, but the NLO correction is negative and undershoots the exact result. The N 2 LO correction is then positive and comes extremely close to the exact energies. In fact, by a close inspection of the figures, it can be seen that the red crosses, corresponding to the perturbative N 2 LO corrections, fit the exact energies (black solid lines) better than the N 2 LO fits (dark-blue dashed lines). Although negligibly small, we expect the remaining discrepancy between the N 2 LO corrected energies and the exact energies to be due to the approximation of linearizing the equation of motion for the perturbation and for truncating the expansion at the third order.
The fits shown in fig. 3 contain no linear term in for the NLO corrections and no linear and quadratic terms in for the N 2 LO corrections. Nevertheless they describe the perturbative calculations very well and this is thus the a posteriori confirmation that δφ is of order (and possibly containing higher-order corrections too). It should be fairly convincing that the LO corrections, which are linear in are orders of magnitude larger than the NLO and N 2 LO corrections, that thus do not contain a linear term.
Because we now have a dependence on m 1 in this model, compared with the compacton case of ref. [17], we have calculated the perturbative corrections for two different values of m 1 , see fig. 3. Since the physical case has a nonvanishing pion mass, we have chosen two nonvanishing values, namely m 1 = 0.5 and m 1 = 1.
We will now use the fits to the NLO and N 2 LO results in fig. 3 to write an approximate expression for the N 2 LO energy: where we have assumed that the dependence on m 1 is quadratic. This fit should work reasonably well, at least for m 1 ∈ [0.5, 1]. The reason for including a fourth-order term in in the fit is to allow for some higher-order (residual) behavior of the perturbation. The higher-order result for the energies in turn gives a correction to N of eq. (3.6), i.e.:   The perturbatively corrected N is shown in fig. 4(a) as a function of for m 1 = 0.5, 1. It is seen from the figure that for sufficiently small (i.e. < 0.1) the perturbative corrections to N are too small to decisively make another topological charge sector the stable one. In order to check this explicitly, we plot in fig. 4(b) the ratio of the energy per topological charge for N = 1, 3 to the same quantity with N = 2; that is E( ,N,m 1 ) N 2 E( ,2,m 1 ) , for m 1 = 0.5, 1 as functions of . We can see that there is only a mild dependency on m 1 and all values in the plot are greater than unity: Hence, the N = 2 axially symmetric baby Skyrmion remains the stable solution. Increasing in the range shown in the figure actually makes the 2-Skyrmion more stable, as can be seen from fig. 4(b).

Numerical calculations
We will now present numerical solutions to the full nonlinear equations of motion. It is computationally a very expensive task and because the static baby Skyrmions require only 2-dimensional PDEs, we are able to use refined enough square grids for the computations. The numerical method used is the so-called arrested Newton flow on grids with about 4096 2 lattice points and lattice spacing down to about ∼ 0.0038. The numerical derivatives are approximated using a finite-difference method on a fourth-order 5-point stencil. The numerical accuracy is about 10 −6 or better. The code is an adaptation of the CUDA C code used in ref. [17], which is executed on a GPU cluster.
In this section, we explore numerical solutions in the topological sector Q = 4 and Q = 2 and include the axially symmetric solutions, which can be directly compared to the very precise (exact) numerical ODE calculations. This gives another estimate on our numerical precision. It also confirms in which part of the parameter space, the axially symmetric solutions exist and which type of solution has the lowest energy in the given topological sector.
Our first and main aim of this numerical exercise, is to confirm that the lowest energy baby Skyrmion solution -for small but finite 0 < 1 -is made of two N = 2 baby Skyrmions sitting at some unknown separation distance from one another. This is supported by the LO energy correction, which in the axially symmetric case is confirmed to be the dominant correction to the energy in the small-limit.
We begin with two N = 2 baby Skyrmions placed side-by-side, which thus is in the topological charge sector Q = 4. The numerical results are shown in fig. 5 for representative values of . This figure and the following four are organized into four columns displaying the topological charge density, the baby Skyrmion orientation in O(3) space, the total energy density and finally the perturbation part of the energy density (−L 2 + m 2 1 V 1 ). The figures in columns 1, 3 and 4 are shown with a color scheme normalized to the content on each graph. The figures in column 2 show the pion vector of the baby Skyrmion in the following sense: φ 3 = 1 corresponds to the vacuum and is shown with white, φ 3 = −1 is the antivacuum and is shown with black, arg(φ 1 + iφ 2 ) = 0 is shown with red, arg(φ 1 + iφ 2 ) = 2π 3 is shown with green and finally arg(φ 1 + iφ 2 ) = 4π 3 is shown with blue. The two N = 2 baby Skyrmions sitting side-by-side exist for all values of (in our scanned range), even = 1 and = 0. Starting with = 1, the solution is made of two ring-like baby Skyrmions, that however are prolonged along the axis that joins them. Upon lowering , the greatest effect appears to be on the distribution of the energy and topological charge within each (approximately) axially symmetric N = 2 baby Skyrmion. The separation distance between them is affected only mildly. More precisely, the separation distance between the two N = 2 Skyrmions is seen to increase slightly for decreasing from = 1 to ∼ 0.379 and below that, the separation distance appears to be constant in the limit of → 0. This situation is quite different from the compacton case of ref. [17], where the limit of → 0 yields two perfectly axially symmetric compactons touching each other at a mathematical point with the pion vectors anti-aligned, so the point they touch has the same pion orientation on both sides of said point. In this case, it appears that the limit of → 0 yields a nontrivial GRH solution, where neither of the two N = 2 Skyrmions are axially symmetric. This is obviously because of their Gaussian tails in the BPS limit (i.e. in the limit of → 0) that makes it impossible to place two perfectly axially symmetric BPS solutions next to each other at any finite distance. Although we are able to lower all the way down to = 0 with extreme resolution and in the same time taking into account the tails up to radii r ∼ 15, the two N = 2 baby Skyrmions do not tend to perfectly axially symmetric solutions, see that last row of fig. 5 (i.e. for = 0). The BPS solution shown in the last row of the figure is evidently the BPS solution closest to the limiting sequence of solutions for 0 < 1, but we should point out that the moduli space for = 0 is infinite dimensional and hence this is just one solution out of infinitely many BPS solutions that can take any shape as long as the volume is preserved. 3 A slight polarization of the constituent N = 2 Skyrmions -stretching the baby Skyrmions along the axis that joins them -persists in the limit of → 0, see fig. 5. This configuration in the topological charge sector Q = 4, consisting of two N = 2 Skyrmions is the stable solution (i.e. with the smallest energy per Q) for all values of studied in this paper.
Apart from the global minimum of the energy functional, which for Q = 4 sector Figure 6. The Q = N = 4 baby Skyrmion for various values of . The columns show the topological charge density, the pion vector orientation using a color scheme described in the text, the total energy density and the perturbation energy density (−L 2 + m 2 1 V 1 ). The last row (i.e. for = 0.379) shows that the solution has decayed into a lower-energy configuration with almosttriangular symmetry instead of axial symmetry. In this figure m 1 = 0.5. we claim is attained by two N = 2 baby Skyrmions with a small separation distance, there exists also local minima -viz. metastable states. For this reason, we have tried with several good guesses as initial conditions and observed what they flow to under the numerical minimization of the energy by means of the arrested Newton flow algorithm.
As N ∼ 1.5 and the energy contains a term that grows quadratically with N , it is expected that the N = 4 axially symmetric baby Skyrmion is just a metastable state, if it exists at all. An explicit check shows that the N = 4 solution exists for = 1, but it  decays to a(n almost) triangularly symmetric arrangement of four N = 1 baby Skyrmions for ≈ 0.4, see fig. 6. This is analogous to what happens in the compacton case [17]. This solution is also a metastable state and in particular for large ∼ 1, two N = 1 baby Skyrmions are unstable and quickly combine into an axially symmetric N = 2 solution, see below. However, the almost triangularly symmetric configuration enjoys metastability due to the fact that the center N = 1 Skyrmion is confused by 3 attractive channels and cannot decide which other N = 1 to combine into an N = 2 with.
We increase for the (almost) triangularly symmetric configuration of four N = 1 Skyrmions from = 0.379 all the way up to = 1, see fig. 7. The arrangement, although metastable, turns out to remain, even for = 1 and in fact it becomes perfectly triangularly symmetric for the larger values of due to larger attraction between the four N = 1 constituents. We also decrease all the way down to = 0 and find that it becomes less triangularly symmetric for the smaller values of -although we do not know exactly why. The attractive force between the N = 1's become very weak and we can see from the second-last row of fig. 7 that it would be fairly easy to knock off the N = 1 baby Skyrmion on the right-hand side of the solution. In that sense, the model will be able to describe weakly bound nuclear clusters. Of course, in two dimensions it makes little sense to try and compare with actual nuclear physics knowledge. In the last row of the figure, where = 0, the position of the constituents are just moduli and can be moved freely, of course.
We could in principle search for further arrangements of four N = 1 baby Skyrmions, but we already know from the leading order calculation of the energy in , that the energetically preferred constituents are the axially symmetric N = 2 baby Skyrmions. We will leave such a search for more clusters to future work as they will not be relevant in the small-limit.
We will now determine the phase diagram in the Q = 4 topological charge sector. To this end, we plot the energies of the baby Skyrmions displayed in figs. 5, 6 and 7 in fig. 10 (with exception of the = 0 ones, due to the logarithmic scale of the ordinate). Starting with the least stable solution, the axially symmetric N = 4 baby Skyrmion (magenta pluses) only exists for 0.4 and then decays to the (almost) triangularly symmetric solution (blue dotted squares), which in turn is only slightly higher in energy compared with the lowest-energy state -the two N = 2 baby Skyrmions side-by-side. Indeed it is the globally stable solution in the entire range of considered here. We also note that the LO energy is actually an extremely good approximation to that of the solutions made of two N = 2 Skyrmions or four N = 1 Skyrmions for 0.1.
We now turn to the Q = 2 topological charge sector, for which the combinatoric possibilities obviously are more limited. Of course, from sections 3.1 and 3.3 we know that the N = 1 baby Skyrmion has a larger mass per N than the N = 2 one, both of them with axial symmetry. Nevertheless for a full understanding of the phase diagram in the Q = 2 sector, we need to study where the two N = 1 baby Skyrmions side-by-side are metastable and where they are unstable to collapse into an axially symmetric N = 2 baby Skyrmion, see fig. 8. With the intuition from the compacton case of ref. [17], we expect the two N = 1's side-by-side to exist only for small values of . The critical value of above which the two N = 1's side-by-side are unstable is roughly crit ≈ 0.15. In the small-limit, the two N = 1 baby Skyrmions side-by-side appear to be two nearly axially symmetric baby Skyrmions vaguely connected in the attractive channel. Because of the Gaussian tail in the BPS limit, the two N = 1 baby Skyrmions cannot recover axial symmetry in the → 0 limit -in sharp contradistinction to the compacton case of ref. [17].
We do not need to solve the full PDEs to obtain the energies for the stable (in the Q = 2 sector) axially symmetric N = 2 baby Skyrmions. Nevertheless, we perform these calculations as a check of our numerical accuracy. The solutions are shown in fig. 9.
Finally we plot the energies of the two types of solution in the topological charge Q = 2 sector in fig. 11. In the left-hand side panel of the figure, we can see that both types of solution are very close in energy and are actually described quite well by the LO correction to the energy for 0.1. In the right-hand side panel, we show the same energies but with the BPS and LO corrections subtracted off. First of all, we can see that the full numerical PDE solutions of the N = 2 baby Skyrmions (magenta pluses) match extremely well with the exact energies from the ODE calculations (black solid line). Second of all, we confirm that the two N = 1 baby Skyrmions sitting side-by-side always have a slightly higher energy than the N = 2 one. The energy difference, however, is very small and hence the metastability is possible even with a small energy barrier between the two kinds of solution.  Table 1. Energies of BPS solutions (i.e. with = 0) compared with the theoretical BPS mass for given Q.
In sharp contrast to the compacton case, we have in this model been able to send all the way to zero, obtaining numerical BPS solutions, 3 of which do not possess axial symmetry, see the last row of figs. 5, 7, 8 and 9. Since = 0 cannot be plotted on the logarithmic ordinate of figs. 10 and 11, and they are all of the same value (the BPS mass), we present the numerically calculated energies in tab. 1 as a handle on the numerical accuracy of our solutions. For the Q = 4 sector, the accuracy (discrepancy) of the numerically calculated energy is about 6 × 10 −5 , whereas for the Q = 2 sector it is 1 × 10 −5 and 1 × 10 −4 , for the two solutions.  We evaluate the binding energy of the Q = 2 + 2 bound state solution of fig. 5 in fig. 12(a) and of the Q = 1 + 1 bound state solution of fig. 8 in fig. 12(b)  and these binding energies are fitted for solutions with m 1 = 0.5. The binding energies are tiny and the slight oscillation of the numerically evaluated binding energies around the fit shows that we are close to the edge of our numerical precision. In principle, we could also numerically evaluate the energy of the bond between two axially symmetric N = 2 Skyrmions, that is 2E 2 − E 4 , but this quantity turns out to be at the level of 10 −4 and beyond our numerical precision, hence we do not attempt at plotting it.
We notice that there seems to be a linear contribution to the binding energy, in contradistinction with the compacton case of ref. [17]. We will discuss this point in later sections.

Long-range interaction
In the case of the compactons, discussed in ref. [17], the restricted harmonic condition was not sufficient to fix the right background solution in the case of two interacting baby Skyrmions. In particular, no restrictions in the choice of the relative orientation or separation distance between the non-overlapping compactons emerge from such a condition.
In this paper, we have turned on a more physical potential and now have long-range forces. This means that at any separation distance, two baby Skyrmions feel each other and the force between them is dependent on the relative orientation, as was shown by Piette-Schroers-Zakrzewski [20]. Since we are only interested in baby Skyrmion bound states, among all the possible GRH maps we must choose among those that are in the attractive channel. This uniquely fixes the relative orientation and hence the choice of the GRH map. In the following, we briefly review the calculation of ref. [20], adapting it to the case of our Lagrangian (2.1).
Let u and v be baby Skyrmion solutions of charge N and M and let φ u and φ v be their representation in 3-vector coordinates. Each field φ is a solution of the full equations of motion: where n ≡ (0, 0, 1) is the vacuum of φ and When the two baby Skyrmions are well separated, the composite solution of total charge Q = N + M can be written as φ w ≡ φ u+v . In order to evaluate the interaction potential V between the two solitons, it is necessary to decompose the total energy in the form To this end, it is useful to separate the coordinate space R 2 into two regions such that φ u ≈ n is close to the vacuum n in region 2 and φ v ≈ n is close to the vacuum n in region 1. We can thus decompose the field into the vacuum n and an infinitesimal correction: where δφ u · n = 0. Since φ u solves the Euler-Lagrange equation (5.1), δφ u satisfies the linearized equation (∆ − m 2 1 ) δφ u = 0, (5.4) and we have that in region 1: where v is linear in δφ v . The expansion for the field φ w in terms of δφ u in the region 2 is analogous. The asymptotic result (5.4) is independent of as we expected and equivalent to that of ref. [20]. Once we have the form of φ w in the regions 1 and 2 (see eq. (5.5)), it is possible to evaluate the energy E[φ w ] as where u is linear in δφ u and is defined by and v is equivalent to eq. (5.7) with the superscripts u and v exchanged. In the end, using the equations of motion (5.1) and Gauss's law, we can finally write the long-range interacting potential V as where Γ is a curve without self-intersections, separating region 1 from region 2 and dS i = ε ijγj dt for any parametrization γ(t) of Γ. Given eq. (5.4), the form of this potential is equivalent to the one obtained in ref. [20] albeit with an overall factor of .
With the result (5.8), it is now possible to calculate the long-range interaction for different multisoliton configurations to determine the Skyrmions' relative orientation maximizing their attraction. For the multisoliton case of charge Q = 1 + 1, the potential takes the form where d is the relative distance, α and β are the two respective phases (orientations) of the baby Skyrmions and we have assumed that δφ ∼ O( ), see the discussion in sec. 3.2. From this expression, we recognize that the maximally attractive channel is obtained for α − β = π, i.e. when the two baby Skyrmions have opposite orientation. Following this indication, we choose the GRH map of this topological sector with the solitons oppositely oriented. In the case of the baby Skyrmion pair of charge Q = 2 + 2, the asymptotic potential is In the previous section, we have seen that two N -Skyrmions in the attractive channel attract each other at asymptotic distances. This attraction presumably breaks down due to nonlinearities at some finite distance and for > 0 there will be a bound state with two N -Skyrmions separated by a finite distance 2a, see fig. 13.
In order to calculate the binding energies using the perturbative method of theexpansion that we have put forward in ref. [17] and this paper, we place two N -Skyrmions side-by-side in the attractive channel, viz. with the pion field matching on the gluing line (x 1 = 0), see fig. 13. In order to be able to glue the two baby Skyrmions together, we need appropriate boundary conditions at x 1 = 0 (see fig. 13), which we shall call gluing conditions: The reason for the odd condition on φ 2 , is that even boundary conditions on all φ a will turn the mirror Skyrmion into an anti-Skyrmion for odd N , which is not what we want to glue the solution with. The relatively simple-looking conditions above, become nonlinear Robin-type boundary conditions for the fluctuation field ∆, once we use that φ = ϕ + δφ and δφ from eq. (3.14) : 3) It turns out to be a rather tricky boundary condition to implement numerically. For this reason we will use the Ansatz (3.42) for ∆ which reduces the above nonlinear Robin gluing conditions to We will now use a finite difference approximation for the x-derivative 4 of δf and δθ: inserting these into the algebraic conditions and linearizing with respect to δδf 0 and δδθ 0 . This yields expressions for δδf 0 and δδθ 0 which are not particularly illuminating, so we will not display them here. In principle the method is simple, δf 0 is the previous value of δf 0 and δf 0 is updated by adding the solution for δδf 0 to δf 0 at each step of the algorithm (and similarly for δδθ 0 ). Notice that δδf 0 vanishes when δf 0 is a solution to the full nonlinear Robin type gluing condition. Unfortunately, this simplest form of Newton iteration does not have particularly good convergence properties [44]. For this reason we had to implement a line search algorithm following ref. [44], which updates δf 0 and δθ 0 using δf 0 = δf 0 + tδδf 0 , δθ 0 = δθ 0 + tδδθ 0 , t ∈ (0, 2), (6.8) where t is a parameter that should be optimized for each step in the iteration. Note that t = 1 is just the simplest version of Newton iteration. We implement a rather crude line search algorithm that tries out 20 points of t in its interval and refines the search one time.
It turns out that often the optimal value of t is around t ∼ 0.5 or t ∼ 0.6. Although this algorithm is rather computationally expensive, ref. [44] showed the line search is one of the cheapest algorithms with improved convergence properties for this type of algebraic Riccati equation.
We now turn to solving the coupled PDEs (3.20) with the boundary conditions for all x 0 being centers of N -Skyrmions, as well as the gluing conditions described above.
The "boundary condition" at all x 0 is to prevent the perturbation from trying to unwrap the soliton solution. Since the calculation is computationally very expensive, we have implemented the code in CUDA C and run it on an NVIDIA GPU cluster.
In order to calculate the binding energy of two charge-N baby Skyrmions, we first need to determine their optimal separation distance; that is the distance where the total energy of the bound state is at a local or global minimum and then read off the energy at that point.
In an attempt to understand which separation distance (2a) is preferred by two N = 2 baby Skyrmions in the attractive channel, we calculate the energy density as a function of a for various values of . The result is shown for = 0.01, 0.0207, 0.0428, 0.0886 in fig. 14. Let us first see what happens at each order in our expansion. At the leading order (LO), there is no fluctuation field and the BPS background solution ϕ is used to calculate the energy from the Lagrangian (2.1), including the BPS energy. Of course, we have stitched together two baby Skyrmions side-by-side, so the LO energy is evaluated without overlap: this means we integrate the energy over the background BPS solution up to the gluing line (x 1 = 0) and then multiply by two. Clearly, this order is flawed in the following sense: The further we push the two baby Skyrmions together, the less the LO energy is. This cannot be the correct physical picture. For that, the perturbation field ∆ is necessary. In particular, the gluing conditions are crucial as should become clear momentarily. For completeness, we show both the energy calculated to order 2 (NLO) and to order 3 (N 2 LO). As we can see from the figure, what happens is that the LO energy goes down as a decreases, and so does the NLO energy, but the N 2 LO energy goes up. This is simply the gluing condition twitching the field at the gluing line and building up energy localized around the gluing condition. However, neither the N 2 LO nor the NLO semi-analytic energies capture a minimum at the separation distance a ∼ 0.8 that was found in sec. 4, and we conclude that the perturbative method has failed to calculate the bound state for the baby-Skyrmions with (Gaussian) tails, as opposed to the case of the compactons, for which the method  was very successful [17]. We can also see the minimum (which only appears for the N 2 LO energies) is extremely shallow and orders of magnitude smaller than expected.
Before, we discuss the sources of the failure of the method to capture the bound states, we consider the bound state of two N = 1 baby Skyrmions, which are only metastable (see fig. 1), since the axially symmetric N = 2 solution has lower energy than the bound state of two N = 1 Skyrmions. The bound state does, nevertheless, exist for 0.15, see fig. 8.  In order to attempt at calculating the binding energy of the bound state in this case, we again perform a large number of PDE calculations in this semianalytic approach for many values of the separation distance 2a, and = 0.01, 0.0207, 0.0428, 0.0886. The result is shown in fig. 15. Now one slight improvement over the previous case, is that we find a minimum of the energy at order N 2 LO for all values of smaller than ∼ 0.08. 5 The minimum, however, again appears at quite large values of a and the perturbative method has thus failed also in the case of two N = 1 baby Skyrmions. To summarize, we show in fig. 16 a comparison of the separation distances found by the perturbative approach (green crosses) and those found by precise full numerical calculations (red crosses) of sec. 4. We also display the would-be compacton radius of a compacton that corresponds to taking the limit → 0 with √ m 1 =:m 1 held fixed.
What this limit amounts to, is to kill off the kinetic term (L 2 ) while keeping the pion mass term V 1 . In this limit, there is a BPS solution with compacton radius R compacton given by eq. (A.10). None of the approaches agree and the only precise method used for calculating the bound state is the numerical method used in sec. 4.
The reasons why the semi-analytic or perturbative approach fails for baby Skyrmions with tails, as compared to compactons, are: • The two axially symmetric BPS solutions glued together, used as a background solution for the bound state calculation, is not a BPS solution (although it was in the compacton case [17]).
• The axially symmetric BPS solution is no longer close to the exact solution at long distances (see fig. 2), as was the case for the compactons [17], which in turn has implications for using it for the bound state calculation by gluing two of them together.
• The nonlinearities become extremely important for the bound states, because of the above reason.
• In the case of the compactons, the cusp condition gave a crucial and exact condition for the perturbation field at the boundary of the compacton, outside of which the calculations fails to capture it in fig. 15(d), but not by much.
perturbation field was governed by a free theory [17]. Without this crucial nonperturbative "bootstrap", the perturbative method becomes imprecise and fails to capture bound states.
• Without the above-mentioned "bootstrap", the equations of motions to be solved must include the nonlinearities, making the problem as difficult as the original numerical problem -the perturbative approach thus fails in this case as it does not simplify the problem to be solved.
• Finally, the gluing condition in the case of the compactons worked well because the perturbation fields were very small near the gluing boundary [17]; this is not the case here for the baby Skyrmions with (Gaussian) tails, in particular near the optimal separation distance, and hence the gluing becomes unreliable without taking into account nonlinearities.

Conclusion and discussion
In this paper, we have studied the near-BPS regime of the baby Skyrme model with a physical pion mass (i.e. one that does not diverge in the BPS limit) and a potential that is necessary for the BPS solution. The latter potential is chosen such that the BPS solutions have Gaussian tails (i.e. e −r 2 ), which means that we avoid a typical complication of the BPS solutions being compactons -as was the case studied in ref. [17]. The problem with the compactons, is that the fields deviate from the vacuum only on a compact region of space and the derivatives exhibit a jump from a negative value to zero at the border of said region. The true solution in the near-BPS regime does not possess such a discontinuity in the first derivative of the fields and therefore we had invented in ref. [17] a cusp condition to make the total field smooth, by inducing a counter-cusp to the perturbation field. This has all been avoided in the present paper. The motivation as stated above, is to get a pion mass that is of the same order as the kinetic term and both are proportional to a small parameter, , hence giving rise to physical pions in the theory and also to avoid the technical difficulties coming along with the mentioned cusps. The leading order correction to the energy in of an axially symmetric charge-N baby Skyrmion is also in this case obtained by evaluating the Dirichlet (kinetic) energy of the BPS solution, which luckily is finite. This correction is linear in . Bound states, however, cannot be studied unless we go to a higher order in and introduce a perturbation field along with suitable boundary conditions. This is because the tails of the BPS solutions fall off to spatial infinity and would just overlap each other if nothing extra is introduced in the theory. The key to the perturbation field was again a transverse field used in ref. [17], inspired by the work in ref. [20]. The perturbation field in this paper does not require a so-called cusp condition, mentioned above, which is a relieving simplification. First we studied the energies of the baby Skyrmions in the axially symmetric case up to order O( 3 ) which we call N 2 LO and find excellent agreement between the perturbative scheme and full numerical ODE calculations.
Then we move on to performing large-scale full numerical brute-force computations of Skyrmions with charges Q = 4 and Q = 2, which include the simplest bound states of the most stable baby Skyrmions in the theory. For the Q = 4 sector, the most stable solution is a bound state of two N = 2 approximately axially symmetric baby Skyrmions quite well separated, but still having their tales intertwined in a bound state. The axially symmetric Q = N = 4 solution is unstable for small , but another metastable solution composed of four N = 1 baby Skyrmion in a nearly triangular arrangement was found. Interestingly, this solution was completely triangularly symmetric in the compacton case studied in ref. [17], whereas in this case two of the solitons move close together and repel the other two, one more than the other. This is probably explained by a quadrupole force, which is incompatible with the triangular symmetry. In the Q = 2 topological sector, the most stable solution is simply the Q = N = 2 axially symmetric solution, but for small , there is also a bound state of two N = 1 baby Skyrmions. Finally, we would like to stress that we have been able to find numerical BPS solutions for all the mentioned cases with = 0 set strictly to zero. That is, these solutions are nontrivial BPS solutions for which we do not yet know a suitable analytic Ansatz.
The last part of the paper is an attempt at calculating the binding energies perturbatively in our semianalytic expansion scheme, which was very successful for the compacton case of ref. [17]. Unfortunately, it turns out that this perturbative scheme is unreliable for baby Skyrmions with (Gaussian) tails, basically because the background field configuration is no longer close to the true solution and nonlinearities become crucial, thus invalidating the linearized approach used in the perturbative scheme. For more details, see the list at the end of the previous section.
Interestingly, it appears that the binding energy in this model, where the BPS solutions have a Gaussian tail, has a leading term which is linear in -in stark contradistinction to the case of the compactons of ref. [17], where the binding energy appears only at the order O( 2 ).
Although the compactons require painful cusp conditions for allowing one to study their near-BPS limit, the solitons with tails turn out to be even more difficult and can probably be well understood only by using nonlinear techniques or simply full numerical computations.
It would be interesting to consider other potentials and as we found in this paper, compared with ref. [17], everything depends strongly on the choice of the potential. Potentials in the BPS sector, giving rise to a power-law tail [24], would probably provide very different properties with respect to the expansion as well as to the solutions in general. We will leave the investigation of such cases to future studies.
A BPS solution with potentialm 2 1 V 1 + V 2 Writing the BPS equation for the Skyrme term, but including bothm 2 1 V 1 + V 2 instead of only V 2 (as in eq. (2.13)), we have Switching to stereographic coordinates yields ij ∂ i ω∂ jω Notice that the solution (A.9) reduces to eq. (2.22) in the limit ofm 1 → 0. The compacton radius is finite for any finite value ofm 1 > 0, but in the limit ofm 1 → 0 the compacton radius R compacton tends to infinity. This is consistent with the solution tending to eq. (2.22), which indeed has a tail that falls off exponentially (or rather like a Gaussian). Finally, the BPS mass of the solution (A.9) is given in eq. (2.29).