Generalized skyrmion crystals with applications to neutron stars

In this article we study properties of isospin asymmetric nuclear matter in the generalized Skyrme model. This is achieved by canonically quantizing the isospin collective degrees of freedom of the recently found multi-wall skyrmion crystal. We obtain, for the first time, an equation of state from the Skyrme model which interpolates between infinite isospin asymmetric nuclear matter and finite isospin symmetric atomic nuclei. This enables us to describe neutron stars with crusts within the Skyrme framework. Furthermore, we observe that the symmetry energy tends to a constant value at zero density, which can be identified with the asymmetry coefficient in the semi-empirical mass formula for atomic nuclei. The symmetry energy also reveals a cusp in its structure below the nuclear saturation point $n_0$ at $n_*\sim 3n_0/4$. This cusp density point $n_*$ can be interpreted as the nuclear density whereby the infinite crystalline multi-wall configuration undergoes a phase transition to a finite isolated multi-wall configuration. Both of these observations are observed to be generic features of skyrmion crystals that tend asymptotically to somewhat isolated skyrmion configurations in the zero density limit. We find that the resulting neutron stars from our study agree quite well with recent NICER/LIGO observational data.


I. INTRODUCTION
The Skyrme model [1] offers a unique, unified framework in which one can study baryonic matter at all scales -from single baryons and atomic nuclei to infinite nuclear matter which, after coupling the model to gravity, gives rise to neutron stars [2].All of this emerges from an elegantly simple Lagrangian containing a limited number of terms and, in consequence, a few free coupling constants, where the fundamental degrees of freedom (d.o.f.) are the lightest mesons disguised into a matrix valued field.In the minimal version, which is used in this work, they are pions forming an SU(2)-valued field.The attractiveness of this approach originates not only in a very small number of parameters but also in the manifestation of baryons.Namely, they are realized as non-perturbative excitations of the mesonic field, that is, as topological solitons, called skyrmions.Importantly, the topological degree of skyrmions has been identified with the baryon charge in a rigorous way [3,4].
The Skyrme model has been very extensively studied in the context of nucleons [5,6], and light atomic nuclei [7][8][9][10][11][12][13][14][15][16] with many spectacular results.In particular, let us mention the description of the ground and Hoyle states in 12 C [17] and excitation bands of 16 O [18,19] as well as the emergence of α-cluster structure [20] which is expected for not too heavy atomic nuclei.This recent progress to large extent relies on an improved quantization procedure where, contrary to the usual rigid-rotor approach, one takes into account both the zero modes and the softest massive vibrations [21].Also, the long standing problem of binding energies has found a resolution by inclusion of additional terms [22][23][24][25][26] or additional mesonic degrees of freedom [27][28][29][30], both physically well motivated.Finally, it is now clear how to extract nuclear forces from the Skyrme model [31,32], which ultimately may provided a much better contact with more traditional nuclear models.
Obviously, a natural field of application of the Skyrme framework is nuclear matter and neutron stars.However, a correct description of this regime is still a serious challenge for the solitonic Skyrme model.
The problem of infinite nuclear matter at non-zero density can be approached if one considers the model on a finite volume unit cell with periodic boundary conditions [33], which results in an infinite but periodic Skyrme crystal.Varying the volume of the unit cell (while keeping the baryon number fixed) allows one to study skyrmionic matter at finite densities and, inter alia, to obtain an equation of state (EoS).Taking the advantage of the Tolman-Oppenheimer-Volkoff construction, one obtains neutron stars.This approach meets some difficulties both at the mathematical and physical level.
First of all, in the traditional approach to determining the EoS, the geometry of the unit cell was fixed to be cubic and the unit cell volume was varied by homothety about the cell center.This rendered the energy minimizing crystalline solutions to inherit the symmetry group of their corresponding initial configuration.In a consequence one obtained not true energy minimizers but solutions with imposed geometrical structure.
Secondly, even if one accepts this constrained approach, three important physical issues have been reported: 1. High density issue: the EoS is too soft, giving rise to neutron stars that are too light.
2. Low density issue: the presence of a minimum at saturation in the EoS yields negative pressure, which represents a thermodynamically unstable phase at low density.
3. Saturation density issue: nuclear binding energies are too large, which in turn means the compression modulus is too large.
The too softness of the standard Skyrme model EoS was found an elegant resolution by extension to the generalized L 0246 -Skyrme model where the so-called sextic term has been included.Indeed, this component of the L 0246 -Skyrme Lagrangian was essential not only to significantly increases the value of the maximal mass of neutron star (from 1.7M ⊙ [40] to above 2M ⊙ [41]), but also to render nuclear matter more like a perfect fluid, especially at higher densities, which corresponds very well to the standard picture of a (super-)fluid core of neutron star.These results are deeply anchored in the mathematical properties of the sextic term.Namely, if treated together with the (pion mass) potential term, the corresponding energy-momentum tensor has a perfect fluid form [42].In addition it enjoys a volume preserving diffeomorphism symmetry which means that the energy of a solution is degenerate up all deformations which do not change its volume [22].On the contrary, deformations that change the volume are strongly penalized as the corresponding EoS has a maximally stiff form [43,44].This agrees with a physical interpretation of the sextic term as a part of the action which effectively arises after integration of ω-mesons.Indeed, EoS' obtained in the Walecka model at large densities tend to maximally stiff EoS' due to the ω-meson repulsion.
At low density the situation is much less clear due to the appearance of thermodynamically unstable regions.This is directly related to the use of the fixed geometry approach mentioned above.For each fixed classical crystal solution (which, in a natural way, is identified with symmetric nuclear matter), the energy E per unit cell possesses a minimum for a certain volume V * , which may be consistently identified with the saturation point.
Obviously, for V > V * , the solution is thermodynamically unstable as it formally corresponds to negative pressure.However, taking into account the isospin quantum corrections and some further contributions, the classical minimum should disappear, thereby providing a thermodynamically stable description even in the low density regime.This periodic crystal was then expected to be replaced by non-homogeneous solutions in this regime [45][46][47][48].Here [49] for example, a crystal of α-particles and B = 32 skyrmions have been considered.Although these configurations lowered the classical energy per cell, they did not cure the instability issue [50].In conclusion, the Skyrme model provided an EoS, but only above the saturation point, leaving the low density regime rather completely unexplored.This is obviously a serious problem of the Skyrme framework as the lower density regime is typically identified with the crust of neutron star.In fact in these densities many geometrically involved phases are expected, see, e.g., lasagna or pasta phases.
Of course, at low densities the isospin and, especially, Coulomb effects should be taken into account.Although it is conceptually clear how it should be done (e.g., taking into account semi-classical quantization of the isospin d.o.f.[51]), and some interesting outcomes have been recently reported [50], the results where obtained for the fixed SC 1/2 crystal of half-skyrmions, which is not the global minimizer at any density [52].
In any case, a thermodynamically stable phase(s) at low density is the first necessary step in approaching of the problem of the crust and nuclear pasta phases within the solitonic Skyrme model.
Finally, there is a very famous problem of the computation of the compression modulus.It is the quadratic term in expansion of EoS of infinite symmetric nuclear matter at the saturation density n 0 .The widely accept value, based on the vibrating frequency extracted from the Isoscalar Giant Monopole Resonance, is K 0 = 240 ± 20MeV.Depending on a version of the Skyrme model with massive pions, as well as on a particular choice of the coupling constants, one gets value a few times bigger than expected.Namely, K 0 ∼ 1350 − 2300 MeV.More importantly, this result was derived for fixed geometry crystals.It has been therefore advocated that the actual value should be lower if a non-homogeneous solution would be the true minimizer.This again brings us back to necessity of solving the generalized L 0246 -Skyrme model at finite density without any geometric constraints.
It is the aim of the current paper to simultaneously face all these problems by constructing generalized skyrmion crystals at finite density without any symmetry assumptions on the skyrmion nor its fundamental period lattice.This is possible due to the recently developed method of obtaining crystalline solutions by not only considering the variation of the Skyrme field φ : R 3 /Λ → SU(2) but also by allowing non-cubic variations of the unit cell period lattice Λ [52].
The main idea is the identification between all 3-tori (R 3 /Λ, d), with d the Euclidean metric, and the unit 3-torus T 3 = R 3 /Z 3 , where T 3 is equipped with the flat pullback metric g = F * d via a diffeomorphism F : T 3 → R 3 /Λ.Varying the metric g on T 3 is equivalent to considering variations of the period lattice Λ.
It is convenient to think of the metric g as a constant symmetric-positive-definite matrix (g ij ).Then one can address this variational problem by identifying the gradient of the energy with respect to the metric (g ij ) with the stress-energy tensor (S ij ) of the field φ.Auckly and Kapitanski [53] showed that, for fixed metric g, the energy functional E(φ, g) attains a minimum.In [52] they proved that, for fixed field configuration φ, any critical metric g of the energy functional E(φ, g) is in fact a unique local minimum.Hence the period lattice Λ, for which the Skyrme field φ has minimum energy, is unique (up to automorphism).This means that the resulting periodic crystalline configuration is indeed a true energy minimizer with respect to both variations of the period lattice Λ and the Skyrme field φ.This slightly improves on the known crystalline solutions at medium and large densities, i.e. above the nuclear saturation point, but it has a tremendous impact on the low density regime where non-homogeneous solutions are expected to exist.
Here, we apply this method to the generalized L 0246 -Skyrme model and obtain the lattice ground state of the generalized model at all densities, that is, above and below the nuclear saturation point n 0 .In our model, infinite nuclear matter is not necessarily treated as being homogeneous.At saturation n 0 , it appears as an almost homogeneous multi-wall configuration with near cubic symmetry.At low densities (n < n 0 ) then it is considered inhomogeneous, with distinct somewhat isolated multi-wall configurations present.Whereas, at high densities (n > n 0 ), e.g. in the core, it appears even more homogeneous and as a simple cubic crystal of fractional half-skyrmions, i.e. it merges with the SC 1/2 crystal.This allows us, for the first time, to obtain an EoS of the skyrmionic matter which interpolates between low and high density regimes.We use this EoS, with an addition of the isospin quantum contribution and with the assumption of β-equilibrium, and investigate its usefulness and consequences for nuclear physics.

