Skyrme crystals, nuclear matter and compact stars

A general review of the crystalline solutions of the generalized Skyrme model and their application to the study of cold nuclear matter at finite density and the Equation of State (EOS) of neutron stars is presented. For the relevant range of densities, the ground state of the Skyrme model in the three torus is shown to correspond to configurations with different symmetries, with a sequence of phase transitions between such configurations. The effects of nonzero finite isospin asymmetry are taken into account by the canonical quantization of isospin collective coordinates, and some thermodynamical and nuclear observables (such as the symmetry energy) are computed as a function of the density. We also explore the extension of the model to accommodate strange degrees of freedom, and find a first order transition for the condensation of kaons in the Skyrme crystal background in a thermodynamically consistent, non-perturbative way. Finally, an approximate EOS of dense matter is constructed by fitting the free parameters of the model to some nuclear observables close to saturation density, which are particularly relevant for the description of nuclear matter. The resulting neutron star mass-radius curves already reasonably satisfy current astrophysical constraints.


Introduction
With the advent of gravitational wave and multi-messenger astronomy, the available constraints on the equation of state (EOS) of neutron stars, namely, strongly interacting matter at finite density, have improved significantly in the last decade. Furthermore, additional gravitational wave measurements from neutron star mergers with improved precision from LIGO/Virgo are expected to impose even stricter constraints on the EOS in the near future.
The physics of strong interactions is described by a nonabelian gauge theory, Quantum Chromodynamics (QCD), whose nonperturbative nature at small and intermediate energy or density regimes makes the computation of nuclear matter properties from first principles extremely challenging. Indeed, a perturbative treatment is only available at asymptotically high densities, while the main nonperturbative computational tool, lattice QCD, can be used for arbitrarily low temperatures but only small densities, due to the fermion sign problem. Hence, for sufficiently high densities (but not asymptotically high), the phase diagram of QCD is poorly understood, and there is a big theoretical uncertainity on basic observables such as the relevant degrees of freedom or the EOS in this regime, which, on the other hand, is precisely the relevant region for the matter in the interior of neutron stars. Indeed, as depicted in the schematic phase diagram of fig. 1, matter inside neutron stars is expected to be both at finite baryonic and isospin densities, and almost zero temperature (when compared with the relevant density scale). As first principle approaches are currently not suitable for this task, a plethora of different effective theories and phenomenological models have been proposed to either qualitatively or quantitatively study the QCD phase diagram in the intermediate density regime.
Standard nuclear physics methods like relativistic mean field theory [1] or chiral perturbation theory [2,3] can be used for this purpose (for a recent review we refer to [4]). These theories are effective field theories (EFT) in the sense that they introduce nucleons and mesons as their basic fields, instead of the quarks and gluons of QCD. Further, they depend on a rather large number of a priori unknown parameters which are usually determined by fitting to nuclear forces and scattering data in the (non-relativistic) low energy regime. Also, nucleons (and baryons in general) are treated as point particles. The extrapolation of these models to baryon densities beyond nuclear saturation might, therefore, be affected by rather large uncertainties [5].
The Skyrme model represents a slightly different type of EFT. This model, proposed more than 60 years ago by Tony Skyrme [6,7], is an EFT of chiral mesons only, whereas nucleons and baryons are realized as topological solitons ("skyrmions") of the mesonic field, and their interactions are described by the same mesonic lagrangian. In addition, the topological degree of these solitonic solutions can be identified with the baryon number. The argument for the identification of skyrmions with baryons and nuclei was further strengthened within the large N c limit of QCD [8,9], for which the theory becomes an effective weakly interacting model of mesons, where baryons satisfy the usual properties of solitons. The Skyrme model incorporates in a completely natural fashion several important features of strong interaction physics, like chiral symmetry and its breaking, the conservation of baryon number, or the extended character of nucleons. As a consequence of the latter, it also avoids short-distance singularities in nucleon-nucleon interactions. In addition, it naturally incorporates the spin-statistics theorem in that skyrmions quantized with a half-odd integer spin are fermions, whereas for integer spin they are bosons [10]. For a recent review we refer to [11].
Skyrme s idea of baryons as topological solitons is shared with a number of other approaches. First of all, holographic models based on the conjectured gauge/gravity duality have proven to be useful as tools for studying the nonperturbative regimes of QCD-like theories. For a recent review on the application of holographic models to the description of neutron stars, see [12]. The Skyrme model has, in fact, been shown to appear as a holographic boundary theory in several holographic QCD models such as the Sakai-Sugimoto model. A flat space version of these holographic models was proposed in [13], such that the holographic dual of the flat space instantons provides the Skyrme field coupled to an infinite tower of vector mesons. The full holographic dual maintains the conformal symmetry inherited from the instanton, whereas any truncation to a finite number of vector mesons breaks conformal invariance. Secondly, in a different but related approach a class of generalized nuclear effective theories have been developed in the last two decades which maintain the field contents and topological structure of the Skyrme model at low energies while flowing towards an effective theory compatible with the symmetries of QCD at higher energies (for a recent review see [14] in this special issue). This flow is achieved by a coupling of the Skyrme field to the dilaton [15] and an infinite tower of vector mesons via hidden local symmetry [16,17], where these fields flow from a spontaneously broken to an unbroken phase. As a consequence, these nuclear effective theories still share some topological structures with the Skyrme model, while their dynamical contents is, in general, different and much harder to calculate. For that reason, instead of the prohibitively difficult full numerical simulations of those field theories, sometimes it is assumed that they inherit some topological structures from the Skyrme model, and the consequences of these assumptions are worked out using more standard effective field theory methods of nuclear physics. We shall comment on some of these assumptions in our conclusions. A detailed discussion of all the issues mentioned in this last paragraph can be found in [18].

The Skyrme model
As said, Skyrme s original motivation was to reproduce baryons as topological soliton solutions from a purely mesonic field theory. In its simplest version, the basic fields of the Skyrme model are just the pion fields combined in an SU (2) valued field U as U = σI + iπ a τ a , U † U = I −→ σ 2 + π a π a = 1, a = 1, 2, 3 (1.1) where τ a are the Pauli matrices. Left (L) and right (R) chiral transformations act on this field like U → LU R. The two simplest terms of the resulting effective (Skyrme model) lagrangian can be easily guessed from low-energy considerations, namely the so-called nonlinear sigma model (or Dirichlet) term providing a kinetic term for the pions, and the pion mass potential Here, f π is the pion decay constant, whose physical value in the conventions used here is f π = 186 MeV. Further, L µ are the components of the left-invariant Maurer-Cartan form of the SU (2) group, m π is the pion mass, and I is the 2 × 2 identity matrix. These two terms alone, however, cannot support the existence of static (soliton) solutions, because they are unstable under a spatial rescaling x → Λx. To stabilize them, Skyrme added the term where e is a dimensionless parameter. Here the notation L n means that the corresponding term contains n first derivatives. L 4 is the only Poincare invariant four derivative term that leads to a positive hamiltonian which is quadratic in momenta.
In the original version [6,7], Skyrme only considered the model L 24 = L 2 + L 4 . The resulting static solutions of the Skyrme field are maps from real space R 3 to the SU (2) group manifold, which is S 3 . Besides, in order to have finite energy solutions, the field must tend to the vacuum of the theory at spatial infinity. We take U (x → ∞) = I, such that the infinity of R 3 is compactified into a point. Then the base manifold has the topology of the S 3 , so the Skyrme field maps the S 3 onto itself and we may conclude from homotopy theory that the Skyrme model allows for the existence of topologically nontrivial configurations [19]. These solitonic solutions can be classified according to their different topologies via the topological number B, and the idea of Skyrme was to identify this integer number with the baryon number, and these topological solutions, which are called skyrmions, with baryons [20,21]. An explicit expression for B in terms of the Skyrme field can be found from the topological current (1.5) It is straightforward to see that the divergence of the topological current vanishes identically, (∂ µ B µ = 0), which implies the conservation of the topological number, which is an integer. Furthermore, the choice of a vacuum value for the Skyrme field at spatial infinite represents the spontaneous chiral symmetry breaking of QCD. These analogies between the Skyrme model and QCD at low energies support the idea that the fundamental fields of the Skyrme model, π a , may be identified with the physical pions. In our approach we will also add a term which is of sixth order in first derivatives, where λ is a parameter. This term is, again, singled out as being the only Poincare invariant term of sixth order which leads to a standard hamiltonian quadratic in momenta. We will argue that this term is, in a certain sense, the most important one for our purposes. The generalized Skyrme model that we will consider is, therefore, L = L 2 + L 4 + L 6 + L 0 (1.8) and we will refer to the individual terms as the quadratic, quartic, sextic and potential terms, respectively. For the simplest model L 24 , the collective coordinate quantization of the spin and isospin degrees of freedom of the B = 1 skyrmion allowed to describe the proton, the neutron and some higher excitations (e.g., the delta resonance) and calculate some of their observables with a precision of about 30% [22], as could be expected naively from the large N c arguments with an expansion parameter N −1 c for N c = 3. A better precision, therefore, requires the inclusion of further terms and further (meson) fields into the effective model. The simplest Skyrme model L 24 also has some problems in the description of higher B skyrmions which should correspond to atomic nuclei with weight number A = B. While the model is partially successful in describing some nuclear spectra in terms of spin and isospin excitations, one major problem is that the binding energies of skyrmions of baryon charge B against their decomposition into B nucleons are up to ten times higher than the binding energies of the corresponding physical nuclei. Also, the resulting skyrmions for large B are rather hollow structures, at variance with the quite constant baryon densities inside physical nuclei. The inclusion of the mass term for the pions, apart from reproducing the explicit chiral symmetry breaking, already improves some of these shortcomings. Indeed, this term induces, e.g., the α-clustering for larger values of B [23], which is a known property of some nuclei. As a consequence, for some nuclei which are known to have α particle subclusters, the Skyrme model with pion mass term already provides an excellent description of nuclear spectra [24], particularly when vibrational degrees of freedom are taken into account in addition to spin and isospin [25]. However, the binding energy problem remains.
Finally, the sextic term was first considered in [26]. In [27] it was shown that combined with a potential term it leads to a BPS model which reproduces many important features of physical nuclei, like small binding energies and the spherical and compact shapes [28]. This potential should be chosen different from the pion mass term (e.g. L 0,k Tr{U −I} k , k ∈ N), because otherwise the BPS submodel would correspond to the unphysical limit of infinite pion mass. Based on additional BPS bounds for generalized Skyrme models discovered in [29,30], it was later found that these additional potentials L 0,k serve to reduce binding energies already on their own, both without [31] and with [32,33] the sextic term.
A further improvement was achieved by the inclusion of the rho mesons [34], based on the instanton-inspired approach to vector meson coupling in the Skyrme model of Ref. [13]. The cases B ≤ 8 were investigated numerically, and it was found that realistic cluster structures emerged for all skyrmions, even for substructures different from α particles. In addition, the binding energies were reduced significantly. The upshot of all this is that i) there has been significant progress in the last years in the Skyrme model as a model for nuclei and nuclear matter, where several shortcomings of the model have been improved significantly, and clear strategies for their resolution have been found; and ii) progress is still slower than one would like, not because of a fundamental problem of the model, but because calculations are hard. Already the starting point for any investigation, i.e., the skyrmion solution for a given B, is difficult to find, especially for the extended versions of the model where more terms and/or more fields are included. A parameter scan of the extended models with the aim of fitting the parameters and coupling constants of a model to physical observables in order to identify promising regions in parameter space is a challenging numerical problem with current techniques.
We shall, therefore, use a simpler version of the Skyrme model and fit only those observables which are most relevant for nuclear matter at sufficiently high densities because, as we will see, the most widely used configuration for Skyrme matter (the Skyrme crystal) allows to model nuclear matter only above nuclear saturation density 1 . At such high densities potential terms like L 0,k are irrelevant, as follows from simple scaling arguments, and may be ignored. For reasons of simplicity, we will also ignore vector mesons, although the adequacy of this approximation is more difficult to gauge 2 . On the other hand, the sextic term (1.7) will provide the leading contribution to the Skyrme crystal energy at high density, as follows from scaling arguments, again.

