Symmetry energy in holographic QCD

We study the symmetry energy (SE), an important quantity in nuclear physics, in the Witten-Sakai-Sugimoto model and in a much simpler hard-wall model of holographic QCD. The SE is the energy contribution to the nucleus due to having an unequal number of neutrons and protons. Using a homogeneous Ansatz representing smeared instantons and quantizing their isospin, we extract the SE and the proton fraction assuming charge neutrality and beta-equilibrium, using quantization of the isospin zeromode. We also show the equivalence between our method adapted from solitons and the usual way of the isospin controlled by a chemical potential at the holographic boundary. We find that the SE can be well described in the WSS model if we allow for a larger 't Hooft coupling and lower Kaluza-Klein scale than is normally used in phenomenological fits.


Introduction
The equation of state (EOS) of nuclear matter is central to nuclear physics from neutron stars to heavy ion collisions, and an important feature is the symmetry energy (SE) as a function of the density.The symmetry energy is the symmetric increase in energy as one moves away from the isospin symmetric point, that is, the point where the number of protons equals that of neutrons, i.e.E(ρ) = E 0 (ρ)+S(ρ)β 2 +• • • , with β = N −Z A being the difference between the number of neutrons N and the number of protons Z, normalized by the atomic mass number, A = Z + N -for a nice review see Ref. [1].The symmetry energy is experimentally well constrained around saturation density, ρ 0 ∼ 0.16fm −3 , to be near S(ρ 0 ) ∼ 30 MeV -both from astrophysical observations as well as heavy ion collision data -but much less so at larger densities.The symmetry energy around saturation density is conventionally expanded as with ε := (ρ −ρ 0 )/ρ 0 , whereas L and K sym are proportional to the slope and second derivative of the SE with respect to the density.Expectedly, the constraints on L and K sym are less tight than those on S 0 .Traditionally, the symmetry energy was defined for nuclear matter, which can be thought of as an infinitely large nucleus at density ρ, and so surface effects are absent.
The symmetry energy can equally well be defined for a fixed, but finite, atomic number A.
Current experimental bounds on the first 3 observables of the symmetry energy as in the expansion of the density come from mass, radius and tidal deformation of neutron stars, excitation energies of isobaric analog states, neutron skin in Sn isotopes and 208 Pb as well as heavy ion collision data [2][3][4][5][5][6][7][8][9].
The equation of state in nuclear physics relates the energy density with the pressure and is the main ingredient in the understanding of neutron stars as well as heavy ion collisions.The problem with obtaining the equation of state for nuclei is that the strong nuclear force is governed by Quantum Chromodynamics (QCD), an inherently strongly coupled theory and hence cannot be tackled by perturbation theory or first-principles calculations.Nuclear physics, in particular, ab initio methods, like the no-core shell model [10], utilize pion scattering data to reconstruct the interaction potential of nuclei and this approach leads to solid predictions for the interaction potential and the chiral effective field theory can accurately determine the EOS [11], albeit only at relatively small densities.QCD at high energies is perturbative due to its asymptotically free nature, and hence can be used to make solid predictions for the EOS [12], unfortunately at pressures far larger than those of a neutron star -the most compact object known, not collapsed into a black hole (BH).
A new paradigm of studying QCD and attempting to extract observables for nuclear physics and hadronic physics, was envisioned by Maldacena at the end of the '90-ies [13] and further elaborated by Witten [14].After a couple of decades, the mentioned framework known as holography or AdS/CFT, has been coined holographic QCD (HQCD) when applied to the strong nuclear force [15][16][17].There are two main approaches to HQCD, top-down and bottom-up; the top-down approach is based on string-theory constructions and the most prominent example is the Witten-Sakai-Sugimoto model (WSS) [14,18,19].For the bottom-up construction, which shares similar theoretical ingredients, two main types of models are known as soft-wall (SW) (e.g.Improved HQCD [20,21] and V-QCD [22][23][24]) and hard-wall (HW) models [25][26][27][28][29][30][31][32][33].Especially, the top-down type of HQCD have quite some predictive power, in the sense that the models have very few adjustable parameters [15][16][17].For the WSS, there is the mass scale and the 't Hooft coupling, where the mass scale is normally fitted to the mass of the ρ meson and the 't Hooft coupling is determined from the pion decay constant [18].
Attempts have already been made at extracting the SE from various HQCD models, including top-down approaches as in the D4/D6 model [34] and even the WSS model, see Ref. [35], and NSs have been constructed using the holographically extracted EOS to solve the governing Tolman-Oppenheimer-Volkov (TOV) equations [36,37].The SE, however, comes out too large in the WSS [35].In HQCD in contrast to traditional nuclear physics, the proton and the neutron are not point particles, but are described by a topological soliton, the Sakai-Sugimoto soliton [18,[38][39][40], which is initially isospin symmetric -that is, the proton is equal to the neutron.In order to compute the SE, we must distinguish the proton from the neutron and this can be done by the introduction of an isospin chemical potential [35].
HQCD at finite isospin chemical potential has been object of inspection in the context of many models: top-down approaches include the D3/D7 model [41][42][43] and the WSS model [44,45] in the case of the DBI action and without employing the homogeneous Ansatz (see below).
In this paper, we propose using the homogeneous Ansatz in the WSS, but quantizing the isospin symmetry -a technique known from the Skyrme model [46][47][48], which is the leadingorder low-energy effective theory of the WSS model [18].The homogeneous Ansatz represents an approximation to describe densely packed nucleons that form nuclear matter above saturation density.It relies on the assumption that nuclear matter forms a spatially homogeneous distribution, in which nucleons lose their individual properties: despite being shown in Ref. [49] that such a configuration is not admitted in holographic models under assumptions of regularity of the gauge fields, the Ansatz can still be employed with modifications, such as formulating it at the level of the field strengths [50] or (as we will do) introducing a discontinuity that acts as a source of baryon number [51].The quantization of the isospin symmetry introduces the isospin quantum number, which makes it possible to extract the SE as the coefficient of the square of the difference between the number of protons and neutrons.The quantization method of introducing isospin is also shown to be equivalent to using a chemical potential, see Appendix A. We find a lower SE compared to previous attempts in the WSS [35], since we include all the needed fields in our Ansatz and because we choose a different N c scaling such that the nucleons are states with the minimal isospin quantum number; however, there is no difference coming from using either the chemical potential or the quantization method -they are equivalent as shown explicitly in Appendix A. In particular, we find a phenomenologically viable value of the constant S(ρ 0 ) at saturation density and the first two coefficients L and K sym are compatible with current experimental bounds from astrophysics and heavy ion collision data for a certain choice of the model parameters.

Model
We will treat the WSS and the HW model on equal footing in the following.The model at low energies is described by the Yang-Mills (YM) and Chern-Simons (CS) actions in 5-dimensional AdS 5 (or AdS 5 -like, for the WSS model) spacetime (M , g): with S = S YM + S CS being the total action, the constant κ = α, β = 0, 1, 2, 3, 4 with x 4 = z the holographic coordinate, and the power of forms is understood with the wedge product.The metric is for the HW model, the index µ = 0, 1, 2, 3 is summed over in the metric and µ is raised with the Minkowski metric here.It is convenient to work in dimensionless units.Both models have a free mass scale M KK (WSS) and L −1 (HW): the WSS model is already presented with the choice M KK = 1, to which we add the choice L = 1 for the HW model.The correct powers of the energy scales M KK , L −1 can then easily be restored via dimensional analysis.
The gauge field can be decomposed for later convenience in the Abelian and non-Abelian parts as where the generators of SU(2) T a are chosen as T a = 1 2 τ a so that Tr T a T b = 1 2 δ ab , and the spacetime indices follow the convention: In writing Eq. ( 2), we performed dimensional reduction in the WSS, integrating out S 4 from the original nine-dimensional action for the stack of D8−Branes, while in the HW we do not explicitly include an action for the scalar field encoding chiral symmetry breaking, since we set the scalar field to zero, which is appropriate in the homogeneous baryonic phase, following Ref.[37].Despite not appearing explicitly in our computation, the scalar field plays an important role: its vacuum energy, determines the density of nuclear matter at the baryonic onset, hence defining saturation density within this model.For details on how the scalar field defines the saturation density, but otherwise vanishes in the baryonic phase, see Appendix D.
Here we will utilize the fit found in Ref. [37] and only adjust the overall energy scale.Two further steps were employed in order to write the action and equations of motion for the two models in a compact way.For the HW model we assumed the symmetry properties for the fundamental fields L M , R M as follows: For the WSS model, we similarly assumed parity properties of the fields with respect to z: With these procedures, we halve the number of fields in the HW model (from L, R to A) and the integration interval in the WSS model (from (−∞, +∞) to [0, +∞)) generating an overall factor of 2 in the action in both cases.As a last step, we introduce generic symbols z IR , z UV to indicate the infrared and ultraviolet boundary values of the holographic coordinate, 1 which in the two models assume the values The classical homogeneous Ansatz for isospin-symmetric matter, reasonable for large-density computations, is defined as where { a 0 , H} = { a 0 , H}(z) are functions of the holographic coordinate z.We have suppressed the unit 2-by-2 matrices in the terms without a Pauli matrix τ.
The function H(z) encodes the baryonic density through its value at z = z IR : if either H(z IR ) or H ′ (z IR ) vanish, then the baryon number would also vanish as noted in Ref. [49], so we will assume that H(z) obeys a Dirichlet boundary condition H(z IR ) = H 0 , with the value of H 0 to be determined by minimization of the action.This defines the baryon density ρ (assuming H(z) → 0 for z → z UV ) as follows: so that the infrared boundary condition for the numerical integration of the function H(z) is directly related to the baryon number density as: where for convenience of putting the two models on same footing, we have defined the integral in the holographic direction as which we will use throughout the paper and ε assumes a different sign depending on the model: Thus the integral is defined in such a way to take into account the different orientation in the integration along z, dictated by the choice of coordinates for the two models.Note that this choice of boundary condition for H(z IR ) means that in the WSS, once we restore the full domain of integration z ∈ (−∞, ∞), the function H(z) will be discontinuous.This still leads to a continuous field strength, since both H ′ and H2 are continuous functions.For the HW model instead, this choice just means that we cannot enforce the standard boundary condition L µ (z IR ) − R µ (z IR ) = 0, which has to be replaced with the one above, implying

Time-dependent configurations
We wish to include the effects of isospin asymmetry in the system.To do so, we follow a method inspired by the single-soliton analysis: we know that for the single baryon, the proton and the neutron are described as degenerate (in absence of quark mass terms 2 ) quantum states of the effective Hamiltonian obtained by considering a slow rotation in SU (2).The homogeneous Ansatz shares a similar structure with the single-soliton configuration, made easier by the absence of translational moduli3 X M and ρ (but with the minor complication of not having an analytical configuration to approximate our static Ansatz ( 10)), so we can attempt to follow steps similar to the ones in Refs.[38] and [39], in order to obtain a timedependent configuration -yet to be quantized.We start by assuming a configuration of the form which implies the following transformations in the field strengths: where V (z, t) encodes the time-dependent rotation in SU(2), and Φ is defined as Notice, this is not a gauge transformation since the field A 0 is not transformed along with the rest.The function V (z, t) needs to depend on z in order to allow us to satisfy the equation of motion The function V (z, t) is holographically dual to the SU(2)-valued collective coordinate a(t), as we choose it such that which in turn implies where χ is the boundary angular velocity.The presence of a nonvanishing F 0z will also enable a source term for the fields A i via the Chern-Simons action, so we will have to complete the field content by turning on A i = − 1 2 Lχ i : here we already guessed that the vector field will be proportional to the angular velocity χ i , and we can do so without loss of generality, since in the homogeneous case this is the only three-vector available to the Abelian field.
At this stage the problem is well posed and the function Φ(z, t) can be found by solving Eq. ( 22), but it is more convenient to perform a gauge transformation to make the system easier to treat.
We perform the gauge SU(2) transformation where the superscript "S" stands for "singular", because this is reminiscent of the transformation changing from the regular gauge to the singular gauge in the single-soliton case.With this choice (dropping the superscript "S" for convenience, since we will use this gauge henceforth) the field content becomes Now we can factorize the function Φ(z, t) as and since we imposed Eq. ( 24), we see that We then conclude that the field A 0 in this gauge vanishes at the UV boundary, and can be expressed as We notice that this result is exactly what one would expect by allowing for the most general field configuration respecting spherical symmetry, homogeneity in three-dimensional flat space, and the gauge choice A z = 0. Taking the functions H, a 0 , G, L to be independent of χ amounts to considering a slow rotation, thus including only linear terms in χ in the Ansatz.Whereas 1 2 χ • τ is the matrix form of the boundary angular velocity, 1  2 aχ • τa −1 = −iȧa −1 is the matrix form of the boundary angular isospin velocity (i.e.describing rotations in SU(2) instead of in space).Thus, although one may think we are spinning the fields in space, this is really an isospin action on the homogeneous fields.
The final form of our time-dependent homogeneous Ansatz is then summarized in compact notation as: with the mandatory boundary condition G(z → z UV ) = 0.This Ansatz leads to the action which gives rise to the equations of motion where the first two equations of motion are truncated to order |χ| 0 , whereas the latter two only appear at quadratic order in χ.Including the subleading |χ| 2 corrections to the solutions of H and a 0 has a negligible impact, which we checked explicitly.Moreover, in the limit of small χ, quadratic corrections in χ to the functions H(z), a 0 (z) would generate terms in the on-shell action at order |χ| 4 , not contributing to the symmetry energy, and also being subleading in the small χ expansion.This set of equations is composed by ODEs in the holographic coordinate z and can be solved with standard off-the-peg solvers in packages like MATHEMATICA or MATLAB, once we specify all the boundary conditions: and all fields are vanishing at z = z UV .These boundary conditions are obtained by imposing the vanishing of also the total derivative that comes about when deriving the Euler-Lagrange field equations; for more details, see Ref. [58]. 4 We recall that that in the chosen coordinates for the WSS (HW) model.
Another possible approach would be to keep the fields in a static configuration, hence keeping the freedom to set the standard orientation of Eq. ( 10), and introduce an external isospin chemical potential, which holographically amounts to introducing a finite UV boundary value for the field A 0 : in Appendix A we show that this approach is related to ours by a gauge transformation, hence leading to the same physics.This formalism is the one employed in [35,55]: the two calculations, however, differ in that in the present work we have turned on the Abelian field A i , which turns out to be linear in χ, and we are effectively truncating the χ dependence of the gauge fields at linear order.The inclusion of the Abelian component is necessary to have a self-consistent Ansatz, as the equations of motion cannot be solved by setting L(z) = 0 (the Chern-Simons term provides a source for L(z)).Moreover, it turns out that the field L(z) dominates the small-λ behavior of the symmetry energy: despite the holographic model being developed with the large-λ limit in mind, for the practical application of extracting a value for the symmetry energy, we need to extrapolate to a finite-λ, and the most popular fit of the model employs the value of λ = 16.63, which does not realize the large-λ nor the small-λ regimes (see Appendix E).On top of the difference at the level of the Ansatz, another difference with respect to Refs.[35,55] lies in the implicit definition of the isospin number of nucleon states in the large-N c limit.We choose as proton (neutron) state the lowest-lying isospin state, which can be thought of as being composed of 1  2 (N c + 1) up (down) and 1  2 (N c − 1) down (up) quarks.With this definition, the angular velocity χ of a nucleon state is of order N −1 c , and so are the isospin chemical potential and the symmetry energy. 5A different choice that still reduces to the familiar N c = 3 case is that in which the proton (neutron) is composed of N c −1 up (down) and one down (up) quarks: in this scenario the isospin number is of order N c , and so is the symmetry energy.We find appropriate the former definition for nucleon states, in that it keeps the nucleons as the ground state baryons, and preserves the familiar electric charge following the Gell-Mann-Nishijima formula with Q, I 3 , N B being the electric charge, the third component of the isospin, and the baryon number, respectively.
The truncation of the χ dependence (and so of the dependence on µ I ) is an approximation that does not affect the computation, since the symmetry energy is by definition obtained by evaluating the first nonvanishing term in the expansion of the energy per nucleon around an isospin symmetric configuration.
Despite this technique being equivalent to the usual introduction of a boundary chemical potential µ I , it has a series of advantages, particularly manifest in the small µ I limit.In this limit, the complicated picture of the isospin asymmetric homogeneous matter becomes similar to the well understood one of a slowly rotating bulk instanton, and the problem of finding the symmetry energy becomes the computation of a moment of inertia.
It is also built-in in this formalism what the smallest value of the isospin is.Identifying this smallest unit of isospin with that of a single quark being flipped from down to up fixes the isospin quantum number without N c -scaling ambiguities.
On top of this simplification, our alternative technique also helps in identifying the solution to the problem of the ambiguity of the Chern-Simons term when dealing with homogeneous nuclear matter.As pointed out in Ref. [58], boundary terms arising from the Chern-Simons term can contribute with an IR effective action because of the discontinuity of H(z).In this case, different choices for the Chern-Simons action that differ by a boundary term are not physically equivalent (as opposed to the case of a smooth instanton), as they enforce different IR boundary conditions on the flavor fields, with the consistent case giving rise to the boundary conditions (39).
A way of solving the ambiguity is to require that the holographic currents on the boundary match with their sources in the bulk: in Ref. [58] this was done for the baryonic density ρ (requiring that the topological charge matches with the baryonic current in the tail of a 0 ) and for the isospin density ρ I (requiring that the angular velocity description matches with the isospin chemical potential one).
This kind of problem arises in every holographic model containing a Chern-Simons term in the action, hence this result obtained with the aid of the new quantization technique is very generalizable and will prove useful in identifying the correct Chern-Simons action for future works exploring isospin asymmetry in holographic models.The correct choice of the action is crucial to obtain the correct thermodynamic quantities, so an improvement in this regard directly translates to a more precise equation of state, and more reliable predictions for properties of neutron stars.

Symmetry energy
The terms quadratic in χ exactly produce the SE upon Hamiltonian quantization: where canonical quantization of a m , m = 0, 1, 2, 3, a coordinate on the 3-sphere (a 2 m = 1), leads to the momentum conjugate and hence to π 2 m = ℓ(ℓ + 2) being the spherical harmonics and ℓ = 2I, with I the isospin quantum number [46]. 6The identification of V χ 2 and I(I + 1)/V coming from Hamiltonian quantization is also justified by the holographic dictionary, since it can be obtained by computing the third component of the isovector charge density, see Appendix B for a detailed computation.The functionals Λ and U are defined as where V denotes the spatial 3-volume.Using now the relation between isospin and the number of protons and neutrons: with Z the proton number and N the neutron number, as well as the atomic number being the product of the 3-volume and the baryonic density.β is defined as the normalized difference between the number of neutrons and protons, β = (N − Z)/A, hence we have where S(ρ) is the symmetry energy as a function of the density.
Using the standard phenomenological fit for the WSS model of Ref. [19], we set λ = 16.63 and find the first SE expansion parameters as S 0 = 74.9 M KK 949 MeV MeV , MeV , which are somewhat larger than values typically obtained from phenomenological models [1], but much smaller than obtained in the WSS previously [35].For the HW model, we fix 12π 2 using the leading OPE coefficient of the vector current correlator [30], for which the first few SE expansion parameters are where L −1 is the mass scale of the HW model, which we set as L −1 = 150 MeV following Ref. [37], which provides phenomenologically good results for neutron stars in terms of massradius data.The saturation density, ρ 0 , in HQCD is defined to be at the onset of baryonic matter, obtained by minimization of the free energy with the baryonic chemical potential as the boundary condition a 0 (z UV ) = µ B .The value of ρ 0 obtained this way in the WSS model with λ = 16.63 is which is about 2.9 times too large with respect to the phenomenological value -an overestimate by the same order of magnitude as the other baryonic quantities.The HW model, however, yields a more realistic saturation density which is only about 22% too large.We explore a larger range of densities for both the WSS and the HW model in Fig. 1a).For the WSS model, we have used the phenomenological value of the 't Hooft coupling (λ = 16.63) and shown the range of M KK ∈ [300, 1200] MeV with a red shaded area, which includes M KK = 949 MeV [19] (red solid line), whereas for the HW model the range of L −1 ∈ [110, 320] MeV is shown with a green shaded area, which includes L −1 = 150 MeV (green solid line) that is chosen from neutron star phenomenology [37] and the highest mass scale is from meson physics [31,59].Up-to-date constraints from astrophysics and heavy-ion collision data are shown with gray and cyan shaded areas near and above saturation density and constraints using nuclear excitation energies from isobaric analog states (IAS) are shown with a purple shaded area below saturation density.As can be seen from the figure, the phenomenologically fitted value of the mass scale M KK at 949 MeV [38] overestimates the SE with about a factor of 2.4; however, the fit is made using mesonic observables and is known to overestimate baryonic observables; for instance, the baryon mass is typically overestimated by a factor of 1.7-1.8[38,60,61] using the mesonic fit.
Although both models come in the ball park of the experimental constraints above saturation density and nuclear physics predictions below saturation density if we allow ourselves to adjust the energy scale (M KK or L −1 ), the shape of the SE does not quite satisfy all the constraints.In the WSS model, however, we can dial the 't Hooft coupling to see whether we can fit in the allowed regions and indeed it is possible by raising both the 't Hooft coupling from the phenomenological value to λ = 60, as well as lowering the KK scale from 949 MeV to 390 MeV, see Fig. 1b); this corresponds to the lower black curve of the red shaded area.With the larger 't Hooft coupling, the SE of the WSS has a compatible shape to pass all the constraints, but the ρ meson is too light and the baryon mass is too heavy -often a problem in HQCD; the baryon mass is reduced from about 1600 MeV to 1191 MeV.If we keep the pion decay constant at its phenomenological value, the KK scale, however, is lowered too and is shown in Fig. 1b) with a solid red curve -not too far from a viable solution.Recomputing the symmetry energy expansion parameters, we obtain MeV , which are compatible with phenomenological constraints.The saturation density for this fit now also has improved as which is only about 10% from the phenomenological value.Lowering the KK scale from 390 to 380 MeV will improve both S 0 and ρ 0 (S 0 = 32.2MeV and ρ 0 = 0.153 fm −3 ), but will create a bit more tension with the constraint coming from the neutron skin thickness of 208 Pb, shown with a gray-shaded area in Fig. 1b).A similar choice of fit is employed in Ref. [62], where we choose M KK and λ in order to reproduce the correct saturation density and SE.The resulting equation of state (hybridized with phenomenological low-density EOSs) then is found to reproduce viable properties of neutron stars.It is not surprising that a somewhat larger value of λ is needed in order to more reliably reproduce the physics of baryonic matter.We can understand it by considering the single baryon in the WSS model: with the BPST instanton approximation it is possible to compute both the classical mass and its quantum corrections [38].In particular, since we are interested in the symmetry energy, we want to consider the (iso)spin quantum correction, given by integer values of l in the formula [38] where n ρ , n Z are quantum numbers for the size and bulk position excitations.
It is well known that the usual mesonic fit with M KK = 949 MeV, λ = 16.63 largely overestimates the nucleon masses, but a major contribution in this result comes from the fact that the quantum corrections are close in magnitude to the classical mass.We can see that the classical The gray, cyan and purple shaded areas are the same as in Fig. 1.The vertical gray line marks the phenomenological 't Hooft coupling λ = 16.63 [19].
mass is of order O(λ), while the quantum corrections are of order O(1).By increasing λ it is then possible to reduce the relative magnitude of the quantum corrections as compared to the classical result, and by reducing M KK it is possible to obtain an overall more realistic nucleon mass.
The same kind of mechanism is inherited by the homogeneous system, where increasing λ, the symmetry energy contribution becomes smaller, moving towards the real world value, provided that we also adjust M KK accordingly.While it is fairly easy to expect that some values of λ, M KK exist that both fit the phenomenology of saturation density and symmetry energy, it is less trivial that once these values are employed, other baryonic observables are then improved, as it happens in our case with the BPST instanton mass and with the expansion parameters L, K sym .
What this analysis suggests is that when describing baryons within the WSS model, the errors introduced by the many approximations employed to make the system approachable seems to be partially mitigated by an alternative choice of the values of the parameters λ, M KK .
The dependence on the 't Hooft coupling for the WSS model is shown in Fig. 2 for the KK scale in the interval 300-1200MeV at saturation density.

