Date: ]

# Quantum Skyrmion crystals and the symmetry energy of dense matter

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

([)

###### Abstract

The canonical quantization method for collective coordinates in crystalline
configurations of the generalized Skyrme model is applied in order to find the
quantum ground state of Skyrmion crystals and study the quantum corrections to
the binding energy resulting from the isospin degrees of freedom. This leads
to a consistent description of asymmetric nuclear matter within the Skyrme
framework and allows us to compute the symmetry energy of the Skyrmionic
crystal as a function of the baryon density, and to compare with recent
observational constraints.

###### Contents

  1. I Introduction
  2. II Generalized Skyrme model and classical crystals
    1. II.1 The model
    2. II.2 The classical crystal of Skyrmions
    3. II.3 The Skyrme crystal EoS for symmetric nuclear matter
  3. III Quantization of the Skyrme crystal
    1. III.1 The internal symmetry of Skyrme crystals and the isospin group
    2. III.2 Quantum isospin states and Hilbert space
    3. III.3 Isospin correction to the energy per baryon
  4. IV Quantum skyrmion crystals and the symmetry energy of $npe(\mu)$ matter
    1. IV.1 The electric charge density of pure Skyrmion matter
    2. IV.2 The symmetry energy
    3. IV.3 Particle fractions of $npe\mu$ matter in $\beta$-equilibrium
  5. V Conclusions
  6. A Calculation of matrix elements
  7. B Isospin inertia tensor of Skyrme crystals

##  I Introduction

Although the strong interactions are, in principle, completely described by
Quantum Chromodynamics (QCD), a derivation of the properties of nucleons,
atomic nuclei and nuclear matter from this fundamental theory is still one of
the biggest challenges of theoretical physics. This is related to the
nonperturbative nature of the low energy regime where the standard expansion
in terms of Feynman diagrams loses its applicability. As a consequence, other
approaches like lattice QCD or effective phenomenological models, usually
tailored for certain regimes, must be applied.

The Skyrme model Skyrme (1962) provides a very attractive field theoretical
framework for this purpose. On the one hand, it is very well anchored in the
underlying fundamental quantum theory Witten (1979) and, therefore, provides a
rare opportunity to study the non-perturbative features of strongly
interacting matter at all scales, from single baryons and nuclei to neutron
stars. Importantly, the model contains a rather limited number of degrees of
freedom which, in the simplest version, are simply pions disguised into an
$SU(2)$ valued matrix field $U$, and a small number of terms in the action,
which translates into an equally small number of free parameters. Then
nucleons and atomic nuclei emerge as collective excitations in such a mesonic
fluid. Mathematically, they are topological solitons called Skyrmions, whose
topological degree can be rigorously identified with the baryon number
Balachandran  _et al._ (1982). Moreover, the relevant quantum numbers are
introduced via an appropriate quantization procedure of the corresponding
collective coordinates.

In the last years, there has been significant progress in the application of
the Skyrme model to the description of atomic nuclei. The first source of this
development is the use of the vibrational quantization. Here the Hilbert state
is built not only on the zero modes, as in the standard rigid rotor
quantization Adkins  _et al._ (1983), but also the lightest massive
deformations are take into account. This approach elevated the Skyrme model to
a quantitative tool for the understanding of excitation bands of light nuclei
Battye  _et al._ (2009); Lau and Manton (2014); Halcrow  _et al._ (2017). In
addition, it has very recently been shown how the spin-orbit interaction
leading to a phenomenologically consistent nucleon-nucleon force emerges in
the Skyrme model Halcrow and Harland (2020). Secondly, it is now well
understood how to reduce the typically quite large binding energies of
Skyrmions, which translate into unphysically large binding energies of atomic
nuclei. This requires the addition of new, physically motivated terms, e.g.,
the so-called sextic term Adam  _et al._ (2010), generalized potentials
Gillard  _et al._ (2015); Gudnason (2016) or new degrees of freedom, e.g.,
vector mesons Naya and Sutcliffe (2018).

In a parallel development, the Skyrme model coupled to gravity has been
applied to the investigation of neutron stars and their properties. Indeed, in
the last few years, astrophysical observations Miller  _et al._ (2021); Riley
_et al._ (2021) and the analysis of gravitational wave measurements Abbott
_et al._ (2017, 2020a, 2020b) have led to a big improvement in our
understanding of the nuclear matter equation of state (EoS), especially above
the saturation density. From the Skyrme model perspective, again the role of
the sextic term is crucial, as this part of the action governs the high
pressure and density regime Adam  _et al._ (2015a). Indeed, it makes
Skyrmionic matter much stiffer at extreme conditions, which results in
physically acceptable values of the maximal masses of neutron stars. This has
been first found in the case where the EoS was motivated by certain limits of
the Skyrme model Adam  _et al._ (2020). Later on, it has been confirmed in the
full Skyrme model computation Adam  _et al._ (2022). Remarkably, some Skyrme
model-based EoS lead to observables which pass the current available
observational data, concerning, e.g., mass-radius curves, tidal
deformabilities and quasi-universal relations between the moment of inertia,
Love numbers and the quadrupole moment of slowly rotating stars Adam  _et al._
(2021).

The starting point in these computations consists in the numerical derivation
of the lowest energy periodic Skyrmion solution representing infinite nuclear
matter, at a given density, using a variational approach Adam  _et al._
(2022). Next, one can find the corresponding EoS by studying how the energy of
such solutions depends on the density. Once the EoS is known, the standard
Tolman-Oppenheimer-Volkoff (TOV) approach Tolman (1939), Oppenheimer and
Volkoff (1939) allows us to study the resulting NS. However, all these results
have been obtained only in the purely classical limit, where no quantum
corrections are taken into account and, importantly, the proton and the
neutron are described by the same classical solution, which preserves isospin
symmetry. Hence, classical skyrmion crystals would correspond to the so-called
symmetric nuclear matter. However, realistic models of nuclear matter inside
neutron stars do not predict an isospin-symmetric state at all. Instead, only
a small fraction of protons over total nucleons is generally allowed. Indeed,
for a nuclear system with total number of baryons $A=N+Z$, where $N$ is the
number of neutrons and $Z$ the number of protons, and defining the isospin
asymmetry parameter as $\delta=(N-Z)/A=(1-2\gamma)$, with $\gamma$ the proton
fraction or ratio between $Z$ and $A$, the binding energy is usually defined
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}},$
(1)

where the term $E_{N}(n_{B})$ would be the binding energy of isospin-symmetric
matter, and $S_{N}(n_{B})$ is the so-called _symmetry energy_. The symmetry
energy is thus a measure of 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. However, its dependence on the density has proven
difficult to measure experimentally, and usually it is 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$ (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)

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 Fiorella Burgio and Fantina (2018), 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 Essick  _et al._ (2021); Tang
_et al._ (2021); de Tovar  _et al._ (2021); Gil  _et al._ (2021); Li  _et al._
(2021).

The main aim of the current work is to carry the description of Skyrmionic
matter beyond the classical realm by semiclassically quantizing the isospin
DoF of the Skyrme crystal. This allows us to calculate the isospin moment of
inertia (hence the symmetry energy) of a given classical solution, taking into
account the electric charge neutrality of the crystal, for which the
introduction of a neutralising charged lepton background is crucial. The
effect of this background is to introduce a certain isospin asymmetry due to
$\beta$ equilibrium with the baryonic matter. Finally, we find the EoS of the
full Skyrmion crystal plus lepton background system, which models $npe\mu$
matter, and, in particular, the particle fraction of each component as a
function of density.

##  II Generalized Skyrme model and classical crystals

###  II.1 The model

The generalized Skyrme model is defined by the following Lagrangian density

$\displaystyle\mathcal{L}=-\frac{f^{2}_{\pi}}{16}$
$\displaystyle\Tr\\{L_{\mu}L^{\mu}\\}+\frac{1}{32e^{2}}\Tr\\{\left[L_{\mu},L_{\nu}\right]^{2}\\}$
$\displaystyle-\lambda^{2}\pi^{4}\mathcal{B}_{\mu}\mathcal{B}^{\mu}+\frac{m^{2}_{\pi}f^{2}_{\pi}}{8}\Tr\left(U-I\right),$
(4)

where $L_{\mu}=U^{\dagger}\partial_{\mu}U$ is the left invariant Maurer-Cartan
current and the matrix field decomposes as

$U=\sigma+i\pi_{k}\tau_{k}.$ (5)

Here, $\pi_{k}$ ($k$ = 1, 2, 3) are the pions and $\tau_{k}$ are the Pauli
matrices. The unitarity of the matrix field puts a condition on the fields.
Namely, $\sigma^{2}+\pi_{i}\pi_{i}=1$. Furthermore, $\mathcal{B}^{\mu}$ is the
conserved topological current which, in the standard manner, defines the
topological index of maps $U$, i.e., the baryon charge $B$

$B=\int d^{3}x\mathcal{B}^{0},\hskip
5.69054pt\mathcal{B}^{\mu}=\frac{1}{24\pi^{2}}\epsilon^{\mu\nu\alpha\beta}\Tr\left\\{L_{\nu}L_{\alpha}L_{\beta}\right\\}.$
(6)

It is a remarkable property of the generalized Skyrme effective model that it
contains only four terms, which translate into four coupling constants
$f_{\pi},m_{\pi},e,\lambda$. Two of them have a direct phenomenological
interpretation as the pion decay constant and the pion mass. In addition,
$\lambda$ can be related to a ratio between the mass and the coupling constant
of the $\omega$ meson. We assume the physical mass of the pions, $m_{\pi}=140$
MeV, in the whole paper. The parameters of the standard Skyrme model are set
to

$f_{\pi}=129\>\text{MeV}\,,\quad e=5.45.$ (7)

These values of $f_{\pi}$ and $e$ are often used in the Skyrme model
literature Castillejo  _et al._ (1989); Baskerville (1996), because they were
obtained in Adkins  _et al._ (1983) by fitting the $B=1$ solution of the
standard Skyrme model to the masses of the proton and the $\Delta$ excitation.
Our main reason to choose them is that this choice facilitates the comparison
of our calculations with previous results. Finally, we will consider different
values for the coupling constant $\lambda^{2}$ multiplying the third term in
the action, containing six powers of derivatives. Although often omitted in
the context of light nuclei, this term is essential when one studies the
properties of nuclear matter at high density, which is a natural environment
in the core of neutron stars. Indeed, this sextic term governs the equation of
state at this regime and asymptotically leads to the maximally stiff EoS Adam
_et al._ (2015b).

###  II.2 The classical crystal of Skyrmions

In the Skyrme model, infinite skyrmionic matter is described by a periodic
minimizer of the static energy

