# Skyrme crystals, nuclear matter and compact stars

Christoph Adam<EMAIL_ADDRESS>Departamento de Física de Partículas,
Universidad de Santiago de Compostela, Spain, and Instituto Galego de Física
de Altas Enerxias (IGFAE) E-15782 Santiago de Compostela, Spain  Alberto
García Martín-Caro<EMAIL_ADDRESS>Departamento de Física de
Partículas, Universidad de Santiago de Compostela, Spain, and Instituto Galego
de Física de Altas Enerxias (IGFAE) E-15782 Santiago de Compostela, Spain
Miguel Huidobro<EMAIL_ADDRESS>Departamento de Física de
Partículas, Universidad de Santiago de Compostela, Spain, and Instituto Galego
de Física de Altas Enerxias (IGFAE) E-15782 Santiago de Compostela, Spain
Andrzej Wereszczynski<EMAIL_ADDRESS>Institute of Physics,
Jagiellonian University, Lojasiewicza 11, Kraków, Poland

###### Abstract

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.

Dedicated to the memory of our unforgettable friend and colleague Ricardo
Vázquez.

Invited contribution to the Special Issue: Symmetries and Ultra Dense Matter
of Compact Stars, Symmetry 2023, 15(4), 899

###### Contents

  1. 1 Introduction
    1. 1.1 The Skyrme model
    2. 1.2 Skyrme matter
  2. 2 Skyrme crystals
    1. 2.1 Crystal ansatz
      1. 2.1.1 Face centered cubic crystal of skyrmions
      2. 2.1.2 Body centered cubic crystal of half-skyrmions
      3. 2.1.3 Enhanced face centered cubic crystal of skyrmions
    2. 2.2 Energy curves
    3. 2.3 Phase transitions
      1. 2.3.1 Low density phase transition
      2. 2.3.2 High density phase transition
      3. 2.3.3 Fluid-like transition
    4. 2.4 New lattice solutions
  3. 3 Isospin quantization and symmetry energy
  4. 4 Kaon condensate in skyrmion crystals
    1. 4.1  The kaon condensate effective potential
    2. 4.2  Quantum corrections to the effective potential
  5. 5 Neutron Stars
    1. 5.1 TOV system of equations
    2. 5.2 A generalized Skyrme EOS
    3. 5.3 Skyrme crystal EOS
      1. 5.3.1 Scan of the parameters
    4. 5.4 Including kaons
  6. 6 Conclusions and Outlook

##  1 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).

Figure 1: Schematic phase diagram of QCD

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].

###  1.1 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=\sigma\mathbbm{I}+i\pi_{a}\tau_{a},\quad
U^{\dagger}U=\mathbbm{I}\longrightarrow\sigma^{2}+\pi_{a}\pi_{a}=1,\quad
a=1,2,3$ (1.1)

where $\tau_{a}$ are the Pauli matrices. Left ($L$) and right ($R$) chiral
transformations act on this field like $U\to LUR$. 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

$\mathcal{L}_{2}=-\frac{f_{\pi}^{2}}{16}\Tr\\{L_{\mu}L^{\mu}\\}\;,\quad
L_{\mu}=U^{\dagger}\partial_{\mu}U$ (1.2)

providing a kinetic term for the pions, and the pion mass potential

$\mathcal{L}_{0}=\frac{m_{\pi}^{2}f_{\pi}^{2}}{8}\Tr\\{U-\mathbbm{I}\\}.$
(1.3)

Here, $f_{\pi}$ is the pion decay constant, whose physical value in the
conventions used here is $f_{\pi}=186$ MeV. Further, $L_{\mu}$ are the
components of the left-invariant Maurer-Cartan form of the $SU(2)$ group,
$m_{\pi}$ is the pion mass, and $\mathbbm{I}$ is the $2\times 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\rightarrow\Lambda x$. To stabilize them, Skyrme added the term

$\mathcal{L}_{4}=\frac{1}{32e^{2}}\Tr\\{[L_{\mu},L_{\nu}]^{2}\\}$ (1.4)