Proton fraction
We will now consider the proton fraction at β-equilibrium with charged leptons, imposing charge neutrality.Using the Gell-Mann-Nishijima formula, we can relate the baryon density, ρ, and isospin density, ρ I , with the proton/neutron densities: where the upper sign is for protons and the lower for neutrons.Charge neutrality is imposed by with ℓ = e, µ being a sum over the charged leptons and the β-equilibrium (from the decay MeV, the solid magenta curve corresponds to the phenomenological pion decay constant, and the dashed magenta curve corresponds to M KK → 0, thus eliminating the muons.The gray shaded area is the result chiral EFT [63]. where µ X is the chemical potential of the particle species X .The lepton density is calculated assuming it to be a (massive) Fermi gas as [35] with Θ H being the Heaviside step function, µ ℓ the chemical potential and m ℓ being the mass of the lepton ℓ.Using the definition of the isospin chemical potential as being the conjugate variable of the isospin density, we get Inserting Eq. ( 58) into the charge neutrality condition (56) and using the β-equilibrium condition (57), we obtain an implicit solution for the isospin density, ρ I , as a function of the density ρ: where we have set the electron mass to zero and the dimensionless muon mass parameter is the ratio of the physical mass (105.7 MeV) to the mass scale M KK and L −1 , for the WSS and the HW models, respectively.In Fig. 3 we show the numerical results for both the WSS and the HW model for the proton fraction at various densities around saturation density.We find that the WSS model phenomenologically fitted to mesons gives more realistic proton fractions (red shaded area) than the HW model (green shaded area) and yields even better proton fractions below saturation density if we use the calibration from Fig. 1b), i.e. λ = 60 and M KK = 390 MeV.Since we take the electrons to be massless, the mass scale of the model only enters in the muon mass parameter.The dashed curves correspond to the muon being infinitely heavy (or the mass scale of the model being sent to zero).