II. SKYRME CRYSTALS AND PHASES OF SKYRMION MATTER
A. The generalized L0246-Skyrme model The generalized L 0246 -Skyrme model consists of a single scalar field φ : Σ → SU(2) where spacetime is given by the (3 + 1)-dimensional Lorentzian manifold Σ = R × M with the product metric g = −dt 2 + h, and (M, h) is an oriented 3-dimensional Riemannian manifold with Riemannian metric h.Let us introduce oriented local coordinates (x 0 , x 1 , x 2 , x 3 ) on the domain Σ and let {∂ 0 , ∂ 1 , ∂ 2 , ∂ 3 } be a local basis for the tan-gent space T x Σ at x ∈ Σ, where we have denoted ∂/∂x µ ≡ ∂ µ .We equip su (2) 2) be an su(2)-valued two-form on SU(2) and ω ∈ Ω 1 (SU(2)) ⊗ su(2) be the left Maurer-Cartan form.Then, for any left invariant vector fields X, Y ∈ T φ(x) SU(2), where x ∈ Σ, we define where Let us write the pullback as Then the curvature can be expressed in terms of the su(2)valued left current as Consider the trivial foliation of spacetime Σ = R × M into spacelike hypersurfaces M and let M be compact and without boundary.This is the case if, for example, M is a 3-torus or the usual vacuum boundary condition φ(x → ∞) = I 2 is imposed on M = R 3 .Then Hopf's degree theorem ensures that such mappings φ : M → SU(2) ∼ = S 3 , for M = R 3 ∪ {∞} ∼ = S 3 and T 3 , are characterized by a homotopy invariant: the topological degree B ∈ Z, since π 3 (S 3 ) = H 3 (T 3 ) = Z.This topological degree is identified with the physical baryon number upon quantization, so we often to refer to B as the baryon number, which may be computed using where the topological current is given by We consider the generalization of the massive Skyrme Lagrangian which yields an ω-meson-like repulsion on short distances, while also allowing the quartic Skyrme term to describe scalar meson effects.This generalized Skyrme Lagrangian is composed of four terms and is given by where the the index i denotes the degree of each term as a polynomial in spatial derivatives.The four terms appearing in the energy functional are the potential, Dirichlet, Skyrme and sextic terms, respectively.It is conventional to label the models by terms used in the energy functional, e.g. the generalized model is labeled L 0246 , the standard massive model is denoted L 024 , the massless Skyrme model L 24 and the BPS model L 06 .The first term is the potential which provides a mass for the pionic fields, The Dirichlet, or kinetic, term is given by and the Skyrme term, corresponding to the four pion interaction, is Finally, we include the sextic term, defined by [54] where B µ is the topological Chern-Simons current defined in (5).The c i are coupling constants and, for the usual L 0246 -Skyrme model, take the values c 0 = c 2 = 1 and c 4 = 1/4.The pion mass is fixed to take its physical value of m π = 140 MeV.So, the free parameters of the model are the pion decay constant F π , the dimensionless Skyrme parameter e, and λ which is related to the mass m ω and coupling constant g ω of the ω meson via λ 2 = g 2 ω ℏ 3 /(2π 4 m 2 ω ) [55].The reduced Planck constant is ℏ = 197.33MeV fm.Throughout we will use the values Qualitatively, the parameters (11) don't have much affect on the ground state configuration.However, quantitatively this is not true.We fit the parameters of the model to give us approximately the binding energy at saturation and the nuclear density, while also allowing the symmetry energy and the pion decay constant not to deviate too much from their experimental values.Other studies have done similar fittings to, e.g., the symmetry energy, but there is always a trade-off where if you fix one parameter accurately then another physical quantities will suffer in consequence.In other studies [50], the symmetry energy and saturation energy can be fitted correctly, but the saturation density can not also be simultaneously fitted correctly.That is the caveat of using the Skyrme model alone to model nuclear matter.For example, in our model, the symmetry energy at saturation is lower than expected but accurately predicts the asymmetry coefficient in the SEMF.If the model is tuned to give the correct symmetry energy value at saturation then the asymmetry coefficient would be off.For a more general review of the quantitative effects of the free parameters on a ground state configuration, see [49,50].
We are interested in static solutions and adopt the usual Skyrme units of length and energy.The classical energy scale is Ẽ = F π /4e (MeV) and the length scale is L = 2ℏ/eF π (fm).Thus the quantum energy scale is defined by ħ = 2e 2 .In these dimensionless Skyrme units, the rescaled pion mass for our studies is and the dimensionless sextic coupling constant is It will prove useful throughout to introduce the Hilbert energy-momentum tensor (in dimensionless Skyrme units): The static energy functional can be obtained from the timelike part of the energy-momentum tensor, T 00 = E stat + E kin , and is given by where A field configuration φ which minimizes the static energy functional (15), for some choice of domain metric g, is referred to as a skyrmion and the static energy M B is often interpreted as the classical mass of the skyrmion.The associated Euler-Lagrange field equations can be approximately solved by discretizing the static energy (15) and employing a 4th order central finite-difference method.This is carried out using the quaternionic formulation detailed below.We can then regard the static energy as a function M B : C → R, where the discretised configuration space is the manifold C = (S 3 ) N1 N2 N3 ⊂ R 4 N1 N2 N3 .To solve the Euler-Lagrange field equations we use arrested Newton flow: an accelerated gradient descent method with flow arresting, with some appropriate initial configuration.That is, we are solving the system of 2nd order ODEs with initial velocity φ(0) = 0. Setting ψ := φ as the velocity with ψ(0) = φ(0) = 0 reduces the problem to a coupled system of 1st order ODEs.We implement a 4th order Runge-Kutta method to solve this coupled system.In general, the initial configuration φ 0 is not a minimizer and so it swaps its potential energy for kinetic energy as it evolves.During the evolution we check to see if the energy is increasing.If the energy is indeed increasing, we take out all the kinetic energy in the system by setting ψ(t) = φ(t) = 0 and restart the flow (this is the arresting criteria).Naturally the field will relax to a local, or global, minimum in some potential well.The evolution then terminates when every component of the energy gradient δM B δφ is zero within some specified tolerance, e.g.tol = 10 −5 .

