BPS Skyrmions as neutron stars

The BPS Skyrme model has been demonstrated already to provide a physically intriguing and quantitatively reliable description of nuclear matter. Indeed, the model has both the symmetries and the energy-momentum tensor of a perfect fluid, and thus represents a field theoretic realization of the"liquid droplet"model of nuclear matter. In addition, the classical soliton solutions together with some obvious corrections (spin-isospin quantization, Coulomb energy, proton-neutron mass difference) provide an accurate modeling of nuclear binding energies for heavier nuclei. These results lead to the rather natural proposal to try to describe also neutron stars by the BPS Skyrme model coupled to gravity. We find that the resulting self-gravitating BPS Skyrmions provide excellent results as well as some new perspectives for the description of bulk properties of neutron stars when the parameter values of the model are extracted from nuclear physics. Specifically, the maximum possible mass of a neutron star before black-hole formation sets in is a few solar masses, the precise value depending on the precise values of the model parameters, and the resulting neutron star radius is of the order of 10 km.


Introduction
The calculation of physical observables of strongly interacting matter at low energies -relevant, e.g., to nuclear physics -directly from QCD is a notoriously difficult problem, which led to the introduction of low-energy effective field theories (EFTs) as a more tractable alternative. The Skyrme model is a well-known example of such a low-energy EFT. It was introduced originally by Skyrme [1] as a purely mesonic nonlinear field theory for the description of nuclei. Skyrme's idea was that nucleons should be described as a kind of "vorticity" in a mesonic "fluid" or, in a more modern language, as topological solitons of the underlying mesonic nonlinear field theory. And, indeed, the Skyrme model is known to possess topological solitons ("Skyrmions") whose topological index is identified with the baryon number. The original idea of Skyrme gained further support when it was observed that QCD in the limit of a large number of colors (large N c ) becomes a theory of weakly interacting mesons (interaction strength ∼ N −1 c ) [2]. Such weakly interacting nonlinear field theories frequently possess solitonic solutions with soliton masses proportional to the inverse of the (weak) coupling, which in the present case of the large N c mesonic model of QCD are identified with baryons and nuclei, recovering thereby the proposal of Skyrme.
The Skyrme model has been applied to the description of nuclei with notable success, e.g., in the description of rotational excitation bands of some light nuclei [3,4]. The version of the model originally proposed by Skyrme, however, has some drawbacks in the description of physical nuclei. First of all, Skyrmions with higher baryon number B have rather high binding energies (i.e., masses significantly below B times the B = 1 Skyrmion mass, see, e.g., [5]), which is in striking contrast to the low binding energies of physical nuclei. Also, Skyrmions for large baryon number tend to form crystals of lower B substructures [6,7], which is at odds with the liquid-type behavior of physical heavy nuclei. These problems recently led to propose several "near BPS" Skyrme models, that is, generalizations of the original Skyrme model which are close to BPS models [8,9]. Here by a BPS model we understand a field theory which has both an energy bound for static field configurations which is exactly linear in the baryon charge B and solutions saturating the bound for all values of B (we shall assume B ≥ 0 in the sequel, i.e., consider only matter not antimatter). The original Skyrme model is not BPS. It has a lower topological energy bound, but it may be shown easily that this bound cannot be saturated. Specifically, we consider the following near BPS Skyrme model [8]  (for the moment in flat Minkowski space; we use the "mostly minus" metric sign convention diag(g μν ) = (+, −, −, −)), and L 0 = −λ 0 U (tr U ), Here U : Maurer-Cartan current and U is a potential. The λ n are dimensionful, non-negative coupling constants, and B μ is the topological or baryon number current giving rise to the topological degree (baryon number) B ∈ Z, is BPS. That is to say, the static energy functional E 06 [U ] has an energy bound linear in B and (in fact, infinitely many) minimizing field configurations saturating the bound for each B, [8]. Further, this energy functional is invariant under volume-preserving diffeomorphisms (VPDs) on physical space, which are the symmetries of a perfect fluid. The energy-momentum tensor of the model L 06 is, in fact, the energy-momentum tensor of a perfect fluid, as we shall see below. These findings lead to the intriguing hypothesis that the near-BPS Skyrme model (1) might be the correct low-energy EFT for the description of nuclear matter, as the BPS submodel L 06 already provides a rather good description of some of its static properties. Indeed, the BPS Skyrme model allows for a very accurate description of nuclear binding energies [10,11], especially for heavy nuclei. It is the purpose of the present letter to couple the BPS Skyrme model to gravity and to use the resulting self-gravitating BPS Skyrmions for the description of neutron stars. We remark that there already exist several attempts to describe neutron stars using the original Skyrme model. In [12] the hedgehog ansatz for higher B was coupled to gravity but it turned out that -as in the non-gravitating case -higher B hedgehogs are not stable. In [13,14] approximate Skyrmion configurations based on rational maps were used. Probably the most promising attempt within this context is using Skyrmion crystals [15,16] because Skyrmion crystals are the true minimizers of the original Skyrme model for large B. The crystal structure is, however, at odds with the fact that, most likely, the core of neutron stars is in a superfluid phase. Also, full numerical calculations are not possible in this case such that certain assumptions about the right equation of state of Skyrme crystals under strong gravitational fields must be made. An accessible review can be found in [17].