$\displaystyle E=\frac{1}{24\pi^{2}}$ $\displaystyle\int
d^{3}x\left[-\frac{1}{2}\Tr\\{L_{i}L_{i}\\}-\frac{1}{4}\Tr\\{\left[L_{i},L_{j}\right]^{2}\\}\right.$
$\displaystyle\left.+8\lambda^{2}\pi^{4}\frac{f^{2}_{\pi}e^{4}}{\hbar^{3}}(\mathcal{B}^{0})^{2}+\frac{m^{2}_{\pi}}{f^{2}_{\pi}e^{2}}\Tr\left(I-U\right)\right]$
(8)

and, therefore, it is usually referred to as the Skyrme crystal. For numerical
purposes, we have adopted Skyrme units, in which energy and length are
measured in units of $3\pi^{2}f_{\pi}/e$ and $\hbar/(f_{\pi}e)$, respectively.
Obviously, while the total energy of the crystal is infinite, the energy per
baryon number remains finite

$\frac{E}{B}=\frac{N_{\text{cells}}\>E_{\text{cell}}}{N_{\text{cells}}\>B_{\text{cell}}}=\frac{E_{\text{cell}}}{B_{\text{cell}}}.$
(9)

Here, $N_{\text{cells}}$ is the number of cells and $E_{\text{cell}}$,
$B_{\text{cell}}$ are the energy and baryon charge in a single, periodic cell.
The energy of the unit cell strongly depends on the assumed geometry and its
size, characterized by the length parameter $L$. Concretely, $L$ is the
distance between nearest-neighbor skyrmions in the maximally attractive
channel. As a consequence, the resulting field configuration is not periodic
in $L$. The period length and the size of the unit cell is $2L$, instead.

Comparing various geometries, i.e., types of crystals, at a particular volume
of the cell, one can find the ground state crystalline solution at a given
density. Note that a particular geometry of the crystal translates into
particular symmetries of the Skyrme field. In practice, it is impossible to
check all possible geometries and one constrains considerations to well
motivated cubic crystals. These are the simple cubic (SC) and face centered
(FCC) crystals of Skyrmions as well as the body centered (BCC) and face
centered (FCC) crystals of half-Skyrmions.

As the cubic symmetries are shared by all these geometries, the following
relations must be obeyed by these crystals,

$\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}),$
(10) $\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}).$
(11)

In this work we will focus on the FCC crystal of half-skyrmions, which has 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}),$
(12) $\displaystyle\text{D}_{4}$ $\displaystyle:(x,y,z)\rightarrow(x+L,y,z),$
(13)
$\displaystyle(\sigma,\pi_{1},\pi_{2},\pi_{3})\rightarrow(-\sigma,-\pi_{1},\pi_{2},\pi_{3}).$
(14)

A more detailed description of the construction of the Skyrme crystal and the
comparison of different symmetries can be found in Adam  _et al._ (2022), and
we have kept the same notation for this work. As in that previous work, the
unit cell has size $2L$ and a baryon content of $B_{\text{cell}}=4$. Then, for
each value of $L$ we obtain the minimum of energy as explained in Adam  _et
al._ (2022). It turns out that the energy-size curve, $E_{\text{cell}}(L)$, is
a convex function which has a minimum at a certain $L_{*}$.

###  II.3 The Skyrme crystal EoS for symmetric nuclear matter

For the crystal solutions one can define the relevant thermodynamical
quantities in the usual way, that is, energy density $\rho$, pressure $p$ and
baryon charge density $n_{B}$

$\displaystyle\rho$
$\displaystyle=\frac{E}{V}=\frac{E_{\text{cell}}}{V_{\text{cell}}},$ (15)
$\displaystyle p$ $\displaystyle=-\frac{\partial E}{\partial
V}=-\frac{\partial E_{\text{cell}}}{\partial V_{\text{cell}}},$ (16)
$\displaystyle n_{B}$
$\displaystyle=\frac{B}{V}=\frac{B_{\text{cell}}}{V_{\text{cell}}}.$ (17)

Again they are functions of $L$ or, in other words, the volume of the unit
cell $V_{\text{cell}}$.

At the point where $L=L_{*}$, the given crystal solution describes skyrmionic
matter at equilibrium, i.e., at zero pressure. In the region where $L<L_{*}$
the volume of the cell decreases, which corresponds to a squeezed crystal.
This translates into growing pressure and density. The remaining region
$L>L_{*}$, where the volume increases in comparison to the equilibrium, is
thermodynamically unstable. Indeed, it formally gives negative pressure. Due
to that, the low density regime cannot be described by any of the previously
mentioned crystals. On the contrary, it is expected that the crystal is
replaced by an inhomogeneous phase, where lumps of nontrivial energy density
are surrounded by regions of void.

In the absence of the sextic term, the landscape of the crystal energy
minimizers was very well understood. Namely, for $L\leq L_{*}$ the ground
state is formed by the FCC half-Skyrmion phase, which at extremely high
density is replaced by the BCC half-Skyrmions. This phase transition occurs at
densities much beyond the values expected at the cores of neutron stars.

It has been recently shown Adam  _et al._ (2022) that this picture is
significantly modified if the sextic term is added. First of all, the FCC to
BCC phase transition is moved to much smaller densities, approximately 4-5
saturation densities, which can be easily found in the center of heavy NS. In
addition, the appearance of the sextic term introduces a fluidity into the
model which mathematically results from the volume diffeomorphism invariance
of this part of the action. Such a fluidity is visible in a more and more
homogenous distribution of the energy density as the pressure increases. This
does not happen if this term is absent.

However, perhaps the most crucial result is that the equation of state, which
relates energy density and pressure, stiffens in the generalized Skyrme model.
Indeed, the sextic term alone leads to the maximally stiff EoS, $\rho=p$. In
the full, generalized model it occurs asymptotically at high density. In any
case, this stiffening is responsible for a significant rising of the maximal
masses of NS to values which are in accordance with current observations.

These very encouraging findings have been obtained in the classical regime of
the Skyrme model, which, as we already underlined, corresponds to symmetric
nuclear matter. Therefore, it is of vital importance to semiclassically
quantize the isospin DoF of the skyrmionic crystal. This is the first step to
describe neutron matter, which is the ground state of the matter in the core
of NS.

##  III Quantization of the Skyrme crystal

It is well known that the largest quantum correction to the classical energy
of Skyrmion configurations comes from the contribution of the isospin degrees
of freedom, which are usually quantized as zero-modes via some collective
coordinate parametrization. To add the contribution of the quantization of the
global isospin zero modes to the total energy, we need to know the quantum
isospin state of the full crystal. This task, however, becomes impossible in
the thermodynamic limit in which the number of particles forming the crystal
goes to infinity. Instead, we can make the following assumptions on the
quantum wavefuntion of the full crystal:

  * •

The isospin wavefuntion of the total crystal $\ket{\Psi}$ can be written as a
superposition of states constructed from the (infinite) product of isospin
wavefuntions of individual unit cells, $\ket{\Psi}=\bigotimes\limits_{\rm
cells}\ket{\psi}$. In other words, as a first approximation we will not
consider the quantum correlation on the isospin state between cells.

  * •

The symmetry of the classical crystal configuration is inherited by the total
wavefunction, and shared with the wavefunction of each of the unit cells, i.e.
both $\ket{\Psi}$ and $\ket{\psi}$ share the same point symmetry group.

These two assumptions imply that finding the quantum states of the total
crystal is equivalent to finding the state of each unit cell. The latter is in
fact a more plausible task as we may use the tools developed for the
quantization of multi-skyrmion configurations.

###  III.1 The internal symmetry of Skyrme crystals and the isospin group

The full symmetry group of the massless Skyrme Lagrangian is given by the
direct product of the Poincaré and chiral groups. However, we are interested
in solutions that minimize the energy functional (measured on a given
reference frame). The internal symmetry group of such functional is the same
as the Lagrangian, but the Poincaré symmetry is broken to the Euclidean
subgroup corresponding to spatial rotations and translations,
$E_{3}=SO(3)\times\mathbb{R}^{3}$. Thus, the symmetry group of the energy
functional is $\widetilde{G}=E_{3}\times SU(2)_{L}\times SU(2)_{R}\simeq
E_{3}\times SO(4)_{chiral}$. The action of an element of such group on the
Skyrme field is given by

$U(\vec{x})\rightarrow g_{L}U(R_{S}\cdot\vec{x}+\vec{a})g_{R}^{\dagger}$ (18)

where $\vec{a}\in\mathbb{R}^{3}$, $R_{S}\in SO(3)$ represents the spatial
rotations and $g_{L/R}\in SU(2)_{L/R}$ are the left and right-handed chiral
transformations, respectively.

Moreover, the presence of nontrivial boundary conditions imposed on the
relevant field configurations may further reduce their symmetry. Indeed, in
the case of finite energy field configurations a boundary condition of the
form $U(\vec{x})\xrightarrow{|\vec{x}|\rightarrow\infty}\mathbb{1}$ must be
imposed (i.e. the Skyrme field must decay to its vacuum value far from the
center of the soliton). Such boundary condition not only allows us to classify
the field configurations into different topological sectors labeled by their
homotopy class within the third homotopy group of the target space
$\pi_{3}(SU(2))=\mathbb{Z}$, but also reduces the symmetry of such
configurations since the vacuum is only preserved under the subgroup
$G=E_{3}\times\rm{diag}[SU(2)_{L}\times SU(2)_{R}]\simeq E_{3}\times
SU(2)_{I}$, i.e. transformations of the form (18) with $g_{L}=g_{R}=g\in
SU(2)_{I}$. The remaining internal symmetry group is called the _isospin
group_ , since it corresponds to the isospin degrees of freedom when the
solitons (Skyrmions) are identified with baryons and nucleons of low-energy
QCD and nuclear physics.

A further reduction of the symmetry group $G$ may occur on individual
configurations minimizing the energy functional for each topological sector.
For instance, the $B=1$ Skyrmion does present the full group $G$ as a symmetry
of its energy density isosurfaces, whereas the spatial rotations are broken to
$O(2)=SO(2)\times\mathbb{Z}_{2}$ for the $B=2$ Skyrmion, which presents a
toroidal shape. As $B$ increases, the symmetry of the configurations
minimizing the static energy becomes more complicated, and for $B\geq 3$ it is
given by a point group, the $B=3$ Skyrmion presenting tetrahedral symmetry,
the $B=4$ cubic symmetry and so on.

On the other hand, for crystalline configurations, the boundary conditions
imposed on the Skyrme field are of different nature: instead of imposing the
vacuum at large distances, one must impose periodicity conditions on the
boundary of a compact space region- the _unit cell_ \- and hence these
solutions do not yield a finite value when its energy density is integrated
over all space. Instead, this requirement is relaxed to yield a finite value
of the energy over the unit cell, or, equivalently, a finite energy per
baryon.