B. Metric independent integral formulation
For numerical purposes, it is convenient to utilize the quaternionic representation of the target group SU(2), which is topologically isomorphic to S 3 .Let us parameterize the unit quaternion φ ∈ H by the mesonic fields (φ 0 , φ 1 , φ 2 , φ 3 ): with the unitary condition σ 2 + ⃗ π • ⃗ π = 1, where ⃗ π = (φ 1 , φ 2 , φ 3 ) is normally identified with the triplet of pion fields and σ = φ 0 with the σ-field.Then the Maurer-Cartan left current can be expressed as the vector quaternion: where τ a are the isospin Pauli matrices and, similarly, the curvature in the quaternionic representation is given by From this we get the following contractions, The baryon number density in contraction form is For numerical simulations involving the minimization of the energy functional with respect to variations of the metric, it will be convenient to define the metric inde-pendent integrals: In terms of these metric independent integrals, the static energy can be compactly written as C. Skyrme crystals Our aim is to study Skyrme fields φ : R 3 → SU(2) that are periodic with respect to some 3-dimensional period lattice Λ, i.e. we impose the condition φ(x + X) = φ(x) for all x ∈ R 3 and X ∈ Λ.We can equivalently interpret the field as a map φ : R 3 /Λ → SU (2), where (R 3 /Λ, d) is a 3-torus equipped with the standard Euclidean metric d.
In particular, we define a Skyrme crystal to be an energy minimizing map where R 3 /Λ ⋄ is some fixed 3-torus such that the field φ is also critical and stable with respect to variations of the lattice Λ about Λ ⋄ .The problem of determining Skyrme crystals was addressed by Harland et al. [52].They prove that, for a fixed field configuration φ, there is a unique period lattice Λ ⋄ (up to automorphism) that minimizes the static energy M B .Therefore, the problem of determining skyrmion crystals is one of finding critical points of the static energy functional (15) with respect to variations of both the field φ and the period lattice Λ ⋄ .For massless L 24 -skyrmions, the period lattice can be determined explicitly.However, only a numerical approach seems possible for generalized L 0246 -skyrmions.For some initial period lattice Λ 0 , the static energy can minimized with respect to variations of the period lattice using the method detailed in §II E. In tandem, with some appropriate initial field configuration φ 0 , the static energy functional can also be minimized with respect to variations of the field by using arrested Newton flow (ANF), which is detailed in §II A.
Skyrme crystals have been studied extensively in the literature, with it being previously accepted that the SC 1/2 crystal found independently by Kugler & Shtrikmann [35] and Castillejo et al. [37] is the minimum energy Skyrme crystal.However, in the massless L 24 -Skyrme model, this SC 1/2 crystal is just one point on an SO(4) orbit of solutions, i.e. it is not an isolated critical point and all of these solutions are all energy degenerate.When the pion mass is turned on, there is no reason to expect these degenerate L 24 critical points to extend to L 0246 critical points upon perturbation.However, there are four critical points which survive perturbation as argued by [52].These are the SC 1/2 , α, chain and multiwall crystals.Each crystal has baryon number B cell = 4 per unit cell, with three of the crystals having lower energy classically than the SC 1/2 crystal for non-zero pion mass and non-cubic (trigonal) lattice geometry.
The SC 1/2 crystal can be obtained from the Fourier series-like expansion of the fields as an initial configuration [37], and cyclic, where s i = sin(2πx i /L) and c i = cos(2πx i /L).From the SC 1/2 crystal, the other three crystals can be constructed by applying a chiral SO( 4) and minimizing the energy with respect to variations of the field and the lattice.These chiral transformations Q ∈ SO(4) can be determined by considering a deformed energy functional on the moduli space of critical points of the Skyrme energy functional, and are found to be [52] The other three rows of the chiral transformations Q α , Q multi-wall and Q chain , labeled by the asterisk, can be obtained by using the Gram-Schmidt process.
Out of the four crystal configurations, the most of interest to astrophysicists are the α-crystal, chain-crystal and the multi-wall crystal; these resemble non-uniform phases of nuclear matter, known as nuclear "pasta".The iron rich crust of a neutron star could be modeled by B = 56 chunks of α-particle crystals, such as those modeled by Feist et al. [56], describing the "gnocchi" phase.As we descend deeper towards the outer core, the pressure due to gravity increases and nuclei are squeezed together into long thin tubes of "spaghetti".This spaghetti phase can be modeled using the chain-crystal.Deeper still and the spaghetti flattens into parallel multi-walls, resembling "lasagna", of which the multi-wall crystal is reminiscent of.Of course, for realistic applications the Coulomb interaction must be added.This is because of the fact that different crust phases arise due to a balance between the strong and electrostatic forces.Nevertheless, the Skyrme model has a built-in ability to model such phases.
The multi-wall-crystal is the lowest energy solution at all baryon densities and also yields a lower compression modulus than the other three crystals.This makes it an ideal candidate for nuclear matter and an equation of state (EoS) at high and low densities.With φ 0 = Q multi-wall φ SC 1/2 as an initial configuration and by considering fixed baryon density variations, as laid out in §II F, the energy-volume curve can be computed and an EoS obtained.

D. The stress-energy tensor
To determine Skyrme crystal solutions, we identify every 3-torus (R 3 /Λ, d), equipped with the standard Euclidean metric d, with the unit 3-torus (T 3 , g) where g is a Riemannian metric and Let the Skyrme field be the map φ • F : We vary the metric g s on T 3 with g 0 = F * d which is equivalent to varying the lattice Λ s with Λ 0 = Λ.The energy minimized over variations g s of the domain metric is equivalent to determining the energy minimizing period lattice Λ ⋄ .Now let the static Skyrme field be the smooth map φ : T 3 → SU(2).Let (x 1 , x 2 , x 3 ) be oriented local coordinates on T 3 and {∂ 1 , ∂ 2 , ∂ 3 } be a local frame for the tangent space T x T 3 at x ∈ T 3 .Let g s be a smooth oneparameter family of metrics on T 3 with g ), a symmetric 2-covariant tensor field on T 3 .Denote the inner product on the space of 2-covariant tensor fields of the tangent space T x T 3 to T 3 at x ∈ T 3 by ⟨•, •⟩.Then for any pair of symmetric bilinear forms A, B we have In particular, we have the following result: Let us consider the rate of change of the energy of the Skyrme field φ with respect to varying the domain metric g.The first variation of the energy with respect to the variation g(s) of the metric on T 3 is given by dM B (φ, g s ) where S(φ, g) ∈ Γ(⊙ 2 T * T 3 ) is a symmetric 2-covariant tensor field on T 3 , known as the stress-energy tensor, defined by This stress-energy tensor is related to the spatial part of the (static) energy-momentum tensor,

E. Numerical optimization of the lattice geometry
Let us fix the field φ : T 3 → SU(2) and think of the energy M B as a function of the metric g on T 3 .That is, we define a map E φ : SPD 3 → R such that E φ := M B ( φ| fixed , g), where SPD 3 is the space of symmetric positive-definite 3 × 3-matrices.To minimize the energy functional E φ with respect to variations of the metric g s , we use arrested Newton flow on SPD 3 .The essence of the algorithm is as follows: we solve Newton's equations of motion for a particle on SPD 3 with potential energy E φ .Now let g s be a smooth one-parameter curve in SPD 3 with g 0 = F * d.Explicitly, we are solving the system of 2nd order ODEs with initial condition (g ij ) 0 = ⃗ X i • ⃗ X j , and where S φ = S(g) is the stress-energy tensor for fixed field configuration φ.Setting δg s = ∂ s g s as the velocity with initial velocity δg 0 = ∂ s g s | s=0 = 0 reduces the problem to a coupled system of 1st order ODEs.We implement a 4th order Runge-Kutta method to solve this coupled system.The components of the stress-energy tensor for fixed field φ, given in the metric independent integral formulation, reads In general, the dimension of SPD n is dim(SPD n ) = n(n + 1)/2.In our case, we are working with SPD 3 and consider the energy as a function E φ : SPD 3 → R.So we are implementing arrested Newton flow on a 6 dimensional manifold.After each time step t → t+δt, we check to see if the energy is increasing.If E φ (t + δt) > E φ (t), we take out all the kinetic energy in the system by setting δg(t + δt) = 0 and restart the flow.The flow then terminates when every component of the stress-energy tensor S φ is zero to within a given tolerance (we have used 10 −6 ).

