Near-BPS Skyrmions: Constant baryon density

Although it provides a relatively good picture of the nucleons, the Skyrme Model is unable to reproduce the small binding energy in nuclei. This suggests that Skyrme-like models that nearly saturate the Bogomol'nyi bound may be more appropriate since their mass is roughly proportional to the baryon number A. For that purpose, we propose a near-BPS Skyrme Model. It consists of terms up to order six in derivatives of the pion fields, including the nonlinear and Skyrme terms which are assumed to be relatively small. For our special choice of mass term, we obtain well-behaved analytical BPS-type solutions with constant baryon density configurations, as opposed to the more complex shell-like configurations found in most extensions of the Skyrme Model. Fitting the four model parameters, we find a remarkable agreement for the binding energy per nucleon B/A with respect to experimental data. These results support the idea that nuclei could be near-BPS Skyrmions.


I. INTRODUCTION
One of the most original and successful attempts to describe the low-energy regime of the theory of strong interactions comes from an idea suggested by Skyrme [1] that baryons (and nuclei) are topological soliton solutions arising from an effective Lagrangian of mesons. The proposal is supported by the work of Witten [2] who realized that the large N c limit of QCD points towards such an interpretation. More recently, an analysis of the low energy hadron physics in holographic QCD [3] have led to a similar picture, i.e. the Skyrme Model. The model, in its original form, succeeds in predicting the properties of the nucleon within a precision of 30% [4]. This is considered a rather good agreement for model which involves only two parameters. Some attempts to improve the model have given birth to a number of extensions or generalizations. Most of them rely, to some extent, on our ignorance of the exact form of the low-energy effective Lagrangian of QCD namely, the structure of the mass term [5][6][7], the contribution of other vector mesons [8,9] or simply the addition of higher-order terms in derivatives of the pion fields [5].
Unfortunately, one of the recurring problems of Skyrme-like Lagrangians is that they almost inevitably give nuclei binding energy that are too large by at least an order of magnitude. Perhaps a better approach would be to construct an effective Lagrangians with soliton solutions that nearly saturate the Bogomol'nyi bound. If this indeed the case, then the classical static energy of such BPS-Skyrmions (Bogomol'nyi-Prasad-Sommerfeld) grows linearly with the baryon number A (or atomic number) much like the nuclear mass. Support for this idea comes from a recent result from Sutcliffe [10] who found that BPS-type Skyrmions seem to emerge for the original Skyrme Model when a large number of vector mesons are added. The additional degrees of freedom bring the mass of the soliton down to the saturation of the Bogomol'nyi bound. A more direct approach to construct BPS-Skyrmions was also proposed by Adam, Sanchez-Guillen, and Wereszczynski (ASW) [11]. Their prototype model consists of only two terms: one of order six in derivatives of the pion fields [12] and a second term, called the potential, which is chosen to be the customary mass term for pions in the Skyrme Model [13]. The model leads to BPS-type compacton solutions with size and mass growing as A with L 0 = −µ 2 V (U ) (2) where L µ = U † ∂ µ U is the left-handed current and we write for simplicity, the commutators as f µν = [L µ , L ν ] . Here the pion fields are represented by the SU (2) matrix U = φ 0 + iτ i φ i and obey the nonlinear condition φ 2 0 + φ 2 i = 1. The subscript i in L i denotes to the number of derivatives of the pion fields which determines how each term changes with respect to a scale transformation.
In the original Skyrme Model, only the nonlinear σ term, L 2 , and the Skyrme term, L 4 , contribute. This implies that α, β > 0 otherwise the static solution would not be stable against scale transformations. A mass term -or potential term -L 0 , is often added to take into account chiral symmetry breaking so as to generate a pion mass term for small fluctuations of the chiral field in V (U ). We shall analyze this term in more details in the coming sections but, as it turns out, the choice of potential V (U ) will have a direct bearing on the form of the solutions and on the predictions of our model. Finally, the term of order six in derivatives of the pion fields, L 6 , is equivalent to L J6 = −ε J6 B µ B µ with ε J6 = 9π 4 λ 2 /4 that was first proposed by Jackson et al. [12] to take into account ω-meson interactions. Here, B µ stands for the topological current density The constants µ, α, β, and λ are left as free parameters although we shall focus on the regime where α and β are relatively small, i.e. in the limit where the solutions remain close to that of the BPS-solitons. It is well known that setting the boundary condition for U at infinity to a constant in order to get finite energy solutions for the Skyrme fields also characterizes such solutions by a conserved topological charge which Skyrme identified as the baryon number B (or mass number A in the context of nuclei) Note that the static energy arising from L 6 corresponds to the square of the baryon density It is often associated to the energy that would emerge if the Skyrme field is couple to the ω−meson [19] E ω = 1 2 where instead of following the e −mω|r−r | / |r − r | law, the interaction is replaced by a δ−function δ 3 (r − r ). Historically, L 0 and L 6 were introduced to provide a more general effective Lagrangian than the original Skyrme Model and indeed, the Lagrangian in (1) represents the most general SU (2) model with at most two time derivatives. Since one generally relies on the standard Hamiltonian interpretation for the quantization procedure, higher-order time derivatives are usually avoided. On the other hand, it should be kept in mind that an effective theory based on the 1/N c expansion of QCD should, in principle, include terms with higher-order derivatives of the fields.
The model (1) has been studied rather extensively in the sector where the values of parameters µ, α, β, and λ close to that of the original Skyrme Model [12,20]. Clearly these choices were made so that L 2 and L 4 would continue to have a significant contribution to the mass of the baryons and thereby preserve the relative successes of the Skyrme Model in predicting nucleon properties and their link to soft-pion theorems (α is proportional to the pion decay constant F π ). Yet this sector of the theory fails to provide an accurate description of the binding energy of heavy nuclei.
Noting that this caveat may come from the fact that the solitons of the Skyrme Model do not saturate the Bogomol'nyi bound, ASW proposed a toy model [11] (equivalent to setting α = β = 0) whose solutions are just BPS solitons. In principle however, the model cannot lead to stable nuclei since BPS-soliton masses are exactly proportional to the topological number, so B > 1 solutions have no binding energies. A more realistic approach was proposed in Refs. [14,15] where the Lagrangian (1) is assumed to be in the sector where α and β are relatively small, treating these two terms as perturbations. The solutions almost saturate without reaching the Bogomol'nyi bound so that it allows for small but non-zero binding energies. However, in spite of a very good agreement with experimental nuclear masses, there remain a few obstacles to the acceptance of such model. For instance, nuclear matter is believed to be uniformly distributed inside a nucleus whereas the solutions of the aforementioned models [11,14,15] display either compact, shell-like or gaussian-like baryon and energy densities respectively. The main purpose of this work is to demonstrate that it is possible to construct an effective Lagrangian which leads to a uniform densities and still preserves the agreement with nuclear mass data.
Let us consider the static solution for U . It can be written in the general form wheren is the unit vectorn = (sin Θ cos Φ, sin Θ sin Φ, cos Θ) (9) and F, Θ, and Φ depend in general on the spherical coordinates r, θ, and φ.
We first consider the model in (1) in the limit where α and β are small. For that purpose, we introduce the axial solutions for the α = β = 0 case, where A is an integer that correspond to the baryon number or mass number of a nucleus. A word of warning is in order here. The solution (10) is only one of an infinite dimensional families of solutions of the BPS model and, is not expected to be the true minimizing solution of the static energy of the model or, for that matter, of the total energy which includes also the (iso)rotational energy, the Coulomb energy and an isospin symmetry breaking term. Since α and β are assumed to be small, the nonlinear σ and Skyrme terms are not expected to a determining factor in minimizing the total energy. In fact, the dominant effect should come from the repulsive Coulomb energy which would have a tendency to favor a most symmetric configuration. Which form is the true minimizer remains an open question only to be answered by heavy numerical calculations. In the absence of such an analysis and for the sake of simplicity, we chose to consider ansatz (10) which allows to easily estimate all the contributions to the mass of the nuclei.
From hereon, we shall use whenever possible the dimensionless variable x = ar where a = (µ/18Aλ) 1/3 in order to factor out the explicit dependence on the model parameters µ, α, β, and, λ and baryon number A. In fact, most of the relevant quantities can be written in terms of three fundamental objects The total static energy E s gets a contribution from each term in (1), respectively, where I k lmn are parameter-free integrals given by and write I k lmn = I k lmn (∞) for simplicity. Note that some of these integrals are related in our case since I k lmn = A n I k l,m+n,0 . In the last equality, we assume that one can recast V (F ) as a power series of sin F, i.e. V (F ) = m C V m sin m F as suggested in Ref. [5]. The terms E 0 and E 6 are proportional to the baryon number A as one expects from solutions that saturate the Bogomol'nyi bound whereas the small perturbations E 2 = A 1/3 (a 2 + b 2 A 2 ) and E 4 = A −1/3 (a 4 + b 4 A 2 ) have a more complex dependence. Part of this behavior, the overall factor A ±1/3 , is due to the scaling. The additional factor of A 2 comes from the axial symmetry of the solution (10) that can be factored out from I k lm2 = A 2 I k l,m+2,0 . The topological charge also simplifies to The root mean square radius of the baryon density is given by which is consistent with experimental observation for the charge distribution of nuclei r 2 . The minimization of the static energy for α = β = 0 leads to the differential equation for F : Multiplying by ∂ x F, this expression can be integrated which leads to where x 0 is an integration constant. Finally, the expression for F (x) can be found analytically provided the integral on the left-hand side is an invertible function of F. For example, assuming that the potential may be written in the form where u = cos (F/2) and, g (u) = ∂g/∂u, equation (19) leads to Such solutions saturate the Bogomol'nyi bound [11], so their static energy is proportional to the baryon number A. One would like ultimately to reproduce the observed structure of nuclei, i.e. a roughly constant baryon density becoming diffuse at the nuclear surface which is characterized by a skin constant thickness parameter. Unfortunately the chiral angle F in (21) cannot reproduce this last feature since F can only be a function of the ratio r/A 1/3 . So the resulting thickness parameter is not constant and should scale like A 1/3 . It is interesting to note that (18) implies that for the minimum energy solutions where the last equality arises from Derrick scaling. Furthermore according to (7) and (22), the square root of the potential where B 0 (x) corresponds to the radial baryon density B 0 (x) = dΩ B 0 (x). Thus, in order to obtain a nonshell baryon density, it suffices to construct a potential V that does not vanish at small x or, equivalently, a solution such that ∂ x F (0) = 0. Expression (20) must be used with caution: it only applies for potentials V which turn out to be function of u alone or, in other words, for potentials that depends on the real part of the pion field matrix U or TrU. On the other hand, L 0 in (1) needs to be explicitly written in terms of the fields U . A simple but not unique approach to construct such potential is to identify u = cos(F/2) to the expression where U ± = (2I ± U ± U † )/8 and I is the 2 × 2 identity matrix. Then, a convenient expression for V (U ) is given by In the context of the BPS-Skyrme Model, not only the potential V appears as one of the dominant term in the static energy but it is also a key ingredient in the determination of the solution. In principle, the full effective theory including the potential should emerge from the low-energy limit of QCD, but apart from a few symmetry arguments, little is known on the exact form of V . A most simple expression for V that reads was first proposed by Adkins et al. [13] and served as an additional term to the original Skyrme Lagrangian. Its main purpose was to recover the chiral symmetry breaking pion mass term − 1 2 m 2 π π · π in the limit of small pion field fluctuations U = exp(2iτ a π a /F π ). It is sometimes useful to recast the potential in the form [5] Taking the limit of small pion field fluctuations, this allows fixing the parameter µ in terms of the pion mass m π through the relation ∞ k=1 The choice of potential (23) corresponds to the choice g(u) = u 3 /3 in (21) and solving for F leads to the BPScompacton solution of ASW [11]: Note here that ∂ x F (x) diverges as x → 3 −1/3 which implies that E 2 and E 4 are not well defined. Unfortunately, this solution as well as those arising from other similar models [21] saturate the Bogomol'nyi bound and as such, they give no binding energies for the classical solitons with B > 1. Several alternatives to (23) have also been proposed [5,7] but recently, the major role played by the potential in the predictions for BPS-Skyrme Models was realized and it has led to a few interesting cases: • One such example is a potential based on Ref. [14] V BoM (U ) = −8Tr U + U 3 − which correspond to the choice −C 1 = C 2 = C 3 = 4C 4 = µ 2 /128 and C k>4 = 0 in (24). It leads to well-behaved solutions where ∂ x F remains negative and finite for all x. In order to set the baryon number to A, the boundary conditions are chose to be F (0) = ±π and F (∞) = 0 for positive and negative baryon number respectively. Note that the exponential fall off of F at large x prevents some quantities such as the moments of inertia from becoming infinite. However, ∂ x F (x) vanishes at x = 0 and so does the baryon density, leading to an unsatisfactory shell-like configuration.
• In that regard, a solution similar to that proposed in Ref. [15] seems more appropriate since it possesses the kind of non-shell like baryon density configurations observed in nature. It emerges from the potential of the form These models display compact, shell-like or gaussian-like baryon and energy densities (see Figs. 1 and 2). However here, we shall demonstrate that it is possible to construct an effective Lagrangian which leads to a uniform baryon density and still preserves and even improves the agreement with nuclear mass data.  If we assume for now that the observed baryon density can be appropriately approximated by the parametrization ρ B (r, A) then, one is looking for a solution for F (r) such that Separating variables and integrating both sides of the equation we get the expression of the form where In order to be consistent, the boundary conditions for Z must obey Z(∞) − Z(0) = −π/2. Matching expressions (21) and (29) then provides an approach to construct a model, i.e. to choose a potential V , that reproduces the empirical baryon density ρ B . Again we stress that our model leads to BPS-Skyrmions with a profile F that must be a function of the ratio r/A 1/3 . Unfortunately, this excludes most parametrizations in the literature, for example, densities such as the 2-parameter Fermi or Wood-Saxon form since they tend to reproduce two empirical observations: (a) a baryon density that is roughly constant for all nuclei up to their boundary where (b) it is suppressed within a thickness t ≈ 4.4τ that is practically constant. The last behavior is inconsistent with the r/A 1/3 dependence of F.
Let us instead construct our model by modifying the gaussian-like profile F BHM (x) in such a way that baryon density B 0 (x) is approximately constant. The solution F BHM (x) leads to a nonshell baryon density but it falls off too rapidly. In order to suppress this behavior we propose a solution of the form (see Figs. 1 and 2) and fix the coefficient a 4 = 7/5 by setting to zero the first coefficient of the series expansion of B 0 (x) near x = 0.
(Note that we could, in principle, extend this procedure by changing the argument of the exponential to a truncated series X(x) = x 2 + N i=2 a 2i x 2i . Imposing that the density remains constant further from the core would require to set a 6 = 1384/525, a 8 = 6302/1125, and so on.). It is easy to find a potential that would allow such a solution Note that in the limit of small pion field fluctuations U = exp(2iτ a π a /F π ), the potential has no quadratic term in the pion field i.e. the pion mass remains zero in this model. where the last result is obtained assuming the axial solution (10).
Using the profile F in (30), the static energy in (12) can be calculated. Recalling that I k lmn = A n I k l,m+n,0 for the form of axial solution at hand, we need to evaluate numerically only four parameter-free integrals: In order to represent physical nuclei, we have taken into account their rotational and isorotational degrees of freedom and quantize the solitons. The standard procedure is to use the semiclassical quantization which is described in the next section.