Therefore, the symmetry group of a crystalline configuration is reduced from
$\widetilde{G}$ due to these periodicity conditions, which in turn result in a
particular point group symmetry for the full crystal. Indeed, in section II.2
we have studied different crystalline structures that present different
symmetries, although all of them are based on a simple cubic symmetry.
However, it is not straightforward to define the isospin subgroup $SU(2)_{I}$
in such configurations, since the vacuum value is no longer imposed at the
boundaries. Instead, one should in principle consider internal rotations of
the full chiral group $SO(4)_{\rm chiral}$, since there is not a natural way
to select the diagonal subgroup corresponding to isospin.

The problem of defining the isospin group in crystal configurations is treated
in Baskerville (1996). As explained there, the procedure of defining the
isospin subgroup in crystalline configurations is subject to some ambiguities,
but the energy per baryon spectrum in the infinite crystal limit is unique.
Indeed, it turns out that the 4-dimensional representation of the cubic point
group of minimal energy crystals is _reducible_ into the trivial 1-D irrep and
a 3-D irrep, which singles out one direction in isospin space. We may then
choose the $\sigma$ field to transform in the trivial irrep, and then to
define the isospin group as the subgroup of isorotations within the 3-D irrep,
i.e. rotations between the three pion fields.

In the rest of this work, we will use this definition of isospin in crystals,
as it is also the natural choice if a symmetry-breaking potential (such a mass
term for the pions) is added to the Skyrme Lagrangian.

###  III.2 Quantum isospin states and Hilbert space

In general, the classical field configuration of a $B>1$ Skyrmion presents a
non-trivial symmetry characterized by a finite point group $H$. Hence, the
corresponding allowed quantum states of spin and isospin will correspond to
different linear combinations of the basis vectors $\ket{ii_{3}k_{3}}$ (see
appendix A). To find these combinations explicitly, the group $H$ has to be
known in terms of a set of generators, $\\{H_{k}\\}$ that can be written as a
product of a rotation $\mathcal{R}_{k}$ and an isorotation
$\mathcal{R^{\prime}}_{k}$. Then, the Finkelstein-Rubinstein (F-R) constraints
may be written as

$\exp{i\alpha\mathbf{n}\mathbf{L}}\exp{i\beta\mathbf{N}\mathbf{K}}\ket{\psi}=\chi_{FR}\ket{\psi},$
(19)

where $\mathbf{L}$ and $\mathbf{K}$ are the (body-fixed) angular and isospin
angular momentum operators, respectively, and $\chi_{FR}=\pm 1$. To
consistently do this, it is important to know as well whether each of the
symmetries $H_{k}$ corresponds to a contractible or a non-contractible loop in
configuration space. Once all the generators $H_{k}$ have been parametrized in
terms of products of rotations and isorotations, one needs to find a solution
of the Finkelstein-Rubinstein system of constraints on each subspace of fixed
$i,j$. Solving the system of F-R constraints is equivalent to finding a set of
common eigenvectors in this subspace of a particular set of matrices.

Indeed, the general transformation law of angular momentum states under a
rotation $\mathcal{R}(\alpha,\beta,\gamma)$ parametrized in terms of the Euler
angles $\alpha,\beta,\gamma$ (in the $Z-Y-Z$ convention) is

$\mathcal{R}(\alpha,\beta,\gamma)\ket{j,m}=\sum\limits_{m^{\prime}}D^{j}_{mm^{\prime}}(\alpha,\beta,\gamma)\ket{jm^{\prime}},$
(20)

where $D^{j}_{mm^{\prime}}(\alpha,\beta,\gamma)$ is the Wigner D-matrix,
corresponding to the irreducible, spin-$j$ representation of the rotation
group.

Hence, for each $H_{k}$, a solution of the corresponding F-R constraint is
given by a state $\ket{\Psi}$ which is both an eigenstate of
$\mathcal{R}_{k}(\alpha,\beta,\gamma)$ with eigenvalue $\lambda_{k}$ and of
$\mathcal{R}^{\prime}_{k}(\alpha^{\prime},\beta^{\prime},\gamma^{\prime})$
with eigenvalue $\lambda^{\prime}_{k}$, and such that
$\lambda_{k}\times\lambda^{\prime}_{k}=\chi_{\rm{FR}}$. On each subspace given
by a fixed value of total spin and isospin, this is equivalent to finding a
common eigenvector of the corresponding Wigner D-matrices. The allowed
(physical) states will be those which satisfy all the F-R constraints.

In the case of a Skyrmion crystal, the unit cell presents a concrete set of
symmetries, some of which relate rotations both in space and isospace, and
hence we should consider as physical states only those that are compatible
with such symmetries. Let us now proceed to calculate the Finkelstein-
Rubinstein (FR) constraints in order to obtain the corresponding quantum
states of the Skyrme crystal. The relevant symmetries of the FCC half-skyrmion
crystal linking rotations and isorotations are $A_{2}$ and $C_{3}$, which are
represented by the following operators,

$\displaystyle\exp{i\frac{\pi}{2}\frac{1}{\sqrt{3}}\left(K_{1}+K_{2}+K_{3}\right)}=\mathcal{R}^{\prime}\left(0,-\pi/2,-\pi/2\right),$
(21)
$\displaystyle\exp{i\frac{\pi}{2}K_{1}}=\mathcal{R}^{\prime}\left(\pi/2,-\pi/2,\pi/2\right).$
(22)

We only write the operators that correspond to isospin transformations, since
the rotations in the real space are the same. Recall that we are using the
Euler angles representation for the rotation and isorotation operators in the
$ZYZ$ convention,

$\mathcal{R}\left(\alpha,\beta,\gamma\right)=R_{z}\left(\alpha\right)R_{y}\left(\beta\right)R_{z}\left(\gamma\right).$
(23)

From (20) we know how these operators act on a state
$\ket{j,l_{3}}\otimes\ket{i,k_{3}}$ (again we will only consider the isospin
part, the quantum numbers resulting from the spin quantization are the same),

$\mathcal{R}\left(\alpha,\beta,\gamma\right)\ket{i,k_{3}}=\sum_{k^{\prime}_{3}}D^{i}_{k_{3},k^{\prime}_{3}}\left(\alpha,\beta,\gamma\right)\ket{i,k^{\prime}_{3}},$
(24)

where $D^{I}_{k_{3},k^{\prime}_{3}}$ are the Wigner D-matrices. Then we can
consider this as a problem of finding the eigenvalues and eigenvectors of the
Wigner D-matrices, and the quantum states will be the combination of
$J,L_{3},I,K_{3}$ that satisfy the FR constraints (19). In our case, we will
consider the possible quantum states of a unit cell, which carries a baryon
number $B_{\rm cell}=4$. We will not prove here that the transformations (19)
with the corresponding spatial rotations correspond to contractible loops on
the configuration space. Instead, we will assume that this is indeed the case,
i.e. $\chi_{FR}=+1$, as tends to be the case for $B=0\mod 4$ Skyrmions Krusch
(2006). Also, we will consider all the possible values of $i$, namely
$i=0,1,2$ (eigenvalue of the isospin moment of inertia) and show the
(unnormalized) eigenvectors for each symmetry.

####  $A_{2}$ symmetry: $\mathcal{R}(0,-\pi/2,-\pi/2)$

  * •

For $i=0$, $D^{0}_{00}=1$ and there is only one state $\ket{i=0,i_{3}=0}$.

  * •

For $i=1$, we have

$D^{1}_{k_{3},k^{\prime}_{3}}=\begin{pmatrix}i/2&-1/\sqrt{2}&-i/2\\\
i/\sqrt{2}&0&i/\sqrt{2}\\\ i/2&1/\sqrt{2}&-i/2\end{pmatrix}.$ (25)

This matrix has three different eigenvalues and corresponding eigenstates

$\displaystyle\lambda_{1}=1,$
$\displaystyle\ket{\psi^{1}_{1}}=-i\ket{1,-1}+(1+i)/\sqrt{2}\ket{1,0}+\ket{1,1}$
$\displaystyle\lambda_{2}=-1/2-\sqrt{3}/2i,$
$\displaystyle\ket{\psi^{1}_{2}}=(2-\sqrt{3})i\ket{1,-1}-\sqrt{2-\sqrt{3}}(1+i)\ket{1,0}+\ket{1,1}$
$\displaystyle\lambda_{3}=-1/2+\sqrt{3}/2i,$
$\displaystyle\ket{\psi^{1}_{3}}=(2+\sqrt{3})i\ket{1,-1}+\sqrt{2+\sqrt{3}}(1+i)\ket{1,0}+\ket{1,1}$
(26)

  * •

For $i=2$, we have

$D^{2}_{k_{3},k^{\prime}_{3}}=\begin{pmatrix}-1/4&-1/2i&\sqrt{6}/4&1/2i&-1/4\\\
-1/2&-1/2i&0&-1/2i&1/2\\\ -\sqrt{6}/4&0&-1/2&0&-\sqrt{6}/4\\\
-1/2&1/2i&0&1/2i&1/2\\\ -1/4&1/2i&\sqrt{6}/4&-1/2i&-1/4\end{pmatrix}$ (27)

which shares the same eigenvalues as the corresponding $i=1$ case, but this
time $\lambda_{2}$ and $\lambda_{3}$ present multiplicity two. The
corresponding eigenstates are

$\displaystyle\lambda_{1}=1,$
$\displaystyle\ket{\psi^{2}_{1}}=-\ket{2,-2}+(1-i)\ket{2,-1}+(1+i)\ket{2,1}+\ket{2,2},$
(28) $\displaystyle\lambda_{2}=-1/2-\sqrt{3}/2i,$
$\displaystyle\ket{\psi^{2}_{2a}}=\ket{2,-2}-\sqrt{2}i\ket{2,0}+\ket{2,2},$
(29)
$\displaystyle\ket{\psi^{2}_{2b}}=-\ket{2,-2}-\frac{(1+\sqrt{3})}{2}(1-i)\ket{2,-1}+$
$\displaystyle\hskip
36.98866pt+\frac{(\sqrt{3}-1)}{2}(1+i)\ket{2,1}+\ket{2,2},$
$\displaystyle\lambda_{3}=-1/2+\sqrt{3}/2i,$
$\displaystyle\ket{\psi^{2}_{3a}}=\ket{2,-2}+\sqrt{2}i\ket{2,0}+\ket{2,2},$
(30)
$\displaystyle\ket{\psi^{2}_{3b}}=-\ket{2,-2}+\frac{(\sqrt{3}-1)}{2}(1-i)\ket{2,1}-$
$\displaystyle\hskip
36.98866pt-\frac{(1+\sqrt{3})}{2}(1+i)\ket{2,1}+\ket{2,2}.$

####  $C_{3}$ symmetry: $\mathcal{R}^{\prime}\left(\pi/2,-\pi/2,\pi/2\right)$

  * •

Again, for $i=0$, the only state is $\ket{0,0}$.

  * •

For $i=1$, the corresponding Wigner matrix