F. Phases of skyrmion matter
Determining phases of nuclear matter and phase transitions in the Skyrme model is a difficult task, and is important if one wants to understand symmetric and asymmetric nuclear matter in high/low density regimes.To study phases of matter at various densities, we consider fixed density variations of the energy functional, i.e. we allow the lattice to vary but keep its volume fixed.Then the volume form vol g is required to be invariant under variations g s of the metric, viz.
That is, δg has to be an element of the space of traceless parallel symmetric bilinear forms E 0 .
In terms of the energy, we are dealing with a constrained minimization problem: minimize the energy functional for fixed field configuration φ = φ| fixed subject to the constraint that det(g) = constant.We can approach this using the method of Lagrange multipliers.This leads to modifying the stress-energy tensor in (34) via the mapping and our convergence criterion becomes max( Sφ ) < tol.Likewise, to ensure numerically that δg is traceless, we need to project out the component of variation vector δg parallel to the metric g via the mapping By employing this process at various volumes it enables us to determine an energy-volume curve or, equivalently, an energy-density curve.This is key to obtaining an EoS within our framework, as the EoS is directly related to the E − V curve.

G. The results
The first main result of this section is the observation that, as it is for the massive L 024 -Skyrme model, the multi-wall crystal is also the ground state crystalline solution for the generalized L 0246 -Skyrme model at all densities.In the low density regime the solution clearly exhibits a two-layer structure, extending parallel to the xy-plane with the vacuum (σ = 1) occupying the regions above and below the multi-wall.This can be seen in Fig. 1.Inside the multi-wall center the σ-field is approximately the anti-vacuum (σ ≈ −1).Therefore, the multi-wall crystal is similar to that of a domain wall crystal, hence the name convention.As the density increases, the regions occupied by the vacuum reduces and the noncubic period lattice becomes more cubic, tending asymptotically to the SC 1/2 crystal in the zero volume limit.In Fig. 2 we plot the classical static energy per baryon E = M B /B of the multi-wall crystal as a function of the baryon density n B .This is interpreted as EoS of the symmetric nuclear matter since the classical Skyrme model does not distinguish between protons and neutrons.
Expansion of the energy function E(n B ) around the minimum n 0 gives As always the local minimum which is identified with the nuclear saturation point with saturation energy E 0 .The curvature of the energy curve is controlled by the compression modulus K 0 and determines the increase in energy due to compression.For our choice of the coupling constants (11) the saturation energy per baryon and saturation density are respectively E 0 = 912 MeV and n 0 = 0.160 fm −3 , which almost perfectly corresponds to the physical values of the saturation energy and density.An important observation is that the difference between the energy at nuclear saturation and the classical energy at zero density is much smaller than in previous works.Indeed, the difference is now roughly ∆E ≈ 7 MeV, which is a 0.8% difference with respect to the total energy.Whereas, for a B = 32 or B = 108 α-crystal the difference is found to be approximately 3% and 1.7%, respectively.This small difference in energy between the nuclear saturation and low-density asymptotic solutions is crucial for the existence of a purely skyrmion generated EoS at all densities.Unfortunately, the compression modulus still exceeds the experimental value by a factor of 6 ∼ 7. Although, in comparison with studies involving the SC 1/2 crystal, where K 0 ∼ 1700 MeV, we do observe a significant improvement in the (in)compressibility by approximately 500 MeV (K 0 = 1169), the non-homogeneous solution alone cannot solve the compression modulus problem in the Skyrme model.Nevertheless, this negative result is very important as it shows that the purely pionic Skyrme model cannot lead to a physically acceptable compression modulus.Consequently, it seems to be necessary to include other mesonic d.o.f. which may soften the skyrmionic matter at the saturation point.

III. QUANTUM SKYRMION CRYSTALS AND THE SYMMETRY ENERGY
In general, the full symmetry group of the generalized L 0246 -Lagrangian ( 6) is the direct product of the Poincaré group and chiral group: G = O(3) ⋉ R 3 × SO(4) chiral .However, static energy minimizers break the Poincaré symmetry group O(3) ⋉ R 3 to the Euclidean subgroup E 3 = SO(3) × R 3 , corresponding to spatial translations and rotations.The resulting symmetry group of the static energy functional ( 15) is thus The action of this group on the Skyrme field is given by where For skyrmions on M = R 3 , one must impose finite boundary conditions φ(x → ∞) = I 2 .This allows for the compactification of the domain R 3 {∞} ∼ = S 3 and further reduces the symmetry group G to the subgroup where SU(2) I is the isospin internal symmetry group.The corresponding action of the subgroup H on the Skyrme field is given by the transformation (40) with When considering crystals on M = R 3 /Λ, one must be careful when defining the isospin subgroup SU(2) I ; the vacuum boundary condition is no longer imposed and there is not a natural way to select the diagonal isospin subgroup SU(2) I .This problem was addressed by Baskerville [51] in the context of the SC 1/2 crystal in the L 24 -model, wherein she considered full SO(4) chiral rotations.She deduced that there are two cubic point groups that can define the SC 1/2 crystal, one of which is related to the centers of the half-skyrmions.The cubic point group symmetry corresponding to the half-skyrmion centers is reducible into the trivial 1-dimensional irreducible representation and a 3-dimensional irrep.We choose the σ = φ 0 field to transform in the 1-dimensional irrep.Then the isospin group SU(2) I can be defined as the group of isorotations of the pion fields ⃗ π = (φ 1 , φ 2 , φ 3 ), corresponding to transformations in the 3-dimensional irrep.If the pion mass potential term L 0 is included then this is a natural choice of isospin group SU(2) I .

A. Isospin quantization
As a field theory, the Skyrme model is nonrenormalizable.One must quantize it semi-classically.It is well-known that the classical dynamics of slowly moving solitons corresponds to geodesic motion on the moduli space of static configurations [57].Minimal energy configurations in the Skyrme model are unique, for a given baryon number B, up to actions of the symmetry group H = E 3 × SU(2) I .The classical configuration space Q of skyrmions is split into connected components, labeled by the baryon number B, Q = B∈Z Q B .The covering space QB of each component is a double-cover with a two-to-one map π Q : QB → Q B [58].It was argued by Finkelstein and Rubinstein [59] that the wave functions Ψ ∈ H must be defined on the covering space of the configuration space Q, where H is a formal Hilbert space such that Ψ is normalizable and square integrable.That is, the wave functions are defined by the map Ψ : Q → C. We make a simple approximation and require the wave function Ψ to be non-vanishing only on minimal energy configurations and their symmetry orbits.This quantization procedure is known as rigid-body, or zero mode, quantization.
In the zero mode quantization method, a skyrmion is treated as a rigid body that is free to translate and rotate in physical space and also rotate in isospace, with the action defined by (40).These solutions are all degenerate in energy and this classical degeneracy is removed when one quantizes the theory.We wish to quantize the isorotational degrees of freedom and work in the zeromomentum frame, ignoring the translational and spin degrees of freedom.The action of the group of isorotations SU(2) I on the Skyrme field φ is defined by the mapping φ(x) → Aφ(x)A † .Semi-classical quantization is performed by promoting the the collective coordinate A ∈ SU(2) to a dynamical degree of freedom A(t) [8].The dynamical ansatz for the Skyrme field is then given by the transformation Define the isorotational angular velocity ⃗ ω to be A † Ȧ = i 2 ω j τ j such that ω j = −i Tr(τ j A † Ȧ).Then, under the dynamical ansatz (41), the Maurer-Cartan left current transforms as where The dynamical ansatz (41) induces a rotational kinetic term in the energy functional, which is given by where the Chern-Simons current transforms as (44) The restriction of the kinetic energy functional of the model to the isospin orbit of a given static solution defines a left invariant metric on SO(3) called the isospin inertia tensor, which is the symmetric 3 × 3-matrix given by where the isospin inertia tensor density is Therefore, the effective Lagrangian on this restricted space of configurations is L eff = L rot − M B , where M B is the static mass of the skyrmion defined by ( 15) and L rot is the induced isorotational part of the Lagrangian given by The rigid-body wavefunctions will be on SU(2) with isospin half-integer if B is odd and integer if B is even.The isorotation angular momentum operator canonically conjugate to ⃗ ω is the body-fixed isospin angular momentum operator ⃗ K, defined by This is related to the usual space-fixed isospin angular momentum ⃗ I by the relation where A ∈ SU(2) has been recast in the SO(3) form via the map These two classical momenta are promoted to quantum operators ⃗ K and ⃗ Î, both satisfying the su(2) commutation relations, and the Casimir invariants satisfy ⃗ Î2 = ⃗ K2 .On the double cover of the group of isorotations SU(2) I , there is a basis of rigid-body wavefunctions where I is the total isospin quantum number, K 3 is the third component of ⃗ K and I 3 is the third component of isospin relative to the spacefixed axes (in units of ℏ) as defined in nuclear physics.