BPS Skyrme model and parameter values
Conveniently redefining the coupling constants λ 6 = λ 2 /(24) 2 and λ 0 = μ 2 , the static energy functional of the theory is Its BPS bound (where √ U S 3 is the average value of √ U on the target space SU(2) ∼ S 3 ) is saturated by infinitely many BPS solutions [8,18,19], and the corresponding BPS equation is We now have to determine the values of the parameters λ and μ to be used in our calculations. The product m ≡ λμ has the dimensions of mass (energy; we use units where the speed of light c = 1). Further, l ≡ (λ/μ) 1/3 has the dimensions of length.
We fit m by requiring that the BPS Skyrmion mass is B times one-fourth of the mass of the helium nucleus, E 06 = Bm N where m N = m He /4 = 931.75 MeV. We use m N instead of the nucleon mass m N ∼ 940 MeV because the latter will receive contributions from (iso)spin excitations in a Skyrme model description, but these are absent for helium. Even helium receives small (e.g., coulombic) contributions in addition to the Skyrmion mass, but the uncertainty will be at most a few MeV. To fix l, we use the fact that BPS Skyrmions for many potentials (in particular, for the potentials considered in this letter) are compactons with a strictly finite volume V , and this volume is the same for all solutions with a given baryon number B and is exactly linear in B. This permits to define a Skyrmion radius R via V = (4π /3)R 3 . We now require that this radius coincides with the nucleon radius r N = 1.25 fm for B = 1, i.e., R = r N B 1/3 . We think that the fit for the mass parameter m is quite accurate, because by far the biggest contribution to the nuclear masses must always come from the Skyrmion mass.
On the other hand, the fit for the length parameter l is probably less precise. Firstly, although the compacton radius is quite natural, there are additional definitions for radii (diverse charge radii) which could be used. For compactons these charge radii are always smaller than the compacton radius, indicating that the latter could be slightly bigger than the nucleon radius. Secondly, going beyond the BPS submodel by including, e.g., the Dirichlet term L 2 , the effects of the pion cloud will tend to increase the radius, indicating that the compacton radius R without pion cloud could be somewhat smaller. To summarize, although our simple fit for l certainly provides a reasonable value, the true best fit value could easily deviate about 20% or 30% in either direction. Determining this true value, however, requires the knowledge of the complete low-energy EFT with all terms (also the non-BPS ones) included, which is beyond the scope of this letter.
Concretely, we shall consider the pion mass potential U π = 1 − cos ξ and the pion mass potential squared U 4 = U 2 π with a quartic behavior near the vacuum (here U = exp(ξ n · τ ), n 2 = 1 and τ are the Pauli matrices) with energies, compacton radii and fit values Bλμ, These expressions for the energies and compacton radii may be calculated directly from the potentials, see [20] (knowledge of the Skyrmion solutions is not required).