$D^{1}_{k_{3},k^{\prime}_{3}}=\begin{pmatrix}-1/2&i/\sqrt{2}&1/2\\\
-i/\sqrt{2}&0&-i/\sqrt{2}\\\ 1/2&i/\sqrt{2}&-1/2\end{pmatrix}$ (31)

has two eigenvalues, with multiplicity $2$ and $1$, respectively. The
associated eigenstates are

$\displaystyle\lambda_{1}=-1,$
$\displaystyle\ket{\phi^{1}_{1a}}=-\ket{1,-1}+\ket{1,1},$ (32)
$\displaystyle\ket{\phi^{1}_{1b}}=\ket{1,-1}+\sqrt{2}i\ket{1,0}+\ket{1,1}$
$\displaystyle\lambda_{2}=1,$
$\displaystyle\ket{\phi^{1}_{2}}=\ket{1,-1}-\sqrt{2}i\ket{1,0}+\ket{1,1}$ (33)

  * •

$i=2$,

$D^{2}_{k_{3},k^{\prime}_{3}}=\begin{pmatrix}1/4&-i/2&-\sqrt{6}/4&i/2&1/4\\\
1/2i&1/2&0&1/2&-i/2\\\ -\sqrt{6}/4&0&-1/2&0&-\sqrt{6}/4\\\
-i/2&1/2&0&1/2&i/2\\\ 1/4&i/2&-\sqrt{6}/4&-i/2&1/4\end{pmatrix}$ (34)
$\displaystyle\lambda_{1}=-1,$ (35)
$\displaystyle\ket{\phi^{2}_{1a}}=-\ket{2,-2}+i\ket{2,-1}-i\ket{2,1}+\ket{2,2},$
$\displaystyle\ket{\phi^{2}_{1b}}=\ket{2,-2}+\sqrt{6}\ket{2,0}+\ket{2,2}$
$\displaystyle\lambda_{2}=1,$ (36)
$\displaystyle\ket{\phi^{2}_{2a}}=\ket{2,-1}+\ket{2,1},$
$\displaystyle\ket{\phi^{2}_{2b}}=-\ket{2,-2}-2i\ket{2,-1}+\ket{2,2},$
$\displaystyle\ket{\phi^{2}_{2c}}=\ket{2,-2}-\sqrt{2/3}\ket{2,0}+\ket{2,2}.$

Physical states with fixed $i$ will correspond to mutual eigenstates of both
spin-$i$ Wigner matrices corresponding to the two symmetries, $A_{2}$ and
$C_{3}$. For each value of total isospin, we have been able to find a unique
state that satisfies this property,

$\begin{split}i=0&\rightarrow\ket{0,0},\\\
i=1&\rightarrow\ket{\psi^{1}_{1}}=(1+i)/2\ket{\phi^{1}_{1a}}+(1-i)/2\ket{\phi^{1}_{1b}},\\\
i=2&\rightarrow\ket{\psi^{2}_{1}}=\ket{\phi^{2}_{2b}}+(1+i)\ket{\phi^{2}_{2a}}\end{split}$
(37)

and the corresponding normalized states will be denoted by $\ket{\psi^{i}}$.
Note that each $\ket{\psi^{i}}$ corresponds to an isospin multiplet with
degeneracy $i$, since the $i_{3}$ quantum number is not restricted by the F-R
constraints (and the same happens to the corresponding spin states, in which
$j_{3}$ is not constrained either).

Then, we just find the quantum ground state of the crystal as the state with
lower total energy. Naively, one would think that it corresponds to the
$\ket{\psi^{0}}=\ket{0,0,0}$ state, as the isospin contribution vanishes.
However, this would result in a non-vanishing total electric charge of each
unit cell, so the crystal will be unstable due to an infinite contribution of
the Coulomb energy, as was already noticed in Klebanov (1985). Thus,
electrical neutrality implies that the true ground state of pure skyrmion
matter corresponds to the $i_{3}=-2$ state in the multiplet $\ket{\psi^{2}}$.

Following the previous reasoning we have obtained the allowed quantum states
for a single unit cell. On the other hand, we have assumed that a basis for
the Hilbert space of the total crystal isospin state can be constructed from
the direct product of states of individual cells. We are now in a position to
be more specific about this statement. Indeed, consider the quantum state of
two unit cells with $i=\\{i_{1},i_{2}\\}$ and $i_{3}=\\{m_{1},m_{2}\\}$. The
allowed values of $\\{i_{1},i_{2}\\}$ are $\\{0,1,2\\}$ , and $-i_{a}\leq
m_{a}\leq i_{a}$, $a=1,2$. The total quantum state will be an eigenstate of
the total isospin, so it will be better described in the coupled angular
momenta basis. Indeed, from representation theory, the tensor product of two
spin $j_{a}$ representations may be decomposed as a direct sum as

$D^{j_{1}}\otimes
D^{j_{1}}=\bigoplus\limits_{k=\absolutevalue{j_{2}-j_{1}}}^{j_{2}+j_{1}}D^{k}.$
(38)

Thus, the basis of states for the two-cell system will be
$\ket{I,I_{3},i_{1},i_{2}}$, with $\absolutevalue{j_{1}-j_{2}}\leq I\leq
j_{1}+j_{2}$ the total isospin number and $-I\leq I_{3}\leq I$ the third
component of total isospin. Therefore, to find a general basis state for an
arbitrary number of unit cells, we should just need to generalize the previous
construction to the coupling of an arbitrary number of different angular
momenta. There are arbitrarily many ways to do it, which should be equivalent
up to a unitary transformation (at least, for an arbitrary, but finite, number
of unit cells). An important final remark is that, following such
construction, the total isospin and third component of isospin of the full
crystal will remain good quantum numbers, independently of the quantum state,
so they will correspond to well defined observables in the quantum theory, as
opposed to the isospin of each of the individual cells.

###  III.3 Isospin correction to the energy per baryon

Let us rewrite the Skyrme Lagrangian (4) as

$\displaystyle\mathcal{L}=\frac{1}{24\pi^{2}}$
$\displaystyle\big{[}a\Tr\\{L_{\mu}L^{\mu}\\}+b\Tr\\{\left[L_{\mu},L_{\nu}\right]^{2}\\}$
$\displaystyle+c\,\mathcal{B}_{\mu}\mathcal{B}^{\mu}+d\Tr(U-I)\big{]}.$ (39)

The values of $a,b,c,d$ are easily obtained from (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}}.$ (40)

We now 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).$ (41)

The Maurer-Cartan (M-C) form transforms as ($\dot{g}=dg/dt$)

$\tilde{U}^{\dagger}\partial_{\mu}\tilde{U}=\left\\{\begin{array}[]{ll}gU^{\dagger}\partial_{i}Ug^{\dagger},&(\mu=i=1,2,3)\\\
g(U^{\dagger}[g^{\dagger}\dot{g},U])g^{\dagger},&(\mu=0).\end{array}\right.$
(42)

We now define the isospin angular velocity $\vec{\omega}$ as
$g^{\dagger}\dot{g}=\tfrac{i}{2}\omega_{a}\tau_{a}$. Then, we may write the
time component of the Maurer-Cartan current as
$\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},$ (43)

where we have made use of the parametrization (5). Moreover, the spatial
components of the M-C form can be written in terms of the sigma and pion
fields as well,

$\begin{split}L_{k}=&(\sigma-i\pi_{a}\tau_{a})(\partial_{k}\sigma+i\partial_{k}\pi_{b}\tau_{b})=\\\
=&i(\sigma\partial_{k}\pi_{c}-\pi_{c}\partial_{k}\sigma+\pi_{a}\partial_{k}\pi_{b}\varepsilon_{abc})\tau_{c}\equiv
L_{k}^{c}\tau_{c}.\end{split}$ (44)

The time dependence of the new Skyrme field induces a kinetic term in the
energy functional, given by 111Remember 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}]\\}$
$\displaystyle-c\,\,\mathcal{B}^{i}\mathcal{B}_{i}\\}d^{3}x,$ (45)

with $\mathcal{B}^{i}$ the spatial components of the topological current:

$\mathcal{B}^{i}=\frac{1}{24\pi^{2}}\varepsilon^{i\alpha\beta\gamma}\Tr\\{L_{\alpha}L_{\beta}L_{\gamma}\\}=\frac{3}{24\pi^{2}}\varepsilon^{ijk}\Tr\\{L_{0}L_{j}L_{k}\\}.$
(46)

We may rewrite the kinetic isorotational energy in the standard way as a
quadratic form acting on the components of the isospin angular velocity,

$T=\frac{1}{2}\omega_{i}\Lambda_{ij}\omega_{j}$ (47)

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

$\begin{split}\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\\}\,.\end{split}$
(48)

As shown in appendix B, the complete isospin inertia tensor for the unit cell
of a cubic crystal will be proportional to the identity, and its eigenvalue
(the isospin moment of inertia) will be given by

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

The numerical results for $\Lambda$ for the $\mathcal{L}_{240}$, and the full
$\mathcal{L}_{2460}$ cases are plotted as a function of the lattice parameter
in fig. 1, for different values of the sextic term coupling constant
$\lambda^{2}$. The $\Lambda$ curve for $\mathcal{L}_{24}$ was obtained by
Baskerville in Baskerville (1996) so it has been helpful to check our results.
The value of $\Lambda$ becomes smaller at high densities both without sextic
or with a small value of $\lambda^{2}$, hence the isospin correction to the
energy will grow with $n_{B}$. However the increase of the sextic coupling
constant produces a final increase of $\Lambda$ at high densities (at which
the sextic term becomes more relevant), so a non-trivial behavior will be
found in the symmetry energy. Therefore, the result for $\mathcal{L}_{2460}$
diverges from that for $\mathcal{L}_{240}$ at high densities, whereas both
join in the opposite, low density limit.

Figure 1: The eigenvalue $\Lambda$ of the isospin inertia tensor is displayed
against the lattice length parameter in Skyrme units.

Note that the pion mass potential term does not directly contribute to the
value of $\Lambda$, but it does so indirectly because it modifies the
classical solution.

The kinetic term in the Lagrangian of an isospinning cubic crystal with a
number $N_{\text{cells}}\equiv N$ of unit cells can thus be written in terms
of the isospin moment of inertia $\Lambda$ (49) as

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

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}.$ (51)

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)$ (52)

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) 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}.$ (53)

For example, for the unit cell of the FCC half-Skyrmion lattice, the isospin
contribution to the energy per baryon is

$\frac{E_{\rm{iso}}}{B_{\text{cell}}}\Bigg{\rvert}_{\rm{FCC}}=\frac{\hbar^{2}}{2\Lambda}.$
(54)

The value of $\hbar$ in Skyrme units is important to calculate the
contribution of the isospin energy, however for our choice of the parameters
$\hbar=e^{2}/(3\pi^{2})\approx 1$.