The operator
⃗ Î2 has eigenvalue I(I + 1) and I 3 the eigenvalue for the operator Î3 .The rigid-body Hamiltonian takes the general form For Skyrme crystals, we can set the principal axes of inertia to be the usual orthogonal axes such that U ij = 0 for i ̸ = j.Let us label the eigenvalues U i = U ii , then the quantum Hamiltonian takes the form The energy eigenstates of the Hamiltonian ( 52) can be classified by I and I 3 .To determine bound states with definite energy one must solve the static Schrödinger equation corresponding to this Hamiltonian, H |Ψ⟩ = E |Ψ⟩.The Schrödinger equation can be expressed more explicitly within a particular (I, I 3 ) sector by expanding the quantum state |Ψ⟩ in terms of the total wavefunctions Ψ as with q ∈ Q and substituting this into the Hamiltonian (52).

B. Symmetry energy and the cusp structure
So far we have only considered symmetric nuclear matter, which we have described by using the classical multiwall skyrmion crystal.In order to study nuclear matter in neutron stars we must consider isospin asymmetric nuclear matter, whereby a small fraction of protons are permitted.Now let us consider asymmetric nuclear matter with baryon number B = N + Z, where N is the number of neutrons and Z the number of protons.The asymmetry of such matter is determined by the isospin asymmetry parameter δ = (N − Z)/(N + Z) = 1 − 2γ, where γ is the proton fraction.We define the nuclear density to be n B = B/V , with the nuclear saturation density n 0 defined to be the nuclear density such that (∂M B )/(∂n B )| n B =n0 = 0. Then the binding energy per baryon number of asymmetric nuclear matter is given by The two terms appearing in the asymmetric binding energy (54) are the binding energy of isospin-symmetric matter E N and the symmetry energy S N .In terms of our model, the symmetric binding energy is defined by E N = (M B − BM 1 )/B.The symmetry energy S N dictates how the binding energy changes when going from symmetric (δ = 0) to asymmetric (δ ̸ = 0) nuclear matter.We can expand the isospin symmetric binding energy E N and the symmetry energy S N around the saturation density n 0 for symmetric matter [60], where ϵ = (n B − n 0 )/n 0 , K 0 is the incompressibility at the saturation point and S 0 = S N (n 0 ) is the symmetry energy coefficient at saturation.We remind ourselves that, for our choice of coupling constants (11), the nuclear saturation point is characterized by the density n 0 = 0.160 fm −3 and energy (per baryon) M B /B = 912 MeV.The higher-order coefficients, L sym and K sym , appearing in the symmetry energy S N are defined as The precise values of these coefficients are not known, but are predicted to be L sym = 57.7 ± 19 MeV and K sym = −107 ± 88 MeV [61].
Consider an infinitely extended and rigidly iso-spinning Skyrme crystal with each unit cell containing baryon number B cell .In order to calculate the isospin correction to the energy of the crystal we would need to know the quantum state of the whole crystal.This is obviously a very difficult computation since the crystal is infinitely extended and is therefore composed of an infinite number of baryons.However we may impose the following restrictions to solve this problem: • The total isospin quantum state of the crystal |Ψ⟩ is written as the superposition of each individual unit cell state |ψ⟩.That is |Ψ⟩ = ⊗ N cell |ψ⟩, where N cell → ∞ in the thermodynamic limit.
• The symmetry of the classical configuration in each unit cell is extended to the whole crystal, so both wavefunctions share the same point symmetry group.
Under these assumptions, and since we have B cell = 4 within our unit cell, there are a finite number of possible quantum states with allowed quantum numbers I = 0, 1, 2 [50].The I 3 = 0 case, which corresponds to symmetric nuclear matter, would be the one with the lowest energy since it has no isospin energy compared to the other cases.This is obviously the most symmetric state possible.However, it is known that inside neutron stars there is a huge asymmetry between protons and neutrons.Baskerville investigated the charge neutral case I 3 = −2, corresponding to a pure neutron crystal, and computed the quantum isospin corrections to the energy [51].However, a realistic description of neutron stars would require the presence of protons.Although the concrete value is still unknown, simulations yield values around γ ∼ 10 −2 −10 −1 [62,63].Therefore, following the arguments in [50] we perform a mean-field approximation considering a larger chunk of crystal, enclosing an arbitrary number of unit cells N cell , which is in a generic quantum state with fixed eigenvalue, Note that in this case the nuclear density of the crystal chunk can be directly interpreted as the nuclear density of the unit cell, since In previous applications of skyrmion crystals to model neutron stars (see, for example, [49,50,[64][65][66]), the SC 1/2 crystal was considered.This crystal has an isotropic inertia tensor with eigenvalue U i = λ, with λ some constant.However, the multi-wall crystal considered in this paper is not isotropic and the isospin inertia tensor generically has the eigenvalues The Schrödinger equation corresponding to such a rigidly iso-spinning crystal with N cell unit cells can be written as where the isospin correction to the energy of the crystal is given by It should be noted that in addition to the quantum numbers I, I 3 being density n B dependent, the inertia tensor is also density dependent, that is The eigenvalue I 3 is already fixed from the mean-field approximation (58), and the value of I = I 3 is the one which minimizes the isospin energy, since by definition I 2 ≥ I 2 3 .In the thermodynamic limit N cell → ∞ we obtain a final expression for the quantum correction (per unit cell) to the energy due to the isospin degrees of freedom, This quantum isospin energy is explicitly related to the proton fraction γ, and so we will need to include leptons if we are to allow the crystal to have a non-zero proton fraction.This is required in order for the system to remain electrically neutral.Thus the proton fraction, and hence the quantum state of the crystal, will be obtained by imposing β-equilibrium for each value of the density.
From the quantum isospin energy (62), we can determine the nuclear symmetry energy of the multi-wall crystal, which in general plays a crucial role in the structure of neutron-rich nuclei and, of more interest to us, in neutron stars.For general skyrmion crystals the symmetry energy is given by where the eigenvalue U 3 of the isospin inertia tensor ( 45) is implicitly dependent on the nuclear density n B .We determine the symmetry energy at at saturation to be S 0 = 22.7 MeV, which is in okay agreement with the experimentally observed value S 0 ∼ 30 MeV [67].The resulting symmetry energy curve S N (n B ) for the multiwall crystal is plotted in Fig. 3. Having obtained the symmetry energy curve we can determine its slope and curvature, which are computed at the nuclear saturation point.We find that they are, respectively, L sym = 36.6 MeV and K sym = −15.1 MeV.
Let us now summarize the results obtained for the multi-wall crystal.First of all, we find that at lower densities the isospin moment of inertia, and specifically its eigenvalue U 3 , tends to a constant value.This is an obvious consequence of the inhomogeneous nature of the solution which, in the limit V cell → ∞, tends to an "isolated" multi-wall configuration on M = S 1 ×S 1 ×R.This simple fact has an important consequence.Namely, it leads to a non-zero value of the symmetry energy at zero density, S N (0) = 23.8MeV.At a first glance, this seems to be in contradiction with the standard description of nuclear matter where the symmetry energy vanishes at zero density.However, we want to argue that this is a desirable property of the Skyrme model as it indicates a smooth transition between infinite nuclear matter and finite atomic nuclei.Indeed, the asymmetry energy in the Bethe-Weizsäcker SEMF reads where a A ≈ 23.7 MeV.Thus, our symmetry energy at zero density can be directly identified with a A with excellent agreement.
We remark that the assumed identification here between the zero density symmetry energy and the asymmetry energy in the Bethe-Weizsäcker formula is not a unique possibility.In fact, in the seminal paper by Natowitz el.al. [68] they computed the symmetry energy of the low density, warm nuclear matter using a quantumstatistical approach.Their results agree amazingly well with values extracted from heavy-ion collisions [69].The symmetry energy, still taking a non-zero value at zero density, is approximately only one fourth of its value at saturation n 0 .It would definitely be very desirable to investigate whether the Skyrme model may lead to similar results or not.
Moving away from zero nuclear density towards n * ∼ 3n 0 /4, the isospin energy and consequently the symmetry energy slowly decreases, as can be seen in Fig. 3.This again is not an unexpected result in the Skyrme model.It was noticed by Kopeliovich et al. [70] that the careful analysis of mass splittings of nuclear isotopes leads to the symmetry energy decreasing with increasing baryon number B. Here, we reproduce this result, however, using a completely different setup, i.e. the collective coordinate quantization of the crystal ground state.
Below the nuclear saturation point n 0 at the density n * ∼ 3n 0 /4, the symmetry energy exhibits a cusp structure.This cusp also seems to be a generic feature of the Skyrme model, independent of the choice of values for the coupling parameters (11) but rather can be interpreted as the point where the multi-wall crystal begins transitioning to an "isolated" multi-wall.On the other hand, its position with respect to the saturation point certainly may be affected by a choice of the model parameters.One can also expect such a cusp to be present where a crystalline configuration transitions to an isolated configuration at zero nuclear density, e.g. for the α and chain crystals.It is interesting to remark that such a cusp, albeit above the saturation density n B > n 0 , has been advocated in [71,72] as an effect of an assumed topological phase transition from the FCC crystal of B = 1 hedgehogs to the SC 1/2 crystal of fractional skyrmions as the nuclear density grows.Although, in reality such a transition does not occur in the Skyrme model as it is found to occur in the thermodynamically unstable regime n B < n 0 [48].To conclude our findings on the symmetry energy cusp, we propose that the origin of the cusp can be associated with a phase transition between an infinite crystalline state and a somewhat isolated state that is non-homogeneous and nucleated.

IV. PARTICLE FRACTIONS OF npeµ MATTER IN β-EQUILIBRIUM
For a more realistic description of cold nuclear matter inside neutron stars we need to consider not completely asymmetric nuclear matter.As was shown in the previous section, this can be achieved by allowing a small fraction of protons over neutrons.The presence of protons gives the crystal positive electric charge, so we need to include a background of negatively charged leptons to neutralize the system.To determine the proton fraction γ at a prescribed nuclear density n B we impose charge neutrality and β-equilibrium conditions, and then we solve the underlying equilibrium equation.Additionally, the presence of protons would require the inclusion of Coulomb interaction within the unit cell and between neighbouring cells.It has been argued [33] that the contribution of this energy diverges in the crystal due to infinitely many interactions between the cells.However, including a background of negatively charged particles in the system suppresses the Coulomb interaction between neighbouring cells and hence has a negligible contribution to the energy [50].
In the neutron star interior, the interaction between leptons and nuclear matter is mediated by the weak force.We can describe the exchange of leptons and nucleons by electron capture and β-decay processes, respectively, These processes take place simultaneously at the same rate, assuming that the charge neutrality, and the β-equilibrium conditions [73], are satisfied.Here µ I is the isospin chemical potential given by Leptons inside a neutron star are treated as a noninteracting, relativistic, highly degenerate Fermi gas.The corresponding chemical potential for such a type of lepton is given by [64] where k F = (3π 2 n l ) 1/3 is the associated Fermi momentum and m l the lepton mass.For electrons we take the ultra-relativistic approximation µ e ≈ ℏk F,e .From the charge neutrality condition (66), the electron number density is The β-equilibrium condition (67) for electrons yields the following relation and for muons gives In the low density regime the electron chemical potential will be smaller than the muon mass, µ e < m µ .So we can solve (71) in the low density regime considering only electrons, by setting n µ = 0 until µ e ≥ m µ .Once the electron chemical potential µ e reaches the muon mass m µ = 105.66MeV at high densities, it will be energetically favourable for muons to appear.Then we solve ( 71) and ( 72) simultaneously [64], and construct the proton fraction curve γ = γ(n B ).
In Fig. 4 we plot the particle fractions of npeµ matter in β-equilibrium for the multi-wall crystal.Note that the cusp structure present in the symmetry energy, or equivalently in the isospin energy, results in an appearance of a similar structure in the particle fractions.This reinforces the proposition that the cusp density point n * is the density at which a phase transition between isospin asymmetric infinite nuclear matter and symmetric finite matter begins.Furthermore, the fact that the symmetry energy S N tends to a constant value at zero density leads to a similar behavior for the proton, neutron and electron particle fractions.Namely, they take their minimal/maximal value at n * then they increase/decrease as zero density is approached.This is once again a direct consequence of a non-zero value of the isospin moment of inertia at this limit and, therefore, a generic feature of the Skyrme model.We remark that at zero density n B = 0, which, in the Skyrme model framework, can be interpreted as a limit where we find nuclei in the vacuum, the nuclear matter becomes totally isospin symmetric with γ p (0) = 0.5.This corresponds quite well to the proton fraction in 56 Fe, γ p = 0.46, which is the element expected to be present in the crust of neutron stars [74].Further, it appears that there is a phase transition at (n/n B = 0.91, p = 0.023MeV fm − 3).The n * density occurs in this region of constant pressure, so it could very well be related to the liquid-gas phase transition.
We remark that at the high density, which corresponds to the core of neutron star, the proton fraction is quite small.This agrees with previous computations in the Skyrme model with the SC 1/2 crystal [50].Fortunately, inclusion of strange d.o.f.resolves this issue and brings the proton fraction to the widely accepted ∼ 0.4 value, see [64].We expect that the same mechanism applies for the multi-wall crystal.Especially considering this ground state crystalline solution and the SC 1/2 crystal are basically identical at high density.On the other hand, inclusion of Kaon condensate does not have any impact on the low density regime.
We now summarize our findings and compute the total energy per unit cell in a β-equilibrated multi-wall skyrmion crystal, that is where the isospin energy for a β-equilibrated crystal is given by The lepton energies are the energies of a relativistic Fermi gas at zero temperature, The crucial observation is that, in the case of the multiwall skyrmion crystal, the inclusion of the β-equilibrated isospin energy and lepton energies does not completely erase the small minimum in the classical energy M B .Strictly speaking there is still a very shallow minimum at a density smaller than the saturation density, n B = 0.146 fm −3 .For smaller densities the total energy grows, until a small maximum is reached.After that the total energy decreases as the nuclear density approaches the zero density limit, n B → 0. Importantly, the asymptotic value of the total energy per unit cell is smaller than the energy at the minimum.This means that, although the total energy per unit cell still possesses a thermodynamically unstable region, we can take advantage of the Maxwell construction and derive an EoS which is valid at all densities.This is a valid construction and has a minute affect on the EoS since the difference in energy between the asymptotic solution and the minimum is ∆E ∼ 0.1 MeV.The formulation of the Maxwell construction is detailed below and the resulting β-equilibrated asymmetric nuclear matter is plotted in Fig. 5, alongside the classical isospin symmetric matter and the pure neutron matter.The pure neutron matter is obtained for the entirely isospin asymmetric case δ = 1 with I 3 = −2.Unlike the β-equilibrated matter, the pure neutron matter approaches the asymptotic solution from below.This is due to the non-vanishing of the quantum isospin energy contributions E iso (n B ) in the zero density limit n B → 0. Consequently, the Maxwell construction cannot be used on the pure neutron matter EoS.
We remark that for the α-crystal the total energy in the zero density limit is greater than the energy at the minimum, so the Maxwell construction is not possible.On the other hand, for B = 32 and B = 108 crystals constructed from α-particles, such a construction is possible but it extends over a non-physical range of densities and occurs for relatively high values of the pressure.For example, the neutron stars obtained from these crystals would almost be entirely made from the Maxwell construction phase.
The Maxwell construction (MC), or equal area rule, is implemented as follows.We find three points that have the same gradient/pressure, i.e. p(V i ) =: p MC .These three points are chosen such that the area enclosed between p([V 1 , V int ]) and p MC is equal to the area enclosed between p([V int , V 2 ]) and p MC , where p([V 1 , V int ]) ≤ p MC and p([V int , V 2 ]) ≥ p MC .This ensures that the total energy of the thermodynamic system remains the same while implementing this construction.Then, in the corresponding MC density regime V 1 < V cell < V 2 , the total energy function is replaced by a straight line connecting E(V 1 ) and E(V 2 ).The resulting total energy per unit cell function can be summarized as (76) Now we are in a position to determine the EoS for the multi-wall configuration.The multi-wall crystal EoS for isospin asymmetric nuclear matter can be obtained by defining the energy density ρ and pressure p as, respectively, This EoS ρ = ρ(p), generated purely from the generalized multi-wall skyrmion crystal, is valid at all densities.In our case, the pressure at which the Maxwell construction is applied is quite small, p MX = 0.023 MeV fm −3 , which corresponds to an energy difference of ≈ 0.1 MeV over a large density range (0.91n 0 to 0.36n 0 ).The resulting EoS is shown in Fig. 7, alongside the EoS without the Maxwell construction applied.
Although the obtained equation of state covers the full range of densities one has to be aware that the multiwall crystal does not describe the low density regime in its entirety.As we have already mentioned, to get a more realistic description of the crust the electrostatic interaction should be included.This can have an impact on the structure and symmetry of the skyrmions, which could potentially lead to the appearance of other non-homogeneous solutions with different baryon numbers per unit cell.

V. NEUTRON STARS FROM QUANTUM SKYRMION CRYSTALS COUPLED TO GRAVITY
In order to describe neutrons stars within the Skyrme framework, we need to couple the generalized Skyrme model to gravity.We do this by introducing the Einstein-Hilbert-Skyrme action [75] where G = 1.3238094 × 10 −42 fm MeV −1 is the gravitational constant and R the Ricci scalar.The matter part of the Einstein-Skyrme action, S matter , describes matter in the interior of the neutron star.It is well known that the interior of a neutron star is well described as a perfect fluid of nearly free neutrons and a very degenerate gas of electrons.We exploit this and use a perfect fluid model such that the energy-momentum tensor takes the form where the energy density ρ and the pressure p are related by the multi-wall crystal EoS ρ = ρ(p).
A. The Tolman-Oppenheimer-Volkoff system Our aim is to calculate the maximum permitted mass and radius for a neutron star described by our system, and obtain the mass-radius curve.Therefore we have to solve the resulting Einstein equations for some particular choice of metric ansatz.The simplest case is that of a static non-rotating neutron star.We use a spherically symmetric ansatz of the spacetime metric, which in Schwarzschild coordinates reads [43] The mass and radius of the neutron star can be calculated by inserting this spherical metric ansatz into the Einstein equations where G µν = R µν − 1 2 Rg µν is the Einstein tensor, and solving the resulting Tolman-Oppenheimer-Volkoff (TOV) equations, 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 center of the neutron star (p(0) = p 0 ) until the condition p(R NS ) = 0 is achieved.The radial point R NS at which the pressure vanishes defines the radius of the neutron star, and the mass M is obtained from the Schwarzschild metric definition outside the star, In order for the metric function B(r) to be non-singular at r = R NS , the pressure p(r) must obey p ′ (R NS ) = 0.The TOV system (83) is solved via a central shooting method from some initial central pressure p 0 at r = 0 until the edge of the star has been reached (corresponding to p(R NS ) = 0).The amount of matter contained at r = 0 should be zero, which gives the boundary conditions B(0) = A(0) = 1.That is, the spacetime metric should approach the Minkowski metric towards the neutron star core.We can simultaneously apply a 4 th order Runge-Kutta method to the system of IVPs (83b), (83c), for the initial conditions B(0) = 1 and p(0) = p 0 , until the condition p(R NS ) = 0 is achieved.This yields the metric function B(r) and the pressure profile p(r) satisfying the necessary boundary conditions.Then the metric function A(r) can be easily obtained by numerically integrating (83a).The corresponding radius R and the stellar mass M = M (R NS ) can be extracted from the Schwarzschild definition (84).Increasing the central pressure p 0 in succession corresponds to determining a sequence of neutron stars of increasing mass, until the mass limit is reached [73].The observational mass limit is approximately 2.5M ⊙ [76], where the solar mass is M ⊙ = 1.116 × 10 60 MeV.

B. Neutron star properties and the mass-radius curve
Now we solve the TOV equations using the EoS obtained from the isospin asymmetric multi-wall crystal solution in the generalized L 0246 -Skyrme model.In Fig. 6 we present the mass-radius curve for the MC crystal (blue line) together with recent astrophysical observations.It can be seen clearly that the obtained mass-radius curve passes through many observational constraints.For our FIG.7: Plots at M max = 2.0971M ⊙ of the pressure p, energy density ρ, metric function B(r) and equations of state ρ = ρ(p).The blue curve is for the crystal EoS with the Maxwell construction applied, removing any negative pressure from the system, whereas the red curve is for the "true" crystal EoS.
choice of coupling constants (11), the Skyrme model generates an EoS which supports rather heavy neutron stars, M > 2M ⊙ .Indeed, the maximum mass is predicted to be M max = 2.0971M ⊙ , occurring for a neutron star of radius R = 13.12 km.For this solution the central energy density is ρ(0) = 784 MeV fm −3 , while the central pressure is p(0) = 155.7 MeV fm −3 .The associated plots as a function of the maximal neutron star radius is shown in Fig. 7.We find that the speed of sound in the core is approximately half of the speed of light, c s = 0.491c.
The maximal mass can be further increased if we assume higher value of the sextic term coupling constant λ, at the cost of increasing the corresponding radius.
The main improvement presented by the generalized multi-wall crystal, in comparison to previous studies involving the SC 1/2 crystal, is in the low density regime.In previous attempts, except the pure BPS Skyrme case, neutron stars obtained from Skyrme models did not have crusts, i.e. the EoS was only defined up to the nuclear saturation point n B ≥ n 0 , and not in the low density region n B < n 0 .In order to obtain a crust, the SC 1/2 crystal EoS can be smoothly joined with an EoS that well describes the low density regime, e.g. the BCPM EoS, as in [41].In the resulting hybrid EoS, the high density region is still described by the SC 1/2 crystal.This typically increases the radius of neutron star by 1-2 km, depending on the mass of the neutron star.However, such a construction is not required here as the EoS from the multi-wall crystal with the Maxwell construction is valid at both high and low densities, naturally giving the neuron star a crust.

VI. CONCLUSION
In the present paper, for the first time, we have obtained a ground state crystalline configuration for the generalized L 0246 -Skyrme model at finite densities.In contrast to previous studies on the generalized model, it has been carried out without imposing any constraints on the geometry.The only limiting assumption is the amount of the baryon charge hosted by the unit cell, which is B cell = 4.For that, we had to solve a variational problem which involves both the matter Skyrme field φ and the metric g of the unit 3-torus T 3 .
For our choice of the values for the coupling constants (11), we determine the ground state solution in the L 0246 -model to be the multi-wall crystal, as was recently observed by Harland et al. [52] in the context of the L 024model.At low densities this solution takes the form of an isolated and planar two-wall layer of skyrmionic matter.As the baryon density grows n B > n 0 then there appears to be a restoration of chiral symmetry, and the solution tends to the cubic SC 1/2 crystal.
We have used this multi-wall crystal to investigate the three most outstanding issues of the Skyrme model in its application to dense nuclear matter and neutron stars.Namely, (i) the problem of the thermodynamic instability at low densities; (ii) the maximal mass problem; and (iii) the compression modulus problem.
Firstly, in comparison with the SC 1/2 crystal or nonhomogeneous crystals (e.g.B = 32 or B = 108 crystals composed of α-particles), the use of the true ground state solution allowed to resolve the issue of thermodynamically instability at low densities.Namely, the classical energy per baryon (of the unit cell) again reveals a minimum identified with the nuclear saturation point, but now the difference between the energy at this point and at zero density is less than one percent.After inclusion of the quantum corrections to the total energy, due to the isospin d.o.f., and the lepton energy contributions for a β-equilibrated crystal, the total energy E cell of the isospin asymmetric multi-wall crystal as a function of the nuclear density n B was obtained.This minimum still existed but had reduced significantly and is practically negligible.The energy difference used in the Maxwell construction is so small that it is difficult to tell if the minimum truly exists or if it is just an artifact of our numerical algorithm.Nevertheless, it was still present so we had to use the Maxwell construction, which allowed us to obtain an EoS valid at all densities within the Skyrme model.
We remark that the Maxwell construction was required to avoid a thermodynamically unstable region which formally has negative pressure.Similar regions were found in previous studies where α, B = 32 or B = 108-crystals were studied.However, it is worth underlining that in these cases the Maxwell construction was impossible (c.f. the α-crystal) or extended to unacceptably large pressure/density regions (e.g. the corresponding neutron stars would possess cores mainly filled up by such regions).In the current work, the pressure at which the Maxwell construction is applied is only p MX = 0.022 MeV fm −3 and it extends to densities below the saturation point.Consequently, our neutron stars are mainly governed by the part of EoS above p MX , which is described by the multi-wall crystal EoS.
Of course, it is premature to identify the nonhomogeneous low density solution found here with nuclear pasta or lasagna phases in the crust of neutron stars.This is due to the fact that such phases emerge due to a balance between the nuclear and electrostatic forces.However, in our study, the Coulomb interaction has not been taken into account.In particular, we emphasize that, while our crystal qualitatively looks like nuclear pasta, it does not model nuclear pasta.Be that as it may, our result shows that the Skyrme model itself has a tendency to form complicated, geometrically nontrivial and non-homogeneous structures at low density.It should be again underlined that, on the contrary to all previous studies, we did not impose any geometry restrictions on the solutions, e.g. by assuming particular boundary conditions as in [77,78].
However, already at this stage of research, the multiwall crystal in the density regime below saturation, n B < n 0 , leads to novel and intriguing observations.The first is the symmetry energy's disclosure of its cusp structure below the nuclear saturation density, n * ∼ 3n 0 /4 < n 0 , and, secondly, the finite value of the symmetry energy in the zero density limit, n B → 0. A cusp in the symmetry energy has previously been advocated for in [71], wherein they attributed the presence of this cusp to a change in topology due to a transition between the FCC crystal of hedgehog skyrmions and the SC 1/2 crystal.A key component of their argument relies on this transition occurring in the high density regime n B > n 0 , however, this transition is believed to take place in the low density regime n B < n 0 [48].However, we have argued that these two features are generic of the Skyrme model and should occur for any infinite nuclear matter that undergoes a phase transition to somewhat isolated and finite matter in the zero density limit.This asymptotic transition to finite matter in the zero density limit is essential as the isolated solution will have a finite isospin moment of inertia tensor.A prime example of a crystalline solution in which such a transition occurs is that of the α-crystal, which tends to the isolated α-particle solution as n B → 0. Therefore, both the presence of the cusp and the non-zero value of the symmetry energy at the vacuum can be attributed as generic properties of the Skyrme model.
In fact, we have observed a further key feature of the symmetry energy.That is, a direct correspondence between the value of the symmetry energy at the vacuum and the asymmetry energy in the Bethe-Weizsäcker SEMF for nuclear binding energies.This strengthens our suggestion that the Skyrme model can be interpreted as a natural interpolation between infinite isospin asymmetric nuclear matter and finite (almost) symmetric atomic nuclei.This is further supported by the observation that the proton fraction γ p → 0.5 in the zero density limit n B → 0, which describes almost totally isospin symmetric nuclear matter, and then, for small densities, decreases yielding asymmetric matter.In this pattern one may again recognize finite nuclei.Indeed, the proton number and neutron number are approximately equivalent (δ ≈ 0) for smaller atomic nuclei while for larger nuclei there is an asymmetry (δ ̸ = 0) caused by a surplus of neutrons.
The second big issue is also resolved since the inclusion of the sextic term makes the EoS sufficiently stiff at large densities.Using this EoS we were able to compute the mass-radius curve for the resulting neutron stars.The maximal mass was found to be M max = 2.0971M ⊙ , which is a acceptable large mass and the mass-radius curve fits very well to known astrophysical data.
Finally, we shown that the problem of the compression modulus cannot be solved solely by consideration of the newly discovered non-homogeneous ground state crystalline configuration.Although reduced by approximately 200 MeV, the compression modulus is still a few times larger than the experimental value.We underline that this negative result is of high importance for the Skyrme model.It simply shows that the the solitonic model based entirely on the lightest, pionic d.o.f. is not able to correctly describe this quantity.Fundamentally, the compression modulus is related to nuclear binding energies, which is also a problem within the Skyrme model.If a variant of the Skyrme model has low binding energies then, naturally, the compression modulus will closer to its accepted value.Therefore, inclusion of more massive mesons, which are known to soften the EoS at the saturation point, seems to be unavoidable.Interestingly, this coincides with the role playing by ρ mesons in reducing of the binding energies of the Skyrmions.
It should be underlined that, if compared with other effective nuclear models, the generalized Skyrme model has an extremely small number of free parameters.It has only four coupling constants {F π , m π , e, λ}, of which the pion mass m π and the pion decay constant F π are, from the onset, fixed to their physical values, or as close to them as possible.The two other parameters e and λ, which, respectively, multiply the quartic (Skyrme) and sextic terms can be treated as free parameters in this model.They can be constrained by fitting the multiwall crystal to nuclear observables, i.e. they can be chosen such that the symmetric energy M B (n B ) and nuclear density n B at saturation n 0 are close to the experimentally determined values.
There are several directions in which our study can be continued.First of all, it is widely known that the lower density phases of nuclear matter are governed by a balance between nuclear and Coulomb forces, which leads to a plethora of geometrically different structures.The fact that the generalized Skyrme model, even without the inclusion of electrostatic interactions, gives rise to the multi-wall crystal (a lasagna like structure) can be viewed as an intrinsic ability of the model to provide such solutions.Other non-homogeneous configurations have been observed in the Skyrme model [77,78], however they were an effect of the imposed boundary conditions and therefore their applications to nuclear physics remain to be clarified.Undoubtedly, inclusion of the Coulomb interaction seems mandatory, see e.g.[79].It seems likely that including Coulomb interactions will not only give insight into such geometric phases but could also allow one to avoid use of the Maxwell construction.Thus it could possibly provide a complete description of the crust in neutron star within the Skyrme model framework.
More importantly, the inclusion of other d.o.f., like for example ρ or ω mesons, seems inevitable to resolve the issue of the compressibility at nuclear saturation.This, combined with the inhomogeneous multi-wall crystal detailed in the paper, may possibly lead to the correct value of the compression modulus.
where the isospin inertia tensor density contribution from the Skyrme field φ is given by (46).In the quaternionic formulation, the su(2) current T i is expressed by the vector quaternion The corresponding contractions are found to be Therefore, the Dirichlet contribution can be written as Then, finally, we need to consider the term Putting all of this together by using the quaternion representation (18), the isospin inertia tensor density takes the form

FIG. 1 :
FIG. 1: L 0246 -Skyrme multi-wall crystal at a fixed baryon density n B < n 0 .The isobaryon density is depicted in (a) and isosurface plots of the σ field, where the vacuum (σ = +0.9) is colored red and the anti-vacuum (σ = −0.9)blue, are shown in (b).

FIG. 2 :
FIG.2:The classical static energy per baryon M B /B as a function of the nuclear density n B .The nuclear density at which at the cusp in the symmetry energy appears is labeled by n * .This corresponds to the density at which the infinite crystalline multi-wall solution begins transitioning to an isolated multi-wall configuration.

38 FIG. 3 :
FIG. 3: The nuclear symmetry energy S N as a function of the baryon density n B , exhibiting the cusp structure detailed in the text at n * ∼ 3n 0 /4.

FIG. 4 :
FIG.4: Plot of the particle number densities n i as functions of the baryon density n B .The particle number densities are normalized such that the total number density is i n i = 1.The transition between isospin asymmetric infinite matter and symmetric finite matter at the cusp density n * is now manifest.

FIG. 5 :
FIG.5: Comparison between the classical isospin symmetric crystal, the pure neutron crystal, and the β-equilibrated asymmetric crystal with the Maxwell construction applied.

5 FIG. 6 :
FIG. 6: Mass-radius curves for neutron stars obtained from the multi-wall crystal EoS with (blue curve) and without (red curve) the Maxwell construction.The maximal mass M max obtained from the MC multi-wall crystal EoS is also shown.