BPS Skyrmions coupled to gravity
The action of the BPS Skyrme model in a general metric g ρσ leads to the energy-momentum tensor (B ρ is defined in Eq. (4)) (12) which is the energy-momentum tensor of a perfect fluid (the perfect-fluid property of the term L 6 alone, as well as its coupling to gravity, have already been discussed in [21]), (13) where the four-velocity u ρ , energy density ρ and pressure p are In the static case, and for a diagonal metric (which is sufficient for our purposes) we have u ρ = ( g 00 , 0, 0, 0) and T 00 = ρ g 00 , In the flat space case, e.g., this implies that the pressure must be constant (zero for BPS solutions, nonzero for non-BPS static solutions [20]), as a consequence of energy-momentum conservation, whereas ρ will be a nontrivial function of the space coordinates. In general, ρ and p will be quite arbitrary functions of the spacetime coordinates, so there does not exist a universal equation of state (EoS) p = p(ρ) which would be valid for all solutions.
We now want to couple the BPS Skyrme model to gravity and solve the resulting Einstein equations for a static, spherically symmetric metric which in standard Schwarzschild coordinates reads For us the following observation is crucial. The above ansatz for the metric together with the axially symmetric ansatz for the Skyrme field with baryon number B ξ = ξ(r), n = (sin θ cos Bφ, sin θ sin Bφ, cos θ) (19) leads to a baryon density B 0 , energy density ρ and pressure p which are functions of r only. The ansatz is, thus, compatible with the Einstein equations (here G ρσ is the Einstein tensor and κ 2 = 16π G = 6.654 · 10 −41 fm MeV −1 ) and the static Euler-Lagrange equations for the Skyrme field, and reduces these equations to a system of ordinary differential equations (ODEs) in the variable r for the three unknown functions A(r), B(r) and ξ(r). Before presenting this system of ODEs and the results of a numerical integration, we want to make some comments. Firstly, in flat Minkowski space the same axially symmetric ansatz (19) (but referring to spherical polar coordinates in that case) was used in the calculations of nuclear binding energies in [10]. As said, the resulting binding energies are very accurate for heavier nuclei, but, nevertheless, once additional terms (like, e.g., the Dirichlet term E 2 ) are taken into account, there are strong arguments indicating that the axially symmetric BPS Skyrmions are not the adequate ones (they do not minimize E 2 among all BPS Skyrmions) [22]. An improved calculation using the true minimizers of E 2 and taking the contribution of E 2 into account should lead to even better results for the binding energies.
Here we just want to emphasize that in the case of self-gravitating BPS Skyrmions the axially symmetric ansatz leading to a spherically symmetric metric, energy density and pressure is the correct one, essentially because gravity straightens out all deviations from spherical symmetry. Secondly, in the subspace of spherically symmetric solutions we may define a kind of EoS p = p(ρ), because both ρ and p are functions of r. We find numerically that a simple power law p = aρ b (21) reproduces this EoS with a high precision. Here, however, a and b are not universal constants. Instead, they depend on the baryon number B. In particular, for "small" baryon number (small compared, e.g., to the solar baryon number B ), where the effect of gravity may be neglected, the constant a vanishes, lim B→0 a = 0 (the pressure is zero like in the case without gravity). If we treated the gravitational coupling κ as a parameter which may vary then, of course, it would also hold that lim κ→0 a = 0.