III. QUANTIZATION
Skyrmions are not pointlike particles so we resort to a semiclassical quantization method which consists in adding an explicit time dependence to the zero modes of the Skyrmions and applying a time-dependent (iso)rotations on the Skyrme fields by SU (2) matrix A 1 (t) and A 2 (t) where R(A 2 (t)) is the associated SO(3) rotation matrix. The approach assumes that the Skyrmion behaves as a rigid rotator. Upon insertion of this ansatz in the time-dependent part of the full Lagrangian (1), we can write the (iso)rotational Lagrangian as where The moment of inertia tensors U ij are given by where T i = iU † τi 2 , U . The expressions for W ij and V ij are similar except that the isorotational operator T i is replaced by a rotational analog S i = − ikl x k L l as follows: Following the calculations in [14] for axial solution of the form (10), we find that all off-diagonal elements of the inertia tensors vanish. Furthermore, one can show that U 11 = U 22 and U 33 can be obtained by setting A = 1 in the expression for U 11 . Similar identities hold for V ij and W ij tensors. Finally the general expressions for the moments of inertia coming from each pieces of the Lagrangian read (37) where due to the axial form of our solution, we can extract an explicit dependence on A through the relation I k lmn = A n I k l,m+n,0 . The axial symmetry of the solution imposes the constraint L 3 + AK 3 = 0 which is simply the statement that a spatial rotation by an angle θ about the axis of symmetry can be compensated by an isorotation of −Aθ about the τ 3 axis. It follows from expressions (33)-(35) that W 11 = W 22 = 0 for |A| ≥ 2 and A 2 U 33 = AW 33 = V 33 . Otherwise, for |A| = 1, the solution have spherical symmetry and 8α a 2 I 2 020 + 64β I 2 220 + I 2 040 + 9λ 2 a 2 4 where here A = 1 in a as well.
The general form of the rotational Hamiltonian is given by [22] where (K i ) L i the body-fixed (iso)rotation momentum canonically conjugate to (a i ) b i . It is also easy to calculate the rotational energies for nuclei with winding number |A| ≥ 2 These momenta are related to the usual space-fixed isospin (I) and spin (J) by the orthogonal transformations According to (41) and (42), we see that the Casimir invariants satisfy K 2 = I 2 and L 2 = J 2 so the rotational Hamiltonian is given by We are looking for the lowest eigenvalue of H r which depends on the dimension of the spin and isospin representation of the nucleus eigenstate |N ≡ |i, i 3 , k 3 |j, j 3 , l 3 . For α = β = 0, we can show that ξ is negative and we shall assume that this remains true for small values of α and β. Then, for a given spin j and isospin i, κ must take the largest possible eigenvalue k 3 . Note that K 2 = I 2 and L 2 = J 2 , so the state with highest weight is characterized by k 3 = i and l 3 = j. Furthermore, since nuclei are build out of A fermions, the eigenvalues k 3 are limited to k 3 ≤ i ≤ A/2. On the other hand, the axial symmetry of the static solution (10) implies that k 3 = −l 3 /A ≤ j/A where j ≤ A/2 as well.
In order to minimize H r , we need the largest possible eigenvalue k 3 , so for even A nuclei, κ must be an integer such that Similarly for odd nuclei, |k 3 | must be a positive half-integer so the only possible value is This last relation only holds for the largest possible spin eigenstate j = A/2 which is not the most stable in general and so it signals that the ansatz (10) may not be the most appropriate for odd nuclei. The axial symmetry may however be only marginally broken if we consider the odd nucleus as a combination of an additional nucleon with an even nucleus especially for large nuclei. Nonetheless, we shall retain the ansatz (10) for both even and odd nuclei and choose the largest possible eigenvalue k 3 for the most stable isotopes as κ = 0 for A = even The lowest eigenvalue of the rotational Hamiltonian H r for a nucleus is then given by [14] E r = 1 2 The spins of the most abundant isotopes are well known. This is not the case for the isospins so we resort to the usual assumption that the most abundant isotopes correspond to states with lowest isorotational energy. Since i ≥ |i 3 |, the lowest value that i can take is simply |i 3 | where i 3 = Z − A/2. For example, the nucleon and deuteron rotational energy reduces respectively to The explicit calculations of the rotational energy of each nucleus then require the numerical evaluation of the following four parameter-free integrals in (36), (37) and (38) which, in our model, turn out to be I 2 020 = 0.142868 I 2 220 = 1.43364 I 2 040 = 0.352712 I 2 240 = 3.94598.
So far, both contributions to the mass of the nucleus, E s and E r , are charge invariant. Since this is a symmetry of the strong interaction, it is reflected in the construction of the Lagrangian (1) and one expects that the two terms form the dominant portion of the mass. However, isotope masses differ by a few percent so this symmetry is broken for physical nuclei. In the next section, we consider two additional contributions to the mass, the Coulomb energy associated with the charge distribution inside the Skyrmion and an isospin breaking term that may be attributed to the up and down quark mass difference.