where $e$ is a dimensionless parameter. Here the notation $\mathcal{L}_{n}$
means that the corresponding term contains $n$ first derivatives.
$\mathcal{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
$\mathcal{L}_{24}=\mathcal{L}_{2}+\mathcal{L}_{4}$. The resulting static
solutions of the Skyrme field are maps from real space $\mathbbm{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\rightarrow\infty)=\mathbbm{I}$, such that the infinity
of $\mathbbm{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

$B^{\mu}=\frac{\epsilon^{\mu\nu\rho\sigma}}{24\pi^{2}}\Tr\\{L_{\nu}L_{\rho}L_{\sigma}\\}.$
(1.5)

It is straightforward to see that the divergence of the topological current
vanishes identically, ($\partial_{\mu}B^{\mu}=0$), which implies the
conservation of the topological number,

$B=\int d^{3}xB^{0},$ (1.6)

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, $\pi_{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,

$\mathcal{L}_{6}=\lambda^{2}\pi^{4}B_{\mu}B^{\mu}$ (1.7)

where $\lambda$ 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,

$\mathcal{L}=\mathcal{L}_{2}+\mathcal{L}_{4}+\mathcal{L}_{6}+\mathcal{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 $\mathcal{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_{c}^{-1}$ 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
$\mathcal{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 $\alpha$-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 $\alpha$ 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.
$\mathcal{L}_{0,k}\simeq\Tr\\{U-\mathbbm{I}\\}^{k}$, $k\in\mathbbm{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 $\mathcal{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\leq 8$ were investigated numerically,
and it was found that realistic cluster structures emerged for all skyrmions,
even for substructures different from $\alpha$ 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 density111In section 2.4 we will discuss non-crystalline
solutions consisting of regularly arranged clusters of skyrmions surrounded by
empty space. These solutions have a much better behavior at low density, but
we will not attempt to model nuclear matter below saturation with these
configurations in the present paper.. At such high densities potential terms
like $\mathcal{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 gauge222We shall, however, take into account the effect of kaon
condensation in section 4, because a kaon condensate can be studied in the
background of an unmodified Skyrme crystal in leading approximation.. 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.

###  1.2 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_{\odot}$ [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_{\odot}$. 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_{s}^{2}<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}\sim n_{0}=0.16$ fm-3, to
the cores of neutron stars $n_{B}\sim 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_{\pi}$, the Skyrme parameter $e$, the pion mass $m_{\pi}$, and
$\lambda^{2}$. From these parameters, we fix the pion mass to its physical
value, $m_{\pi}=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_{\pi}$ 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_{\pi}$ partly takes into account the effect of neglected
terms, whereas in more complete, general versions of the model the optimum
value of $f_{\pi}$ 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,

$E_{s}=\frac{3\pi^{2}f_{\pi}}{e},\quad x_{s}=\frac{1}{f_{\pi}e}.$ (1.9)

Then the static energy functional in these units becomes

$\displaystyle E=\frac{1}{24\pi^{2}}\int d^{3}x$
$\displaystyle\left[-\frac{1}{2}\Tr\left\\{L_{i}^{2}\right\\}-\frac{1}{4}\Tr\left\\{\left[L_{i},L_{j}\right]^{2}\right\\}+8\lambda^{2}\pi^{4}f^{2}_{\pi}e^{4}(\mathcal{B}^{0})^{2}+\left(\frac{m_{\pi}}{f_{\pi}e}\right)^{2}\Tr(\mathbbm{I}-U)\right]=$
$\displaystyle\frac{1}{24\pi^{2}}\int d^{3}x$
$\displaystyle\left[(\partial_{i}n_{A})^{2}+\left(\partial_{i}n_{A}\partial_{j}n_{B}-\partial_{i}n_{B}\partial_{j}n_{A}\right)^{2}+c_{6}\left(\epsilon_{ABCD}n_{A}\partial_{1}n_{B}\partial_{2}n_{C}\partial_{3}n_{D}\right)^{2}+c_{0}(1-\sigma)\right],$
(1.10)

where we have defined $c_{6}=2\lambda^{2}f_{\pi}^{2}e^{4}$ and
$c_{0}=2m_{\pi}^{2}/(f_{\pi}e)^{2}$. In the last expression, we have
introduced the field configuration eq. 1.1 and adopted the vectorial notation
$n_{A}=(\sigma,\pi_{a})$. Recall that $U\in 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 ($\lambda^{2}=m_{\pi}=0$) becomes independent
of the value of the parameters $f_{\pi}$ 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_{\rm cell}$ and side length $l_{\rm cell}$ (i.e.,
volume $V_{\rm cell}=l_{\rm cell}^{3}$), such that the Skyrme field is
periodic in all three Cartesian directions under the translation $x^{i}\to
x^{i}+l_{\rm 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_{\rm cell}(l_{\rm cell})$ is a convex function with a
minimum for a certain $l_{0,{\rm cell}}$. For some crystals, the resulting
energy at the minimum $l_{0,{\rm cell}}$ is only slightly above the
topological energy bound. For the FCC crystal for the original Skyrme model,
e.g., $E_{\rm cell}(l_{0,{\rm 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_{\rm cell}\not=l_{0,{\rm cell}}$, the energy grows rather quickly.
Here, the region $l_{\rm cell}>l_{0,{\rm 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,{\rm 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_{{\rm cell}}>l_{0,{\rm cell}}$ still is slightly above
$E_{\rm cell}(l_{0,{\rm cell}})$, but the difference between $E_{\rm
cell}(l_{0,{\rm cell}})$ and, say, $\lim_{l_{{\rm cell}}\to\infty}E_{\rm
cell}(l_{{\rm cell}})$ can be made smaller than 5%.

The nuclear pasta phase between $l_{0,{\rm 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 $\alpha$-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_{\rm
cell}/V_{\rm cell}$, corresponding to the nuclear saturation density
$n_{0}=B_{\rm cell}/V_{{0,\rm 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}\leq 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_{\rm cell}<l_{0,{\rm 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.

##  2 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_{\odot}/m_{p}\sim 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
$\mathbbm{T}^{3}$, however since the $\mathbbm{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,

$\frac{E_{\rm crystal}}{B_{\rm
crystal}}=\frac{N_{\rm{cell}s}E_{\rm{cell}}}{N_{\rm{cell}s}B_{\rm{cell}}}=\frac{E_{\rm{cell}}}{B_{\rm{cell}}}.$
(2.1)

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.

###  2.1 Crystal ansatz

The starting point in the construction of the Skyrme crystal proposed in [40]
is the expansion of the fields in Fourier series,

$\displaystyle\sigma=\sum^{\infty}_{a,b,c}\beta_{abc}\cos\left(\frac{a\pi
x}{L}\right)\cos\left(\frac{b\pi y}{L}\right)\cos\left(\frac{c\pi
z}{L}\right)$ (2.2)
$\displaystyle\pi_{1}=\sum^{\infty}_{h,k,l}\alpha_{hkl}\sin\left(\frac{h\pi
x}{L}\right)\cos\left(\frac{k\pi y}{L}\right)\cos\left(\frac{l\pi
z}{L}\right)$
$\displaystyle\pi_{2}=\sum^{\infty}_{h,k,l}\alpha_{hkl}\cos\left(\frac{l\pi
x}{L}\right)\sin\left(\frac{h\pi y}{L}\right)\cos\left(\frac{k\pi
z}{L}\right)$
$\displaystyle\pi_{3}=\sum^{\infty}_{h,k,l}\alpha_{hkl}\cos\left(\frac{k\pi
x}{L}\right)\cos\left(\frac{l\pi y}{L}\right)\sin\left(\frac{h\pi
z}{L}\right).$

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 $\beta_{abc}$ and $\alpha_{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 $\sim 1$, we have calculated that the next
orders decay to the $\sim 4\%$, $\sim 0.3\%$ and $\sim 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

$n_{A}\longrightarrow\frac{n_{A}}{\sqrt{n_{A}n_{A}}}.$ (2.3)

The crystal considered by Klebanov has the simplest unit cell. It is invariant
under cubic symmetry transformations:

$\displaystyle\text{A}_{1}$ $\displaystyle:(x,y,z)\rightarrow(-x,y,z),$
$\displaystyle(\sigma,\pi_{1},\pi_{2},\pi_{3})\rightarrow(\sigma,-\pi_{1},\pi_{2},\pi_{3}),$
(2.4) $\displaystyle\text{A}_{2}$ $\displaystyle:(x,y,z)\rightarrow(y,z,x),$
$\displaystyle(\sigma,\pi_{1},\pi_{2},\pi_{3})\rightarrow(\sigma,\pi_{2},\pi_{3},\pi_{1}),$
(2.5)

and it has an additional periodicity symmetry on the side length of the unit
cell,

$\displaystyle\text{A}_{3}$ $\displaystyle:(x,y,z)\rightarrow(x+L,y,z),$
$\displaystyle(\sigma,\pi_{1},\pi_{2},\pi_{3})\rightarrow(\sigma,-\pi_{1},\pi_{2},-\pi_{3}).$
(2.6)

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.

####  2.1.1 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,

$\displaystyle\text{C}_{3}$ $\displaystyle:(x,y,z)\rightarrow(x,z,-y),$
$\displaystyle(\sigma,\pi_{1},\pi_{2},\pi_{3})\rightarrow(\sigma,-\pi_{1},\pi_{3},-\pi_{2}),$
(2.7) $\displaystyle\text{C}_{4}$
$\displaystyle:(x,y,z)\rightarrow(x+L,y+L,z),$
$\displaystyle(\sigma,\pi_{1},\pi_{2},\pi_{3})\rightarrow(\sigma,-\pi_{1},-\pi_{2},\pi_{3}).$
(2.8)

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 $\pi$ 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_{\rm{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.

Figure 2: Energy density contour plots of the face-centered cubic (FCC) unit
cell of skyrmions for a large value of $L=10$. Each plot shows different
heights of the unit cell, which are $z=0$, $L/2$ and $L$ because in this case
also the energy density has periodicity in $2L$. Figure 3: The energy density
of the FCC unit cell is shown in 3D. The FCC symmetry is clearly visible as
well as in the contours fig. 2 plot. We adopt the Runge colouring convention
[52] in this figure to represent the pion fields.

####  2.1.2 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,

$\displaystyle\text{B}_{4}$
$\displaystyle:(x,y,z)\rightarrow(L/2-z,L/2-y,L/2-x),$
$\displaystyle(\sigma,\pi_{1},\pi_{2},\pi_{3})\rightarrow(-\sigma,\pi_{2},\pi_{1},\pi_{3}).$
(2.9)

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 $\sigma=-1$ at the centre of the
cube of side length $L$ until some radius $r_{0}$ for which $\sigma=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
$(\sigma,\pi_{a})\rightarrow(-\sigma,-\pi_{a})$, these new solutions are
located in the corners of the cube, connected to the $\sigma=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 $\sigma$ field in this cube is exactly 0, so
the energy coming from the potential term $\mathcal{L}_{0}\sim(\sigma-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_{\rm{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.

  * •

$\beta_{abc}=\beta_{bca}=\beta_{cab}$.

  * •

$\alpha_{hkl}=-(-1)^{\frac{h+k+l}{2}}\alpha_{khl}$.

  * •

$\beta_{abc}=-(-1)^{\frac{a+b+c}{2}}\beta_{bac}$.

We show the resulting energy density plots in fig. 4 and fig. 5.

Figure 4: Energy density contour plots for the unit cell of the body-centered-
cubic (BCC) crystal at the minimum of energy. In this case, the energy density
has period $L$ so we show the cuts at $z=0$, $z=L/4$ and $z=L/2$. 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.

####  2.1.3 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,

$\displaystyle\text{D}_{4}$ $\displaystyle:(x,y,z)\rightarrow(x+L,y,z),$
$\displaystyle(\sigma,\pi_{1},\pi_{2},\pi_{3})\rightarrow(-\sigma,-\pi_{1},\pi_{2},\pi_{3}),$
(2.10)

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 $\sigma=-1$ in their centre are located at the corners and
faces of the unit cell. Further, the opposite half-skyrmions with $\sigma=1$
occupy the body center and the link centers of the unit cube. As a
consequence, the mean value of the $\sigma$ 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_{\rm{cell}}=4$. We show the resulting energy density plots in
fig. 6 and in fig. 7.

Figure 6: Energy density contour plots for the unit cell of the enhanced face-
centered cubic crystal (FCC+) at the minimum of energy. It can be appreciated
in the figure that this crystal can also be viewed as a simple cubic crystal
of half-skyrmions. The energy density has the same period as the BCC, so the
cuts are the same. Figure 7: The energy density in three dimensions for the
FCC+ crystal. Again the Runge colouring convention is used as in the last two
cases.

###  2.2 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}\sim L^{3-i}$, where $i$ represents the number of derivatives in the
$i-$term, and the three comes from the integration over $3$D space. Hence we
use the following fit for the energy curves in each case,

$E(L)=k+k_{2}L+\frac{k_{4}}{L}+c_{6}\frac{k_{6}}{L^{3}}+c_{0}k_{0}L^{3}.$
(2.11)

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_{\pi}=129\>\text{MeV},\quad e=5.45,\quad\lambda^{2}=5\>\text{MeV
fm}^{3},\quad m_{\pi}=140\>\rm{MeV}.$ (2.12)

Figure 8:  Energy curves of the different Skyrme crystals considering the
different terms in the generalized Skyrme lagrangian. We also show the BPS
bound in each case to see how close it is to the minima of the crystals.

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_{\pi}$ or $\lambda^{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
$\mathcal{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 $\mathcal{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 [45, 44], 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,

Model | $L_{0}$ | $E_{0}/B$ | Minimum ($\%$)
---|---|---|---
$\mathcal{L}_{24}$ | 4.7 | 1.04 | 4
$\mathcal{L}_{240}$ | 4.1 | 1.09 | 5
$\mathcal{L}_{246}$ | 6.2 | 1.18 | 5
$\mathcal{L}_{2460}$ | 5.3 | 1.30 | 7
Table 1: FCC+ crystal.

Model | $L_{0}$ | $E_{0}/B$ | Minimum ($\%$)
---|---|---|---
$\mathcal{L}_{24}$ | 5.5 | 1.08 | 8
$\mathcal{L}_{240}$ | 4.9 | 1.13 | 9
$\mathcal{L}_{246}$ | 7.3 | 1.23 | 9
$\mathcal{L}_{2460}$ | 6.4 | 1.34 | 10
Table 2: BCC crystal

The side length and the value at which the $E(L)$ finds the minimum. In the
last column we show $[(E_{\rm min}-E_{\rm bound})/E_{\rm bound}]\times 100$,
i.e., the percentage deviation of the minimum crystal energy from the bound
for the lattice.

The values of the coefficients in eq. 2.11 are obtained in each case using the
Python optimization library GEKKO.

Model | $k$ | $k_{2}$ | $k_{4}$ | $k_{6}$ | $k_{0}$
---|---|---|---|---|---
$\mathcal{L}_{24}$ | 0.047 | 0.105 | 2.344 | 0 | 0
$\mathcal{L}_{240}$ | 0.022 | 0.109 | 2.384 | 0 | 0.008
$\mathcal{L}_{246}$ | 0.006 | 0.106 | 2.750 | 0.905 | 0
$\mathcal{L}_{2460}$ | 0.334 | 0.074 | 1.747 | 1.062 | 0.010
Table 3: Fitting constants for the numerically obtained $E(L)$ curves for the FCC+ crystal. Model | $k$ | $k_{2}$ | $k_{4}$ | $k_{6}$ | $k_{0}$
---|---|---|---|---|---
$\mathcal{L}_{24}$ | 0.017 | 0.096 | 2.988 | 0 | 0
$\mathcal{L}_{240}$ | -0.022 | 0.101 | 3.061 | 0 | 0.004
$\mathcal{L}_{246}$ | 0.040 | 0.093 | 3.150 | 1.710 | 0
$\mathcal{L}_{2460}$ | 0.172 | 0.078 | 2.798 | 1.751 | 0.005
Table 4: Fitting constants for the numerically obtained $E(L)$ curves for the
BCC crystal.

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.

###  2.3 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.

####  2.3.1 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 $\sigma\rightarrow-\sigma$. 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 $\sigma$
field was proposed to be the order parameter of the transition.

Figure 9:  The mean value of the $\sigma$ 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 $\sigma$ field in the unit cell for
the cases $\mathcal{L}_{24}$ and $\mathcal{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.

####  2.3.2 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,

$n_{B}=\frac{B_{\rm{cell}}}{V_{\rm{cell}}}.$ (2.13)

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.

Figure 10:  The comparison between the energies of the BCC and FCC crystals at
the same baryon density shows that the first one becomes more favourable at
some point, denoted by the red cross in the plots. The sextic term is crucial
to have this transition at realistic densities. Both $n$ and $E$ are shown in
Skyrme units.

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.

Model | $n_{\rm PT}/n_{0}$
---|---
$\mathcal{L}_{246}$ | 3.7
$\mathcal{L}_{2460}$ | 2.4
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,

$\displaystyle p$ $\displaystyle=-\frac{\partial E}{\partial
V}=-\frac{1}{24L^{2}}\frac{\partial E_{\rm cell}}{\partial L},$ (2.14)
$\displaystyle\rho$ $\displaystyle=\frac{E}{V}=\frac{E_{\rm cell}}{8L^{3}}.$
(2.15)

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,

$p_{1}=p_{2},\quad\mu_{B}^{1}=\mu_{B}^{2}.$ (2.16)

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
($\mu_{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_{PT}$, and join
them with the curve which has the same value of $\mu_{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_{PT}$).

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,

$n_{\rm FCC}=n_{\rm BCC}\longrightarrow V_{\rm FCC}=\frac{V_{\rm BCC}}{2}.$
(2.17)

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,

$E(V)=\left\\{\begin{array}[]{ll}E_{\rm FCC}(V),&\quad V\geq V_{1}\\\ E_{\rm
FCC}(V_{1})-p_{PT}(V-V_{1}),&V_{1}\geq V\geq V_{2}\\\ E_{\rm BCC}(V),&\quad
V\leq V_{2}\end{array}\right.$ (2.18)

Figure 11: Maxwell construction in the FCC-to-BCC phase transition. The red
cross denotes again the transition point, but now the black crosses denote
where the mixed phase starts. At low densities we are in the FCC phase (green
line) until we reach the mixed phase (black line), which joins to the BCC
phase (blue line).

see 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 $\mathcal{L}_{246}$ case and the difference between them in the
$\mathcal{L}_{2460}$ case is negligible.

####  2.3.3 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 ($\mathcal{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, $\rho_{\rm
mean}=E_{\rm cell}/V_{\rm 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 $\mathcal{L}_{240}$ and $\mathcal{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$,

$E(r)=\int_{0}^{r}d^{3}x\>\rho,$ (2.19)

where $\rho$ is the $\rho_{\rm 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 $\chi=E(r)/E_{\rm mean}(r)$ of the two REP for each case.

Figure 12: REP for the $\mathcal{L}_{240}$ (blue) and $\mathcal{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 $\mathcal{L}_{240}$, whereas the energy density becomes much more
homogeneous with increasing baryon density for $\mathcal{L}_{2460}$.

###  2.4 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\geq 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
$\alpha$-lattice.

We calculated the energy for the $\alpha$-lattice at different values of $L$
for the $\mathcal{L}_{240}$ and $\mathcal{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\rightarrow\infty$. 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
$\alpha-$lattice, the $B=32$ lattice.

Figure 13: Energy per baryon number for the different lattices that we
consider with and without sextic term. These $\alpha$ 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 $\alpha$ 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 $\mathcal{L}_{240}$ model.

Figure 14: Energy density contour of the $\alpha$ and $B=32$ lattices at
different sizes of the unit cell. The shape of the isolated skyrmion is
visible at very large values of $L$, whereas the FCC+ crystal is recovered at
high densities. The legend of the colours is as in fig. 6

##  3 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 $\delta=(N-Z)/A=(1-2\gamma)$, with $\gamma=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,

$\frac{E}{A}(n_{B},\delta)=E_{N}(n_{B})+S_{N}(n_{B})\delta^{2}+\order{\delta^{3}},$
(3.1)

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}$,

$S_{N}(n_{B})=S_{0}+\frac{1}{3}L\epsilon+\frac{1}{18}K_{\rm
sym}\epsilon^{2}+\cdots$ (3.2)

with $\epsilon=(n-n_{0})/n_{0}$, and

$L=3n_{0}\partialderivative{S_{N}}{n}\evaluated{}_{n=n_{0}},\quad K_{\rm
sym}=9n_{0}^{2}\partialderivative[2]{S_{N}}{n}\evaluated{}_{n=n_{0}}$ (3.3)

the slope and curvature of the symmetry energy at saturation, respectively.The
symmetry energy at saturation is well constrained ($S_{0}\sim 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, 62, 63, 24]).

Following this approach, we consider a (time-dependent) isospin transformation
of a static Skyrme field configuration,

$U(\vec{x})\rightarrow\tilde{U}(\vec{x},t)\equiv
g(t)U(\vec{x})g^{\dagger}(t).$ (3.4)

Then, the time component of the left invariant form $L_{\mu}$ becomes
$\tilde{U}^{\dagger}\partial_{0}\tilde{U}=gT_{a}g^{\dagger}\omega_{a}$, where
$T_{a}$ is the $\mathfrak{su}(2)$-valued current,

$T_{a}=\frac{i}{2}U^{\dagger}[\tau_{a},U]=i(\pi_{a}\pi_{b}-\pi_{c}\pi_{c}\delta_{ab}+\sigma\pi_{c}\epsilon_{abc})\tau_{b}\equiv
iT_{b}^{a}\tau_{b},$ (3.5)

and $\vec{\omega}\doteq g^{\dagger}\dot{g}=\tfrac{i}{2}\omega_{a}\tau_{a}$ is
the isospin angular velocity. The time dependence of the new Skyrme field
induces a kinetic term in the energy functional, given by 333Remember that we
are using the mostly minus convention for the metric signature.

$\displaystyle T=\frac{1}{24\pi^{2}}\int\\{$ $\displaystyle
a\Tr\\{L_{0}L_{0}\\}-2b\Tr\\{[L_{0},L_{k}][L_{0},L_{k}]\\}-c\,\,\mathcal{B}^{i}\mathcal{B}_{i}\\}d^{3}x=\frac{1}{2}\omega_{i}\Lambda_{ij}\omega_{j}$

where $\Lambda_{ij}$ is the isospin inertia tensor, given by

$\Lambda_{ij}=\frac{1}{24\pi^{2}}\int
d^{3}x\left\\{2a\Tr\\{T_{i}T_{j}\\}-4b\Tr\\{[T_{i},L_{k}][T_{j},L_{k}]\\}\right.-\left.\frac{c}{32\pi^{4}}\varepsilon^{abc}\Tr\\{T_{i}L_{b}L_{c}\\}\varepsilon_{ars}\Tr\\{T_{j}L_{r}L_{s}\\}\right\\}\,,$
(3.6)

and the values of $a,b,c,d$ are easily obtained from (1.8),

$a=-\frac{1}{2},\quad b=\frac{1}{4},\quad
c=-8\lambda^{2}\pi^{4}\frac{f_{\pi}^{2}e^{4}}{\hbar^{3}},\quad
d=\frac{m^{2}_{\pi}}{f^{2}_{\pi}e^{2}}.$ (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

$\Lambda=\frac{1}{24\pi^{2}}\left[2a\Lambda^{(2)}-4b\Lambda^{(4)}-\frac{c}{32\pi^{4}}\Lambda^{(6)}\right],$
(3.8)

where the contribution from each term in the lagrangian, denoted by
$\Lambda^{(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_{\text{cells}}\equiv N$ of unit cells, which is reduced to

$T=\frac{1}{2}\omega_{i}\Lambda_{ij}\omega_{j}=N\frac{1}{2}\Lambda\omega_{a}\omega^{a},$
(3.9)

and, by defining the corresponding canonical momentum $J_{a}=\partial
L/\partial\omega^{a}=N\Lambda\omega_{a}$, we may write it in Hamiltonian form,

$H=\frac{1}{2N\Lambda}J_{a}J^{a}.$ (3.10)

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,

$H=\frac{\hbar^{2}}{2N\Lambda}J^{\rm{tot}}(J^{\rm{tot}}+1)$ (3.11)

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_{\text{cell}}$, the total isospin will be $\frac{1}{2}B_{\text{cell}}$, and
hence the total isospin of the full crystal will be
$J^{\rm{tot}}=\frac{1}{2}NB_{\text{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\rightarrow\infty$)

$E_{\rm{iso}}=\frac{\hbar^{2}}{8\Lambda}B_{\text{cell}}^{2}.$ (3.12)

Such correction could also have been obtained directly by introducing an
’external’ isospin chemical potential $\mu_{I}$ , and promoting the regular
derivatives in the Skyrme Lagrangian to covariant derivatives of the form [64]

$\partial_{\mu}U\rightarrow
D_{\mu}U=\partial_{\mu}U-\frac{i\mu_{I}}{2}\delta_{\mu 0}[\tau^{3},U],$ (3.13)

so that, if $U$ is a static configuration, the time component of the Maurer-
Cartan form becomes

$L_{0}=-\frac{i}{2}\mu_{I}U^{\dagger}[\tau^{3},U]=-\mu_{I}T_{3}.$ (3.14)

This expression is equivalent to that of an iso-rotating field with angular
velocity $\omega_{a}=-\mu_{I}\delta_{3a}$. Thus, it is straightforward to
obtain the isospin chemical potential for the Skyrmion crystal using its
thermodynamical definition $\mu_{I}=-\partialderivative{E}{n_{I}}$, where
$n_{I}$ is the (third component of) the isospin number density. Given that
$(J^{\rm{tot}})^{2}=J_{1}^{2}+J_{2}^{2}+J_{3}^{2}$ and $n_{I}=J_{3}/N$, we may
write the isospin energy per unit cell as

$E_{\rm{iso}}=\frac{\hbar^{2}}{2\Lambda}\quantity(n_{I}^{2}+\frac{J_{2}^{2}}{N^{2}}+\frac{J_{1}^{2}}{N^{2}})$
(3.15)

and then

$\mu_{I}=-\partialderivative{E_{\rm{iso}}}{n_{I}}=-\frac{\hbar^{2}}{\Lambda}n_{I}.$
(3.16)

For $B<\infty$ 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\leq 12$ [63,
24]. 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

$\expectationvalue{I^{0}_{3}}=\frac{\expectationvalue{I_{3}}}{\int
d^{3}x}=\frac{\bra{\Psi}\int I^{0}_{3}d^{3}x\ket{\Psi}}{NV_{\rm
cell}}\doteq\frac{n_{I}}{V_{\rm cell}}$ (3.17)

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_{\text{cell}}$ is
given by

$n_{I}=-\frac{1}{2}(1-2\gamma)B_{\text{cell}}=-\frac{B_{\text{cell}}}{2}\delta.$
(3.18)

Hence, the isospin energy per unit cell of the Skyrmion crystal can be written
in terms of the asymmetry parameter

$E_{\rm iso}=\frac{\hbar^{2}B_{\text{cell}}^{2}}{8\Lambda}\delta^{2},$ (3.19)

from where we can readily read off the symmetry energy for Skyrme crystals

$S_{N}(n_{B})=\frac{\hbar^{2}L^{3}}{\Lambda}n_{B}.$ (3.20)

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 $\Lambda$ (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_

$n_{p}=\frac{Z}{V}=n_{e}+n_{\mu},$ (3.21)

and the _$\beta$ -equilibrium condition_

$\mu_{n}=\mu_{p}+\mu_{l}\implies\mu_{I}=\mu_{l},\qquad l=e,\mu,$ (3.22)

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\leftrightarrow p+l+\bar{\nu}_{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

$\mu_{l}=\sqrt{(\hbar k_{F})^{2}+m_{l}^{2}}$ (3.23)

where $k_{F}=(3\pi^{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, $\mu_{e}\geq m_{\mu}$, 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}=\gamma
B_{\text{cell}}/(2L)^{3}$. The $\beta$ equilibrium condition then provides an
equation which defines $\gamma$ implicitly as a function of the lattice length
parameter,

$\frac{\hbar
L}{\Lambda}(1-2\gamma)=\quantity(\frac{3\pi^{2}}{B_{\text{cell}}^{2}})^{1/3}\gamma^{1/3}$
(3.24)

where we also assumed the ultrarelativistic electron approximation, i.e.
$m_{l}/k_{F}\simeq 0$.

Including the muon contribution to the charge density yields a slightly more
complicated expression for the $\beta$-equilibrium condition, given by

$\frac{\hbar
B_{\text{cell}}}{2\Lambda}(1-2\gamma)=\quantity[3\pi^{2}\quantity(\frac{\gamma
B_{\text{cell}}}{8L^{3}}-n_{\mu})]^{\tfrac{1}{3}},$ (3.25)

where

$n_{\mu}=\frac{1}{3\pi^{2}}\quantity[\quantity(\frac{\hbar
B_{\text{cell}}(1-2\gamma)}{2\Lambda})^{2}-\left(\frac{m_{\mu}}{\hbar}\right)^{2}]^{\tfrac{3}{2}}.$
(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\rightarrow
p+e+\bar{\nu}_{e}$. This process is expected to occur if the proton fraction
reaches a critical value, $\gamma_{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 ultra-high densities. It is, however, not clear
whether this enhanced cooling actually occurs, although there is recent
evidence that supports it [69].

In $npe\mu$ matter, the DU threshold is given by [67]

$x_{DU}=\frac{1}{1+(1+(\frac{n_{e}}{n_{e}+n_{\mu}})^{1/3})^{3}}.$ (3.27)

Figure 15: Fraction density $\gamma_{i}$ for each particle as a function of
the baryon density for $\lambda^{2}=0$ (solid) and $\lambda^{2}=1.5$ (dashed).
The corresponding DU threshold is also shown in black. This figure was
originally published in [48].

We show the particle populations $\gamma_{i}$ in beta-equilibrated skyrmion
matter in fig. 15 for the cases $\lambda^{2}=0,1.5$ MeV fm3.

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_{\pi},e,\lambda^{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-3$ times the nuclear saturation density,
additional degrees of freedom (strange baryons) appear and become important
for the description of nuclear matter which, in particular, may affect the
proton fraction at these densities.

##  4 Kaon condensate in skyrmion crystals

Up to this point, we have shown how to describe $\beta$-equilibrated
$npe\mu$-matter at finite density in terms of a semi-classical skyrmion
crystal and a leptonic sea of electrons and muons, and we have determined the
different particle fractions as a function of density in the mean field
approximation. However, for densities about twice nuclear saturation and
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].

###  4.1  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)$ skyrmion-like 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]:

$U=\sqrt{U_{K}}U_{\pi}\sqrt{U_{K}}.$ (4.1)

In this ansatz $U_{\pi}$ 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_{\pi}$ 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
$\mathfrak{su}$(3) -valued matrix $\mathcal{D}$ which is non trivial in the
off-diagonal elements:

$U_{\pi}=\begin{pmatrix}u&0\\\ 0&1\end{pmatrix},\hskip
8.53581ptU_{K}=e^{i\frac{2\sqrt{2}}{f_{\pi}}\mathcal{D}},\quad
u=\sigma+i\pi_{a}\tau_{a},\quad\mathcal{D}=\matrixquantity(0&K\\\
K^{\dagger}&0)$ (4.2)

where $K$ consists of a scalar doublet of complex fields representing charged
and neutral kaons:

$K=\matrixquantity(K^{+}\\\ K^{0}),\quad K^{\dagger}=(K^{-},\bar{K}^{0}).$
(4.3)

We now extend the Generalized Skyrme Lagrangian (1.8) to include strange
degrees of freedom. First, we need to replace the potential term
$\mathcal{L}_{0}$ in order to correctly take into account the mass of the kaon
fields. The new potential term has the form [72]:

$\displaystyle\mathcal{L}_{0}^{\rm
new}=\frac{f^{2}_{\pi}}{48}\left(m^{2}_{\pi}+2m^{2}_{K}\right)\Tr{U+U^{\dagger}-2}+\frac{\sqrt{3}}{24}f^{2}_{\pi}\left(m^{2}_{\pi}-m^{2}_{K}\right)\Tr{\lambda_{8}\left(U+U^{\dagger}\right)},$

where $\lambda_{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:

$S_{WZ}=-i\frac{N_{c}}{240\pi^{2}}\int
d^{5}x\>\epsilon^{\mu\nu\alpha\beta\gamma}\Tr{L_{\mu}L_{\nu}L_{\alpha}L_{\beta}L_{\gamma}}.$
(4.4)

The onset of kaon condensation in the Skyrme model takes place at a critical
density $n_{\rm cond}$ at which $\mu_{e}$ becomes greater than the energy of
the kaon zero-momentum mode (s-wave condensate). Thus, for baryon densities
$n\geq n_{\rm 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), $\expectationvalue{K^{\pm}}$ ) as a homogeneous field whose time
dependence is given by:

$\expectationvalue{K^{\mp}}=\phi e^{\mp i\mu_{K}t}$ (4.5)

The real constant $\phi$ 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 $\mu_{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

$\tilde{\mathcal{D}}=\matrixquantity(0&0&\phi e^{i\mu_{K}t}\\\ 0&0&0\\\ \phi
e^{-i\mu_{k}t}&0&0)$ (4.6)

which results from substituting the kaon fields in $\mathcal{D}$ as defined in
(4.2) by their corresponding vev in the kaon condensed phase. Also, taking
advantage of the property $\mathcal{D}^{3}=\phi^{2}\mathcal{D}$, we may write
the $SU(3)$ element generated by $\tilde{\mathcal{D}}$ explicitly in matrix
form:

$\Sigma=e^{i\tfrac{\sqrt{2}}{f_{\pi}}\tilde{\mathcal{D}}}=\matrixquantity(\cos\tilde{\phi}&0&ie^{i\mu_{K}t}\sin\tilde{\phi}\\\
0&1&0\\\ ie^{-i\mu_{K}t}\sin\tilde{\phi}&0&\cos\tilde{\phi})$ (4.7)

where $\tilde{\phi}=\tfrac{\sqrt{2}}{f_{\pi}}\phi$ 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=\Sigma
U_{\pi}\Sigma$, where $U_{\pi}$ 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:

$S_{Sk}(U)+S_{\rm WZW}(U)=S_{Sk}(U_{\pi})-\int dtV_{K}(\tilde{\phi}),$ (4.8)

where

$V_{K}=\frac{1}{24\pi^{2}}\int
d^{3}x\Big{[}V^{(2)}_{K}+V^{(4)}_{K}+V^{(6)}_{K}+V^{(0)}_{K}\Big{]}+V^{(WZW)}_{K},$
(4.9)

with

$\displaystyle
V_{K}^{(0)}=2\frac{m^{2}_{K}}{f^{2}_{\pi}e^{2}}(1+\sigma)\sin^{2}\tilde{\phi}$
(4.10) $\displaystyle
V_{K}^{(2)}=\mu_{K}^{2}\sin^{2}\tilde{\phi}[(1+\sigma^{2}+\pi^{2}_{3})\sin^{2}\tilde{\phi}-2(1+\sigma\cos^{2}\tilde{\phi})].$
(4.11) $\displaystyle
V_{K}^{(4)}=-2\mu_{K}^{2}\sin^{2}{\tilde{\phi}}\big{\\{}(1+\sigma)\partial_{i}n^{2}\cos^{2}{\tilde{\phi}}+2[\partial_{i}\sigma^{2}(1-\pi_{3}^{2})+\partial_{i}\pi_{3}^{2}(1-\sigma^{2})+2\sigma\pi_{3}\partial_{i}\sigma\partial_{i}\pi_{3}]\sin^{2}\tilde{\phi}\big{\\}}$
(4.12) $\displaystyle
V_{K}^{(6)}=-\lambda^{2}f^{2}_{\pi}e^{4}\mu_{K}^{2}\sin^{4}(\tilde{\phi})(\partial_{i}\pi_{3}\partial_{j}\sigma-\partial_{i}\sigma\partial_{j}\pi_{3})^{2},$
(4.13) $\displaystyle V_{K}^{\rm(WZW)}=-\mu_{K}N_{c}B_{\rm
cell}\sin^{2}{\tilde{\phi}}.$ (4.14)

###  4.2  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 $\beta$ 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=\Sigma U_{\pi}\Sigma$:

$U\rightarrow\tilde{U}\equiv A(t)UA^{\dagger}(t),$ (4.15)

where $A$ is an element of $SU(3)$ modelling an isospin rotation,

$A=\matrixquantity(a&0\\\ 0&1),\quad a\in SU(2).$ (4.16)

in analogy with the $SU(2)$ case, we define the isospin angular velocity
$\vec{\omega}$ as $A^{\dagger}\dot{A}=\tfrac{i}{2}\omega_{a}\lambda_{a}$
($a=1,2,3$), with $\lambda_{\textsc{a}}$ the Gell-Mann matrices generating
$SU(3)$ for $\textsc{a}=1,\cdots 8$. Notice that $\vec{\omega}$ is a three-
vector, since $A^{\dagger}\dot{A}$ belongs to the isospin $\mathfrak{su}(2)$
subalgebra of $\mathfrak{su}(3)$. Then, we may write the time component of the
Maurer-Cartan current as
$\tilde{U}^{\dagger}\partial_{0}\tilde{U}=AL_{0}A^{\dagger}+AT_{a}A^{\dagger}\omega_{a}$,
where $T_{a}$ is the $\mathfrak{su}(3)$-valued current:

$T_{a}=\frac{i}{2}U^{\dagger}[\lambda_{a},U]\equiv
iT_{a}^{\textsc{a}}\lambda_{\textsc{a}},$ (4.17)

where we have made use of the parametrization (4.2).

Now the kinetic isorotational energy in this case takes the form

$T=\frac{1}{2}\omega_{a}\Lambda_{ab}\omega_{b}+\Delta_{a}\omega_{a}-V_{K}$
(4.18)

where $\Lambda_{ab}$ is the isospin inertia tensor and $\Delta_{a}$ is the
kaon condensate isospin current, given by

$\displaystyle\Lambda_{ab}$
$\displaystyle=\int\left\\{2a\Tr\\{T_{a}T_{b}\\}-4b\Tr\\{[T_{a},L_{k}][T_{b},L_{k}]\\}-\frac{c}{32\pi^{4}}\varepsilon^{lmn}\Tr\\{T_{a}L_{m}L_{n}\\}\varepsilon_{lrs}\Tr\\{T_{j}L_{r}L_{s}\\}\right\\}\,d^{3}x,$
(4.19) $\displaystyle\Delta_{a}$
$\displaystyle=\int\left\\{2a\Tr\\{L_{0}T_{a}\\}-4b\Tr\\{[T_{a},L_{k}][L_{0},L_{k}]\\}-\frac{c}{32\pi^{4}}\varepsilon^{lmn}\Tr\\{L_{0}L_{m}L_{n}\\}\varepsilon_{lrs}\Tr\\{T_{a}L_{r}L_{s}\\}\right\\}\,d^{3}x,$
(4.20)

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.
$\Lambda^{\rm{crystal}}_{ab}=\Lambda\delta_{ab}$. However, the presence of a
kaon condensate breaks this symmetry to a $U(1)$ subgroup, so that
$\Lambda_{ab}$ presents two different eigenvalues in the condensate phase,
$\Lambda_{\rm cond}=\text{diag}(\Lambda,\Lambda,\Lambda_{3})$. Similarly,
$\Delta_{a}=0$ in the purely barionic phase, and its third component acquires
a non-zero value in the condensate phase, $\Delta_{\rm cond}=(0,0,\Delta)$.
The explicit expressions for $\Lambda_{3}$ and $\Delta$ in the condensed phase
can be found in [49]. One can check that in the non-condensed phase, $\phi=0$
and the results of the previous section are recovered, namely,
$\Lambda_{3}=\Lambda$, $\Delta=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}=\Lambda_{3}\omega_{3}+\Delta$.

Thus, after a Legendre transformation to rewrite (4.18) in Hamiltonian form,
and making the $N\rightarrow\infty$ approximation, one can write the quantum
energy correction per unit cell of the crystal in the kaon condensed phase as

$E_{\rm{quant}}=\frac{1}{2\Lambda_{3}}(I_{3}^{2}-\Delta^{2}).$ (4.21)

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
$\Lambda_{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:

$n\leftrightarrow p+K^{-},\qquad l\leftrightarrow K^{-}+\nu_{l}$ (4.22)

such that the chemical equilibrium conditions
$\mu_{n}=\mu_{p}+\mu_{K},\quad\mu_{l}=\mu_{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:

$E=E_{\rm class}+E_{\rm
iso}(\gamma,\tilde{\phi})+E_{K}(\mu_{e},\tilde{\phi})+E_{e}(\mu_{e})+\Theta(\mu_{e}^{2}-m_{\mu}^{2})E_{\mu}(\mu_{e})$
(4.23)

The kaon contribution is the effective potential energy

$E_{K}(\mu_{e},\tilde{\phi})=V_{K}-\frac{\Delta^{2}}{2\Lambda_{3}},$ (4.24)

which depends on the condensate $\tilde{\phi}$ and on the lepton chemical
potential through the explicit dependence on $\mu_{K}$ of both $V_{K}$ and
$\Delta$, and $\mu_{K}=\mu_{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

$\Omega=E-\mu_{e}(N_{e}+\Theta(\mu_{e}^{2}-m_{\mu}^{2})N_{\mu}-\gamma B)$
(4.25)

with respect to $\gamma$, $\tilde{\phi}$ and $\mu_{e}$, i.e.

$\partialderivative{\Omega}{\gamma}\evaluated{}_{n_{B}}\\!\\!\\!(\gamma,\tilde{\phi},\mu_{e})=\partialderivative{\Omega}{\tilde{\phi}}\evaluated{}_{n_{B}}\\!\\!\\!(\gamma,\tilde{\phi},\mu_{e})=\partialderivative{\Omega}{\mu_{e}}\evaluated{}_{n_{B}}\\!\\!\\!(\gamma,\tilde{\phi},\mu_{e})=0.$
(4.26)

The first equation imposes the expected condition
$\mu_{e}=\mu_{I}=2\hbar^{2}(1-2\gamma)/\Lambda_{3}$. Then, after substituting
into the other two conditions we get:

$\displaystyle\gamma n_{B}-$
$\displaystyle\frac{(\mu^{2}_{I}-m^{2}_{e})^{3/2}+(\mu^{2}_{I}-m^{2}_{\mu})^{3/2}}{3\pi^{2}\hbar^{3}}+\frac{n_{B}}{4}\partialderivative{E_{K}}{\mu_{e}}\evaluated{}_{\mu_{e}=\mu_{I}}=0,$
(4.27)
$\displaystyle\partialderivative{V_{K}}{\tilde{\phi}}-\frac{\Delta}{\Lambda_{3}}\partialderivative{\Delta}{\tilde{\phi}}+\partialderivative{\Lambda_{3}}{\tilde{\phi}}\left(\frac{\Delta^{2}}{2\Lambda_{3}^{2}}-\frac{\mu_{I}^{2}}{2\hbar^{2}}\right)=0,$
(4.28)

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 $\gamma$ and $\tilde{\phi}$ 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\mu$ and $npe\mu\overline{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_{\rm
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\tilde{\phi}$ from the second equation and setting $\tilde{\phi}=0$. Then
we may see the system as a pair of equations to obtain the values of
$\gamma_{\text{cond}}$ and $L_{\text{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.

label | $f_{\pi}$ | $e$ | $\lambda^{2}$ | $E_{0}$ | $n_{0}$ | $S_{0}$ | $L_{\text{sym}}$ | $n_{\text{cond}}/n_{0}$
---|---|---|---|---|---|---|---|---
set 1 | 133.71 | 5.72 | 5 | 920 | 0.165 | 23.5 | 29.1 | 2.3
set 2 | 138.11 | 6.34 | 5.78 | 915 | 0.175 | 24.5 | 28.3 | 2.2
set 3 | 120.96 | 5.64 | 2.68 | 783 | 0.175 | 28.7 | 38.7 | 1.6
set 4 | 139.26 | 5.61 | 2.74 | 912 | 0.22 | 28.6 | 38.9 | 1.6
Table 6: Sets of parameter values and observables at nuclear saturation

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.

##  5 Neutron Stars

###  5.1 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 a spherically symmetric ansatz of the spacetime metric,

$ds^{2}=-A(r)dt^{2}+B(r)dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\varphi^{2}),$ (5.1)

in the Einstein equations,

$R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=8\pi GT_{\mu\nu}.$ (5.2)

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,

$T_{\mu\nu}=(\rho+p)u_{\mu}u_{\nu}+pg_{\mu\nu},$ (5.3)

where the pressure $p$ and the energy density $\rho$ 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.

$\displaystyle\frac{dA}{dr}$ $\displaystyle=Ar\left(8\pi
GBp-\frac{1-B}{r^{2}}\right),$ (5.4) $\displaystyle\frac{dB}{dr}$
$\displaystyle=Br\left(8\pi GB\rho+\frac{1-B}{r^{2}}\right),$ (5.5)
$\displaystyle\frac{dp}{dr}$ $\displaystyle=-\frac{p+\rho}{2A}\frac{dA}{dr}.$
(5.6)

We use a 4th order Runge-Kutta method of step $\Delta 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,

$B(r=R)=\frac{1}{(1-\frac{2GM}{R})}.$ (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.

###  5.2 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 $\mathcal{L}_{24}$ model lead to
too small maximum NS masses [36], whereas the EOS based on the
$\mathcal{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
$\mathcal{L}_{24}$ crystal EOS at intermediate densities and the
$\mathcal{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 ($\mathcal{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

$E_{24}(l)=E_{0}\left[a\left(\frac{l}{l_{0}}+\frac{l_{0}}{l}\right)+b\right],$
(5.8)

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_{0}^{-3}=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_{\pi}$ and $e$. From
eq. 5.8 we may obtain the pressure and the energy density using eq. 2.15,

$p_{24}(l)=a\frac{E_{0}}{3l^{2}}\left(\frac{l_{0}}{l^{2}}-\frac{1}{l_{0}}\right),$
(5.9)

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 ($\rho_{6}=p$), with a
speed of sound ($c^{2}_{s}=\partial p/\partial\rho$) equal to 1. The BPS
Skyrme model, even with a general potential term (with no derivatives of $U$),
$\mathcal{L}_{0}=-\mu^{2}\mathcal{U}$, retains this property, besides its
stress-energy tensor is of the perfect fluid form from which we may easily
obtain the EOS.

$\displaystyle T^{\mu\nu}_{\rm BPS}=(p+\rho)u^{\mu}u^{\nu}+pg^{\mu\nu},\quad
u^{\mu}=\frac{B^{\mu}}{\sqrt{g_{\alpha\beta}B^{\alpha}B^{\beta}}},$ (5.10)
$\displaystyle
p=\lambda^{2}\pi^{4}g_{\alpha\beta}B^{\alpha}B^{\beta}-\mu^{2}\mathcal{U},\quad\rho=p+2\mu^{2}\mathcal{U}.$
(5.11)

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$, $\rho$ and $n_{B}$ can be obtained

$p+\rho=2\lambda^{2}\pi^{4}n_{B}^{2}.$ (5.12)

The dependence on the Skyrme field of the potential $\mathcal{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,

$\rho=\rho_{6}+\rho_{0}=p+\rm{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\sim 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_{\rm max}\approx
1.5M_{\odot}$) compared with the most recent measurements of NS ($M_{\rm
max}\geq 2M_{\odot}$). Later, the BPS Skyrme model was also coupled to gravity
[73], using different potentials, and the NS masses were too large ($M_{\rm
max}\approx 3-4M_{\odot}$) 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\sim 2-2.5M_{\odot}$).

In [27] the parameters of the BPS submodel ($\lambda^{2}$, $\mu^{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,

$\rho_{\rm Gen}(p\ll p_{PT})=\rho_{\rm Sk},\quad\rho_{\rm Gen}(p\gg
p_{PT})=p+\text{const.}$ (5.14)

The EOS $\rho_{\rm 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,

$\displaystyle\rho_{\rm Gen}(p)=$
$\displaystyle(1-\alpha(p,p_{PT},\beta))\rho_{\rm
Sk}+\alpha(p,p_{PT},\beta)(p+\rho_{\rm Sk}(p_{PT})),$ (5.15)
$\displaystyle\alpha(p,p_{PT},\beta)=\frac{(p/p_{PT})^{\beta}}{1+(p/p_{PT})^{\beta}}.$

For this interpolation, the value of $\beta$ measures how fast the transition
is. Here we take the value $\beta=0.9$ since depending on the value of
$p_{PT}$ a larger value of $\beta$ 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_{PT}$ 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
$\alpha(p)$ function, taking $\beta=2$ and replacing $p_{PT}$ by $p_{*}$, we
called the combination of the two models the Hybrid EOS,

$\rho_{\rm Hyb}=(1-\alpha(p,p_{*},2))\rho_{\rm
BCPM}+\alpha(p,p_{*},2)\rho_{\rm 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_{PT}=50$ MeV/fm3 and the Hybrid model,
with the same $p_{PT}$ and $p_{*}=2$ MeV/fm3. We also show some usually
considered EOS: APR4, WWF1, BCPM and we fill in orange and blue the regions of
allowed values for $p_{PT}$ and $p_{*}$. This figure was originally published
in [50]

We sample the space of parameters $(p_{PT},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_{\odot}=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 ($\sim 1.4M_{\odot}$) region, also the heavier pulsar observations
PSR J1614 - 2230 ($1.928\pm 0.017M_{\odot}$) [77], PSR J0348 + 0432 ($2.01\pm
0.04M_{\odot}$) [78], PSR J0740 + 6620 ($2.14^{+0.10}_{-0.09}M_{\odot}$) [79],
and some additional constraints from NICER, chiral EFT and multimessenger
observations [80]. As expected, it is found that $p_{PT}$ 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_{PT}\in[25,50]$ MeV/fm3 and $p_{*}\in[0.5,2]$ MeV/fm3.
We even allowed for maximal masses $\sim 2.7M_{\odot}$ 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_{\odot}$) [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.

Figure 17: Mass-Radius curves of the Hybrid model with the values
$p_{PT}=\\{25,40,50\\}$, $p_{*}=\\{0.5,1,2\\}$ MeV/fm3. We compare our results
with those of other EOS and constraints obtained from different analysis. This
figure was originally published in [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 $\widetilde{\Lambda}$.
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$.

Figure 18: The effective tidal deformability $\widetilde{\Lambda}$ of the
Hybrid model against the mass ratio. The orange regions show the PDFs obtained
from [84], and we adopt the notation Hyb_$p_{PT}$_$p_{*}$ for the curves. This
figure was originally published in [50].

The estimated value of $\widetilde{\Lambda}$ 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.

###  5.3 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, 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.

####  5.3.1 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, it means that each of the terms in the energy
functional scales with $L$ independently as $E_{i}\propto 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.

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.

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. 1.

Due to the perfect scaling property, the (adimensional) energy and isospin
moment of inertia approximately satisfy the following expressions,

$\displaystyle E(L,c_{6},c_{0})$
$\displaystyle=K_{2}L+\frac{K_{4}}{L}+c_{6}\frac{K_{6}}{L^{3}}+c_{0}K_{0}L^{3},$
(5.17) $\displaystyle\Lambda(L,c_{6},c_{0})$
$\displaystyle=\Lambda_{2}L^{3}+\Lambda_{4}L+c_{6}\frac{\Lambda_{6}}{L},$
(5.18)

and the values of $K_{i}$ and $\Lambda_{i}$ are ”universal”, i.e., they do not
depend on the parameters or on $L$. Hence, these values are obtained from the
contribution of each term independently, for only one choice of the parameters
and for a single value of the unit cell length $L$. For simplicity, we set
$c_{0}$ and $c_{6}$ equal to one. The values of the universal scaling
constants are shown in table 7.

n | 0 | 2 | 4 | 6
---|---|---|---|---
$K_{\rm n}$ | 0.034 | 0.466 | 9.617 | 4.329
$\Lambda_{\rm n}$ | – | 0.038 | 1.393 | 0.883
Table 7: Perfect scaling parametrization constants

  2. 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_{\pi}$.

  3. 3.

Then, we calculate $L_{0}$ by minimizing (5.17) and the values of
$E_{0},n_{0},S_{0}$ and $L_{\rm sym}$ for different pairs of values
$(e,\lambda^{2})$.

  4. 4.

When we find a set of parameters $(f_{\pi},e,\lambda^{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. 5.

Finally, we accept the sets of values that satisfy the constraints, $M_{\rm
max}\geq 2M_{\odot}$ and $R_{1.4M_{\odot}}\leq 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.

Figure 20: Scan of the parameters $(e,\lambda^{2})$ for a fixed value of
$E_{s}$. Each coloured line corresponds to a fixed value of $\lambda^{2}$ and
the black horizontal lines are the experimental values of $E_{0},n_{0}$ and
$S_{0}$. The yellow dots are the optimal values and they lie in the black
vertical line which is the corresponding value of $e$.

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.

Figure 21: Symmetry energy of Skyrme crystals as a function of the density. We
show 23 different curves from the scanned values. The shaded regions
constrains the symmetry energy at sub-saturation [87] and supra-saturation
[61] densities.

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 $\lambda^{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 $(\rho,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.

Figure 22: 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\mu$
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.

In a more extensive analysis of the parameters, we found $\sim 10,000$
accepted sets of parameters. In order to obtain these values we made the scan
with the steps: $\Delta E_{s}=5$ MeV, $\Delta e=0.01$,
$\Delta\lambda^{2}=0.01$ MeV fm3. As briefly mentioned before, the constraints
on the symmetry energy yield rather stringent upper bounds on the sextic term
coupling constant, we find that $\lambda^{2}\lesssim 3.4$ MeV fm3.
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 $\lambda^{2}\gtrsim 0.29$ MeV fm3.

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\lesssim
2.3M_{\odot}$, 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}\sim 2-2.5M_{\odot}$. 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_{\odot}}\sim 12-13$ km. However a final comment about the radii of
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_{\rm
max}$ and $R_{1.4M_{\odot}}$ 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 [Abbott_2017] 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).

Figure 23: MR curves for the 23 representative sets of parameters considered.
The colours of the lines represent the same as in fig. 22. The shaded regions
correspond to GW (blue and orange) and pulsar (green) constraints.

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 $\sim
2.5M_{\odot}$, the radius increases when we consider the full model. The low-
mass region ($\sim 1.4M_{\odot}$) 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_{\odot}$ 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.

###  5.4 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\mu$
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.

Figure 24: Left plot: energy against the side length of the crystal,
calculated with more points near the condensation values for both branches and
their interpolations. Right plot: pressure against the energy density from
which we conclude that there is a first order phase transition. All magnitudes
are shown in Skyrme units.

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\mu$ matter, the Gibbs
conditions for the phase equilibrium,

$p^{\rm I}=p^{\rm II},\hskip 5.69054pt\mu^{\rm I}_{i}=\mu^{\rm II}_{i},\hskip
4.2679pti=B,q$ (5.19)

cannot be both satisfied in a standard MC. In the last expression $\mu_{B}$
and $\mu_{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:

$\mu_{i}=B_{i}\mu_{B}+q_{i}\mu_{q},$ (5.20)

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 $\chi$ of the kaon condensed phase, charge neutrality is
imposed in the GC as:

$n^{MP}_{q}=(1-\chi)n^{\rm I}_{q}+\chi n^{\rm II}_{q}=0.$ (5.21)

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. 5.19, 5.21 and 4.28.

Figure 25: $E(L)$ curves for the two phases. The different slopes at the point
of phase separation indicate a first-order phase transition. We also show the
curves resulting from a Maxwell construction (MC) and a Gibbs construction
(GC).

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.

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].

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