Numerical results
We find it convenient to introduce the new target space vari- where ρ and p for the axially symmetric ansatz (h r ≡ ∂ r h) read We integrate the system (22)  from Eq. (23) that, for a non-singular metric function B, p at r = R must obey p (R) = 0 which leads to a condition on h r (R), con- In the numerical integration, the free parameter value ρ 0 is varied until this condition is met. Formal solutions which do not obey this condition produce metric functions B which are singular at r = R.
In particular, such a metric function cannot be joined smoothly to the Schwarzschild solution in empty space (for r ≥ R) and is, therefore, physically unacceptable.
We find the following behavior in the numerical integrations.
We remark that neutron star masses up to about M ∼ 2M are firmly established, whereas there are indications for masses up to about 2.5M , see e.g. [24,25] for an overview of recent measurements. The results of our calculations are, therefore, in excellent agreement with these observations, indicating that our model provides a very good description of the bulk properties of nuclear matter also in the presence of the gravitational interaction. Concerning the radii, we remark that the observational results are less precise. Besides, R is the geometric radius which leads to a neutron star surface area of 4π R 2 , whereas when comparing to measurements sometimes other radii are more appropriate, like the proper distance from the origin to the surface, R = R 0 dr √ B(r), or the radiation radius R * = R √ B(R). Both R and R * are somewhat bigger than R because B(r) ≥ 1. In any case, also our values for the radii are in the expected range of about R ∼ 10-20 km. As said already, our fit for the unit of mass m is quite precise (determined by the nuclear mass m N ), but the unit of length l is less so, therefore it is interesting to study the sensitivity of both M max and R max under a change of the length scale, l → l = αl. We find numerically that both M max and R max approximately change by a factor of α (3/2) under this rescaling.
Finally, we show our main numerical results in Figs. 1-4. Concretely, in Fig. 1a we plot the neutron star mass as a function of the non-gravitational Skyrmion mass, both in solar units. We find that for the extremal case M max the gravitational mass loss is about 25%. In Fig. 1b we plot M against the (geometric) neutron star radius R. We find that even in the extremal case the neutron star radius is about a factor of two above the Schwarzschild radius. In Fig. 2 we show the equation of state for different values of B (concretely for n = 1 in Fig. 2a, and for n max in Fig. 2b) together with the fit function p = aρ b for appropriate values of a, b.
In Fig. 3, we plot the metric function B(r) for several values of the baryon number n = B/B close to its maximum value n max .
We find for both potentials that the maximum value which B takes for n = n max is about B max ∼ 2.7. It is interesting to compare this finding with the analogous result for the Skyrmion crystal of Ref. [16]. There the authors calculated the minimum value of B −1 (which was called S in that paper) for different solutions and always found that S min > 0.4, which translates into B max < 2.5.
So the B max we find for the maximum mass case is slightly bigger (i.e., the induced self-gravitation slightly stronger), but still  quite similar to the result of [16]. The position of the maximum of B(r) is quite close to the neutron star surface for the potential U π , whereas it is shifted towards the center for U 2 π . This is related to the fact that, for U 2 π , the energy density is more concentrated about the center (see Fig. 4).
In Fig. 4, we plot the energy densities for several values of the baryon number close to n max . We find that, especially for the potential U 2 π , the energy density is quite sharply concentrated about the center. This may look surprising at first sight, but is simply related to the shape of the potential U 2 π , which is quite peaked about the anti-vacuum (h = 1). Indeed, the BPS equation (8) just states that the baryon density is proportional to the square root of the potential, so peaked potentials lead to peaked baryon density (and energy density) profiles already in the case without gravity. It is perhaps more instructive to compare the central energy density of the case without gravity to the central energy density for n max . The central energy density for the case without gravity does not depend on the baryon number B and is given by . Using the parameter values (9), we find for U π : ρ BPS (r = 0) = 4(m/l 3 ) = 285 MeV fm −3 . The central energy density for n max is, therefore, about 2.7 times the nongravitational energy density ρ BPS (r = 0), see Fig. 4a. Similarly, we get for U 2 π : ρ BPS (r = 0) = 8(m/l 3 ) = 909 MeV fm −3 . In this case, the central energy density for n max is just about 2.2 times the nongravitational energy density ρ BPS (r = 0), see Fig. 4b. These results in both cases indicate a rather high stiffness of the effective (onshell) EoS of strongly self-gravitating BPS Skyrmions, i.e., a nuclear matter which is only weakly compressible in strong gravitational fields. This result, again, compares quite well with the Skyrmion crystal results of Ref. [16], where a compression of the central energy density by not more than a factor of three is observed for all solutions.