The classical Skyrmion crystal configurations can be understood as models for
isospin-symmetric nuclear matter, i.e. nuclear matter with zero total isospin.
Indeed, since the Skyrme Lagrangian is symmetric under isospin rotations of
the chiral fields, in principle there is no distinction between nucleons in
classical configurations. However, as we have seen, the quantum isospin
correction to the crystal energy per baryon does depend on the difference
between protons and neutrons through the total isospin number per unit cell.
Hence, by considering the effect of iso-rotations over classical solutions we
are effectively breaking the isospin symmetry of the static energy functional
by adding a correction of quantum origin that explicitly breaks it. Moreover,
such correction could also have been obtained through the inclusion of an
isospin chemical potential. Indeed, we may introduce a nonzero isospin
chemical potential $\mu_{I}$ in any chiral effective theory (and the Skyrme
model in particular) in terms of a covariant derivative of the chiral fields
of the form Son and Stephanov (2001)

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

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}.$ (56)

Comparing with eq. 42, one sees that 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}})$
(57)

and then

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

As we have seen in the previous subsection, the ground state of the Skyrme
crystal is forced by the charge neutrality condition to have
$(J^{\rm{tot}})^{2}=J_{3}^{2}$, and for a unit cell of $B_{\text{cell}}$
baryon number, the chemical potential is simply
$\mu_{I}(L)=-\hbar^{2}B_{\text{cell}}/(2\Lambda)$. This indeed coincides with
the expression of $\omega_{a}$ in terms of the isospin angular momentum and
the isospin moment of inertia. As a final comment, we remark that the
inclusion of isospin may have non-trivial effects on the geometry of classical
solitonic solutions. In particular, it was shown in Loewe  _et al._ (2006)
that the isospin chemical potential may alter the stability of classical
Skyrmion configurations. Hence, a thorough analysis of dense Skyrmion matter
at finite isospin chemical potential should focus on its effects on the
classical crystal configuration. However, as the contribution to the isospin
energy per unit cell $E_{\rm{iso}}(L)$ is never dominating for any $L$, we may
neglect the backreaction of the isospin term into the background and consider
it simply as a (quantum) correction to the energy.

##  IV Quantum skyrmion crystals and the symmetry energy of $npe(\mu)$ matter

We remark that the charge neutrality constraint directly fixes the total
isospin number density per unit cell of the lattice to be $i=-2$, which
describes a pure neutron crystal. This corresponds to a highly isospin
asymmetric state of matter, similar to what is expected to occur in the core
of neutron stars. However, realistic models of the nuclear matter inside
neutron stars do not predict a totally asymmetric state, but a small fraction
of protons over total nucleons is generally allowed.

In the case of nuclear matter as described by a Skyrmion crystal, the
fundamental degrees of freedom are not protons and neutrons, and in fact the
number of protons and/or neutrons per unit cell is ill-defined in some
crystalline phases in which the baryon charge is fractionalized due to the
presence of half-Skyrmions. However, we can still define the asymmetry
parameter and thus the symmetry energy of a Skyrme crystal Lee  _et al._
(2021), which allows us to calculate the density dependence of $S$, and obtain
the above defined observables (3), within the Skyrme model.

###  IV.1 The electric charge density of pure Skyrmion matter

The Gell-Mann-Nishijima formula tells us how to obtain the charge density of a
Skyrmion field configuration,

$Z=\frac{1}{2}B+I_{3}\hskip 5.69054pt\longrightarrow\hskip
5.69054pt\rho=\frac{1}{2}B^{0}+\expectationvalue{I^{0}_{3}},$ (59)

where $I^{0}_{3}$ is the time-like component of the third isospin Noether
current, and the brackets represent the expectation value on the crystal
quantum state. An explicit expression for (the classical version of) this
current can be obtained from Noether’s theorem, given that the infinitesimal
version of (41) yields the transformation

$U\rightarrow U^{\prime}=U+\epsilon^{k}\delta U_{k},\qquad\delta
U_{k}=\frac{i}{2}[\tau^{k},U],$ (60)

and from the definition of the corresponding Noether current,

$\begin{split}I^{\mu}_{k}=\frac{i}{24\pi^{2}}&\left(2a\Tr\left\\{L^{\mu}T_{k}\right\\}-4b\Tr\left\\{\left[L^{\nu},L_{\mu}\right]\left[L^{\nu},T_{k}\right]\right\\}-\right.\\\
&\left.-c\epsilon^{\mu\nu\alpha\beta}B_{\nu}\Tr\left\\{L_{\alpha}L_{\beta}T_{k}\right\\}\right),\end{split}$
(61)

where $T_{k}$ is the $\mathfrak{su}(2)$ current defined in (43). In
particular, for the third isospin charge we have

$I_{k}^{0}=-R_{kj}\lambda_{ij}\omega_{i},$ (62)

where $R_{ab}$ is the $SO(3)$ matrix representation of the $SU(2)$ element
$g$, $R_{ab}(g)=\Tr{\tau_{a}g^{\dagger}\tau_{b}g}$; $\lambda_{ij}$ is the
inertia tensor density, defined as the integrand of (48), and we have used the
transformation property of the Pauli matrices
$g^{\dagger}\tau_{a}g=R_{ab}\tau_{b}$.

The quantum version of $I_{3}^{0}$ can be obtained by substituting the
classical variables $g$ and $\omega_{i}$ with the corresponding quantum
operators and Weyl ordering the products of two or more non-commuting
operators. In the case of the angular velocity, we will use the quantum (body-
fixed) isospin angular momentum operator, $K_{i}=\omega_{i}/(\Lambda N)$.
Then, the expectation value of $I_{3}^{0}$ in the isospin state $\ket{\psi}$
is given by

$\expectationvalue{I_{3}^{0}}=-\frac{\lambda_{ij}}{2\Lambda
N}\bra{\psi}(K_{i}R_{3j}+R_{3j}K_{i})\ket{\psi},$ (63)

where $\lambda_{ij}$ is the isospin inertia tensor density, i.e the integrand
in (48).

Furthermore, we may write
$R_{3j}K_{i}=K_{i}R_{3j}-[K_{i},R_{3j}]=K_{i}R_{3j}-i\epsilon_{ijk}R_{3k}$,
and since the isospin inertia tensor density is symmetric we just have

$\expectationvalue{I_{3}^{0}}=-\frac{\lambda_{ij}}{\Lambda
N}\bra{\psi}K_{i}R_{3j}\ket{\psi}=-\frac{\lambda_{ij}}{\Lambda
N}\bra{\psi}R_{3j}K_{i}\ket{\psi}.$ (64)

As shown in the previous section, there is only one quantum state that
preserves the symmetry of the unit cell for each value of the isospin quantum
number, that we have denoted by $\ket{\psi^{i}}$, $(i=0,1,2)$ in (37).
Therefore, for a single unit cell, the electric charge density can have three
different profiles, depending on the quantum state of the cell. We are,
however, interested in the case of (infinitely) many crystal cells. To obtain
the charge density of the crystal, we would need to know the full quantum
state. For pure skyrmion matter, we have seen that charge neutrality implies
that the lowest energy state must be the tensor product
$\ket{\Psi}=\bigotimes\limits_{\rm cells}\ket{\psi^{2},i_{3}=-2}$. Then the
charge density becomes simpler, since it is just the same for each unit cell
and given by (59), where $\expectationvalue{I_{0}^{3}}$ is just the quantum
isospin density for each cell, as given explicitly by (64). Performing the
calculation of the matrix element as in appendix A, we may obtain the charge
density of the skyrmion crystal in its ground state, which is represented in
fig. 2. The figure reveals a nontrivial charge fractionalization in the half-
skyrmion phase, in which half-skyrmions are positively charged while the space
between them presents negative charge to keep the charge neutrality of the
unit cell.

Figure 2: Electric charge density of half-Skyrmion matter in the quantum
ground state.

The fact that the charge density presents a non trivial multipolar structure
means that a residual Coulomb interaction between unit cells may exist.
However, we have checked numerically that such a contribution is several
orders of magnitude lower than the isospin correction, therefore unimportant
for the total binding energy of the crystal.

###  IV.2 The symmetry energy

Let us consider a finite Skyrme crystal of $N$ unit cells, and let $B=N\times
B_{\text{cell}}$, where $B_{\text{cell}}$ is the baryon number of a unit cell.
We do not enforce charge neutrality at this step, and further leave unknown
the quantum state $\ket{\Psi}$ of the crystal. We have seen that the total
charge of this system is

$Q=\expectationvalue{e\int
d^{3}x\\{B^{0}/2+I^{0}_{3}\\}}=eN\quantity[\frac{B_{\text{cell}}}{2}+\frac{\expectationvalue{\int
I^{0}_{3}d^{3}x}}{N}]$ (65)

As argued at the end of section III, the total third component of isospin is a
good quantum number for the total quantum state of the crystal, although this
is not true for individual unit cells. In other words, the expectation value

$\expectationvalue{I_{3}}=\bra{\Psi}\int I^{0}_{3}d^{3}x\ket{\Psi}$ (66)

is well defined in an arbitrary quantum state, but
$\int\expectationvalue{I^{0}_{3}}d^{3}x$ is not. Since we are seeking for a
definition of the isospin density in the quantum theory, we may 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{\expectationvalue{I_{3}}}{NV_{\rm
cell}}\doteq\frac{n_{I}}{V_{\rm cell}}$ (67)

where $n_{I}$ is the effective isospin charge per unit cell in this arbitrary
quantum state. We may further consider the effective proton fraction that
would yield such an isospin charge per unit cell with baryon number
$B_{\text{cell}}$ to write

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

Hence, we may write the isospin energy per unit cell of the Skyrmion crystal
in such a state in terms of the asymmetry parameter

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

and thus the symmetry energy for Skyrme crystals is given by

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

Figure 3: Symmetry energy of Skyrme crystals as a function of the density for
different values of $\lambda$. Constraints of the symmetry energy at sub-
saturation densities from isobaric analog states Danielewicz and Lee (2014)
are plotted in red. The grey region corresponds to recent constraints from the
analysis of neutron star observations Li  _et al._ (2021)

In fig. 3 we plot the symmetry energy of Skyrmion crystals for the generalized
Skyrme model with different values of the sextic coupling constant.

We can observe that, for this choice of parameters, the symmetry energy of the
crystal in the $\mathcal{L}_{240}$ submodel ($\lambda^{2}=0$) comes out too
big with respect to the constraints at saturation density $n_{0}$, while the
inclusion the sextic term reduces the value. We found that a value of
$\lambda^{2}\sim 1.5\,{\rm MeV\,fm^{3}}$ fits to the correct value at
saturation. Moreover, we may also compute the slope and curvature parameters
at saturation for all cases, shown in table 1. Interestingly, for values of
$\lambda^{2}$ between $1.5-3$ MeV fm-3 we find that both the saturation
density and the symmetry energy parameters agree quite well with the most up-
to-date experimental values. Here we have defined the saturation density as
the density where the energy per baryon of the skyrmion crystal takes its
minimum value.