Discussion and outlook
In this paper, we have computed the symmetry energy in two holographic QCD models using the method of quantizing the isospin symmetry, namely in the top-down WSS model and the bottom-up HW model.We find fairly good agreement between our model results for the SE and proton fraction in both the HW model, using the fit from neutron stars and in the WSS model with a new calibration (i.e.λ ∼ 60 and M KK ∼ 390 MeV).
We have also shown that the method known from Skyrmions of quantizing the isospin zeromode is equivalent to introducing a chemical potential on the holographic boundary for the gauge fields, see Appendix A. There is mathematically no difference between the two methods.
It would be interesting in future work to take into account the strange quark (3 instead of 2 flavors) or alternatively the kaons, to see at what densities it might have an impact on the SE.Furthermore, there are certain transitions that happen at larger densities, for example the Skyrmion-half-Skyrmion transition [64], which has an analog in holographic instantons [65].Although it is not directly observable in our homogeneous Ansatz, it may have some effect on the SE and proton fractions at large densities.Another approximation we employed in our calculations is that of keeping fixed the position of the discontinuity that sources the baryonic charge.It is expected that the location of the baryon is dynamically determined, moving towards the boundary as the density is increased: this behavior is the (homogeneous version of) the so-called "popcorn transition" that is known to occur in HQCD [66].However, as shown in Ref. [37], the homogeneous Ansatz already reproduces features of the popcorn transition even when keeping the position of the discontinuity fixed in the bulk: it does so by having the baryonic holographic density form a peak at a location in the bulk that depends on the boundary density ρ.Because of this, we expect corrections coming from the dynamical determination of the discontinuity location to be particularly small.density (and µ I ) anyway.For calculations at higher isospin chemical potential, the equivalence still holds, but since the now large angular velocity χ backreacts every field, this framework loses its main advantage of factorizing away the dependence on χ.Moreover, the isospin symmetric Ansatz is not consistent anymore, and the function H(z) has to be substituted with a set of functions H i (z) [55]: where i is not summed over.Then the problem to be solved is that of five (H 1 (z) = H 2 (z) because of residual symmetry) coupled ODEs, with χ dependence in each of them, effectively the same as we would have if we worked with a boundary chemical potential µ I .Despite this not being a necessity for the calculation of the symmetry energy (since by definition it is calculated at vanishing µ I ), when moving to finite isospin density (for example when computing the full phase or properties of neutron stars) the more rigorous approach would be to include the effects above (and possibly non-diagonal terms for A a i when including also quark masses, see Ref. [55] for the validity limits of the diagonal approximation), which is equally difficult with both gauge choices.