IV. COULOMB ENERGY AND ISOSPIN BREAKING
The electromagnetic and isospin breaking contributions to the mass have been thoroughly studied for A = 1, mostly in the context of the computation of the proton-neutron mass difference [23][24][25], but are usually neglected, to a first approximation, for higher A since they are not expected to overcome the large binding energies predicted by the model. There are also practical reasons why they are seldom taken into account. The higher baryon number configurations of the original Skyrme Model are nontrivial (toroidal shape for A = 2, tetrahedral for A = 3, etc.) and finding them exactly either requires heavy numerical calculations (see for example [26]) or some kind of clever approximation like rational maps [27]. In our case however, we are interested in a precise calculation of the nuclei masses and an estimate of the Coulomb energy is desirable, and even more so in our model which generates nonshell configurations. It turns out that the axial symmetry of the solution and the relatively simple form of the chiral angle F (r) in (30) simplify the computation of the Coulomb energy.
Let us first consider the charge density inside Skyrmions. Following Adkins et al. [4], we write the electromagnetic current with B µ the baryon density and J µ3 V the vector current density. The conserved electric charge is given by The vector current is then defined as the sum of the left and right handed currents where U ij and W ij are the moment of inertia densities in (33)-(35). The calculations of the Coulomb energy here follows that in [28]; it differs that from Ref. [15] where only the body-fixed charge density was considered. The anticommutator is introduced to ensure that J 0i V is a Hermitian operator. In the quantized version, a j and b j are expressed in terms of the conjugate operators K i and L i . Here we only need the relation The solution is axially symmetric then the off-diagonal elements of U ij and W ij vanish, W 11 = W 22 = 0 for |A| ≥ 2 and AU 33 = W 33 . Then have Inserting a i in (50), the isovector electric current density reduces to where U ii /U ii may be interpreted here as a normalized moment of inertia density for the i th component of isospin in the body-fixed frame. The expectation value R(A 1 ) 31 K 1 and R(A 1 ) 32 K 2 for eigenstate |N = |i, i 3 , k 3 |j, j 3 , l 3 are equal so that we may simplify where we have used relation (41). The moment of inertia density are given by The expression in brakets in equation (51) integrates to zero so that one recovers the relation Z = A/2+i 3 as expected. But while it does not contribute to the total charge, the charge density is not zero everywhere. Let us examine this contribution in more details. Since the electric charge does not depend on the angular momentum, we can limit our analysis to the isospin wavefunctions. Following Adkins [29] we write the wavefunctions A 1 |i, i 3 , k 3 in terms of the Wigner's functions D n mm : Similarly the matrix R(A 1 ) 33 corresponds to a spin zero and isospin zero transition that can be written The appropriate expectation value is then given by where the last two expressions on the second line are Clebsch-Gordan coefficients. Recalling that we have imposed the condition |k 3 | = κ = 0 or 1/2 for even and odd nuclei respectively and fixed the value of the isospin to i = |i 3 |, we find The last term drops for even nuclei (κ = 0). For odd nuclei, the cancellation in the brackets leads a relatively small contribution which is further suppressed by the factor κ 2 / (i + 1) for large nuclei. It is indicative of the asymmetry in the moments of inertia. The Coulomb energy associated to a given charge distribution ρ(r) takes the usual form Since we have at hand an axially symmetric distribution, it is convenient to expand ρ(r) in terms of normalized spherical harmonics to perform the angular integrations Following the approach described in [30], we define the quantities which, at large distance, are equivalent to a multipole moments of the distribution. The total Coulomb energy is given by The isocalar part to the charge distribution is a spherically symmetric contribution where I k lmn is defined in (13). On the other hand, the isovector contribution in (??) possesses a simple angular dependence so that the summation (56) consists of only two terms in Y 0 * 0 and Y 0 * 2 . The moments Q 00 and Q 20 are then given by where and C = k 2 3 /i(i + 1). Finally, the Coulomb energy then takes the form It is again convenient to regroup the model parameters in the dimensionless quantity such that we may write Here, each element of M ij 00 (M ij 20 ) comes from squaring Q 00 (Q 20 ) in (58) and depend only on the form of the profile F (x) and baryon number A according to The Coulomb energy can explain part of the isotope mass differences, but it is certainly not sufficient. For example for the nucleon, the Coulomb energy would suggest that the neutron mass is smaller than that of the proton. Of course, one can invoke the fact that isospin is not an exact symmetry to improve the predictions. Several attempts have been proposed to parametrize the isospin symmetry breaking term within the Skyrme Model [24,25]. Here we shall assume for simplicity that this results in a contribution proportional to the third component of isospin where the parameter a I is fixed by setting the neutron-proton mass difference to its experimental value ∆M expt n−p = 1.293 MeV. Since both of them have the same static and rotational energies, we find where E n C and E p C are the neutron and proton Coulomb energy, respectively. Summarizing, the mass of a nucleus reads where E s is the total static energy. The prediction depends on the parameters of the model µ, α, β, and λ and the relevant quantum numbers of each nucleus as shown in (62).