$\lambda^{2}$ (MeV fm3) |  $n_{0}$ (fm-3) |  $S_{0}$ (MeV) |  $L$ (MeV) |  $K_{\text{sym}}$ (MeV)
---|---|---|---|---
0 | 0.33 | 44.4 | 72.9 | -143
1.5 | 0.22 | 31.9 | 46.4 | -130
3 | 0.18 | 26.4 | 35.4 | -118
5 | 0.15 | 22.2 | 27.5 | -105
Exp. | 0.16 $\pm$ 0.01 | 31.7 $\pm$ 3.2 | 57.7 $\pm$ 19 | -107$\pm$ 88
Table 1: Symmetry energy coefficients in the Skyrme model with different
sextic couplings. In the last row, we show the most up-to-date fiducial values
of $n_{0},S_{0},L$ and $K_{\text{sym}}$ Li  _et al._ (2021)

Hence, the constraints on the symmetry energy yield rather stringent bounds on
the value of $\lambda^{2}$, where the precise numerical values of these bounds
will, of course, depend on the choices made for the other model parameters
$f_{\pi}$ and $e$. We remark that a lower bound for this constant can also be
obtained from the maximum mass requirement of neutron star EoS Adam  _et al._
(2021).

We end our discussion on the Symmetry energy of Skyrme crystals by pointing
out the absence of the _cusp structure_ predicted by Lee et al. Lee  _et al._
(2011) in our results of $S$. Their argument (recently reviewed in Lee  _et
al._ (2021)) for the appearance of a minimum in the symmetry energy at a given
density $n_{1/2}>n_{0}$ is based on the Skyrmion to half-Skyrmion transition,
which is proposed to happen at $n_{1/2}\sim 2-3n_{0}$. In our pure Skyrme
model setting, however, such a transition has been shown to occur in a
thermodynamically unstable branch of the Skyrmion crystal phase diagram once
the pion mass is taken into account (see, eg. Park and Vento (2010); Adam  _et
al._ (2022)). Therefore, we do not find such a transition, as we consider a
crystal of (nearly) half-Skyrmions to be the correct ground state for
densities $n>n_{0}$. The ground state of the Skyrme model for densities $n\leq
n_{0}$ is still not well understood, and some inhomogeneous configurations
have been proposed Park  _et al._ (2019); Adam  _et al._ (2022) that point
towards a complicated phase structure predicted by the Skyrme model near
saturation. In particular, the transition from regular nuclear matter to a
crystal of half-Skyrmions should take place in such a range of densities.

###  IV.3 Particle fractions of $npe\mu$ matter in $\beta$-equilibrium

As previously argued, any quantum state that deviates from the ground state in
pure skyrme matter would lead to a divergence in the Coulomb energy in the
infinite crystal limit. Indeed, an isolated system of positively charged
matter is unstable due to Coulomb repulsion. Therefore, it is assumed that
there exists a neutralizing background of negatively charged leptons
(electrons and possibly muons), such that this Coulomb repulsion is
compensated. Such a system of nuclear matter plus leptons is characterized in
the equilibrium phase by two equilibrium conditions, namely the _charge
neutrality condition_

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

i.e., the densities of positively charged nucleons (protons) and negatively
charged leptons (electrons and muons) are equal, and the _$\beta$ -equilibrium
condition_

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

i.e., the isospin chemical potential must equal that of charged leptons, such
that the neutron decay and electron capture processes

$n\rightarrow p+l+\bar{\nu}_{l}\quad,\quad p+l\rightarrow n+\nu_{l}$ (73)

take place at the same rate. Moreover, leptons inside a neutron star are
usually 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}}$ (74)

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. Indeed, for sufficiently
large densities the electron chemical potential will be larger than the muon
mass, $\mu_{e}\geq m_{\mu}$, and the appearance of muons in the system will be
energetically favorable. We may now estimate the total proton fraction by
enforcing both charge neutrality and beta equilibrium. Let us start by
neglecting the contribution of muons to the charge density. Then, from the
charge neutrality condition (71), we relate the electron density to the proton
fraction parameter, $n_{e}=\gamma B_{\text{cell}}/(2L)^{3}$, and the $\beta$
equilibrium condition yields an equation that implicitly defines $\gamma$ 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}$
(75)

where we have also made 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}},$ (76)

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}}.$
(77)

On the other hand, the proton fraction inside the beta-equilibrated matter
also determines whether a proto-neutron star will go through a cooling epoch
via neutrino emission through the direct Urca (DU) process $n\rightarrow
p+e+\bar{\nu}_{e}$, which is expected to occur if the proton fraction reaches
a critical value, $\gamma_{p}>x_{DU}$, the so-called DU-threshold Lattimer
_et al._ (1991); Klahn  _et al._ (2006). As 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 merge of binary NS systems Alford  _et al._
(2021) would determine the proton fraction (hence, the symmetry energy) of
matter at ultra-high densities. However, it is not clear whether such enhanced
cooling actually takes place, although there is recent evidence that supports
it Brown  _et al._ (2018).

In $npe\mu$ matter, the DU threshold is given by Klahn  _et al._ (2006)

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

Figure 4: 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.

The particle populations $\gamma_{i}$ in the beta-equilibrated Skyrmion matter
are shown in fig. 4 for the cases $\lambda^{2}=0,1.5$ MeV fm3.

In both cases, a persistent population of protons and leptons with increasing
nucleon density is expected, although in the case with sextic term we see that
the fraction of charged particles is smaller. This is the impact of the sextic
term, since it is much easier to convert protons into neutrons due to the
lower symmetry energy. Finally, the DU-threshold is not reached in any case
for the values of the parameters ($f_{\pi},e,\lambda^{2}$) that we have
chosen. However, one should not take this fact as a prediction of the Skyrme
model, as it strongly depends on the parameter values. Also, it is generally
assumed that around $2-3$ times the nuclear saturation density, additional
degrees of freedom (strange baryons) become important for the description of
nuclear matter, which in particular may affect the proton fraction at such
densities.

##  V Conclusions

In this paper, we have applied standard quantization methods to find the
quantum ground state of classical skyrmion crystal configurations. The
quantization of the isospin degrees of freedom has allowed us to determine how
the isospin contribution to the energy (hence the isospin chemical potential)
depends on the crystal density. Then, we have shown that the effects of a
finite isospin asymmetry can be included once the divergence in Coulomb energy
of the infinite crystal is avoided by the inclusion of a neutralising leptonic
background, yielding a particularly simple model for studying the high density
regime of isospin asymmetric nuclear matter.

In conclusion, the Skyrme model yields a concrete, nonperturbative framework
in which not only the energy of symmetric matter can be calculated at high
densities, but also the corresponding symmetry energy can be computed. This
extends the list of physical, dense nuclear matter observables which can be
described by the model. Further, we can fit the parameters of the model to
physical observables like the symmetry energy at saturation, in order to be
able to make predictions about the nuclear matter EoS at higher densities
within the Skyrme model. Concretely, we found that a simple scan of values for
the sextic term coupling constant allows us to fit the symmetry energy well
within the current uncertainties. However, the saturation energy still comes
out a bit smaller than its experimental value. This is a universal feature of
the Skyrme model, which usually yields too small values for the energy per
baryon of nuclei and nuclear matter, i.e., a too strongly bound nuclear
matter. Our choice for these parameter values is mainly motivated by the
possibility to compare with other Skyrme model calculations, being our main
objective to provide a clear discussion of the most important physical
properties of quantized skyrmion crystals, including the contribution from the
sextic term. The determination of the best set of parameter values to describe
nuclear matter and NS should include the calculation of NS properties and a
careful comparison with the most recent observational data like the ones
inferred from the LIGO and NICER detections, as was done, e.g., in Adam  _et
al._ (2020), Adam  _et al._ (2021), and is out of the scope of this work.

Finally, the results of this paper open the door for a further improvement in
the description of the equation of state for dense Skyrmion matter, including
the contribution of strange degrees of freedom, such as kaon condensates or
hyperonic matter which can be done systematically in the Skyrme model Callan
and Klebanov (1985); Blom  _et al._ (1989). Indeed, hyperons and strange
mesons are believed to constitute a non-negligible fraction of the total
particle content of ultra-dense matter, therefore they may be important for an
accurate determination of the proton fraction at densities $n\gtrsim 2-3n_{0}$
Park  _et al._ (2010).

###### Acknowledgements.

The authors would like to thank C. Naya and N. Manton for helpful discussions
and comments. Further, the authors acknowledge financial support from the
Ministry of Education, Culture, and Sports, Spain (Grant No.
PID2020-119632GB-I00), the Xunta de Galicia (Grant No. INCITE09.296.035PR and
Centro singular de investigación de Galicia accreditation 2019-2022), the
Spanish Consolider-Ingenio 2010 Programme CPAN (CSD2007-00042), and the
European Union ERDF. AW is supported by the Polish National Science Centre,
grant NCN 2020/39/B/ST2/01553. AGMC is grateful to the Spanish Ministry of
Science, Innovation and Universities, and the European Social Fund for the
funding of his predoctoral research activity (_Ayuda para contratos
predoctorales para la formación de doctores_ 2019). MHG is also grateful to
the Xunta de Galicia (Consellería de Cultura, Educación y Universidad) for the
funding of his predoctoral activity through _Programa de ayudas a la etapa
predoctoral_ 2021\.

##  Appendix A Calculation of matrix elements

The physical Hilbert space is spanned by the basis $\ket{X}=\ket{ii_{3}k_{3}}$
of states with fixed total isospin and third components of body-fixed and
space-fixed isospin. We may parametrize the isospin collective coordinates as
an element $A\in SU(2)$ by $A=a_{0}+ia_{k}\tau_{k}$, where $a=(a_{0},\vec{a})$
is subject to the constraint $a_{0}^{2}+a_{k}a_{k}=1$. Thus, the coordinate
vector $a$ can be thought of parametrizing a four-dimensional sphere, and
hence, in the $\ket{A}$ representation, the wave function of each of the basis
states corresponds to the (hyper-)spherical harmonics on $S^{3}$:

$\psi(A)=\innerproduct{A}{ii_{3}k_{3}}=Y_{i_{3},k_{3}}^{2i}(a),$ (79)

with

$Y_{m,m^{\prime}}^{2j}(A)=[(2j+1)/8\pi^{2}]^{1/2}{D}^{j}_{m,m^{\prime}}(A),$
(80)

where ${D}^{j}_{m,m^{\prime}}(A)$ are the Wigner’s $D$-matrices, associated to
an irreducible, spin-$j$ representation of the $SU(2)$ group. Let us write its
explicit expression in terms of the matrix elements of an arbitrary $2\times
2$ matrix $B$ Bander and Itzykson (1966):

(81)