B Isospin density from holographic current
In this section we want to prove that the isospin number density that we defined from the quantized angular momentum coincides with the canonical one obtained through the holographic dictionary via the computation of the isovectorial current.For simplicity, we will show the proof in the WSS model, so that z IR = 0, z UV = +∞, but it holds true in the HW model too, after substitution of the appropriate quantities.As shown in Ref. [39], the vectorial current is obtained as With this quantity we can build the isovectorial charge Q V of which we take the third component to coincide with the isospin operator We plug in this formula the homogeneous Ansatz (32): where we used the fact that G ′ (0) = 0.The angular velocity χ i is related to the angular momentum operator J i by the familiar relation involving the moment of inertia Λ: and we can exploit the useful relationship between angular momentum and isospin operators that holds due to the spherical symmetry of the system: so that we are left with: We see that the validity of this relationship depends on whether the following identification holds true: 4κ To prove this relationship, we first notice that the formula for the current is obtained by differentiating the action with respect to the UV boundary value of the A 0 field, following δ A 0 S = (e.o.m. terms) + 4κV Tr k(z)A ′ 0 δA 0 z=+∞ , (B.9) where the first term means that we neglect contributions that vanish by the equations of motion, we evaluate only boundary terms, and we made use of the boundary condition A ′ 0 (z = 0) = 0. We decide to employ the gauge in which the field A 0 has a finite boundary value, dual to the isospin chemical potential, so we take the field to be as in Eq. (A. 14).With this configuration, the variation of the action assumes the shape We can change this result to our usual "rotating" gauge by noting that we have to rename µ I → χ 3 , and that G ′ = G ′ , so that on-shell we obtain We now look at the definition of Λ: it is nothing but the part of the energy density that is quadratic in the angular velocity, and there is no linear term.In this picture, the system is rotating and there is no chemical potential, so the on-shell action gives the energy of the system, so that we can write