V. RESULTS AND DISCUSSION
The values of the parameters µ, α, β and λ remain to be fixed. Let us first consider the case where α = β = 0. This should provide us with a good estimate for the values of µ, α, β, and λ required in the 4-parameter model (1) and, after all, it corresponds to the limit where the minimization of the static energy leads to the exact analytical BPS solution in (30). For simplicity, we choose the mass of the nucleon and that of a nucleus X with no (iso)rotational energy (i.e. a nucleus with zero spin and isospin) as input parameters. Neglecting for now the Coulomb and isospin breaking energies, the mass of these two states is according to expression (62) For example, if the nucleus X is Calcium-40, a doubly magic number nucleus, with mass E Ca = 37214.7 MeV, then solving for λ and µ, we get the numerical values µ = 12322.3 MeV 2 , α = β = 0 and, λ = 0.00474078 MeV −1 which we shall refer as Set I. The masses of the nuclei are then computed using Eq. (62) which results in predictions that are accurate to at least 0.6%, even for heavier nuclei. This precision is somewhat expected since the static energy of a BPS-type solution is proportional to A so if it dominates, the nuclear masses should follow approximately the same pattern. However, the predictions remain surprisingly good compared to that of the original Skyrme model, another 2-parameter model.
Perhaps even more relevant are the predictions of the binding energy per nucleon B/A = (E − Zm p − (A − Z)m n ) /A, in which case, the calculation simplifies. For example, subtracting from the static energy of a nucleus from that of its constituents we find that the binding energy does not depend on the static energies E 0 or E 6 , whereas the contribution from E I simply cancels out. The dominant contributions come from the (iso)rotational and Coulomb energy differences, respectively, dominated by AE N r for large nuclei and which is, of course, negative due to the repulsive nature of the Coulomb force between nucleons.
The results for B/A are presented in Fig. 3 (dashed line). They are compared to the experimental values (empty circles). We show here only a subset of the table of nuclei in [31] composed of the most abundant 140 isotopes. The parameters of Set I lead to a sharp rise of the binding energy per nucleon at small A followed by a slow linear increase for larger nuclei. The accuracy is found to be roughly within 10% which is relatively good considering the facts that the model involves only two parameters at this point and the calculation involves a mass difference between the nucleus and its constituents.
Experimentally the charge radius of the nucleus is known to behave approximately as r 2 em 1 2 = r 0 A 1 3 with r 0 = 1.23 fm. It is straightforward to calculate the root mean square radius of the baryon density [see Eq (16)] which leads to r 2 1 2 = (2.007 fm) A where ρ(r) is given in expression (56) and r 2 V is given by where for the sake of conciseness we wrote r 2 V in terms of U i.e. the integrals along the radial component in U (2) 11 contains an extra factor of r 2 . Our computation verifies that the charge radius obeys roughly the proportionality relation ∼ r 0 A 1 3 but overestimates the experimental value of r 0 by about 80% with parameter Set I.
Let us now release the constraint α = β = 0, and allow for small perturbations from the nonlinear σ and Skyrme term. In order to estimate the magnitude of the parameters α and β in a real physical case, we perform two fits: the four parameters µ, α, β and λ in Set II optimizes the masses of the nuclei while Set III reaches the best agreement with respect to the binding energy per nucleon, B/A. Both fits are performed with data from the same subset of the most abundant 140 isotopes as before. The best fits on both cases would lead to small negative values for β similar to that of Refs. [14,15]. However, since the classical (static) energy of the model is unbounded below if α, β < 0 we impose the constraint α, β ≥ 0 from hereon to avoid stability problems. (Note that in principle β could take small negative values as long as the Skyrme term is overcome by the repulsive Coulomb energy in which case the physical nuclei would be stable but not the classical soliton.) A summary of the results is presented in Table I   We find that the two new sets of parameters are very close to Set I. In order to make a relevant comparison, we look at the relative importance of the four terms in (1) and how they scale with respect to the parameters of the model, namely for L 0 , L 2 , L 4 , and, L 6 , respectively. So the nonlinear σ and Skyrme terms are found to be very small compared to that of L 0 and L 6 , i.e. by at least five orders of magnitude. This provides support to the assumption that (30) is a good approximation to the exact solution.
The energy scale µλ remain approximately the same for all the sets while the values of µ and λ shows noticeable differences. In particular for the fit involving B/A turns out to be somewhat sensitive to these variations mostly because it involves a mass difference. We also note some variation in the baryonic charge radius r 0 = 1.3982 (λ/µ) 1/3 ; all sets overestimates the experimental value by roughly 80%. Since setting the parameters mainly involves fixing the relevant energy scale µλ, perhaps the process may not be as sensitive to setting a proper length scale for the nucleus so the predicted value of r 0 should probably be taken as an estimate rather than a firm prediction.
Matching the parameters of the model with that of the original Skyrme Model, we identify F π = 4 √ α, e 2 = 1/32β whereas m π = 0 due to the form of the potential. The quantities F π and e 2 take values which are orders of magnitude away for those obtained for the Skyrme Model (see Table I) but this is not surprising since we have assumed from the start that α and β are relatively small. Unfortunately, one of the successes of the original Skyrme Model is that it established a link with soft-pion physics by providing realistic values for F π , m π and baryon masses. Such a link here is more obscure. The departure could come from the fact that the parameters of the model are merely bare parameters and they could differ significantly from their renormalized physical values. In other words, we may have to consider two quite different sets of parameters: a first one, relevant to the perturbative regime for pion physics where F π and m π are closer to their experimental value and, a second set which applies to the nonperturbative regime in the case of solitons. In our model, this remains an open question.
The model clearly improves the prediction of the nuclear masses and binding energies in the regime where α and β are small. Let us look more closely at the results presented in Fig. 3. The experimental data (empty circles) are shown along with predicted values for parametrization Set I , Set II and Set III (dashed, dotted and solid lines, respectively). Setting α = β = 0 (Set I) leads to sharp increase B/A at low baryon number followed by a regular but slow growth in the heavy nuclei sector. This suggests that heavier nuclei should be more stable, in contradiction to observation. However the agreement remains within ∼ 10% in regards to the prediction of the nuclear masses. This is significantly better than what is obtained with the original Skyrme Model which overestimates B/A by an order of magnitude. Since B/A depends on the difference between the mass of a nucleus and that of its constituents, it is sensitive to small variation of the nuclear masses so the results for B/A may be considered as rather good. The second fit (Set II) is optimized for nuclear masses. The behavior at small A is similar to that of Set I (as well as in Set III) while it reproduces almost exactly the remaining experimental values (A 40). Finally, the optimization of B/A (Set III) provide a somewhat better representation for light nuclei at the expense of some of the accuracy found in Set II for A 40. Overall, the binding energy is rather sensitive to the choice of parameters. This is partly because the otherwise dominant contributions of E 0 and E 6 to the total mass of the nucleus simply cancel out in B/A.
The difference of behavior between light and heavy nuclei shown by the model may be partly attributed to the (iso)rotational contribution to the mass. The spin of the most abundant isotopes remains small while isospin can have relatively large values due to the growing disequilibrium between the number of proton and the number of neutron in heavy nuclei. On the other hand, the moments of inertia increase with A, so the total effect leads to a (iso)rotational energy E r < 1 MeV with A > 10 for all sets of parameters considered and its contribution to B/A decreases rapidly as A increases. On the contrary, for A < 10, the rotational energy is responsible for a larger part of the binding energy which means that B/A should be sensitive to the way the rotational energy is computed. So clearly, the variations in shape of the baryon density has some bearing on the predictions for the small A sector not only the values of the parameters.
To summarize, the main purpose this work is to propose a model in a regime where the nuclei are described by near-BPS solitons with approximately constant baryon density configuration. This is acheived with a 4-terms generalization of the Skyrme Model in the regime where the nonlinear σ and Skyrme terms are considered small. The choice of an appropriate potential V allows to build constant baryon density near-BPS solitons, i.e. a more realistic description of nuclei as opposed to the more complex configurations found in most extensions of the Skyrme Model (e.g. A = 2 toroidal , A = 3 tetrahedral, A = 4 cubic,...). Fitting the model parameters, we find a remarkable agreement for the binding energy per nucleon B/A with respect to experimental data. On the other hand, there remain some caveats. First, the Skyrme Model provides a simultaneous description for perturbative pion interactions and nonperturbative baryon physics with realistic values for F π and m π and baryon masses. The connection between the two sectors here seems to be much more intricate. Secondly, there may be place for improvement by proposing more appropriate solutions that would describe equally well the light and heavy nuclei. Finally, the model seems unable to reproduce a constant skin thickness in the baryon or charge density and the experimental size of the nucleus correctly. On the other hand, the concept of BPS-type Skyrmions also arises when one adds a large number of vector mesons to the Skyrme Model as suggested by recent results based on holographic QCD from Sutcliffe [10]. Unfortunately, the emerging large A Skyrmions configurations are rather complicated or simply unknown so that it has yet been impossible to perform an analysis of the nuclear properties comparable to that presented in this work. More recently Adam, Naya, Sanchez-Guillen and Wereszczynski [28] considered the special case of the pure BPS-model (α = β = 0) using the potential V ASW . Although their treatment differ slightly they find a similar agreement for the binding energy per nucleon. Yet, all approaches clearly suggest that nuclei could be treated as near-BPS Skyrmions.