and the constants $n_{i}$ are related through ${n_{1}+n_{2}=j-m}$,
${n_{3}+n_{4}=j+m}$, ${n_{1}+n_{3}=j+m^{\prime}}$, $n_{2}+n_{4}=j-m^{\prime}$.
Indeed, the wavefunctions (80) are correctly normalized under the
corresponding inner product,

$\int_{SU(2)}(Y_{i_{3},k_{3}}^{2j}(A))^{*}Y_{i^{\prime}_{3},k^{\prime}_{3}}^{2j^{\prime}}(A)d\mu(A)=\delta_{jj^{\prime}}\delta_{i_{3}i_{3}^{\prime}}\delta_{k_{3}k_{3}^{\prime}}$
(82)

where $d\mu(A)$ is the Haar measure of $SU(2)$, which can be written in this
case as a volume element in terms of three angular variables parametrizing the
three-sphere $S^{3}$. We may choose these angles to be the standard
hyperspherical coordinates for a unit (Euclidean) 4-vector $a$:

$\displaystyle a_{0}=\cos\chi,\quad a_{1}=\sin\chi\sin\theta\cos\phi,$
$\displaystyle a_{2}=\sin\chi\sin\theta\sin\phi,\quad
a_{3}=\sin\chi\cos\theta$ (83)

where the angles take the values $\chi,\theta\in[0,\pi]$, $\phi\in[0,2\pi]$.
The Haar measure is then

$d\mu(A[\chi,\theta,\phi])=\sin^{2}\chi\sin\theta d\chi d\theta d\phi,$ (84)

and the states $D_{m,m^{\prime}}^{j}(A[\chi,\theta,\phi])$ are obtained from
eq. 81 given that $a=ia_{0}-a_{3}$, $b=ia_{2}-a_{1}$, $c=-b^{*}$, $d=a^{*}$.
Once we have characterized all the basis states, we are ready to calculate the
matrix elements of any observable that can be expressed as an operator
$\hat{F}(A)$ over the Hilbert space $\mathcal{H}$. Indeed, for any pair of
basis sates $\ket{ii_{3}k_{3}},\ket{i^{\prime}i_{3}^{\prime}k_{3}^{\prime}}$,
we may write the matrix element as an integral over $SU(2)$:

$\bra{ii_{3}k_{3}}\hat{F}(A)\ket{i^{\prime}i_{3}^{\prime}k_{3}^{\prime}}=\int_{SU(2)}(Y_{i_{3},k_{3}}^{2i}(A))^{*}F(A)Y_{i^{\prime}_{3},k^{\prime}_{3}}^{2i^{\prime}}(A)d\mu(A).$
(85)

or equivalently as an integral in terms of the hyperspherical angles.

##  Appendix B Isospin inertia tensor of Skyrme crystals

All crystalline configurations that we have considered in section II.2 share
the same basic cubic symmetry group, generated by the $A_{1}$ and $A_{2}$
transformations. This implies that Skyrmion crystals present a cubic symmetry
in isospin space as well as in real space. Therefore, the isospin inertia
tensor of such configurations becomes proportional to the identity, i.e.
$\Lambda^{\rm{crystal}}_{ij}=\Lambda\delta_{ij}$. We will now prove this
statement and obtain an explicit expression for the eigenvalue of the isospin
inertia tensor. To do so, let us consider the three terms in the rhs of (48)
separately.

#### Quadratic term.

We need to calculate the following trace:

$\begin{split}\Tr\\{T_{i}T_{j}\\}&=-2T_{i}^{a}T_{j}^{a}=2(\pi_{i}\pi_{j}-\pi_{b}\pi_{b}\delta_{ij})\end{split}$
(86)

This expression has the form of the standard definition of the inertia tensor.
It is indeed symmetric, but not diagonal. However, we have not used the cubic
symmetry of the configuration yet. Let us see that, as in the case of a cubic
rigid body, the moment of inertia becomes proportional to the identity due to
its symmetry. Consider the integral of $\pi_{i}\pi_{j}$ over a unit cell
$\Omega$ of size $2L$:

$I_{ij}=\int_{-L}^{L}dx\int_{-L}^{L}dy\int_{-L}^{L}dz\pi_{i}\pi_{j}.$ (87)

It is straightforward to see that symmetry $A_{1}$ implies that $I_{ij}=0$ for
$i\neq j$ and that symmetry $A_{2}$ implies $I_{xx}=I_{yy}=I_{zz}$, so that we
may write

$I_{ij}=\delta_{ij}\frac{1}{3}\int_{\Omega}\pi_{a}\pi^{a}d^{3}x,$ (88)

and hence the total contribution of the quadratic term to the inertia tensor
of the crystal can be written

$2a\Lambda^{(2)}_{ij}=2a\Lambda^{(2)}\delta_{ij},\quad\text{with}\quad\Lambda^{(2)}=-\frac{4}{3}\int_{\Omega}\pi_{a}\pi^{a}d^{3}x.$
(89)

#### Quartic term.

The quartic term is proportional to:

$\Tr\\{[T_{i},L_{k}][T_{j},L_{k}]\\}=8[T^{a}_{i}T^{a}_{j}L_{k}^{b}L_{k}^{b}-T_{i}^{a}L_{k}^{a}T^{b}_{j}L_{k}^{b}]$
(90)

Noting that:

$\begin{split}T^{a}_{i}T^{a}_{j}&=\vec{\pi}^{2}\delta_{ij}-\pi_{i}\pi_{j},\\\
T_{i}^{a}L_{k}^{a}&=i[2\sigma(\pi_{i}\vec{\pi}\partial_{k}\vec{\pi}-\vec{\pi}^{2}\partial_{k}\pi_{i})+(\vec{\pi}\times\partial_{k}\vec{\pi})_{i}],\\\
L^{a}_{k}L_{k}^{a}&=-[(\sigma\partial_{k}\pi_{a}-\pi_{a}\partial_{k}\sigma)^{2}+\vec{\pi}^{2}(\partial_{k}\vec{\pi})^{2}-(\vec{\pi}\partial_{k}\vec{\pi})^{2}],\end{split}$
(91)

we have

$\begin{split}\Tr\\{[T_{i},L_{k}][T_{j},L_{k}]\\}=&-8(\vec{\pi}^{2}\delta_{ij}-\pi_{i}\pi_{j})[(\sigma\partial_{k}\pi_{a}-\pi_{a}\partial_{k}\sigma)^{2}+\vec{\pi}^{2}(\partial_{k}\vec{\pi})^{2}-(\vec{\pi}\partial_{k}\vec{\pi})^{2}]+8(\vec{\pi}\times\partial_{k}\vec{\pi})_{i}(\vec{\pi}\times\partial_{k}\vec{\pi})_{j}\\\
+&16\sigma[(\vec{\pi}\times\partial_{k}\vec{\pi})_{i}(\pi_{j}\vec{\pi}\partial_{k}\vec{\pi}-\vec{\pi}^{2}\partial_{k}\pi_{j})+(\vec{\pi}\times\partial_{k}\vec{\pi})_{j}(\pi_{i}\vec{\pi}\partial_{k}\vec{\pi}-\vec{\pi}^{2}\partial_{k}\pi_{i})]+\\\
+&32\sigma^{2}[(\vec{\pi}\partial_{k}\vec{\pi})^{2}\pi_{i}\pi_{j}+\vec{\pi}^{4}\partial_{k}\pi_{i}\partial_{k}\pi_{j}-(\vec{\pi}\partial_{k}\vec{\pi})\vec{\pi}^{2}(\pi_{i}\partial_{k}\pi_{j}+\pi_{j}\partial_{k}\pi_{i})]\end{split}$
(92)

this expression is again symmetric in $i,j$ and can be further reduced when
taking into account the $A_{1}$ and $A_{2}$ symmetries of crystal
configurations. Indeed, following the same reasoning that in the previous
case, one can show, by making use of such symmetries, that for any scalar
function $F(\vec{x})$ and any vector function $G_{k}(\vec{x})$, one has

$\begin{split}&\int_{\Omega}F\pi_{i}\pi_{j}d^{3}x=\delta_{ij}\frac{1}{3}\int_{\Omega}F\vec{\pi}^{2}d^{3}x,\\\
&\int_{\Omega}F\partial_{k}\pi_{i}\partial_{k}\pi_{j}d^{3}x=\delta_{ij}\frac{1}{3}\int_{\Omega}F(\partial_{k}\vec{\pi})^{2}d^{3}x,\\\
&\int_{\Omega}G_{k}\pi_{i}\partial_{k}\pi_{j}d^{3}x=\delta_{ij}\frac{1}{3}\int_{\Omega}G_{k}\vec{\pi}\partial_{k}\vec{\pi}d^{3}x,\\\
&\int_{\Omega}(\vec{\pi}\times\partial_{k}\vec{\pi})_{i}(\vec{\pi}\times\partial_{k}\vec{\pi})_{j}d^{3}x=\delta_{ij}\frac{1}{3}\int_{\Omega}(\vec{\pi}\times\partial_{k}\vec{\pi})^{2}d^{3}x,\\\
&\int_{\Omega}F\partial_{k}\pi_{i}(\vec{\pi}\times\partial_{k}\vec{\pi})_{j}d^{3}x=\int_{\Omega}G_{k}\pi_{i}(\vec{\pi}\times\partial_{k}\vec{\pi})_{j}d^{3}x=0.\end{split}$
(93)

Therefore, the total contribution to the inertia tensor of the quadratic term
is indeed proportional to the identity, and given by

$-4b\Lambda^{(4)}_{ij}=-4b\Lambda^{(4)}\delta_{ij},\quad\rm{with}$ (94)
$\begin{split}\Lambda^{(4)}=-&\frac{8}{3}\int_{\Omega}\\!\\!d^{3}x\left\\{2\vec{\pi}^{2}[(\sigma\partial_{k}\pi_{a}-\pi_{a}\partial_{k}\sigma)^{2}+\vec{\pi}^{2}(\partial_{k}\vec{\pi})^{2}-(\vec{\pi}\partial_{k}\vec{\pi})^{2}]\right.\\\\[5.69054pt]
-&\left.(\vec{\pi}\times\partial_{k}\vec{\pi})^{2}-4\sigma^{2}[(\vec{\pi}\partial_{k}\vec{\pi})\vec{\pi}-\vec{\pi}^{2}\partial_{k}\vec{\pi}]^{2}\right\\}.\end{split}$
(95)

#### Sextic term.

For the sextic term contribution we must calculate the following expression:

$\varepsilon^{abc}\Tr\\{T_{i}L_{b}L_{c}\\}\varepsilon_{ars}\Tr\\{T_{j}L_{r}L_{s}\\}=8\varepsilon_{lmn}\varepsilon_{pqr}T_{i}^{l}L_{a}^{m}L_{b}^{n}T_{j}^{p}L_{a}^{q}L_{b}^{r}.$
(96)

Given that