C The subleading order in N c of the chiral anomaly
Throughout the main body of this work, we have ignored the presence of the chiral anomaly of QCD: we expect on general grounds that the currents dual to the holographic fields A α are conserved, with the exception of the axial U(1) current, since the corresponding symmetry is broken by the chiral anomaly.Analogously, we expect the Goldstone boson associated to the axial symmetry to acquire a finite mass as a consequence of the anomaly.This holds true in the WSS model, where the mechanism is incorporated nontrivially from the top-down construction in string theory.The model includes Ramond-Ramond forms C n of odd rank n: Among these is C 7 , whose action, inclusive of a coupling with the flavor branes, reads where a one-form ω y = δ( y)d y has been introduced to model the distribution of the stack of branes in the y-direction (by definition transverse to z), extending the otherwise 9-dimensional integral to the whole 10-dimensional spacetime.This form is gauge invariant if we allow C 1 to transform with a U(1) transformation of the flavor group: The implication of this fact is that dC 1 is not a gauge invariant form, only F 2 is the correct gauge invariant combination.This is welcome, since in the model C 1 is dual to the θ angle of QCD as Let us now consider a zeromode for the field A z , dual to the η ′ meson, such that and plug it into the action The result is an action that displays a mass term for the η ′ : with the η ′ mass agreeing with the Witten-Veneziano formula The topological susceptibility χ g and the pion decay constant are computed in the model (see Ref. [18]): We now recall that the parameter κ was defined to be κ ≡ λN c 216π 3 : this means that the mass term for the η ′ meson is of order N −1 c (since f 2 π ∝ O(N c )), while the action for the flavor fields that we employed in the main body of this work is of order N c .While true that the angular velocity χ i itself is of order N −1 c , hence pushing the isospin asymmetric action (proportional to χ 2 ) to be of order N −1 c , we have to recall that eventual isospin asymmetric terms will carry factors of N −1 c or higher also in the η ′ mass term, pushing it to even higher order in the N −1 c expansion.Hence it is formally safe to neglect the contribution from the axial anomaly in the large-N c scheme of approximation.