Skyrme matter
The inclusion of the sextic term is, in fact, mandatory for any realistic modeling of nuclear matter by Skyrme matter at sufficiently high densities, and its inclusion should provide already a rather reasonable description there. First of all, this term precisely describes the repulsive nuclear force in the high density regime [35]. Secondly, the Skyrme model without the sextic term leads to a maximum neutron star mass of about 1.5M [36], whereas NS of up to 2 solar masses are firmly established, and there are clear indications that the maximum possible NS mass may, in fact, be as high as 2.3 − 2.5M . Thirdly, a related fact is that the Skyrme model without the sextic term always leads to an EOS with a speed of sound which is below the conformal bound c 2 s < 1/3 for all pressures. But EOS with such a restricted speed of sound are strongly disfavored according to a recent analysis [37].
In principle, the ultimate goal of a sufficiently general Skyrme model would be to describe nuclear matter covering a large range of densities, from individual baryons to isolated nuclei where n B ∼ n 0 = 0.16 fm −3 , to the cores of neutron stars n B ∼ 10n 0 . We will find, however, that the currently available Skyrme matter solutions only allow to model nuclear matter above nuclear saturation density n 0 . We will, therefore, restrict to the model (1.8) which contains four parameters, the pion decay constant f π , the Skyrme parameter e, the pion mass m π , and λ 2 . From these parameters, we fix the pion mass to its physical value, m π = 140 MeV, whereas we will use the other three to fit our solutions to reproduce the nuclear matter properties most relevant for our considerations. In particular, as is common practice in the Skyrme model, we will not fix the pion decay constant f π to its physical value. The idea is that in the still rather restricted version (1.8) of the Skyrme model this modified value of f π partly takes into account the effect of neglected terms, whereas in more complete, general versions of the model the optimum value of f π should flow to its physical value.
We will consider static solutions of the Skyrme field and we adopt the usual Skyrme units of energy and length to work with more manageable expressions, (1.9) Then the static energy functional in these units becomes where we have defined c 6 = 2λ 2 f 2 π e 4 and c 0 = 2m 2 π /(f π e) 2 . In the last expression, we have introduced the field configuration eq. (1.1) and adopted the vectorial notation n A = (σ, π a ).
Recall that U ∈ SU (2) implies that n A is a unit vector, n A n A = 1, with A = 0, 1, 2, 3. It is also possible to find from (1.10) a lower bound for the energy per baryon number, which is known as the BPS (Bogomol'nyi-Prasad-Sommerfield) bound [29,30]. For this choice of units the standard Skyrme model (λ 2 = m π = 0) becomes independent of the value of the parameters f π and e, and its BPS bound is equal to 1.
Unfortunately, trying to numerically calculate the skyrmion solution for B >> 1, corresponding to a macroscopic amount of nuclear matter, is not feasible. Some simplifying assumptions, therefore, must be made for a description of nuclear matter within the Skyrme model. The simplest and most widely used assumption is that skyrmionic matter forms a cubic lattice [38][39][40][41][42][43][44][45]. That is to say, there exists a cubic unit cell with a baryon content B cell and side length l cell (i.e., volume V cell = l 3 cell ), such that the Skyrme field is periodic in all three Cartesian directions under the translation x i → x i + l cell . It is then sufficient to minimize the energy functional within one unit cell, and the extension to an arbitrary amount of Skyrme matter is trivial. In addition to periodicity, usually some further discrete symmetries are assumed, leading to different cubic crystal types like, e.g., simple cubic (SC), face-centered cubic (FCC), or body-centered cubic (BCC), and the resulting solutions are known as Skyrme crystals.
It must be emphasized, however, that the mere existence of solutions for any of these crystal types by itself does not imply their physical relevance. Indeed, the minimization of the Skyrme model energy functional leads to an elliptic PDE, and elliptic systems typically always have a solution for sufficiently regular boundary conditions. The reason that some crystals are considered relevant is that the corresponding solutions are of rather low energy. More precisely, the crystal energy per unit cell as a function of the lattice length, E cell (l cell ) is a convex function with a minimum for a certain l 0,cell . For some crystals, the resulting energy at the minimum l 0,cell is only slightly above the topological energy bound. For the FCC crystal for the original Skyrme model, e.g., E cell (l 0,cell ) is only about 3.7% above the bound, which is very close, because it is well-known that the bound cannot be saturated.
For l cell = l 0,cell , the energy grows rather quickly. Here, the region l cell > l 0,cell defines a thermodynamically unstable regime (formally, the pressure is negative). In this region, it is expected that equilibrium configurations of nuclear matter are, instead, given by larger clusters of nuclear matter (large nuclei, or a "nuclear pasta" phase in an intermediate region close to l 0,cell ) surrounded by almost empty space. We shall find some indications for the formation of larger nuclei in our investigations. This more inhomogeneous phase ameliorates the thermodynamical instability without resolving it completely. That is to say, the energy for l cell > l 0,cell still is slightly above E cell (l 0,cell ), but the difference between E cell (l 0,cell ) and, say, lim l cell →∞ E cell (l cell ) can be made smaller than 5%.
The nuclear pasta phase between l 0,cell and the low density phase where large nuclei appear could, in principle, lead to a thermodynamically stable description for all densities. The formation of nuclear pasta, however, is expected to result from a subtle interplay between the strong and the Coulomb forces. While the coupling of the Skyrme model to the electromagnetic field is known [46], and promising results concerning the computation of Coulomb effect in the case of α-particle skyrmions have been reported [47], a viable treatment of macroscopic amounts of skyrmionic matter coupled to electromagnetism which could give rise to the strong inhomogeneities implied by nuclear pasta is currently not available. As a consequence, all currently existing descriptions of nuclear matter based on Skyrme crystals approach zero pressure already at a finite baryon density n B = B cell /V cell , corresponding to the nuclear saturation density n 0 = B cell /V 0,cell . This implies that compact stars based on Skyrme crystals alone have no outer core and crust. On the other hand, the region n B ≤ n 0 is well described by standard methods of nuclear physics and can be joined with skyrmionic matter at higher densities. We shall make use of this possibility in several occasions.
The region l cell < l 0,cell of the Skyrme crystal is thermodynamically stable, but this does not necessarily imply that it provides the true minimum energy configuration for all baryon densities n B > n 0 . However, other good candidates for these minimizers have not been found within the Skyrme model, therefore we shall simply assume that they are given by Skyrme crystals and work out the consequences of this assumption as far as possible.
The paper is organized as follows: In Section 2, we introduce different types of classical Skyrme crystal solutions. We discuss their properties and possible phase transitions between them in several versions of Skyrme models. We also present evidence for a high density crystal-fluid transition and a low density transition from a crystal to a non-homogeneous phase. This section is partly based on our previous work [45], but most of the results have been newly computed. Section 3 is devoted to the semiclassical quantization of the Skyrme crystals, which allows to pass from symmetric to asymmetric nuclear matter. Here we show how to compute the symmetry energy and particle (proton, neutron and lepton) fractions within the Skyrme theory. This part reviews the material recently published in [48]. Next, in Section 4, we take into account the strangeness degrees of freedom by extending the Skyrme field to a SU (3)-valued matrix field. In particular, we compute the kaon condensation in the semiclassical Skyrme crystal, briefly reviewing the very recent results of [49]. In Section 5 we apply all the previously investigated crystal solutions to the description of neutron stars. Concretely, in section 5.2 we briefly review the construction of an EOS motivated by the generalized Skyrme model [50], which already leads to a very successful description of NS. In Section 5.3 we use the full semiclassical Skyrme crystal of the generalized model to describe both nuclear matter and NS, and we scan the model parameter space to find particularly promising parameter values. In Section 5.4, we include the effects of the kaon condensate on NS properties. Finally, Section 6 contains our conclusions and an outlook to further research.

Skyrme crystals
Skyrmions have been extensively studied, and solutions for finite values of B were found both in the standard Skyrme and the BPS submodels with different shapes and properties. The usual procedure to find a minimal energy solution considers the different possible symmetries for the skyrmion and then the solution is the one with the lowest energy. However, it becomes more difficult to find the minimal energy solutions for increasing B since the number of possible configurations grows quickly [51]. A simple estimation shows that the number of nucleons inside a NS must be of the order of M /m p ∼ 10 57 , then obtaining a solution via the usual procedure becomes an impossible task.
Skyrme crystals are solutions obtained imposing periodic boundary conditions, then they are infinitely spatially extended solutions so they formally have infinite baryon number. For this reason, skyrmion crystals are good candidates to describe infinite nuclear matter and to reproduce the conditions inside NS. To construct these periodic solutions, we split the crystal in finite unit cells where we construct the skyrmion configuration, then the main difference to obtain the crystal with respect to the isolated skyrmions lies in the boundary conditions. Now Skyrme crystals compactify the real space into T 3 , however since the T 3 is still an oriented and compact manifold, the Hopf's degree theorem ensures the existence of topological solitons labelled by an integer number.
From all the possible unit cells in three dimensions that we may use to construct a Skyrme crystal, we will consider cubic unit cells throughout this work, but we will allow for different symmetries within them. Additionally, since the crystal is infinitely extended it has infinite energy and baryon number, however the unit cell is finite in size, hence it carries a finite amount of energy and baryon number. Then the energy per unit cell as well as the energy per baryon number of the crystal are completely well defined and finite, The first Skyrme crystal was proposed in 1985 for the standard Skyrme model by Klebanov [38], motivated by the phenomenological application of crystals to the interior of NS. He considered the simplest possible crystal with a simple cubic (SC) unit cell, in which eight B = 1 skyrmions were located in the corners of the cube in the maximal attractive channel with respect to their nearest neighbours. Then, he computed the minimal energy field configuration respecting these conditions for different values of the unit cell side length and found that the lowest value of the energy was just an 8% above the BPS bound. In the following, we will explain how we construct Skyrme crystals via the procedure given in [40] with the different symmetries that have been proposed, and we will compare them within the generalized Skyrme model.

Crystal ansatz
The starting point in the construction of the Skyrme crystal proposed in [40] is the expansion of the fields in Fourier series,  Here, the unit cell extends from −L to L in all three Cartesian directions, so the side length of the unit cell is l = 2L. Then the symmetries of a given crystal impose some conditions on the Skyrme field which, in the end, are translated into some constraints on the coefficients of the expansions β abc and α hkl . Finally, the constrained coefficients are varied in order to obtain the lowest energy configuration. Although the expansion series of the fields are infinite, even the truncation to the first coefficients provides a good approximation to the solution, then the addition of higher modes produces corrections to the energy which become smaller for higher orders. This conclusion is also seen numerically, while the first coefficients are of order ∼ 1, we have calculated that the next orders decay to the ∼ 4%, ∼ 0.3% and ∼ 0.06% of the first order results. Hence we may safely truncate the series to a finite number of coefficients, we take around 30 coefficients to obtain the solution for each crystal. Finally, these expansions of the fields break the normalization condition of the vector n A , then we have to renormalize it The crystal considered by Klebanov has the simplest unit cell. It is invariant under cubic symmetry transformations: (σ, π 1 , π 2 , π 3 ) → (σ, π 2 , π 3 , π 1 ), (2.5) and it has an additional periodicity symmetry on the side length of the unit cell, Indeed all the symmetries shown in this work are based on the cubic symmetry, hence they will all satisfy symmetries A 1 and A 2 . The last symmetry, A 3 , repeats the location of a skyrmion with period L and performs the mutual isorotation between nearest neighbours. Under these symmetries, the energy is periodic in L but the fields are periodic in 2L, then we take the ranges of the unit cell to be [−L, L] and perform the integrals of the energy and baryon numbers within these limits.
We proceed now to show how to construct other symmetries of interest and show the numerical results.

Face centered cubic crystal of skyrmions
This symmetry was proposed in [41] in order to have a new solution with lower energy for very large values of L. It shares symmetries A 1 and A 2 and also two additional symmetries, C 3 : (x, y, z) → (x, z, −y), (σ, π 1 , π 2 , π 3 ) → (σ, −π 1 , π 3 , −π 2 ), (2.7) In this case, the energy and baryon number are periodic in 2L, and the unit cell has the shape on an FCC lattice of skyrmions. We have eight B = 1 skyrmions in the corners of the cube, symmetry C 4 locates other six skyrmions in centre of the faces and it also isorotates them by π with respect to their nearest neighbours. Hence, this lattice differs from the first in that each skyrmion is surrounded by 12 nearest neighbours all of them in the maximal attractive channel. Since we have the eight skyrmions in the corners and other six in the faces of the cube, the total baryon number in this unit cell is B cell = 4. As we mentioned before, these symmetries impose some constraints on the Fourier coefficients. They can be easily obtained imposing the symmetries on the field ansätze eq. (2.2). In this case the non-vanishing coefficients are obtained from the combination of the following restrictions, • h is odd, k and l are even or h is even, k and l are odd, • a, b, c are all odd or a, b, c are all even.
We show the resulting energy density contours in fig. 2, and the three-dimensional energy density plot in fig. 3. The crystal symmetry is clearly visible in the plots.  We adopt the Runge colouring convention [52] in this figure to represent the pion fields.

Body centered cubic crystal of half-skyrmions
This unit cell was proposed in [39] to be the one with the lowest energy for small values of L. It introduces an additional symmetry, to those of the Klebanov crystal (A 1 , A 2 and A 3 ). The motivation of this new crystal comes from the appearance of an additional symmetry when two B = 1 skyrmions are brought together and form the lower energy B = 2 field configuration, in which the B = 1 skyrmions have lost their individual identity. This new symmetry produces a unit cell which may be seen as a BCC of half-skyrmions, which are solutions for which σ = −1 at the centre of the cube of side length L until some radius r 0 for which σ = 0. This solution carries a half of baryon charge and it is undefined outside r 0 . However, a new half-skyrmion solution can be defined via the transformation (σ, π a ) → (−σ, −π a ), these new solutions are located in the corners of the cube, connected to the σ = 0 value at r = r 0 of the central half-skyrmion, forming a cube of side length L. As a result, the mean value of the σ field in this cube is exactly 0, so the energy coming from the potential term L 0 ∼ (σ − 1) will scale exactly as 8c 0 L 3 . Further, the 8 half-skyrmions in the corners contribute a total baryon number of 1/2, so the cube of side length L contains a baryon charge of B = 1. The energy and baryons densities are also periodic in L but, again, the fields have a 2L periodicity, then we have B cell = 8 within our unit cell of side length 2L. The restrictions imposed on the Fourier coefficients by the last symmetries are • h, k are odd, l is even.
• a, b and c are even.
We show the resulting energy density plots in fig

Enhanced face centered cubic crystal of skyrmions
This new crystal configuration was almost simultaneously found in two different publications [40,42] to be the one with the lowest energy in the standard Skyrme model. It may be seen as the half-skyrmion version of the FCC crystal explained before. Indeed it shares the symmetries A 1 , A 2 , C 3 and an additional symmetry, (σ, π 1 , π 2 , π 3 ) → (−σ, −π 1 , π 2 , π 3 ), (2.10) Figure 5: The tridimensional energy density plot of the BCC unit cell is shown. Again we represent the pion fields using the Runge colouring convention.
then some of the FCC crystal Fourier coefficients are set to 0 in this crystal.

Concretely, this new unit cell only allows the Fourier coefficients which satisfy the conditions
• h is odd, k and l are even, • a, b, c are all odd.
Since this crystal has less free Fourier coefficients than the FCC crystal, it is a particular case of the last one which we shall call FCC + . As a result, the FCC + crystal will always have equal or larger energy than the FCC crystal. This may lead to phase transitions between the crystals at some length of the unit cell, as we will see later. As in the FCC crystal, the half-skyrmion solutions with σ = −1 in their centre are located at the corners and faces of the unit cell. Further, the opposite half-skyrmions with σ = 1 occupy the body center and the link centers of the unit cube. As a consequence, the mean value of the σ field is 0, again, as in the BCC crystal. The energy and baryon densities are periodic in L and they have the appearance of a simple cubic unit cell of half-skyrmions. However, since the fields are periodic in 2L we take that to be the side length of the unit cell, hence the unit cell still has the shape of an FCC crystal with the alternating half-skyrmion solutions. Then, the baryon content within our unit cell is again B cell = 4. We show the resulting energy density plots in fig. 6 and in fig. 7.

Energy curves
Each crystal configuration explained before is constructed imposing the corresponding symmetry transformations, then we fix the value of L and use a Nelder-Mead algorithm [53] from the GSL C++ library to find the optimal values of the remaining free Fourier coefficients that minimize the energy functional. For different values of L we may construct the curve E(L) for all the different crystals and combinations of terms in the lagrangian eq. (1.8). For all cases we always find a convex curve with a minimum located at some value of L which is different in each case.
It will be important for the next sections to obtain an analytical expression of the energy curve. Then we may try to guess a specific fitting curve just by studying the scaling of each term in the energy functional. We find that E i ∼ L 3−i , where i represents  the number of derivatives in the i−term, and the three comes from the integration over 3D space. Hence we use the following fit for the energy curves in each case, In order to show the main properties of Skyrme crystals and the impact that the different terms have in the E(L) curve we fix the value of the physical constants that appear in the generalized Skyrme model to some standard values, f π = 129 MeV, e = 5.45, λ 2 = 5 MeV fm 3 , m π = 140 MeV. (2.12) We may see from fig. 8 how the E(L) curves change turning on and off the different terms. The pion-mass potential is an attractive term, so we would expect that more compact solutions are preferable. On the other hand the sextic term is repulsive so the opposite effect is expected. These behaviors are visible in fig. 8 where the length parameter value at which the minimal energy configuration is achieved, denoted by L 0 , is shifted to smaller values when the potential term is introduced, whereas it increases when the sextic term is present. Larger values of m π or λ 2 will increase these effects, and values of the parameters different from eq. (2.12) lead to completely different results. However the Skyrme units, in which we are computing the curves, produce a universal E(L) curve in the L 24 case in the sense that it does not depend on the value of the parameters. We also note that the expression eq. (2.11) that we are using to fit the energy of the crystals is quite accurate for the FCC + and BCC cases, however the FCC crystal has a non-trivial behaviour for large values of L, invalidating the expression of our fit in this region. For comparison, we also include the topological lower energy bounds in each case in fig. 8. For the L 24 case, this bound coincides with the Skyrme-Faddeev bound originally given already in [7], whereas more stringent bounds can be derived once more terms are included in the lagrangian [29,30]. Here, we always use the most stringent bound.
As in previous work [44,45], we conclude that the FCC + crystal reaches the lowest energy in the standard Skyrme model between the crystals considered here. In that case, the minimum is only a 3.7% above the BPS bound, which also places the Skyrme crystal as the skyrmion solution with the lowest energy ever achieved in the standard Skyrme model. Obviously, when the other terms are included the energy increases, but also the BPS bound of the model which is shown in each plot. We also conclude that for the parameters that we consider here, the lowest energy configuration is also achieved for the FCC + crystal in the generalized model. The location and the value of the energy in the minimum, L 0 and E 0 respectively, for this case, as well as the value that it is above the BPS bound in percentage, are shown in tables 1 and 2 below, The values of the coefficients in eq. (2.11) are obtained in each case using the Python optimization library GEKKO.
We remark that the fitting constants shown in tables 3 and 4 will change if we use different values of c 6 and c 0 . However it will be shown in the next section that a value for the constants of the energy curve fit independent of the parameters may be obtained for the FCC + crystal solution.

Phase transitions
We may anticipate from fig. 8 that even though the FCC + crystal reaches the lowest energy at the minimum, it may not be the crystal with the lowest energy for all values of L. This is clear in the region with large values of L, for which the FCC crystal has lower energy than the FCC + . We will also see that there is a phase transition from the FCC to the BCC crystal at small values of L, however since they do not have the same baryon content within the unit cell, a more careful comparison is necessary. We will study in the following the possible phase transitions that we may have since it may lead to an interesting phenomenology of the Skyrme crystals. For simplicity, since L is a measure of the size of unit cell it is also a measure of the baryon density, then we will also refer to the region of small values of L as the high density regime and for large values of L the low density regime.

Low density phase transition
As we noted in the construction of the crystals, the FCC + crystal will always have larger or equal energy than the FCC. In fig. 8 we see that the energies of both crystals are indistinguishable at high densities, but at some value of the length the curves split and the FCC crystal becomes the ground state. This behaviour in the E(L) curves suggests a phase transition between the crystals, but the presence of the pion-mass potential term is crucial in the understanding of this possible transition. Concretely, this potential term  explicitly breaks chiral invariance, then it is not compatible with the FCC + symmetry σ → −σ. However, the relevance of the potential term in the energy decreases at high densities, so both crystals tend to the same energy in the chiral limit. Hence when this term is not included in the lagrangian both crystals are allowed and we find a FCC to FCC + second order phase transition, but when the potential term is present the FCC crystal is always the ground state and the energy curves approach asymptotically. This phase transition has been extensively studied in [54], where the σ field was proposed to be the order parameter of the transition.  Figure 9: The mean value of the σ field within the unit cell shows the second order transition from the FCC to the FCC + crystal when the pion-mass term is not included, and the asymptotic approach when it is included.
We show in fig. 9 the mean value of the σ field in the unit cell for the cases L 24 and L 240 since they represent the cases without and with pion-mass potential term respectively. The addition of the sextic term does not qualitatively change the curves.
Although the FCC + is not the ground state crystal it is a good approximation to the FCC crystal at large densities. Indeed, for the values of the parameters that we have considered, the transition point (in the case without potential term) always occurs at densities smaller than the minimum of energy, and even with potential term the FCC + crystal is already a good approximation to the FCC crystal.

High density phase transition
Now we want to compare the energy curves between the BCC and the FCC crystals. An important point here is that whilst the FCC unit cell contains 4 baryons, the BCC unit cell has 8 baryon units. If we want to compare both crystals we need to do it at the same baryon density, which may be easily defined, Hence, if we want to compare the energies we may calculate the density of both crystals and find the point at which the BCC crystal is more energetically favourable than the FCC crystal. We find that the different terms that we consider in the lagrangian have an important impact on the transition point. Concretely, the sextic term locates the transition at physically reasonable densities, i.e. the same order of magnitude as the density at the energy minimum. Without the sextic term we find the transition point at very high densities, and the addition of the pion-mass term shifts the transition density to even higher values, therefore we only plot the cases in which we have the sextic term, see fig. 10. Table 5: Ratio between the transition density and the density at which the minimum of the energy is achieved for the FCC + crystal.
Since the energy curves have different slopes at the transition point there is a discontinuity in the derivative of the energy. This implies that the FCC-to-BCC is a first order phase transition and we must perform a Maxwell construction (MC) in order to avoid unphysical regions.
The pressure of a system acquires great relevance when phase transitions are present since it must remain finite and continuous in order to have a physical transition. The pressure, as well as the energy density of the crystal, can be obtained from their thermodynamical definition, From these expressions we may conclude that there is a discontinuity in the pressure of the crystal and this is contradiction with the Gibbs conditions that must be preserved in every phase transition, For this analysis we will identify the FCC crystal as phase 1 and the BCC crystal as phase 2. Besides, in our system the baryon charge is conserved so we must find the mixed phase which has the associated chemical potential (µ B ) common to both phases.
The MC introduces a mixed phase which preserves the Gibbs conditions in this case. The main idea of this construction is to find one point in each of the energy curves which have the same pressure, we denote it by p P T , and join them with the curve which has the same value of µ B for the two phases. Mathematically this means that we have to find the points (V 1 , V 2 ) of each phase that have the same slope in the E(V ) diagram and are both tangent to the straight line with the same slope (which is p P T ).
We must be careful for this calculation since we are dealing now with the volumes of the unit cells. This means that same volumes have different baryon content in each crystal, so we need to rescale them in order to have the same baryon number, The final energy curve with a physical phase transition starts at low densities in the FCC crystal until we find the mixed phase which joins to the BCC crystal, fig. 11. We may use the FCC + energy curve for these calculations since the density at which the transition occurs the FCC and FCC + crystals are the same in the L 246 case and the difference between them in the L 2460 case is negligible.

Fluid-like transition
From the previous calculations we know how to calculate the most important thermodynamical magnitudes (pressure, energy and baryon densities) for the Skyrme crystal. We also know that decreasing the size of the unit cell increases the densities as well as the pressure, hence we may try to calculate how homogeneous the crystal becomes for increasing density and if a transition to a fluid is possible. Indeed, it is known that the BPS Skyrme submodel (L 60 ) describes a perfect fluid due to the properties of the sextic term. Then it seems reasonable that, since the sextic term is the most important one at small values of L, the crystal becomes more homogeneous within the unit cell. To study the degree of homogeneity we will compare the energy density profiles obtained from the numerical minimization with a constant density profile with the value of the mean energy density of the unit cell, At the minimum of the energy we expect to have a highly inhomogeneous crystal, where the skyrmions are surrounded by regions of vacuum, and reducing the size of the unit cell will decrease the inhomogeneity. However, we also expect to increase this effect with the addition of the sextic term, so we will compare the L 240 and L 2460 cases, since we want to consider the more realistic cases in which pions have mass. To compare the energy densities within the unit cell we define the radial energy profile (REP) enclosed within a sphere of radius r, where ρ is the ρ mean in the case of the constant energy density unit cell and the integrand of eq. (1.10) in the real case. We calculate both REP for each case at the baryon density of the minimal energy (n 0 ), and at the higher densities 3n 0 and 7n 0 . We show in fig. 12 the ratio χ = E(r)/E mean (r) of the two REP for each case. Figure 12: REP for the L 240 (blue) and L 2460 (red) cases, for n 0 (continuous lines), 3n 0 (dashed lines), and 7n 0 (dashed-dotted lines). The degree of homogeneity almost does not change for L 240 , whereas the energy density becomes much more homogeneous with increasing baryon density for L 2460 .

New lattice solutions
For the moment, we have mostly focused on the behaviour of Skyrme crystals in the region to the left of the minimum of the E(L) curve. The reason is that, from eq. (2.15) we may observe that the minimum corresponds to the point p = 0, and the region L ≥ L 0 has negative pressure, hence it is unstable. It is the aim of this section to show that there is a new branch of solutions which have different energies in the low density regime and tend to the FCC + crystal at high energies.
The fact that the Skyrme crystal has a minimum is not a bad behaviour, since this is expected to occur in symmetric nuclear matter. However, the energy of the crystals seems to diverge with L, but this is due to the Fourier expansion that we use to construct the Skyrme crystal, in which we impose the skyrmions to be in fixed positions and we do not allow them to move freely within the unit cell to find the lowest energy configuration. This is a correct procedure for small values of L, however if we increase the size of the unit cell the skyrmion can only spread instead of clustering to form a compact configuration surrounded by vacuum.
This motivated us to find new lower energy configurations with a new numerical minimization method which lets the skyrmions move freely within the unit cell. We use a gradient flow method to find the field configurations with minimal energy, locating a B = 4 skyrmion in the centre of the unit cell. The motivation for this new lattice starts with the similarities between the isolated B = 4 skyrmion, which has cubic symmetry, and the FCC + symmetry. Indeed the B = 4 skyrmion is quite similar to the crystal in the sense that it is composed by eight half-skyrmions located in the corners of a cube.
Besides, the study of the B = 4 skyrmion in periodic boundary conditions under different deformations showed the phase transitions it may experiment [55]. Concretely, the phase transition between the FCC + and the B = 4 skyrmion lattice was found at a certain value of L, then the new lattice becomes a more energetically favourable crystal. Since the isolated B = 4 skyrmion aims to describe an alpha particle, we will refer to this configuration as the α-lattice.
We calculated the energy for the α-lattice at different values of L for the L 240 and L 2460 models. This new lattice has lower energy than the crystal at a certain value of L and, furthermore, it tends to a constant value at L → ∞. Indeed it tends to the value of the isolated B = 4 skyrmion, so we may construct other cubic lattices with larger values of B, since it is know that the energy per baryon of skyrmions decreases for increasing values of the baryon charge. We show the energy of the next simplest cubic lattice, which is a multiple of the α−lattice, the B = 32 lattice.  Figure 13: Energy per baryon number for the different lattices that we consider with and without sextic term. These α and B = 32 lattices are important in the low density regime where they have much lower energies than the standard crystal. However, already for densities which are slightly smaller than at the minimum of the energy all these lattices tend to the same FCC crystal.
As expected, both the α and B = 32 lattices have less energy than the FCC crystal at low densities, since they achieve a more compact configuration surrounded by vacuum. Decreasing the size of the unit cell forces the skyrmion within the unit cell to recover the FCC + crystal configuration. The transition for both lattices may be seen in fig. 13. In fig. 14 the corresponding energy density contours are shown for the L 240 model.

Isospin quantization and symmetry energy
The classical cristalline solutions presented in the previous sections can be understood as models for infinite, isospin symmetric nuclear matter, that is, with the same number of protons and neutrons. However, it is well known that nuclear matter in the interior of neutron stars cannot be completely isospin symmetric, with all but a small fraction of the total baryonic degrees of freedom being protons at a given density. The magnitude that determines the fraction of protons over the total baryon number is the so-called symmetry energy, namely, the change in the binding energy of the system as the neutron-to-proton ratio is changed at a fixed value of the total baryon number, and its knowledge is essential to determine the composition of nuclear matter at high densities.
For an finite nuclear system with total baryon number A = N + Z, with N (Z) the number of neutrons (protons), the isospin asymmetry parameter is defined as δ = (N − Z)/A = (1 − 2γ), with γ = Z/A the proton fraction. For an infinite system, the above quantities can be defined with densities instead of total numbers. The binding energy of infinite nuclear matter is thus parametrized as a function of both the baryon density and the asymmetry parameter, with E N (n B ) the binding energy of isospin-symmetric matter, and S N (n B ) the symmetry energy. Although its dependence on the density has proven difficult to measure experimentally, it is usually parametrized as an expansion in powers of the baryon density around nuclear saturation n 0 , with = (n − n 0 )/n 0 , and the slope and curvature of the symmetry energy at saturation, respectively.The symmetry energy at saturation is well constrained (S 0 ∼ 30 MeV) by nuclear experiments [56], but the values of the slope and higher order coefficients are still very uncertain. However, recent efforts on the analysis of up to date combined astrophysical and nuclear observations have allowed to constrain the value of these quantities with reasonable uncertainty above nuclear saturation [57][58][59][60][61]. In the Skyrme model, due to the SU (2) isospin symmetry of the Lagrangian, isospin degrees of freedom represent zero-modes, which are quantized using standard canonical quantization in terms of some collective coordinate parametrization (see, e.g., [22,24,62,63]).
Following this approach, we consider a (time-dependent) isospin transformation of a static Skyrme field configuration, (3.4) and ω . = g †ġ = i 2 ω a τ a is the isospin angular velocity. The time dependence of the new Skyrme field induces a kinetic term in the energy functional, given by 3 where Λ ij is the isospin inertia tensor, given by (3.7) Due to the symmetries of the cristalline phases, the complete isospin inertia tensor for the unit cell of a cubic crystal turns out to be proportional to the identity, and the associated eigenvalue (the isospin moment of inertia) can be written where the contribution from each term in the lagrangian, denoted by Λ (n) , can be found in [48]. This fact enormously simplifies the kinetic term in the Lagrangian of an isospinning cubic crystal with a number N cells ≡ N of unit cells, which is reduced to and, by defining the corresponding canonical momentum J a = ∂L/∂ω a = N Λω a , we may write it in Hamiltonian form, Now, following the standard canonical quantization procedure, we promote the isospin angular momentum variables to operators, so that we may diagonalise the Hamiltonian in a basis of eigenstates with a definite value of the total isospin angular momentum, The total isospin angular momentum of the full crystal will be given by the product of the total number of unit cells times the total isospin of each unit cell, which can be obtained by composing the isospin of each of the cells. In the charge neutral case, all cells will have the highest possible value of isospin angular momentum, so that in each unit cell with baryon number B cell , the total isospin will be 1 2 B cell , and hence the total isospin of the full crystal will be J tot = 1 2 N B cell . Therefore, the quantum correction to the energy (per unit cell) in the charge neutral case (completely asymmetric matter) due to the isospin degrees of freedom will be given by (assuming N → ∞) Such correction could also have been obtained directly by introducing an 'external' isospin chemical potential µ I , and promoting the regular derivatives in the Skyrme Lagrangian to covariant derivatives of the form [64] so that, if U is a static configuration, the time component of the Maurer-Cartan form becomes This expression is equivalent to that of an iso-rotating field with angular velocity ω a = −µ I δ 3a . Thus, it is straightforward to obtain the isospin chemical potential for the Skyrmion crystal using its thermodynamical definition µ I = − ∂E ∂n I , where n I is the (third component of) the isospin number density. Given that (J tot ) 2 = J 2 1 + J 2 2 + J 2 3 and n I = J 3 /N , we may write the isospin energy per unit cell as and then For B < ∞ Skyrmions, isospin quantization (together with spin) induces a set of constraints (Finkelstein-Rubinstein constraints [10,65]) which characterize the allowed states with a given value of total spin and isospin angular momentum, and thus the ground states and lowest energy spin and isospin excitations have been studied for many Skyrmions with B ≤ 12 [24,63]. However, in the case of a crystal, to compute the contribution of the quantization of the global isospin zero modes to the total energy we would need to know, in principle, the quantum isospin state of the whole crystal. This is of course impossible in the thermodynamic limit, and some additional approximation becomes necessary.
Following [48], since the total third component of isospin is a good quantum number for the total quantum state of the crystal, we perform a mean field approximation and consider that the isospin density in an arbitrary skyrmion crystal quantum state is approximately uniform so that where n I is the effective isospin charge per unit cell in this arbitrary quantum state. The (mean field) effective proton fraction corresponding to such an isospin charge per unit cell with baryon number B cell is given by Hence, the isospin energy per unit cell of the Skyrmion crystal can be written in terms of the asymmetry parameter from where we can readily read off the symmetry energy for Skyrme crystals Therefore, the symmetry energy of a Skyrmion crystal is uniquely determined by the (isospin) moment of inertia of each unit cell, which has an intrinsic dependence on the density. On the other hand, the isospin energy correction depends, in the mean-field approximation, both on Λ (hence the density) and the asymmetry parameter (or equivalently, the proton fraction). However, a nonzero proton fraction would rapidly lead to a divergence in the Coulomb energy of the total crystal in the infinite crystal limit. Therefore, a neutralizing background of negatively charged leptons (electrons and possibly muons) must be included in the system, such that the Coulomb forces are screened. The total system is thus characterized at equilibrium by two conditions, namely the charge neutrality condition and the β-equilibrium condition i.e., the isospin chemical potential must equal that of charged leptons, which in turn means that the direct and inverse beta decay processes such as n ↔ p + l +ν l take place at the same rate. Leptons inside a neutron star can be described as a non-interacting, highly degenerate fermi gas, so that the chemical potential for each type of leptons can be written where k F = (3π 2 n l ) 1/3 is the corresponding Fermi momentum, and m l is the mass of the corresponding lepton. For sufficiently high densities, then, the electron chemical potential becomes larger than the mass of the muon, µ e ≥ m µ , and the production of muons is preferred by the system. We can estimate the total proton fraction by enforcing both beta equilibrium and charge neutrality. Neglecting the contribution of muons in a first step, to the charge density, from (3.21) we can relate the electron density to the proton fraction parameter, n e = γB cell /(2L) 3 . The β equilibrium condition then provides an equation which defines γ implicitly as a function of the lattice length parameter, where we also assumed the ultrarelativistic electron approximation, i.e. m l /k F 0. Including the muon contribution to the charge density yields a slightly more complicated expression for the β-equilibrium condition, given by where (3.26) The proton fraction inside beta-equilibrated matter determines, in addition, whether a proto-neutron star will go through a cooling phase via the emission of neutrinos through the direct Urca (DU) process n → p + e +ν e . This process is expected to occur if the proton fraction reaches a critical value, γ p > x DU , the so-called DU-threshold [66,67]. The DU process allows for an enhanced cooling rate of NS. Whether it takes place or not in the hot core of proto-neutron stars or during the merger of binary NS systems [68], therefore, determines the proton fraction (and the symmetry energy) of matter at ultrahigh densities. It is, however, not clear whether this enhanced cooling actually occurs, although there is recent evidence that supports it [69].
In both cases, a persistent population of protons and leptons at higher nucleon density is expected, although we can see that in the case with the sextic term the fraction of charged particles is smaller. This is explained by the lower symmetry energy of the sextic term, which makes it much easier to convert protons into neutrons. Finally, for all values of the parameters (f π , e, λ 2 ) we considered, the DU-threshold is not reached. One should, however, not consider this fact as a prediction of the Skyrme model, because it strongly depends on the parameter values. Besides, it is also generally assumed that around 2 −  Figure 15: Fraction density γ i for each particle as a function of the baryon density for λ 2 = 0 (solid) and λ 2 = 1.5 (dashed). The corresponding DU threshold is also shown in black. This figure was originally published in [48].
above, additional degrees of freedom are expected to become relevant for the physics of dense nuclear matter. In particular, strange mesons (kaons) and hyperons are believed to modify the EOS at sufficiently high density, which motivates the question of how to describe an additional quark flavor using the Skyrme model approach.
Below we shall briefly review the crucial steps to determine whether kaon fields may condense inside a Skyrme crystal for a sufficiently high density, and, if so, whether this critical density value is relevant for the description of matter inside compact stars. A more detailed discussion can be found in [49].

The kaon condensate effective potential
Following the bound-state approach first proposed in [70] we may extend the skyrmion field to a SU (3)-valued field U and consider kaon fluctuations on top of a SU (2) skyrmionlike background u. With the only requirement that unitarity must be preserved, different ansätze have been proposed in the literature for the total SU (3) field describing both pions and kaons. In this work, we choose the ansatz proposed by Blom et al in [71]: In this ansatz U π represents the SU (3) embedding of the purely pionic part u, and the field U K are the fluctuations in the strange directions. It can be shown that this ansatz is equivalent to the one first proposed by Callan and Klebanov in [70] when computing static properties of hyperons, although both may differ in other predictions of the model [72].
In the simplest SU (3) embedding, the SU (2) field u is extended to U π by filling the rest of entries with ones in the diagonal and zeros outside. On the other hand, the kaon ansatz is modelled by a su(3) -valued matrix D which is non trivial in the off-diagonal elements: where K consists of a scalar doublet of complex fields representing charged and neutral kaons: We now extend the Generalized Skyrme Lagrangian (1.8) to include strange degrees of freedom. First, we need to replace the potential term L 0 in order to correctly take into account the mass of the kaon fields. The new potential term has the form [72]: where λ 8 is the eighth Gell-Mann matrix and m K is the vacuum kaon mass. In addition, the effects of the chiral anomaly must be taken into account in the effective theory by means of the Wess-Zumino-Witten (WZW) term, which can be expressed in terms of a 5-dimensional action: The onset of kaon condensation in the Skyrme model takes place at a critical density n cond at which µ e becomes greater than the energy of the kaon zero-momentum mode (s-wave condensate). Thus, for baryon densities n ≥ n cond , the presence of kaons will be more energetically favorable than electrons at the fermi surface, so kaons will be produced at expense of electrons, and the macroscopic contribution of the kaon condensate to the energy must be taken into account when obtaining the EOS. To do so, we write the field condensates (i.e. the non-zero vacuum expectation values (vev), K ± ) as a homogeneous field whose time dependence is given by: The real constant φ corresponds to the zero-momentum component of the fields, which acquires a nonvanishing, macroscopic value after the condensation. Its exact value is determined from the minimization of the corresponding effective potential, which is determined from the Skyrme lagrangian. On the other hand, the phase µ K is nothing but the corresponding kaon chemical potential. First, we will need an explicit form of the SU (3) Skyrme field in the kaon condensed phase. Assuming the charged kaons will be the first mesons to condense, we drop the neutral kaon contribution and define the following matrix which results from substituting the kaon fields in D as defined in (4.2) by their corresponding vev in the kaon condensed phase. Also, taking advantage of the property D 3 = φ 2 D, we may write the SU (3) element generated byD explicitly in matrix form: fπ φ is the dimensionless condensate amplitude. Furthermore, assuming the backreaction from the kaon condensate to the skyrmion crystal is negligible, and thus the classically obtained crystal configuration will be the physically correct background even in the kaon condensed phase, we may write the SU (3) field in this phase as U = ΣU π Σ, where U π is the SU (3) embedding of the SU (2) skyrmion background as in (4.2). Introducing this U in the total action yields the standard Skyrme action for the SU (2) field plus an effective potential term for the kaon condensate: (4.14)

Quantum corrections to the effective potential
In the above calculations, we have taken separately the contributions of a kaon condensate and an isospin angular momentum of the skyrmion crystal. However, since kaons possess an isospin quantum number, the condensate interacts with the skyrmion isospin both indirectly via the charge neutrality and β equilibrium conditions, which relate their corresponding chemical potentials, and also directly, due to the appearance of additional contributions to the total energy when we consider a (time-dependent) isospin transformation of the full Skyrme field + kaon condensate configuration U = ΣU π Σ: where A is an element of SU (3) modelling an isospin rotation, (2). (4.16) in analogy with the SU (2) case, we define the isospin angular velocity ω as A †Ȧ = i 2 ω a λ a (a = 1, 2, 3), with λ a the Gell-Mann matrices generating SU (3) for a = 1, · · · 8. Notice that ω is a three-vector, since A †Ȧ belongs to the isospin su(2) subalgebra of su(3). Then, we may write the time component of the Maurer-Cartan current asŨ † ∂ 0Ũ = AL 0 A † + AT a A † ω a , where T a is the su(3)-valued current: where we have made use of the parametrization (4.2). Now the kinetic isorotational energy in this case takes the form where Λ ab is the isospin inertia tensor and ∆ a is the kaon condensate isospin current, given by where a, b and c are those in eq. (3.7). As argued previously, the symmetries of the crystalline configuration imply that the isospin inertia tensor becomes proportional to the identity, i.e. Λ crystal ab = Λδ ab . However, the presence of a kaon condensate breaks this symmetry to a U (1) subgroup, so that Λ ab presents two different eigenvalues in the condensate phase, Λ cond = diag(Λ, Λ, Λ 3 ). Similarly, ∆ a = 0 in the purely barionic phase, and its third component acquires a non-zero value in the condensate phase, ∆ cond = (0, 0, ∆). The explicit expressions for Λ 3 and ∆ in the condensed phase can be found in [49]. One can check that in the non-condensed phase, φ = 0 and the results of the previous section are recovered, namely, Λ 3 = Λ, ∆ = 0.
The quantization procedure now goes along the same lines as in the previous section. However, the isospin breaking due to the kaon condensate implies that the canonical momentum associated to the third component of the isospin angular velocity will now be different, and given by I 3 = Λ 3 ω 3 + ∆.
Thus, after a Legendre transformation to rewrite (4.18) in Hamiltonian form, and making the N → ∞ approximation, one can write the quantum energy correction per unit cell of the crystal in the kaon condensed phase as The first term on the rhs is just the isospin correction, while now there is an additional second term due to the isospin of the kaons. Indeed, since the kaon field enters also in the expression of the isospin moment of inertia Λ 3 , both terms will depend nontrivially on the kaon vev field. When the kaon field develops a nonzero vev, apart from the neutron decay and lepton capture processes, additional processes involving kaons may occur: such that the chemical equilibrium conditions µ n = µ p + µ K , µ l = µ K are satisfied. These are the extension of eq. (3.22) to the condensate phase. The total energy within the unit cell may be obtained as the sum of the baryon, lepton and kaon contributions: The kaon contribution is the effective potential energy which depends on the condensateφ and on the lepton chemical potential through the explicit dependence on µ K of both V K and ∆, and µ K = µ e due to the equilibrium conditions. Therefore, the energy of the full system depends on the proton fraction, the kaon vev field and the electron chemical potential. Their respective values can be obtained, for fixed n B (or equivalently, fixed L) by minimizing the free energy with respect to γ,φ and µ e , i.e.
The first equation imposes the expected condition µ e = µ I = 2 2 (1 − 2γ)/Λ 3 . Then, after substituting into the other two conditions we get: which are precisely the charge neutrality condition, and the minimization of the grand canonical potential with respect to the kaon field. We note here that we drop the ultrarrelativistic consideration for electrons since the appearance of kaons may decrease hugely the electron fraction. By solving the system of equations 4.27 and(4.28) for γ andφ we obtain all the needed information for the new kaon condensed phase. Then we may compare the particle fractions and energies between both phases, which we will call npeµ and npeµK. Before solving the full system for different values of the lattice length L, we may try to obtain the value of the length at which kaons condense, L cond . This value is indeed important since it will determine whether or not a condensate of kaons will appear at some point in the interior of NS. This is accomplished with the same system of eqs. (4.27) and (4.28) by factoring the sinφ from the second equation and settingφ = 0. Then we may see the system as a pair of equations to obtain the values of γ cond and L cond , the values of the proton fraction and the length parameter for which the kaons condense.
We show in the table below the density at which kaons condense for different values of the parameters as well as the values of some nuclear observables they yield. All the values are given in units of MeV or fm, respectively.
Parameter sets 1 and 2 are chosen so that the energy per baryon and baryon density at saturation are fitted to experimental values, whereas the sets 3 and 4 correctly fit the symmetry energy and slope at saturation. In section 5.4 we shall always use set 1 of parameter values.

TOV system of equations
In order to calculate the mass and radius for a non-rotating NS we have to solve the standard TOV (Tolman-Oppenheimer-Volkoff) system of ODEs. It is obtained inserting label f π e λ 2 E 0 n 0 S 0 L sym n cond /n 0 set 1 133. 71 in the Einstein equations, To describe matter inside the star, in the right-hand side of the equation, it is standard to use the stress-energy tensor of a perfect fluid, where the pressure p and the energy density ρ are not independent but related by the EOS. Hence the EOS describes the nuclear interactions inside the NS and different EOS lead to different observables. The resulting TOV system involves 3 differential equations for A, B and p, which must be solved for a given value of the pressure in the centre of the NS (p(r = 0) = p 0 ) until the condition p(r = R) = 0 is achieved.
We use a 4 th order Runge-Kutta method of step ∆r = 1 m to solve the system in order to obtain the main observables from the solutions. The radial point at which the pressure vanishes defines the radius of the NS, and the mass M is obtained from the Schwarzschild metric definition outside the star, .

(5.7)
For different values of the pressure in the centre we may obtain different couples of values for the mass and radius, then we can represent the Mass-Radius (MR) curve which is the most important result in this static problem. Different EOS may lead to very different MR curves, then the measurement of these observables can constrain the curve and so the EOS. Indeed, the region of low mass values in the MR curve has been tightly constrained from nuclear experiments at low densities.
The TOV system can be generalized via the Hartle-Thorne perturbative formalism to consider rotating NS in a slow rotation approximation. In that case more equations are added to the TOV system hence a new set of observables like the moment of inertia, rotational Love number and quadrupolar moment may be extracted. These observables are typically obtained from isolated pulsars, hence the importance in the extension of the TOV system to the Hartle-Thorne formalism lies in the appearance of a new source of information for NS to restrict the nuclear EOS. Furthermore it is possible to extract from this formalism an additional equation which describes how a compact object is deformed by an external tidal force. The associated observable is the tidal deformability which can be calculated from the GW spectrum emitted in NS collisions.
In this work we will only consider the MR curves and the tidal deformability to compare with recent observations.

A generalized Skyrme EOS
In this section we want to briefly review the construction of a Skyrme-model based EOS first discussed in [50], based on the following two observations. Firstly, neutron stars based on the standard L 24 model lead to too small maximum NS masses [36], whereas the EOS based on the L 60 submodel imply rather large maximum masses [73]. Secondly, the sextic term provides the leading contribution to the EOS at large densities [35], whereas it must be subleading at lower densities because of its scaling properties. This motivates to consider a generalized Skyrme model EOS in the study of NS which interpolates between the standard L 24 crystal EOS at intermediate densities and the L 60 submodel EOS at large densities, as was done in [50]. Then we will compare these results with the full numerical generalized Skyrme crystal solutions.
Even though we have not compared the energy of the crystals with other skyrmion solutions, it is known that the crystalline configurations are, so far, the ones with the lowest energy in the standard (L 24 ) Skyrme model. In one of the works where the FCC + crystal was found [42] a parametrization of the energy in terms of the unit cell side length, similar to eq. (2.11), is shown where a = 0.474 and b = 0.0515 are dimensionless parameters, and E 0 , l 0 are, respectively, the energy and length scales of the minimum in the E(L) curve. The relation between the lattice length used in [42] and the one used in this work is l = 2L. Since the standard Skyrme model is relevant at lower energies, we want to identify the minimum of the crystal with the saturation point of infinite nuclear matter, E 0 = 923.3 MeV, l −3 0 = n 0 = 0.16 fm −3 . The change in the values with respect to the ones originally used may be seen as a redefinition for f π and e. From eq. (5.8) we may obtain the pressure and the energy density using eq. (2.15), which vanishes at the minimum l 0 . On the other hand, for high densities the sextic term is the most relevant one since the EOS yielded by this term is maximally stiff (ρ 6 = p), with a speed of sound (c 2 s = ∂p/∂ρ) equal to 1. The BPS Skyrme model, even with a general potential term (with no derivatives of U ), L 0 = −µ 2 U, retains this property, besides its stress-energy tensor is of the perfect fluid form from which we may easily obtain the EOS.
Additionally, since the topological current B 0 may be identified with the baryon density n B in the Skyrme model, a full thermodynamical relation between p, ρ and n B can be obtained p + ρ = 2λ 2 π 4 n 2 B . (5.12) The dependence on the Skyrme field of the potential U may lead to a non-barotropic EOS, but the idea to finally join both submodels is to take a constant potential which will take into account the residual contributions when we are at high densities. Then, we define the EOS at high densities as, ρ = ρ 6 + ρ 0 = p + const. (5.13) Besides the different properties that each submodel has, some attempts to describe NS also motivate a generalized model from the combination of these two submodels to describe nuclear matter from the saturation point n 0 to very high densities (n ∼ 10n 0 ). Skyrmion crystals described by eq. (5.8) with the original parameters were used to calculate NS observables [74], however the masses turn out to be too small (M max ≈ 1.5M ) compared with the most recent measurements of NS (M max ≥ 2M ). Later, the BPS Skyrme model was also coupled to gravity [73], using different potentials, and the NS masses were too large (M max ≈ 3−4M ) compared to observations. These results are themselves an evident motivation to combine the submodels and study the NS properties, since we intuitively may think that a generalized model would lead to intermedium NS masses, which are precisely in the region of observed masses (M ∼ 2 − 2.5M ).
In [27] the parameters of the BPS submodel (λ 2 , µ 2 ) were fitted to reproduce the nuclear saturation point. However from the previous discussions we use the parameters of the standard Skyrme model to fit the minimum of the crystal to these values and construct a generalized Skyrme EOS from the asymptotic values, ρ Gen (p p P T ) = ρ Sk , ρ Gen (p p P T ) = p + const. (5.14) The EOS ρ Sk (p) may be obtained from E 24 (l)/l 3 and the relation between l and p 24 in eq. (5.9). A simple construction, continuous and with a smooth transition between the two descriptions is an interpolation of the form, ρ Gen (p) =(1 − α(p, p P T , β))ρ Sk + α(p, p P T , β)(p + ρ Sk (p P T )), (5.15) α(p, p P T , β) = (p/p P T ) β 1 + (p/p P T ) β .
For this interpolation, the value of β measures how fast the transition is. Here we take the value β = 0.9 since depending on the value of p P T a larger value of β would lead to a speed of sound larger than 1 in some regions inside the star. We will constrain the range of values for p P T comparing the resulting maximum masses with the observations.
Although we have constructed a general EOS which is able to describe matter inside NS it is still not complete. We have seen that the low density model (the Skyrme crystal) is only valid above nuclear saturation due to the presence of a minimum. For densities below saturation, finite-size and electromagnetic effects become important so we need to take them into account. The modelling of very large B nuclei or the coupling to electromagnetic interactions is possible within the Skyrme model, however, how nuclear matter behaves below saturation is well known using standard methods like many-body calculations. Then for our purpose we may use one of these models to describe the low density regimes inside NS. Indeed, we will use the low-density model to describe nuclear matter until a certain pressure, denoted by p * and then we glue the generalized model eq. (5.15), hence we are adding a crust to our generalized EOS which will mainly affect NS with low masses. We take the BCPM equation of state [75] as the model for low densities and perform a similar interpolation as before between this model and eq. (5.15). For the same α(p) function, taking β = 2 and replacing p P T by p * , we called the combination of the two models the Hybrid EOS, ρ Hyb = (1 − α(p, p * , 2))ρ BCPM + α(p, p * , 2)ρ Gen . (5.16) The resulting EOS, as well as some standard nuclear physics EOS, are shown in fig. 16. Figure 16: Energy and baryon density of the standard Skyrme crystal, the BPS model, the Generalized model with p P T = 50 MeV/fm 3 and the Hybrid model, with the same p P T and p * = 2 MeV/fm 3 . We also show some usually considered EOS: APR4, WWF1, BCPM and we fill in orange and blue the regions of allowed values for p P T and p * . This figure was originally published in [50] We sample the space of parameters (p P T , p * ) and restrict their values solving the TOV equations to compare the mass and radius results with some pulsar and GW measurements. We plot in fig. 17 some representative results of the MR curves to show the accurate agreement between the generalized Skyrme EOS and observations. We use the constraints of maximum mass M/M = 2.16 +0. 17 −0.15 from the combined analysis of the GW170817 event and quasi-universal relations from NS in [76], which constraints the medium mass (∼ 1.4M ) region, also the heavier pulsar observations PSR J1614 -2230 (1.928 ± 0.017M ) [77], PSR J0348 + 0432 (2.01 ± 0.04M ) [78], PSR J0740 + 6620 (2.14 +0. 10 −0.09 M ) [79], and some additional constraints from NICER, chiral EFT and multimessenger observations [80]. As expected, it is found that p P T affects the maximum mass value, while p * mostly determines the radius of the medium mass range. The allowed ranges for the two transition pressures in order to satisfy the constraints are p P T ∈ [25,50] MeV/fm 3 and p * ∈ [0.5, 2] MeV/fm 3 . We even allowed for maximal masses ∼ 2.7M much larger than the observations, mainly motivated by the [81] event in which a black hole and a secondary compact object falls in the mass gap (2.5 − 2.7M ) [82], with the possibility of being the heaviest stable NS ever observed. Besides we want to remark that we may easily accommodate these very large values of NS masses keeping reasonable values of the radii in the medium mass region, which is not an easy task in general for other nuclear models.  [50] It is also possible to extend the calculations and solve an additional equation to obtain the tidal deformability for a non-rotating NS [50]. The importance of this magnitude is that it can be extracted from the waveform of the early phase inspiralling coalescence of two NS. Indeed what can be measured is a mass weighted averaged of the two individual deformabilities which is called effective tidal deformability Λ. Similarly, the individual masses are not directly measured as well, but a magnitude called the chirp mass, M c = m 1 q 3/5 /(1+q) 1/5 , with q = m 1 /m 2 the ratio between the two masses, has been constrained at the 90% confidence level in the GW170817 event [83]. In fig. 18 we show the effective tidal deformability obtained from our Hybrid model with the estimations extracted from [84] of the effective tidal deformability for the measured values M c = 1.188 +0.004 −0.002 and q = 0.7 − 1. The estimated value of Λ from the GW170817 event is rather low (smaller than 800) and it is a new source to eliminate some EOS which may satisfy the mass-radius constraints but do not have a correct value of the deformability.

Skyrme crystal EOS
The results from the last section provide a strong motivation to use the full crystal solutions within the generalized Skyrme model, studied in the first section, to compute NS observables. Besides being solutions of the generalized Skyrme model, we have seen how the crystals allow to include isospin effects which furthermore led to the addition of leptons, leading to a much more realistic description of nuclear matter at high densities. However, Figure 18: The effective tidal deformability Λ of the Hybrid model against the mass ratio. The orange regions show the PDFs obtained from [84], and we adopt the notation Hyb p P T p * for the curves. This figure was originally published in [50]. these solutions and their properties depend on the values of the parameters that appear in the model, which we have set to some standard values. Then, in order to be acceptable candidates to describe realistic nuclear matter they must reproduce observables from nuclear physics and NS. In this section we will study how to fit the parameters to reproduce different nuclear observables using Skyrme crystals and, finally, we will extract the most important NS observables and compare them with the most fiducial measurements.

Scan of the parameters
A remarkable property of the standard Skyrme model without potential term (namely, the quadratic and quartic terms only) is that, with a suitable choice of units, one can factor out all dimensionful constants from the energy functional, so that the constants remaining inside it are just dimensionless numbers. This can be seen, for example, in eq. (2.11). As a consequence, one can just forget about the numerical values of the coupling constants and numerically find the different crystalline solutions at different unit cell lengths. Once the relation between the (adimensional) energy and length E(L) is found, the values of the coupling constants can be adjusted a posteriori in order to fit whatever observable we are interested in. Unfortunately, the addition of the sextic and mass potential terms to the energy functional spoils this property, as there is no choice of units that allows to factor out all coupling constants in the energy functional. This means that, in order to be able to obtain solutions, one needs first to give specific values to (some of) the coupling constants appearing in the problem. However, this is problematic if one wants to fit the values of energy and density to some physical values. For example we want to identify the energy and the density at the minimum of the crystal (L 0 ) with the nuclear saturation point, however, one cannot know the value of L 0 without performing the numerical simulations, but in order to do so you need to fix the values of the parameters.
Hence, the fitting of the skyrmion crystal parameters in the generalized model to values at nuclear saturation is, in principle, a very difficult problem that needs to be solved iteratively until a self consistent solution is achieved. Given the computational cost of simulating a single unit cell at a given length, following this naive approach would make it almost impossible to realize a significant scan of parameters in a reasonable time.
However, from the results of our numerical simulations we have observed that the FCC + crystal displays an almost perfect scaling property [45] with the unit cell length. This property is stronger than the fact that each E(L) curve may be fitted by eq. (2.11),  Figure 19: Energy and isospin moment of inertia of a skyrmion crystal as obtained via the full numerical minimization (dots) and the perfect scaling fit (solid curves), for the models with and without sextic term.
it means that each of the terms in the energy functional scales with L independently as E i ∝ L −i+3 , where i is again the number of spatial derivatives appearing in that particular term. The values of the perfect scaling constants, which we label by K i , are then universal in the sense that they will not change for different values of the parameters. This perfect scaling property is a characteristic of the field configuration, and not only of its energy. For example, in fig. 19 it is observed that also the isospin moment of inertia of the crystal unit cell displays a sufficiently well perfect scaling. This is important in our analysis since we also want to fit the values of the symmetry energy and its derivatives. Although the scaling is not perfect, in general the biggest deviations from the full numerical values of energy start far from the minimum, at which the perfect scaling fit is most precise. Thus, we take advantage of this property in order to fit the magnitudes obtained from generalized Skyrme crystals to their physical values (up to a certain error). We do so following an iterative process based on five main steps: 1. Due to the perfect scaling property, the (adimensional) energy and isospin moment of inertia approximately satisfy the following expressions,  Table 7: Perfect scaling parametrization constants 2. We fix the energy scale E s to a desired value in MeV. This is equivalent to fixing one of the three free parameters of the model, for instance, f π .
3. Then, we calculate L 0 by minimizing (5.17) and the values of E 0 , n 0 , S 0 and L sym for different pairs of values (e, λ 2 ).
4. When we find a set of parameters (f π , e, λ 2 ) that fits the nuclear magnitudes within their respective errors of at most 15% then we calculate the corresponding EOS and solve the TOV system to obtain the mass-radius curve.
5. Finally, we accept the sets of values that satisfy the constraints, M max ≥ 2M and R 1.4M ≤ 12.5 km. These constraints are motivated from pulsar measurements [78,85,86]. We find that there is more than one set of parameters, so there is a residual freedom in the choice of these values that satisfy the nuclear physics magnitudes at saturation and NS observables.
We plot in fig. 20 the result of a scan for a fixed value of the energy scale, for other values of E s new sets of parameters were found.
The scan of parameters underscores the motivation for the generalized Skyrme model. The nuclear physics magnitudes are better fitted for very low or even null values of λ 2 since the sextic term reduces the value of S 0 . However, those sets of parameters are not accepted since they do not satisfy the maximum mass requirement. This reflects the importance of the sextic term in the extension of the Skyrme model to very high densities as inside NS.
In fig. 21 we plot 560 symmetry energy curves obtained from a first quick scan in blue, and in red we plot 23 representative cases from the larger set which have been fully minimized. We also represent at densities larger than n 0 some restrictions obtained from the most recent constraints of the analysis of neutron star observations, and at densities smaller than the saturation point which are more restrictive.
We have obtained the EOS from these 23 sets of parameters and compare them with some constraints obtained from a recent analysis [37]. In that work they build a huge number of physically well motivated EOS and compare the resulting NS with pulsars and GW observations. They conclude that the conformal limit in the speed of sound (c 2 s = 1/3) is expected to be surpassed inside NS. In fig. 22 we show the EOS obtained from our analysis and a good agreement is found with their results. The majority of our EOS exceeds the conformal bound too, and all of them lie inside the constrained region in the (ρ, p) diagram. We have cut the low density region in fig. 22 due to the absence of a crust in our EOS. However it is remarkable how the Skyrme model correctly describes the high density regime, which corresponds to the core of the NS and hence the main responsible for the mass of the star. Also the BCPM [75] EOS is represented as an accepted candidate to compare with in the diagram.
In a more extensive analysis of the parameters, we found ∼ 10, 000 accepted sets of parameters. In order to obtain these values we made the scan with the steps: ∆E s = 5 MeV, ∆e = 0.01, ∆λ 2 = 0.01 MeV fm 3 . As briefly mentioned before, the constraints on the symmetry energy yield rather stringent upper bounds on the sextic term coupling constant, we find that λ 2 3.4 MeV fm 3 . Nevertheless, we remark that a lower bound for this constant can also be obtained from the maximum mass requirement of neutron star EOS [88]. In this analysis we have found a lower bound of λ 2 0.29 MeV fm 3 . We solved the system for the same 23 cases and in fig. 23 we plot the MR curves, again with and without isospin effects for each case. The main conclusion is that the isospin always increases the radii of the stars. On the other hand, the isospin increases the masses of the stars with M 2.3M , for larger values the masses are reduced. This effect is also visible in fig. 22, where the red curves lie below the black lines at lower densities, hence the reds are stiffer, whilst for very high densities the situation is slightly the opposite. Another keypoint is that all the sets of parameters obtained in this analysis allow to have a wide range of maximum masses, M max ∼ 2 − 2.5M . This was an important feature in the NS obtained in the last section, and it is still possible using the fully minimized Skyrme crystals. This is of great important for the Skyrme model since it would be able to describe possible high-mass measurements like [81]. In addition, despite the difference in the maximum masses, the radii of the stars do not change as much when choosing some parameters or others, R 1.4M ∼ 12 − 13 km. However a final comment about the radii of The EOS for the same 23 values shown before. In black we plot the resulting EOS without isospin effects, whilst in red we consider npeµ matter. We find a good agreement between our EOS and the shaded regions obtained from the analysis in [37] at high densities. The purple region is an estimation for the range of the maximum central density inside NS, and the dots represent the maximum central energy densities in our models.
the NS requires the presence of a crust, since it will affect the radii of low mass NS.
As can be seen all the black lines satisfy (in good approximation) our M max and R 1.4M restrictions although they were imposed on the mass-radius curves resulting from the PS approximation. We have checked that the mass-radius curves obtained via the PS approximation are indeed quite similar to those of the fully minimized results for these 23 cases, so it confirms the PS approximation as a powerful and accurate tool for skyrmion crystals.
We also plot in fig. 23 the most likely mass-radius relations for the NS corresponding to GW170817 [?] and GW190425 [81] events. The green regions represent the estimations for the mass and radius values of J0030+0451 (bottom) [89], and a more recent analysis of the PSR J0740+6620 mass and radius from NICER (top) [85]. The purple region constraints the mass-radius curves from the statistical analysis done in [37], besides they also give an estimation for the maximum central energy density that a NS may support. We also do the comparison in fig. 22 between the range of values that they obtain (purple region) and our values (dots).
The greatest difference between the NS obtained using the interpolation between the submodels in the last section and those obtained from the minimization of the generalized Skyrme model using crystal solutions is found in the radii. Although we still have freedom to reach very high masses ∼ 2.5M , the radius increases when we consider the full model. The low-mass region (∼ 1.4M ) in the MR curves of fig. 17 and fig. 23 differ due to the presence of a crust in the first case. We do not include a crust in the EOS shown in this section since it would be more interesting to consider an EOS which already has a crust entirely obtained from the Skyrme model. This is still an open problem due to the behaviour of the Skyrme crystals at low densities, but the study of the new lattices presented in the first section may lead to the correct description of the full EOS within the Skyrme model. Nonetheless we have seen that the inclusion of a crust via the simple construction eq. (5.16) increases the radius of the NS around 1 km.
On the other hand, the 2.5M NS radius in fig. 17 is around 11.5 km, while the radius for the same mass in fig. 23 reaches 13 km. These high-mass NS are hardly affected by the presence of a crust, so the numerical simulation of the Skyrme crystal leads to a stiffer EOS. Nevertheless, the high-mass region of the last plot may be sharply improved with the inclusion of strangeness degrees of freedom in the system. Taking into account this effect provides even more realistic EOS, besides it is known that it will decrease the maximal mass as well as the radius, leading to a softer EOS at high densities.

Including kaons
The main ingredient to obtain the EOS for the Skyrme crystal is the energy dependence on the unit cell size. We have shown the results in the npeµ matter case above. However, once we include kaons, the change in the energy curve leads to a first order phase transition, see fig. 24. This can be clearly seen in the right plot of fig. 24, where we show the resulting EOS. There exists a region of diminishing pressure which implies a first order phase transition which, in principle, can be treated by a Maxwell construction or a Gibbs construction. The Maxwell construction (MC) is the standard procedure to obtain a physical EOS when a first order transition is present. Indeed, the MC has already been used for Skyrme crystals to describe the transition between crystals with different symmetries [45]. The MC is based on a mixed phase of constant pressure which connects the two solutions. However, the MC is only correct if there exists a single conserved charge (in this case, the baryon number) for which the associated chemical potential is enforced to be common for both phases in the mixed phase [90]. If, instead, an additional charge is conserved, like the electric charge in the case of npeµ matter, the Gibbs conditions for the phase equilibrium, cannot be both satisfied in a standard MC. In the last expression µ B and µ q represent the chemical potentials associated to the conserved baryon and electric charges, respectively. Instead, one should perform a Gibbs construction (GC) [90,91]. Indeed, the GC has also been proven useful in the context of a hadron-to-quark phase transition inside NS [92]. We may write the chemical potential of each particle species as a linear combination of the chemical potentials associated to the conserved charges of our system: where B i and q i are the baryon number and electric charge of the particle species i. Then we might identify the baryon and electric charge chemical potentials with the neutron and electron chemical potentials respectively. The main difference between MC and GC is that, in the mixed phase, the first one imposes charge neutrality independently for both phases, whereas in the GC it is only imposed globally in the mixed phase. Considering a volume fraction χ of the kaon condensed phase, charge neutrality is imposed in the GC as: The mixed phase in the GC is calculated by identifying first the contributions to the pressure and charge densities in each phase separately. Then we have to solve the system of equations composed by eqs. (4.28), (5.19) and (5.21). We show our results in fig. 25 and in fig. 26, using in all figures the parameter set 1 from table 6. We find that the mixed phase of the GC and, therefore, the density at which the kaon condensate sets in, is smaller density than the value obtained in table 6. This is also found in [91], for which the GC mixed phase extends to a larger region than the one obtained from the MC, because the mixed phase in the GC no longer is for constant pressure. In our case, even the minimum of E(L) is shifted to slightly lower values, see fig. 25 and, hence, the use of the GC affects the low density regime of the EOS, as can be seen in fig. 26.
We may also calculate the particle fractions in the mixed phase of the GC using an expression equivalent to eq. (5.21) for each particle. We show the new particle fractions in fig. 27. Besides, during the mixed phase, we find that there are more protons in the second phase than in the first one. However, the presence of more kaons than protons in the second phase results in a partial negative charge density. That negative charge is compensated by the overall positive charge density of the first phase. In both phases the number of electrons is much less than that of protons and kaons.   Figure 26: EOS for the three different cases that we have constructed. The jump in the MC due to the first order transition and the different behaviour of the GC at low densities are clearly visible. We also show the standard nuclear physics EOS of [75] (BCPM) and a hybrid EoS obtained by joining the BCPM EOS at low pressure with the GC EOS at high pressure. This figure was originally published in [49].
Once we have obtained the EOS by performing a Gibbs construction, the corresponding NS solutions can be obtained. We show in the left panel of fig. 28 the Mass-radius curves for the 4 different sets of parameters given in table 6. In the right panel of fig. 28, we show the same curves but after the addition of a crust to the NS. We do so by joining the GC equations of state of the Skyrme crystal with the BCPM EOS [75] for low densities, exactly as in Section 5.2. The two EOS join at the pressure p = p * where they coincide, i.e., where ρ BCPM (p * ) = ρ ,crystal (p * ). In terms of the baryon density, the joining occurs at n B, * ∼ 1.1n 0 for the parameter sets 1-3, and for n B, * ∼ 1.2n 0 for the set 4. Again, as in [50], we assume a smooth joining between the two EOS (concretely, described by a quadratic interpolation) in order to avoid an artificial phase transition at p * . We also plot in fig. 28 the most likely mass-radius relations for the NS corresponding to GW170817 [83] and GW190425 [81] events (orange and blue regions). The green regions represent the estimations for the mass and radius values of PSR J0740+6620 (top) [89] and J0030+0451 (bottom) [85]. The purple region constraints the mass-radius curves from the statistical analysis done in [37]. We find that the NS resulting from the addition of a crust to the Skyrme crystal EOS with a nonzero kaon condensate agree very well with these recent constraints. Further, the softening of the EOS due to the presence of kaons and the resulting smaller NS radii are important for this agreement.

Conclusions and Outlook
It was the main purpose of the present paper to report and review recent progress on the modeling of nuclear matter in terms of Skyrme crystals and its application to neutron stars. To put our results in perspective, in the first instant let us underline that, despite its conceptual strengths and some partial successes, the Skyrme model in its current state of exploration is not yet quantitatively competitive with more standard methods of nuclear physics in the description of the large number of available data on nuclei and nuclear matter at low densities. We believe, nevertheless, that the investigation of the Skyrme model as a possible model to describe nuclear matter is a worthy enterprise, for the following reasons. Firstly, once the assumption of periodicity of skyrmionic matter is accepted, the Skyrme model actually simplifies at higher densities and shows a more universal behavior which does not depend on many details (e.g., potential terms) which would be very relevant at low densities. The results reported in the present review are a clear demonstration of this fact. The Skyrme model, therefore, provides a rather simple alternative for the description of nuclear matter at high densities which is based on the qualitative assumptions that i) in the high-pressure region repulsive nuclear forces are more important than, e.g., degeneracy pressure, ii) the extended character of nucleons -which is a built-in property of the model -is important, and iii) the deconfinement transition is irrelevant inside NS. The inclusion of further hadronic degrees of freedom, on the other hand, is in principle straightforward, and we considered the condensation of kaons as a relevant example.
Secondly, the methods and calculations presented in this article can be applied with only minor modifications to more extended versions of the Skyrme model, where additional The effect of adding a standard nuclear physics crust to the Skyrme crystal EoS with kaon condensate obtained from the Gibbs construction (GC). This figure was originally published in [49].
terms and further degrees of freedom (e.g., vector mesons) are included. In other words, once promising candidates among these extended Skyrme models have been identified, their application to the study of nuclear matter at high densities and the resulting NS is relatively straightforward. The difficult part in finding these promising models is related to the calculation of skyrmion solutions, particularly for higher baryon number B, and to the identification and calculation of the most relevant quantum corrections. It is our hope that powerful up-to-date computing resources together with modern machine learning techniques like neural networks or other artificial-intelligence based methods will allow to find such promising models via a systematic exploration of the parameter space of the extended Skyrme models. In any case, the progress achieved in recent years is encouraging. In our investigation, we placed special emphasis on the importance of the sextic term of the Skyrme model for a viable description of strongly interacting matter at high densities. Further, the impact of this term on the resulting EOS and properties of NS was studied in detail. We also included the important effects of isospin quantization, relevant for the modeling of non-symmetric nuclear matter, and of kaon condensation on the skyrmionic matter properties. The overall behavior of the resulting models of nuclear matter at high density and the resulting NS is already in good agreement with observational data for certain choices of the coupling constants of the underlying Skyrme crystal. Still, there remain some minor differences with the EOS and NS that are most favored by recent observations. In particular, there are indications that the Skyrme crystal EOS is slightly too stiff in the intermediate density range (close to saturation), resulting in mass-radius curves with slightly bigger radii than the most likely ones according to a statistical analysis of observational data. On the other hand, these differences typically do not exceed one kilometer, demonstrating that the Skyrme crystal framework in its current state of development already provides a very reasonable description. There are several possible ways to overcome the remaining discrepancies.
There are, in fact, still some open questions remaining in the Skyrme crystal frame-work, and many interesting ways in which the results reviewed here can be extended. For instance, the ground state of the Skyrme model with periodic boundary conditions for values of the lattice parameter L slightly above the value at the minimum is not known in the large baryon number limit. That is to say, skyrmion matter in this region is not correctly described by the Skyrme crystal and will most likely be more inhomogeneous, see, e.g., fig. 13, which would result in a softening of the EOS there. This implies that some observables, such as the compression modulus, which measures the curvature of the energy vs lattice length, will not be well reproduced in the Skyrme crystal at the saturation point. Further investigations on the low density limit of the model are also required in order to be able to describe the physics of neutron star crusts, without recurring to a matching with other nuclear EOS. In particular, a more realistic description of these intermediate and low density regions will require the addition of further terms to the lagrangian and the inclusion of further degrees of freedom (DOF) such as hyperons or additional mesons. Additional DOF, in general, tend to soften the EOS, and some of these DOF may have a significant impact on the EOS. This question should be further investigated. Further, a more meticulous treatment of the quantum corrections to the crystal EOS should include, in addition to the isospin quantum effects, also the effect of quantum fluctuations associated to non-zero modes like, e.g., small fluctuation of the pion fields on top of the crystal background. The excitation of such modes could play a relevant role in the finite temperature case, but they might also be important in the zero temperature limit, in which case they contribute to the one loop correction to the energy of the crystal, or the Casimir energy. However, the computation of Casimir energies for solitons in three dimensions without spherical symmetry is in general extremely difficult, especially in nonrenormalizable theories such as the Skyrme model, and in the case of skyrmion crystals some sort of mean field approximation will probably be necessary. A further natural extension of the EOS presented here is the inclusion of the effects of finite temperature and/or a large magnetic field, as well as the study of transport properties of skyrmion crystals, which are crucial to describe nonequilibrium processes of nuclear matter.
Finally, we want to comment on the relevance of our findings for the generalized nuclear effective field theory (GnEFT) based on the hidden local symmetry [14] which we already mentioned in the introduction. Owing to their different field contents, many results are difficult to compare directly between the two theories. There are, however, some results in the GnEFT which are based on the assumption of a skyrmion-to-half-skyrmion phase transition somewhere above two times the nuclear saturation density [93]. All that we can say about this issue is that in all our investigations we did not find a sign of this phase transition. For zero pion mass, we formally do find a skyrmion-to-half-skyrmion phase transition (from the FCC to the FCC + crystal), but this transition is located on the thermodynamically unstable branch, i.e., at a baryon density below the density n 0 where the energy takes its minimum value, see fig. 8. Because of this instability, the Skyrme crystal and its resulting EOS are physically irrelevant in this region and, further, we know that there exist other, more inhomogeneous Skyrme matter solutions with lower energy, see fig. 13. It was natural for our purposes to identify n 0 with the nuclear saturation density, but the present argument, in fact, does not depend on this identification. In other words, our results for the generalized Skyrme model agree with the standard Skyrme model results in [41], which the authors of that paper summarized as "Thus the phase transition between a crystal of half-Skyrmions to a crystal of Skyrmions that was investigated in Refs. 5 -11 and 14 is not accessible, it appears on a thermodynamically unstable branch of the phase diagram" (the reference numbers are the references of that paper). Of course, even in the IR limit the effective coupling constants in the Skyrme-type model appearing in the GnEFT are different from ours, and at higher densities the impact of additional fields is difficult to gauge. But taking into account that a skyrmion-to-half-skyrmion phase transition in a physically relevant branch of the Skyrme crystal EOS has not been found in the full numerical study of any Skyrme model, this transition is probably rather unlikely to happen.