$\begin{split}\varepsilon_{lmn}T_{i}^{l}L_{a}^{m}L_{b}^{n}=&-\Big{[}A_{ab}\pi_{i}+B_{b}\partial_{a}\pi_{i}-B_{a}\partial_{b}\pi_{i}+C_{b}(\vec{\pi}\times\partial_{a}\vec{\pi})_{i}-\\\
&-C_{a}(\vec{\pi}\times\partial_{b}\vec{\pi})_{i}+D(\partial_{a}\vec{\pi}\times\partial_{b}\vec{\pi})_{i}\Big{]},\end{split}$
(97)

where

$\begin{split}A_{ab}&=\vec{\pi}^{2}(\partial_{b}\sigma\vec{\pi}\cdot\partial_{a}\vec{\pi}-\partial_{a}\sigma\vec{\pi}\cdot\partial_{b}\vec{\pi})-\sigma^{2}\vec{\pi}\cdot(\partial_{a}\vec{\pi}\times\partial_{b}\vec{\pi}),\\\
B_{a}&=\sigma\vec{\pi}^{2}\vec{\pi}\cdot\partial_{a}\vec{\pi}-\vec{\pi}^{4}\partial_{a}\sigma,\\\
C_{a}&=\sigma^{2}\vec{\pi}\cdot\partial_{a}\vec{\pi}-2\vec{\pi}^{2}\sigma\partial_{a}\sigma,\quad
D=\sigma^{2}\vec{\pi}^{2},\end{split}$ (98)

we are ready to write the contribution of the sextic term to the inertia
tensor:

$-\frac{c}{32\pi^{4}}\Lambda^{(6)}_{ij}=-\frac{c}{32\pi^{4}}\Lambda^{(6)}\delta_{ij},$
(99)

with

$\begin{split}\Lambda^{(6)}&\\!=\\!\frac{1}{3}\int_{\Omega}\\!\\!d^{3}x\big{\\{}A_{ab}A_{ab}\vec{\pi}^{2}+2A_{ab}B_{b}\vec{\pi}\partial_{a}\vec{\pi}-2A_{ab}B_{a}\vec{\pi}\partial_{b}\vec{\pi}+2B_{b}B_{b}\partial_{a}\vec{\pi}\partial_{a}\vec{\pi}-2B_{a}B_{b}\partial_{b}\vec{\pi}\partial_{a}\vec{\pi}+\\\\[5.69054pt]
&+2C_{a}C_{a}(\vec{\pi}\times\partial_{b}\vec{\pi})(\vec{\pi}\times\partial_{b}\vec{\pi})-2C_{a}C_{b}(\vec{\pi}\times\partial_{a}\vec{\pi})(\vec{\pi}\times\partial_{b}\vec{\pi})+\\\\[5.69054pt]
&+D^{2}(\partial_{a}\vec{\pi}\times\partial_{b}\vec{\pi})(\partial_{a}\vec{\pi}\times\partial_{b}\vec{\pi})+4C_{a}D(\vec{\pi}\times\partial_{b}\vec{\pi})(\partial_{a}\vec{\pi}\times\partial_{b}\vec{\pi})+A_{ab}D\vec{\pi}(\partial_{a}\vec{\pi}\times\partial_{b}\vec{\pi})\Big{\\}}.\end{split}$
(100)

To obtain this expression, we have used again (93).

## References

  * Skyrme (1962) T. H. R. Skyrme, Nuclear Physics 31, 556 (1962).
  * Witten (1979) E. Witten, Nuclear Physics B 156, 269 (1979).
  * Balachandran  _et al._ (1982) A. P. Balachandran, V. P. Nair, S. G. Rajeev, and A. Stern, Phys. Rev. Lett. 49, 1124 (1982), [Erratum: Phys.Rev.Lett. 50, 1630 (1983)].
  * Adkins  _et al._ (1983) G. S. Adkins, C. R. Nappi, and E. Witten, Nuclear Physics B 228, 552 (1983).
  * Battye  _et al._ (2009) R. A. Battye, N. S. Manton, P. M. Sutcliffe, and S. W. Wood, Physical Review C 80, 034323 (2009).
  * Lau and Manton (2014) P. H. C. Lau and N. S. Manton, Phys. Rev. Lett. 113, 232503 (2014).
  * Halcrow  _et al._ (2017) C. J. Halcrow, C. King, and N. S. Manton, Phys. Rev. C 95, 031303 (2017).
  * Halcrow and Harland (2020) C. Halcrow and D. Harland, Phys. Rev. Lett. 125, 042501 (2020).
  * Adam  _et al._ (2010) C. Adam, J. Sánchez-Guillén, and A. Wereszczyński, Physics Letters B 691, 105–110 (2010).
  * Gillard  _et al._ (2015) M. Gillard, D. Harland, and M. Speight, Nucl. Phys. B 895, 272 (2015).
  * Gudnason (2016) S. B. Gudnason, Phys. Rev. D 93, 065048 (2016).
  * Naya and Sutcliffe (2018) C. Naya and P. Sutcliffe, Phys. Rev. Lett. 121, 232002 (2018).
  * Miller  _et al._ (2021) M. C. Miller _et al._ , Astrophys. J. Lett. 918, L28 (2021).
  * Riley  _et al._ (2021) T. E. Riley _et al._ , Astrophys. J. Lett. 918, L27 (2021).
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (LIGO , VIRGO), Phys. Rev. Lett. 119, 161101 (2017).
  * Abbott  _et al._ (2020a) B. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 892, L3 (2020a).
  * Abbott  _et al._ (2020b) R. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 896, L44 (2020b).
  * Adam  _et al._ (2015a) C. Adam, C. Naya, J. Sanchez-Guillen, R. Vazquez, and A. Wereszczynski, Physics Letters B 742, 136–142 (2015a).
  * Adam  _et al._ (2020) C. Adam, A. G. Martín-Caro, M. Huidobro, R. Vázquez, and A. Wereszczynski, Phys. Lett. B 811, 135928 (2020).
  * Adam  _et al._ (2022) C. Adam, A. G. Martin-Caro, M. Huidobro, R. Vazquez, and A. Wereszczynski, Phys. Rev. D 105, 074019 (2022), arXiv:2109.13946 [hep-th] .
  * Adam  _et al._ (2021) C. Adam, A. G. Martín-Caro, M. Huidobro, R. Vázquez, and A. Wereszczynski, Phys. Rev. D 103, 023022 (2021).
  * Tolman (1939) R. C. Tolman, Physical Review 55, 364 (1939).
  * Oppenheimer and Volkoff (1939) J. R. Oppenheimer and G. M. Volkoff, Physical Review 55, 374 (1939).
  * Fiorella Burgio and Fantina (2018) G. Fiorella Burgio and A. F. Fantina, Astrophys. Space Sci. Libr. 457, 255 (2018).
  * Essick  _et al._ (2021) R. Essick, P. Landry, A. Schwenk, and I. Tews, Phys. Rev. C 104, 065804 (2021).
  * Tang  _et al._ (2021) S.-P. Tang, J.-L. Jiang, M.-Z. Han, Y.-Z. Fan, and D.-M. Wei, Phys. Rev. D 104, 063032 (2021).
  * de Tovar  _et al._ (2021) P. B. de Tovar, M. Ferreira, and C. Providência, Phys. Rev. D 104, 123036 (2021).
  * Gil  _et al._ (2021) H. Gil, P. Papakonstantinou, and C. H. Hyun, (2021), arXiv:2110.09802 [nucl-th] .
  * Li  _et al._ (2021) B.-A. Li, B.-J. Cai, W.-J. Xie, and N.-B. Zhang, Universe 7, 182 (2021).
  * Castillejo  _et al._ (1989) L. Castillejo, P. Jones, A. Jackson, J. Verbaarschot, and A. Jackson, Nuclear Physics A 501, 801 (1989).
  * Baskerville (1996) W. K. Baskerville, Phys. Lett. B 380, 106 (1996).
  * Adam  _et al._ (2015b) C. Adam, M. Haberichter, and A. Wereszczynski, Phys. Rev. C 92, 055807 (2015b).
  * Krusch (2006) S. Krusch, Proc. Roy. Soc. Lond. A 462, 2001 (2006).
  * Klebanov (1985) I. Klebanov, Nuclear Physics B 262, 133 (1985).
  * Note (1) Remember that we are using the mostly minus convention for the metric signature.
  * Son and Stephanov (2001) D. T. Son and M. A. Stephanov, Phys. Rev. Lett. 86, 592 (2001).
  * Loewe  _et al._ (2006) M. Loewe, S. Mendizabal, and J. C. Rojas, Phys. Lett. B 632, 512 (2006).
  * Lee  _et al._ (2021) H. K. Lee, Y.-L. Ma, W.-G. Paeng, and M. Rho, (2021), arXiv:2107.01879 [nucl-th] .
  * Danielewicz and Lee (2014) P. Danielewicz and J. Lee, Nucl. Phys. A 922, 1 (2014), arXiv:1307.4130 [nucl-th] .
  * Lee  _et al._ (2011) H. K. Lee, B.-Y. Park, and M. Rho, Phys. Rev. C 83, 025206 (2011), [Erratum: Phys.Rev.C 84, 059902 (2011)], arXiv:1005.0255 [nucl-th] .
  * Park and Vento (2010) B.-Y. Park and V. Vento, “Skyrmion approach to finite density and temperature,” in  _The multifaceted skyrmion_, edited by G. E. Brown and M. Rho (2010) pp. 115–146, arXiv:0906.3263 [hep-ph] .
  * Park  _et al._ (2019) B.-Y. Park, W.-G. Paeng, and V. Vento, Nuclear Physics A 989, 231 (2019).
  * Lattimer  _et al._ (1991) J. M. Lattimer, C. J. Pethick, M. Prakash, and P. Haensel, Phys. Rev. Lett. 66, 2701 (1991).
  * Klahn  _et al._ (2006) T. Klahn _et al._ , Phys. Rev. C 74, 035802 (2006).
  * Alford  _et al._ (2021) M. G. Alford, A. Haber, S. P. Harris, and Z. Zhang, Universe 7 (2021).
  * Brown  _et al._ (2018) E. F. Brown, A. Cumming, F. J. Fattoyev, C. J. Horowitz, D. Page, and S. Reddy, Phys. Rev. Lett. 120, 182701 (2018).
  * Callan and Klebanov (1985) C. G. Callan and I. Klebanov, Nuclear Physics B 262, 365 (1985).
  * Blom  _et al._ (1989) U. Blom, K. Dannbom, and D. Riska, Nuclear Physics 493, 384 (1989).
  * Park  _et al._ (2010) B.-Y. Park, J.-I. Kim, and M. Rho, Phys. Rev. C 81, 035203 (2010).
  * Bander and Itzykson (1966) M. Bander and C. Itzykson, Rev. Mod. Phys. 38, 330 (1966).