D The vanishing of the hard-wall tachyon in the baryonic phase
Let us quickly review the setup of the hard-wall model of Ref. [37] that utilizes a scalar field with an IR potential to dynamically stabilize it.The metric is given by where L := 1 is the curvature scale of AdS 5 and set equal to one.Upon restoring units, energies are multiplied by a physical scale.
For two flavors, we have left and right U(2) gauge fields, L M , R M and the minimal action [37]: where a(z) = L/z, the U(2) gauge field L M is split into SU(2) and U(1) parts as and similarly for R M , the field strength for L M is and similarly for R M N , the covariant derivative for the scalar field is defined as the boundary condition for the scalar field is which is stabilized by the boundary potential as the minimization of the vacuum solution and is given by the mass of the scalar is M 2 Φ L 2 = −3, the would-be quark mass in the model is switched off, the indices M , N = 0, 1, 2, 3, z run over all AdS 5 , and finally N c is the number of colors and M 5 is a coupling of the theory (playing the role of κ in the WSS model), which we have set as Chiral symmetry breaking is done in Ref. [37] following [68,69] as (L zµ + R zµ ) z=z IR = 0, and hence we choose L z = R z = 0 (gauge choice), L i = −R i , L 0 = R 0 and Φ diagonal, which means that the scalar field only couples to L i = −R i via the covariant derivative  with µ being the (baryonic) chemical potential.
In the phase ρ > 0, the vacuum of the theory is still given by H = a 0 = 0 and Φ given by the vacuum solution (D.14) until the baryonic onset, which corresponds to the nuclear saturation density.At the onset, there are two vacua: a mesonic and baryonic one each with the same value of the action (by definition), see Fig. 4(a).Once ρ > ρ crit one may ask at what configuration the scalar field stabilizes at.It turns out by numerical computations that the scalar turns off, which corresponds to ξ = 0 in the baryonic phase, see Fig. 4(b).This corresponds to Φ = 0 and when studying only the baryonic phase, the impact of the scalar is to set the saturation density of this simplistic hard-wall model.