Discussion
We used the BPS Skyrme model (6) for the description of neutron stars and found that by simply fitting the two model parameters to the nucleon mass and radius we already get very rea- sonable results for the resulting neutron star masses and radii. In particular, for the maximum possible neutron star mass we find M max = 2.44M or M max = 3.73M , respectively, for the two potentials considered. This compares extremely well with the observational constraint M max ∼ 2.5M . We take this, together with the perfect fluid behavior of the model, as a further very strong indication that, indeed, the BPS Skyrme model provides the most important contribution to the static bulk properties of nuclear matter. In a strict sense, our results are not yet final predictions of neutron star properties, because genuine predictions require the knowledge of the full near-BPS Skyrme model (1) together with the values of all its coupling constants, which should follow from an application to nuclear physics and the corresponding detailed fit to nuclear data. The full near-BPS Skyrme model may also lead to a further improvement in the description of neutron stars, in the following sense. Even if the additional (standard Skyrme) terms are quite unimportant in the bulk, this is not true at the surface, because at the surface the Skyrme field is close to its vacuum value, and the term L 6 approaches the vacuum much faster than the standard Skyrme model terms. The standard Skyrme model is known to prefer crystalline structures for large B, so crystalline structures ("neutron star crust") can be expected at the surface of a neutron star described by the near-BPS Skyrme model, whereas the bulk and core remain in a fluid phase. But precisely this structure is expected in current models of neutron stars (see, e.g., [26]).
When compared with other, more traditional methods of nuclear physics, the advantage of the (near-) BPS Skyrme model at this moment is not so much its ability to make quantitative predictions -although this, too, should change with more detailed investigations and with advanced numerical methods, assisted by a rigorous analytical control which follows from the integrability properties of the BPS model. After all, the methods and models of nuclear physics are well developed and lead to very precise descriptions of nuclei and nuclear matter. However, a drawback of many models of nuclear physics is that they are tailor-made to describe rather specific physical phenomena, therefore it is difficult to use them for extrapolations to new phenomena or parameter values where they have not been employed before. We think it is one of the outstanding features of the BPS Skyrme model that it captures a generic property of (bulk) nuclear matter and allows, therefore, for far-reaching extrapolations. Concretely, in the present letter we extrapolated from B = 1 (which provided the parameter fit values) to B ∼ 10 57 (the neutron star) and from a nonrelativistic to a highly relativistic regime, with very accurate results. In other words, the (near) BPS Skyrme model provides a unified description of nuclear matter, reaching from nucleons and atomic nuclei to neutron stars.
There are two particular (related) results of our calculations which are somewhat different from most traditional nuclear physics calculations of neutron stars using the TOV equations (22), (23), although they are completely compatible with all observational data. In the traditional approach, the metric function B(r), the energy density ρ(r) and the pressure p(r) are considered as independent field variables, so the two TOV equations (22), (23) must be closed by a third equation. For this, usually a universal algebraic equation of state (EoS) p = p(ρ) resulting from the thermodynamic limit of a nuclear effective field theory (EFT) (like Quantum Hadron Dynamics (QHD) [27]) is assumed. In our model, on the other hand, we find that already the EFT itself is of the perfect-fluid type defining its own energy density and pressure, both of which depend on the metric in an explicit fashion. It is, therefore, not possible to define a universal, algebraic off-shell EoS, and the true off-shell EoS relating ρ and p is a complicated and metric-dependent differential equation. We remark that our offshell EoS share some features with the "quasi-local" EoS explicitly depending on the geometry (e.g., metric or curvature), which were introduced in [28] for the description of anisotropic stars and further studied in [29] and, in relation with neutron stars, in [30]. It turns out that in stars with anisotropic matter such quasi-local EoS are even required for consistency [28]. In our case, it is still possible to find (numerically) an on-shell algebraic EoS for solutions ρ(r) and p(r), but this on-shell EoS is no longer universal and depends on the neutron star mass or baryon number B. This does not mean that the EoS of nuclear matter depends on the sample size. The EoS for the BPS Skyrme model without gravity is always the same, p = 0 at equilibrium (nuclear saturation), for arbitrary B.
The B dependence of the on-shell EoS for self-gravitating nuclear matter in the BPS Skyrme model is exclusively a consequence of self-gravitation. Due to the nonlinearity of gravity, the effects of self-gravitation are stronger for larger B (larger neutron star mass) and the effective on-shell EoS, therefore, gets stiffer. Concretely, This increasing stiffness has a particular physical effect in the cases we considered, namely a neutron star radius R which grows with the neutron star mass M, i.e., (dM/dR) > 0 (except for stars very close to their maximum mass in the case of the potential U π , see Fig. 1). This behavior is at variance with the results found for solutions of the TOV equations for many (fixed, universal) nuclear physics EoS, where the neutron star radius is either essentially constant for a range of neutron star masses or even shrinks with increasing mass [24]. The reason for this behavior is that for a fixed EoS the increasing strength of self-gravitation for larger masses may collapse the star to much higher densities and, for softer EoS, even to smaller sizes. Only sufficiently stiff universal EoS are compatible with (dM/dR) > 0. We remark that one particular case of an EoS which is sufficiently stiff to support (dM/dR) > 0 for almost all values of M is precisely given by the Skyrme crystal of Ref. [16]. The M(R) curve found there is, in fact, quite similar to the one we find for the pion mass potential U π , see Fig. 1b.
In the BPS Skyrme model, the squeezing effect of nonlinear selfgravitation is balanced by the increasing stiffness of the on-shell EoS. We emphasize that, at present, (dM/dR) > 0 is compatible with observations and that the observational data are not yet sufficiently precise to settle this question. If (dM/dR) > 0 finally turns out to be true, this either rules out a large class of EoS which are well motivated from nuclear physics, because only very stiff fixed EoS are compatible with (dM/dR) > 0. Or it may indicate that in the traditional derivation of the EoS from an EFT like QHD one has to go beyond mean field theory, such that backreaction effects of gravity on the EoS may be taken into account, at least for nuclear matter in sufficiently strong gravitational fields. A detailed discussion of these issues will be given elsewhere. In any case, the qualitative results we found for the EoS within the BPS Skyrme model also point towards possible improvements of the standard nuclear physics approach to neutron stars in strong gravitational fields.
There are many ways in which the present investigation can be deepened and extended. One obvious possibility is to use additional potentials and to study how the shapes of these potentials influence the properties of the resulting neutron stars, e.g., which maximal masses can be reached and for which potentials the relation (dM/dR) > 0 remains true. Another interesting research direction is related to rotating neutron stars and to neutron stars in magnetic fields. In principle, both these tasks are rendered feasible by the fact that it is known how to rotate Skyrmions (for a recent discussion see, e.g., [31]) and what is the correct, QCD induced coupling of Skyrmions to the electromagnetic interaction [32]. Still, the resulting systems are no longer spherically symmetric, so a full system of PDEs has to be solved numerically in these cases. A further step in the analysis would be to use the full near-BPS Skyrme model as a basis for the calculation of neutron star solutions and properties, but here, in a first step, the detailed application of the near-BPS Skyrme model without gravity to nuclei and nuclear matter is required. As the BPS and integrability properties are no longer available in this case, full three-dimensional numerical calculations will be necessary.