E Comparison with large-λ approximation
In the large-λ limit, one may consider ignoring the inclusion of the Abelian field L (in Eq. ( 32)), which however is not a consistent choice for finite values of λ, as the equation of motion for L is not satisfied by L = 0 (it is sourced by H and G, see Eq. ( 38)).This approximation can be seen as consistent in the large-λ limit by considering that the field L is only sourced by the Chern-Simons, which is indeed subleading in λ with respect to the Yang-Mills terms.One may be led to believe that the same is true for the field G, and that they have to appear at the same order in λ, but a closer look at its equation of motion (37) shows that a source is present already in the Yang-Mills terms, coming from the time derivative of the field A i (in our time-dependent gauge, while the same contribution arises from the UV boundary condition of G in the static gauge).
Neglecting the Abelian field L was one of the approximations made in Ref. [35] in addition to choosing a different N c -scaling of the isospin chemical potential (essentially defining the large-N c baryon's proton and neutron by maximally flipping the down and up quarks, as opposed to our definition where only one quark is flipped from down to up).In order to see quantitatively how good the approximation of using L = 0 in the homogeneous Ansatz is, we perform the numerical calculation corresponding to Fig. 2 with and without L taken into account, see Fig. 5. From the figure, we can see that the large-λ approximation works in the sense that the two results asymptote to the same curve for λ ≳ 200.At λ = 16.63 and for M KK = 949 MeV the correct computation including L yields a symmetry energy of 74.9 MeV compared to 109.1 MeV when neglecting L in the Ansatz, which at the meson-fitted value of λ gives an excess in the symmetry energy of 46%.The smaller λ is the worse it gets, as expected from a large-λ approximation.

λN c 216π 3
for the WSS model and κ = M 5 for the HW model, λ = g 2 YM N c the 't Hooft coupling, N c the number of colors of QCD (i.e. 3 in nature), the field strength 2-form is

Figure 1 :
Figure 1: The symmetry energy (SE) calculated in the WSS model with the phenomenological value of the 't Hooft coupling and in the HW model, both using quantization of isospin as functions of the density.a) The red area corresponds to the WSS model with M KK ranging from 300 MeV to 1200 MeV and the red line in the middle is at 949 MeV.The green area corresponds to the HW model with L −1 ranging from 110 MeV to 320 MeV and the green line in the middle is at 150 MeV.The constraints from the PREX-II experiment using the neutron skin thickness of 208 Pb [6] are shown with a gray shaded area, while the extensive 2021 survey of constraints on the symmetry energy of Li et.al.[4] using neutron stars, are shown with a cyan shaded area.Constraints from isobaric analog states below saturation density are shown with a purple shaded area [2].b) The SE calculated in the WSS model with the 't Hooft coupling λ = 60.The red shaded area corresponds to M KK ∈ [390, 949] MeV and the red solid curve is the rescaled phenomenological mass scale, that keeps the pion decay constant at 93 MeV, corresponding to M KK = 500 MeV.

Figure 2 :
Figure 2: The symmetry energy (SE) calculated in the WSS model using quantization of isospin as a function of the 't Hooft coupling λ at saturation density.The red area spans the mass scale M KK from 300 MeV to 1200 MeV and the red line is at 949 MeV.The gray, cyan and purple shaded areas are the same as in Fig.1.The vertical gray line marks the phenomenological 't Hooft coupling λ = 16.63[19].

Figure 3 :
Figure 3: The proton fraction calculated in the WSS and HW model as functions of the density.The red shaded area corresponds to the WSS model with the phenomenological 't Hooft coupling λ = 16.63 and M KK ∈ [300, 1200] MeV, the red line at M KK = 949 MeV and the red dashed line at M KK → 0. The green shaded area corresponds to the HW model with L −1 ∈ [110, 320] MeV, the green line atL −1 = 150 MeV and the green dashed line is the limit L −1 → 0. The magenta shaded area corresponds to the WSS model with the calibration of Fig.1b), i.e. λ = 60 and M KK ∈ [390, 949] MeV, the solid magenta curve corresponds to the phenomenological pion decay constant, and the dashed magenta curve corresponds to M KK → 0, thus eliminating the muons.The gray shaded area is the result chiral EFT[63].

Figure 4 : 2 ,
Figure 4: The action evaluated as a surface for densities ρ and scalar field with coefficient ξ at (a) (for (b) twice) the chemical potential that corresponds to saturation density.In this figure, the IR potential is chosen as m b = 0.657 and λ b = 0.001, giving λ b ξ 4 0 = 1.024.

Figure 5 :
Figure 5: The symmetry energy calculated as function of λ at saturation density in the WSS model with L taken into account (green solid line) and without L (i.e.L = 0 in the Ansatz) (red solid line).In this figure, we have used M KK = 949 MeV.