SU(3) parity doubling in cold neutron star matter

We present a phenomenological model to investigate the chiral phase transition characterized by parity doubling in dense, beta equilibrated, cold matter. Our model incorporates effective interactions constrained by SU(3) relations and considers baryonic degrees of freedom. By constraining the model with astrophysical data and nuclear matter properties, we find a first-order phase transition within realistic values of the slope parameter L. The inclusion of the baryon octet and negative parity partners, along with a chiral-invariant mass $m_{0}$, allows for a non-massless chiral symmetric phase. Through exploration of parameter space, we identify parameter sets satisfying mass and radius constraints without requiring a partonic phase. The appearance of the parity partner of the nucleon, the N(1535) resonance, suppresses strangeness, pushing hyperonization to higher densities. We observe a mild first-order phase transition to the chirally restored phase, governed by $m_{0}$. Our calculations of surface tension highlight its strong dependence on $m_{0}$. The existence of mixed phases is ruled out since they become energetically too costly. We compare stars with metastable and stable cores using both branches of the equation of state. Despite limited lifespans due to low surface tension values, phase conversion and star contraction could impact neutron stars with masses around 1.3 solar masses or more. We discuss some applications of this model in its non-zero temperatures generalization and scenarios beyond beta equilibrium that can provide insights into core-collapse supernovae, proto-neutron star evolution, and neutron star mergers. Core-collapse supernovae dynamics, influenced by chiral symmetry restoration and exotic hadronic states, affect explosion mechanisms and nucleosynthesis.


I. INTRODUCTION
First principle calculations on the lattice [1,2] support the idea that hadronic matter undergoes a transition to a chirally restored phase known as the quark-gluon plasma at high temperatures and low baryon densities.Recently, the existence of yet another type of chirally symmetric phase dominated by chromoelectric interactions has been suggested [3].Despite these expectations, the precise nature and location of this transition in the phase diagram remain uncertain.When dealing with finite chemical potentials, Monte Carlo techniques on the lattice are hindered by the notorious sign problem [4,5].Consequently, adopting effective models becomes the practical approach to capture the microphysics of the relevant degrees of freedom.Perturbative QCD at finite densities is also an option but can only be safely applied for extremely large densities, surpassing the typical central densities of compact stars [6][7][8][9], which are the densest known physical systems in the universe.
In this study, our focus is directed to the chiral phase transition occurring at high baryochemical potentials and zero temperature.We specifically examine how the interplay between parity doubling and hyperonization influences the properties of neutron star matter.Chiral approaches are employed, introducing a nonvanishing nucleon mass that remains finite even during the restoration of chiral symmetry.This is achieved through the mirror assignment of chirality within the parity doublet model [10].
Here, we address the effects of hyperonization and parity doubling within neutron star composition.Our approach is built upon a prior model [25], where higher orders of the scalar potential proved to be necessary for a good description of nuclear physics and also to avoid the Lee-Wick instability [26] without a dynamically generated vector meson mass [27].Our goal is to incorporate a chirally invariant mass enabled by the inclusion of a negative parity baryon octet.This understanding is pivotal for comprehending the properties of compact stars that may potentially contain chirally restored matter.To study dense nuclear matter and investigate the chiral phase transition characterized by parity doubling, we have developed a phenomenological model that incorporates baryonic degrees of freedom and effective interactions, with constraints derived from SU(3) relations.By constraining the model with astrophysical data and properties of nuclear matter, we find that a first-order phase transition for reasonable values of the slope parameter L can occur in the core of neutron stars.
The incorporation of parity doubling in our model allows for stable static configurations of stars with a metastable matter core, enabling stars with masses higher than the expected minimum mass of a neutron star formed via corecollapse supernova [28] and around the value of the less massive observed neutron star [29] which makes metastability related phenomena particularly relevant.A key advantage of using a unified model that describes both the neutron matter phase with fully broken chiral symmetry and the approximately chirally restored phase is the ability to calculate the surface tension and to determine, without external input, the location of the phase transition and also assess the possibility of mixed phases.This information is essential for determining the timescales of phase conversion via nucleation and understanding the potential signals emitted by compact objects that can, potentially, undergo such a phase transition.
This work is organized as follows: in Sec.II, we outline the formalism.In particular, we briefly review the SU(2) parity doublet model in Sec II A, which serves as a warm up for the SU(3) model discussed in Sec.II B. The complete Lagrangian is discussed in Sec.II C, stationary equations in Sec.II D, and a discussion on how to fix the parameters is left to Sec.II E. We present our results in Sec.III.The final section, Sec.IV, summarizes our findings and points to possible future investigations and improvements within this approach.

II. THE MODEL
A. Parity doublets in SU (2) Our approach is based on the "mirror fermion" field formalism [10,11,30], in which negative and positive parity fermion chiral transformation properties are reversed.As an example, take the nucleon fields under SU (2 Then, we assume For these fields, one can construct the following chirally invariant mass terms where Ψ is a doublet in parity space over which the Pauli matrices ρ i act upon.Here, we also have the sigma field and pions.The term proportional to m 0 mixes parity fields, giving rise to a nonzero nucleon mass even when chiral symmetry is restored.At mean field level and defining σ ≡ ⟨σ⟩, the part of the Lagrangian corresponding to the mass of the fields reads It is convenient to rewrite the Lagrangian in terms of mass and parity eigenstates.This can be accomplished by the transformation that diagonalizes the fields in parity space defined by: with δ implicitly defined by sinh δ = g 1 σ/m 0 .In terms of the new defined fields, we have where B. Parity doublets in SU(3) The inclusion of baryon parity doublets can be done, as in Ref. [18], by considering the doubling of each baryon field in Eq. ( 40) of Appendix A, such that the mass term in a mean-field approximation reads s Tr Ψ{Σ, Ψ} + F (1)  s Tr Ψ [Σ, Ψ] + S (1)  s Tr (Σ) Tr ΨΨ (8) + D (2)  s Tr Ψρ 3 {Σ, Ψ} + F (2)  s Tr Ψρ 3 [Σ, Ψ] + S (2)  s Tr (Σ) Tr Ψρ 3 Ψ .
The mass matrix for each Ψ ij has nondiagonal terms that couple opposite parity states.These matrices are diagonalized by rotations such that Ψ ij , and each δ ij is determined during the diagonalization procedure analogously to what was done for the SU(2) case (see Appendix A for details).
The 16 coupling constants thus introduced, {g Σj , g Λj , g Ξj } (i = 1, 2; j = σ,ζ), are given as linear combinations the 6 independent parameters D (1) s present in Eq. (8).If one wants to respect the structure imposed by chiral symmetry, one can choose 6 of them freely.Choosing, as free parameters, g Λσ , the remaining coupling constants are given by

C. Mean-field Lagrangian density
In terms of mass and parity eingenstates, the full baryonic Lagrangian assumes the form where the index i = N, Λ, Σ, Ξ, runs through the baryon octet and the index p sums over parity eigenstates.The baryon masses are denoted by m i,p and read where we also include explicit symmetry breaking terms in the hypercharge direction as in [31] to improve on the vacuum value of baryon masses and hyperon potential depths: where S = diag(0, 0, 1).The effective chemical potentials are affected by the vector meson condensates.Assuming only the condensation of the fields ⟨ρ 0 0 ⟩ ≡ ρ, ⟨ω 0 ⟩ ≡ ω and ⟨ϕ 0 ⟩ ≡ ϕ, they read where only three of these couplings are free parameters while the rest are given by symmetry relations (see Appendix A).The mesonic part of the Lagrangian reads In our approach, higher than the usual fourth-order linear sigma potential contributions are necessary to describe the vacuum and nuclear matter properties (see Sec. II E), and to assure stability.Different relativistic mean field approaches [32,33] also find that higher-order contributions are necessary to describe nuclear matter properties correctly.We make the assumption that ζ ≡ σ, which is reasonable if we take into account that via PCAC we can relate the values of the scalar condensates to the decay constants of the pion and the kaon via: This can fulfill the condition ζ ≡ σ when considering f K around 112 MeV for pion decay constant of f π = 93 MeV, which aligns well with the experimental value of the kaon decay constant f K ≃ 110 MeV [34].There are two practical reasons for using this approximation: first, for all the tested generalized potentials only an abnormal solution with σ ≃ 0 and ζ > ζ V EV was found as the global minima solution for all densities; second, for performing most of the calculations we are concerned about, having only one order parameter is considerably simpler.With this assumption the chirally invariant piece of the vacuum potential, U 0 (σ, ζ = σ) = U 0 (Φ), must be given as a function of the invariant scalar Φ = 1 2 σ 2 .We take the Taylor expansion around the vacuum value Φ, which is done by a simple constant shift in the potential Taking σ vac = f π leads to where the parameter a 1 = m 2 π guarantees the correct pion vacuum mass, m π = 139 MeV, and an explicit symmetry breaking term was added to ensure that U (σ) has a minimum at σ = f π , with the pion decay constant set to f π = 93 MeV.In this approach, the mass of the chiral condensate is given by m 2 σ = m 2 π + f 2 π a 2 which for the range of values that the parameter a 2 assumes in our results yields m σ ∼ (500 − 700) MeV.The vacuum potential in (19) is the same adopted previously for two-flavors [35][36][37][38] and also in Ref. [25], where hyperonization was considered.
The vector meson potential in (16) reads where for the vector meson masses we take m ω = 782 MeV, m ρ = 775 MeV and m ϕ = 1020 MeV.The last term is the self interaction that introduces quartic contributions, ∼ ω 4 , to the vector interactions which are responsible for the high density behavior of c 2 s = 1/3, as was explicitly demonstrated for isospin symmetric matter [25].Moreover, the coupling ∼ ω 2 ρ 2 yields a better agreement with the experimental data on the value of the slope parameter L (for a recent overview of the various estimates for L, see Ref. [39]).The usual parametrization of the matrix containing the vector meson fields is There are two SU(3) symmetric structures for the vector self-interactions at fourth order in the vector fields: where in the last line we kept only the fields assumed to condense.In the following we set d 1 = 0 and denote d = d 2 in line of the discussion above.Finally, noninteracting electrons and muons are added to the system via where lepton masses are taken to be m e = 0.5 MeV and m µ = 106 MeV.In this work we assume that neutrinos have mean free paths larger than the size of the system, which in our case is the typical size of a cold neutron star (∼ 10 km).For this reason, they will be omitted throughout the rest of the text.

D. Free energy and stationary equations
In the no-sea approximation (fermionic vacuum contributions ignored), the standard field theory computation of the zero-temperature free energy density leads to with Ω hom defined as which depends only on the condensates themselves, not their derivatives, being the part of the potential that determines the homogeneous solution.The electric charge is denoted e and, in Heaviside-Lorentz units, has a value of e ≃ 0.3.In Eq. ( 25) the partial pressures have the form with the Fermi momentum k F = µ 2 − m 2 .The chemical potentials for each baryon are not independent, but in three-color QCD they are fundamentally related to the three chemical potentials of u, d and s quarks.The condition of electroweak equilibrium makes it possible to narrow this freedom down to two chemical potentials: µ B , related to the baryon number density, and µ e (µ µ = µ e , in beta equilibrium), related to the lepton number density that will affect only the charged baryons, such that: independent of the parity of the baryonic state.From the free energy density of Eq. ( 24), we derive the following Euler-Lagrange equations where with We can also derive a Euler-Lagrange equation for µ e such that local electric charge neutrality corresponds to ∇ 2 µ e = 0.

E. Parameter fixing
In this section we explain how the free parameters of the model are fixed to reproduce properties of the vacuum and of isospin symmetric nuclear matter at saturation (see Appendix B for more details).The parameters g  Other assignments are also possible since the experimental data shows multiple resonances that could potentially be identified to the chiral partners of the baryon in the baryon octet [34,40].
One still has to fix the values of 12 other parameters: a 2 , a 3 , a 4 , d, g N ω , g N ρ , g N ϕ , g N σ , g Σσ , g Λσ , m 1 , m 2 .This set of parameters will be fixed with the help of 9 properties of isospin-symmetric nuclear matter at saturation and the MeV and U Ξ = −20 MeV.We checked that mild modifications are found if we instead consider other values for hyperon potentials within U Σ ∈ [0, 30] MeV and U Ξ ∈ [−30, 0] MeV.These values have been selected by considering experimental evidence that makes reasonable placing the U Λ value near -28 MeV and the U Ξ values near -20 MeV (albeit with less certainty).Although the U Σ values remain uncertain, it is noteworthy that Σ baryons are absent in hypernuclei bound states.This observation provides support to the notion of a repulsive, positively valued optical potential.For recent research and comprehensive discussions regarding hypernuclei physics and its relevance to neutron star physics, see Refs.[41][42][43][44][45][46][47][48].The rest of the couplings are all determined by the SU(3) relations (11a) and (45).In our approach the masses of the cascade baryon m vac Ξ,+ and its parity counterpart m vac Ξ,− are both determined after the parameter fixing.For the parameters adopted here, m vac Ξ,+ ∼ 1330 MeV (see Table I) which is in good agreement with usually adopted mean masses of 1318 MeV.The value of m vac Ξ,− in our approach is m vac Ξ,− ∼ 1845 MeV.At high densities, the low-energy properties of matter provide limited information, and it is unreasonable to expect that the relevant degrees of freedom are the same as those in confined matter.At some point, partonic degrees of freedom, as predicted by perturbative QCD, become relevant.To overcome the lack of experimental microphysics input at high densities, we incorporate constraints from astronomical observations, which significantly restrict the high-energy portion of the equation of state derived from our model (more on this matter in the next section).Furthermore, we enforce the condition that, at asymptotically high densities, the speed of sound approaches the conformal limit of c 2 s → 1/3.Within our model, this is achieved through the quartic self-coupling of the vector meson.

A. Parameter space
We narrow down our parameter space by selecting values that allow for the existence of stable star configurations with masses of at least 2.0 solar masses.We solve the Tolman-Oppenheimer-Volkoff equations using as input the equation of state (EoS) obtained via the usual thermodynamic relation ϵ = −P + µ B n B for densities higher than 1.1 n 0 .For lower densities, we adopt the EoS proposed by Hebeler et al. [49].Although the use of this low-density chiral effective theory EoS has a negligible impact on the maximum mass of neutron stars, it provides significant corrections to the radii, particularly for the less massive stars [50][51][52][53].Therefore, it is reasonable to impose additional constraints on the parameter space by considering recent constraints on the radii of neutron stars.We specifically select parameters that lead to neutron star families falling within the one-sigma error bars of the measured radii of 1.4 M ⊙ and 2.0 M ⊙ stars consistent with the analyses conducted in Refs.[54][55][56][57] using data from NICER on the pulsars PSR J0030+0451 and PSR J0740+6620.This approach allows us to identify the different qualitative characteristics that emerge from the possible parameter choices within reasonable experimental and observational constraints.In our analysis, we primarily focus on varying key properties, namely the effective mass at saturation M 0 , the invariant mass m 0 and the slope parameter L. These parameters have a significant impact on our results and are responsible for the most pronounced changes.For simplicity, we only vary these parameters as minor modifications within their experimental uncertainty range.The remaining nuclear properties yield negligible changes in our final conclusions.
For a given baryon chemical potential, we solve the stationary equations [ignoring the Laplacian contribution in Eqs.(28)] for homogeneous meson condensates, incorporating the equation of local charge neutrality and enforcing beta equilibrium.Figure 1 displays the M 0 × m 0 parameter space for three values of L = 65, 70 and 80 MeV, respectively.The parameter space is limited from above by the condition that the potential U (σ) remains bounded, which translates to the requirement that parameter a 4 in Eq. ( 19) is greater than zero.The excluded region are shaded in gray.The blue region corresponds to parameters that lead to star masses and radii within the astrophysical constraints.The red region represents parameter values for which the chiral restoration occurs via a crossover.
A crossover can be attained for the highest values of effective mass M 0 and invariant mass parameter m 0 , being restricted to the upper right part of the parameter space in all cases.Its intersection with the astrophysical region becomes smaller for lower values of the slope parameter L. For L < 70 MeV, and therefore toward more realistic values of the slope parameter, the astrophysical constraints become incompatible with a crossover.We also observe that lower values of L tend to drag the blue region toward lower values of M 0 and higher values of m 0 for each given M 0 .
To exemplify the different ways that chiral restoration can take place in the model, we take four points in the allowed parameter space, corresponding to the asterisks in the middle panel of Fig. 1, with L = 70 MeV, and with values of M 0 and m 0 in Table I.We also show the critical chemical potentials for the first-order phase transitions and the  In all panels the blue region corresponds to parameters that satisfy the astrophysical constraints, the red region corresponds to parameters that lead to crossovers.The gray region leads to a unstable vacuum potential with a4 < 0. The asterisks mark the sample sets of parameters of Table I potential for the onset of strangeness, µ c and µ S , respectively.In all cases, the negative parity baryons just onset at the (approximate) chiral restored phase being present already just after the transition, which implies that the density of negative parity partners serves as a order parameter for the transition, being zero at the chirally broken phase and nonzero at the chirally restored phase.We also provide the corresponding m 1 and m 2 parameter values related to explicit symmetry breaking in the hypercharge direction introduced in Eq. ( 14).The significant m 1 values indicate a considerable strange quark bare mass.As our focus is on understanding the broad phenomenological effects of the interplay between parity doubling and hyperonization, we anticipate our outcome to be only qualitatively similar to a more realistic model with a realistic strange quark mass (m s = 95 MeV).In this context, we prioritize the accurate representation of nuclear properties, such as baryon masses and hyperon potential depths, rather than emphasizing the reproduction of bare quark masses, which are not even treated as degrees of freedom in our approach.Figure 2 illustrates the behavior of the scalar condensate as a function of the chemical potential for µ B ≥ 980 MeV, which corresponds to a baryonic density of 1.1 n 0 in all cases.The figure includes stable (solid lines) and spinodal curves (dashed lines).The turning points where the first-order transition occurs are marked with dots.In all cases the high density phase is just approximately chirally symmetric, approaching σ = 0 only for very high chemical potentials, which is expected since there are explicit chiral symmetry breaking terms in the Lagrangian.
A remarkable aspect of the model presented in Ref. [25], which did not consider parity doublets, is its implication for the formation of massive stars.In that context, stars with 2.0 solar masses can only be formed when a robust and early first-order phase transition takes place, resulting in stars predominantly composed of chirally restored matter.However, with the inclusion of negative parity baryons, the chiral phase transition, while remaining firstorder, gains the flexibility to manifest in different regions of the phase diagram.When examining curves that share the same invariant mass and yet exhibit different effective nucleon masses, it becomes evident that the main parameter influencing the critical potential for the phase transition is the magnitude of the chiral invariant mass m 0 , assuming larger values for larger values of m 0 .Furthermore, variations in the nucleon effective Dirac mass M 0 appear to directly impact the jump in the order parameter ∆σ during the phase transition, with ∆σ consistently increasing as M 0 increases.This observation suggests that adjusting the values of the invariant mass within the blue region of Fig. 1 allows for precise tuning of both the position and the magnitude of the phase transition.This "tuning" feature can be used to explore different phase transition scenarios within the same model, all consistent with the same experimental constraints.
The invariant mass, represented by m 0 , plays a pivotal role in quantifying the contribution of chiral symmetry breaking to the generation of baryon mass.A higher invariant mass signifies that the chirally restored phase primarily consists of massive baryons.On the other hand, as m 0 → 0, the chirally restored phase becomes dominated by lighter,   almost massless baryons.In a broader perspective, m 0 may arise through the condensation of scalar fields, such as the glueball condensate [12,58].If the theory exhibits dilatation invariance, we would anticipate m 0 to vary with density, decreasing as density increases.In our present context, our approach can be viewed as a "frozen glueball" description, where we neglect any density dependence of the glueball condensate.The running of m 0 and its connection to the trace anomaly of QCD and deconfinement will be postponed to a later work.With this in mind, our focus in the subsequent section of this study will solely be on the value of M 0 /M N = 0.65, corresponding to the black set of parameters.

Domain wall
The surface tension is determined by the classical configuration that connects different homogeneous phases [25,35,36,[59][60][61][62][63][64][65][66], as specified by Eqs.(28).In the specific case of domain walls, the surface tension is well-defined as the energy difference per unit area between the configuration of the domain wall and the homogeneous configuration of each phase.Essentially, the surface tension represents the energy cost associated with transitioning from the first minimum to the second minimum by following a path in configuration space that satisfies the equations of motion.The domain wall solution can be seen as the critical bubble solution when its radius approaches infinity, particularly applicable in the immediate vicinity of the phase transition.
In the domain wall geometry, the Euler Lagrange equations become one-dimensional.Thus, one can make the simplification ∇ 2 → d 2 /dx 2 in Eqs.(28).To solve the system of differential equations at the phase transition, we employ the numerical method of over-relaxation with boundary conditions b±(x = ±∞), where b = (σ, ω, ρ, ϕ, µe).At first glance, one might consider using the degenerate solutions at the phase transition as boundary values.However, a subtle issue arises due to the disparity in lepton chemical potential between the phases connected via a Maxwell construction, despite having identical pressure.Interpolating between these lepton chemical potential values would imply the presence of a nonhomogeneous lepton background across the interface region, which contradicts expectations during nucleation timescales that should be much smaller than the typical timescales associated with electromagnetic forces.To address this problem, we set µ e to its value in the chiral broken phase µ e,χBroken .This decision is based on our primary focus being the nucleation process of the chirally restored phase within the predominantly chirally broken phase.Subsequently, we search for the chirally restored solution that possesses the same lepton chemical potential.Although this adjustment slightly modifies the critical value of the baryon chemical potential for the phase transition, it ensures that, when calculating the surface tension, we interpolate between condensate solutions embedded in a On the other hand, the region shaded in blue indicates values of m0 for which nuclear matter exists solely as a metastable solution and the chirally broken phase used in the surface tension calculation is the vacuum phase.
Once the numerical solution is obtained, the surface tension is computed as To evaluate the error in the choice of lepton chemical potentials, we also conduct the same calculation using the lepton chemical potential of the chirally restored phase µ e,χRestored .Figure 3 exhibits the surface tension as a function of the invariant mass for both charge chemical potentials.The yellow-shaded region represents the invariant mass values that meet astrophysical constraints, while the blue-shaded region corresponds to values of m 0 that lead to nuclear matter existing solely as a metastable solution.For the case of metastable nuclear matter, the surface tension is computed by connecting the vacuum phase, where all baryons have their vacuum masses, directly to the chirally restored phase.Both calculated surface tensions exhibit similar values when m 0 assumes essentially the same magnitude for m 0 > 400 MeV, with the chirally broken chemical potential yielding slightly higher values for Σ.In the region of 225 MeV < m 0 < 400 MeV, the µ e,χRestored results in a higher surface tension value relative to the values associated with µ e,χBroken .In the bluish region, the behavior is reversed, with the chirally broken chemical potential leading to larger values of Σ after a jump occurring at the point where nuclear matter ceases to be stable.This jump is associated with a corresponding leap in the values of the lepton chemical potential, shifting from a finite value for m 0 > 225 MeV to exactly µ e = 0, which represents the vacuum value.Interestingly, achieving physical masses and radii for neutron stars requires very low values of surface tension.These values are approximately an order of magnitude smaller than those previously calculated in works such as [35,36,[67][68][69][70], where no parity doubling was considered.The reduced surface tension values directly stem from the fact that the energy density discontinuity connecting the phases is much less pronounced when the baryonic species remain massive in the chirally restored phase, which explains the monotonic decreasing behavior of its value with respect to the invariant mass.

Bubbles
To estimate the behavior of the surface tension in the spinodal regions, where the nucleating critical bubble assumes a finite size, we assume spherical symmetry [59][60][61][62][63][64][65][66].The critical bubble profile is obtained by solving Eqs. ( 28 phase, which have a higher free energy.The values of the condensates at the center of the bubble are dynamically determined via a shooting method algorithm. To calculate the bubble profile for a single field, the standard semiclassical approach [62] is mathematically equivalent to the problem of a point particle moving along a potential valley under the influence of a time-dependent viscous force1 .To leverage this approach, we utilize the "one-condensate approximation" described in Ref. [35].This approximation involves neglecting spatial derivatives in Eqs.(28b,28c,28d), allowing us to solve these three equations for (ω(σ), ρ(σ), ϕ(σ)) for a given σ value.Consequently, the Euler-Lagrange equation for the σ condensate simplifies to a straightforward second-order differential equation: where Ω ≡ Ω(σ, ω(σ), ρ(σ), ϕ(σ), µ e ) represents a function that depends solely on the σ condensate.We fix the lepton chemical potential at its values within the chirally broken, metastable phase in the same fashion performed previously, in the case of the domain wall.
In the one-condensate approximation, we define the bubble surface tension as Figure 4 illustrates the behavior of the surface tension in the spinodal region, which has been normalized by the corresponding critical chemical potential µ c in each case.Since the critical bubble size diverges at the phase transition, the bubble computation is initiated at µ c + 2MeV toward the end of the spinodal region, where Σ = 0.The solutions obtained within the domain wall approach, with µ e = µ e,χBroken , are denoted by the red squares, serving as a consistency check between the numerical methods employed.One can see a clear agreement between the results obtained from both approaches, validating their reliability.Higher values of m 0 correspond to narrower spinodal regions and smaller magnitudes of Σ.

Mixed phases
So far, we have described a first-order phase transition using the Maxwell construction.However, due to the relatively low values of surface tension observed in the previous section, there is a possibility that mixed phases, rather than homogeneous matter, could emerge.For m 0 = 460 MeV, the value of surface tension is Σ = 0.32 MeV/fm 2 as seen in Figure 3.To explore this possibility, we adopt the Gibbs construction, relaxing the local charge neutrality condition to global conditions.
In the Gibbs construction, we disregard the Coulomb and surface energies of the crystal structures and focus on calculating the energy for solutions with a nonspecified structure but with a chirally restored volume fraction χ, with χ ∈ [0, 1].The free energy density of the homogeneous solution (red curve) relative to the Gibbs constructed heterogeneous phases (black line) is illustrated in Fig. 5.The dotted line marks the critical baryon chemical potential.We observe that the Gibbs constructed phase is only slightly favored when compared to Maxwell constructed solutions.
Including the Coulomb and surface corrections to the Gibbs solutions and searching for specific structures of various sizes such as cylinders, rods, and bubbles, that minimize the free energy density per unit volume, we find the mixed phase solution represented by the purple curve in Fig. 3. Nevertheless, it becomes evident that despite the low values of surface tension, homogeneous phases are still energetically favored over mixed phases in this model.This outcome is attributed to the minimal energy gain in constructing heterogeneous phases, as demonstrated by the marginal difference in energy density between the Gibbs and Maxwell constructed phase transitions.
This finding should be compared to the previously encountered mixed phases found to be favored in neutron matter near the chiral phase transition [36].The introduction of hyperons is expected to hinder the occurrence of mixed phases, as the additional baryonic hyperon degrees of freedom alleviate the energy cost of isospin asymmetry in charge-neutral homogeneous phases, lowering their energy cost in relation to the Gibbs solutions.The same effect occurs, for example, in hadron-quark mixed phases [71].In this study, we observe a different scenario, where the mixed phases consist of a chirally broken nucleonic phase and chirally restored parity-doubled neutrons (see the left panel of Fig. 6).Hyperons only appear at significantly higher densities, so they cannot explain the disappearance of mixed phases.Apart from these differences, parity doubling and the presence of an invariant mass m 0 also act to prevent the formation of mixed phases.C. Neutron stars

Hyperonization and parity partners
In this section, we present our results regarding the impact of parity doubling on the composition of neutron stars and the stiffness of the equation of state.Figure 6 compares the resulting particle fraction with and without parity doubling for zero-temperature, beta-equilibrated matter using the same parameters of M 0 /M N = 0.65 and m 0 = 460 MeV.On the left panel, we prevent the appearance of negative parity partners by decoupling the negative baryon states by assigning an extremely high mass to these states.As a result, their onset are artificially pushed to very high chemical potentials.On the right panel, we allow the onset of all parity partners by keeping their vacuum masses as described in Sec.II E. We observe that the inclusion of negative parity states induces a first-order phase transition (marked by the blue band), with the N(1535) partner of the neutron emerging as the second most abundant baryonic state immediately after the phase transition.Additionally, the onset of strangeness, indicated by the presence of Σ baryons, is shifted to significantly higher densities, transitioning from approximately n B /n 0 = 3.02 to n B /n 0 = 4.68.Moreover, the particle fraction of strangeness is reduced in relation to the amount of its value in when parity doubling is excluded, going from around 10% to around 7% at the maximum density achieved by the most massive stable stars built using the corresponding equations of state, which corresponds to the red vertical red line in the figure.We also observe that this line is pushed to higher densities when parity doubling is considered, which makes the cores of those stars more densely packed with baryons.
To distinguish between the impact of parity doubling and hyperonization on the equation of state (EoS), we compare these two scenarios in Figure 7. On the left panel, we present the EoS considering only nucleons and allow for the appearance of N(1535) states but no hyperons.Just after the phase transition, when the N(1535) states onsets, the parity doubled EoS has less pressure for the same energy density of the case without parity doublets, which is related to the fact that after the transition the system has more degrees of freedom to fill up as the chemical potential rises.
A significant aspect of parity doubling within this model is the preservation of a relatively consistent slope in the EoS, similar to the scenario involving only nucleons.This uniformity arises due to the identical coupling with the vector meson sector for both parity-based baryonic states.Consequently, the contribution of N(1535) to the repulsive component of the nuclear force becomes comparable to that of the nucleon.Correspondingly, when considering a specific energy density, the speed of sound squared, denoted as c 2 s = ∂P ∂ϵ , remains indistinguishable between the two setups.On the right panel, we observe softening caused by the onset of hyperons with the inclusion of both positive and negative parity baryons.In this case the onset of strangeness is continuous and the softening becomes more pronounced as the fractional amount of strange baryons in the system increases for higher energy densities.This progressive softening stems from the fact that hyperons tend to couple weakly to the vector sector alleviating the repulsion and, consequently, the pressure.This results in a decrease in the value of c 2 s as pressure increases, in contrast to the scenario where strangeness is absent.It is important to note that all equations of state were computed using the same parameters, which leads to the same couplings.Also notice that the softening resulting from hyperonization occurs relatively late.This delay is due to the fact that parity doubling, as explained previously, tends to push hyperonization to higher densities.
The stellar configurations corresponding to the equations of state (EoS) shown in Fig. 7 are depicted in Fig. 8.The gray ellipses represent the results of combined one-sigma confidence analyses extracted from the studies [54][55][56][57].
The black bars indicate the radius values of 1.4 and 2.0 solar masses.On the left panel, the EoS exhibits a noticeable softening effect due to parity doubling, particularly in the context of nucleons.This softening effect has a more pronounced impact on the more massive stars, whereas it has a minimal effect on the 1.4 solar mass star, which has a relatively small chirally restored core.Despite this overall reduction in radii, none of the curves deviate to an extent that would render them incompatible with NICER data.In our analyses, this compatibility is represented by a mass-radius configuration that remains within the elliptical regions.
In contrast, the right panel illustrates the inclusion of hyperonization across the entire baryon octet, incorporating both parities.The inclusion of hyperons does not lead to a significant alteration in the maximum mass when parity doubling is present.Hyperons emerge only at high densities, thereby impacting solely the most massive stars with masses exceeding 2.0 solar masses.Consequently, in this scenario it is anticipated that only the most massive stars will have a small core region where hyperonization can occur.
In Fig. 9, we compare the increased repulsion within the system for a given density.As in this work we assume that negative parity partners couple in the same way with the vector meson sector, these curves are identical for corresponding parity partners.The vertical red line indicates the maximum density achieved in the most massive stable star.Additionally, the vertical dotted line represents the onset of the Σ − baryon, the only hyperon that onsets in our model.For a single condensate, e.g.ω condensate, a simple comparison between g i,ω suffices.However, given   9: Y-axis: gi,ωω + gi,ρρ + g i,ϕ + µ i,charge , with i = n, p, Σ.A vertical red line indicates the maximum density achieved in the most massive stable star.Additionally, a vertical dotted line represents the onset of the Σ − baryon, the sole strange baryon onset in our model.The plotted quantity serves to quantify the increased repulsion within the system as the number of the ith species is elevated for a given density.
the presence of ρ and ϕ condensates, the plotted quantity serves to quantify the increased repulsion within the system as the number of the i th species is increased.This analysis provides insight into the dynamic interplay of the repulsive forces introduced by each species in the system.The figure makes it clear that the addition of Σ − baryons contributes significantly less to the overall repulsion of the system than the addition of nucleons (and their parity partners).Supporting the main conclusion, while both hyperonization and parity doubling soften the equation of state due to the extra Fermi levels available, hyperonization softens relatively more when it comes to repulsive interactions alone.As parity doubling tends to push hyperonization to occur at higher densities, its net effect is to stiffen the EoS when compared to what is expected from hyperonization alone, without parity doubling.

Metastable cores
The advantage of having a single description for both the chirally restored and chirally broken phase is that phenomena associated with metastability and their timescales can be investigated in a unified framework.We see from Figs. 2 and 4 that, in the cases where the chiral restoration occurs via a first-order phase transition, the spinodal region can extend through a considerable range of baryon chemical potential.In Fig. 10 we show the resulting family of stars obtained by solving the TOV equations for the EoS corresponding to the black parameter set of Table I.If we follow the thermodynamically stable branches through the phase transition, we branch off of the neutron star curve by following the chirally restored branch.We obtain hybrid stars, shown by the orange curves in Fig. 10, i.e., stars with a chirally broken mantle and a chirally restored core.If instead we follow the thermodynamically metastable solution after the phase transition we obtain the purple curve, where all the star is formed by matter in its chirally broken phase, but in the core matter is in a metastable phase (see Fig. 10).The orange band highlights the mass range at which both solutions are possible and we can encounter this mass degenerate star configurations.These bands begin when the critical central pressure for phase transition is reached at the center of the star and end where the metastable branch reaches the spinodal point and no metastable solution exists for greater central pressure.
The gravitational and baryonic masses are defined, as usual, by [72] so that the gravitational mass is defined as M G ≡ M (R) and with the baryonic number for each species A i,p given by Since no known physical process violates baryon number conservation in compact stars, M A should be conserved in isolated stars.Thus, the difference between baryonic mass and gravitational mass can be used as a measurement of the binding energy of the star B = M A − M G .We have checked that for the same baryonic mass the change in gravitational mass associated with the core phase conversion is ≲ 0.02% as shown on the left panel of Fig. 11 for the metastable stars and for the completely stable stars corresponding to the orange band of Fig. 10.This energy difference, associated with the change in binding energy, could in principle drive a late heating and late emissions from the star, depending on the timescales for the phase transition.On the right panel of Fig. 11, the difference in radii between metastable and stable stars as a function of baryonic mass is shown.We see that metastable stars are only slightly larger than the correspondingly (same gravitational mass) stable stars, and the order of magnitude of the radius difference is ∆R ≲ 20m.
Given the relatively small values of surface tension, as seen in Sec.III B, it is not unreasonable to anticipate that these metastable stars would have a limited lifespan.Nevertheless, the metastable phases can be relevant in the early development of protoneutron stars.Additionally, the relatively short timescales involved indicate that, in scenarios characterized by dynamical processes, where localized areas can undergo destabilization due to fluctuations in pressure or lepton number, the transition from a state of chirally restored matter to a state of chirally broken matter (or vice-versa) would be constrained primarily by the timescales of weak decay that drives beta equilibration.
To illustrate this point, consider the case of a neutron star undergoing vibrations.These oscillations lead to variations in pressure, with the response differing according to the timescales of phase conversion.Another instance in which such dynamical configurations can manifest is in the aftermath of compact star mergers and in core-collapse supernovae.In these scenarios, regions of lepton-rich matter can stabilize areas of chirally broken symmetry.As the system cools and progresses toward lower leptonic densities and temperatures, these regions can destabilize, and transition to the chirally restored phase in a timeframe determined by weak interactions.

IV. SUMMARY AND OUTLOOK
We have investigated the impact of chiral symmetry restoration through parity doubling in beta-equilibrated cold matter.To achieve this, we built a chiral Lagrangian that incorporates the baryon octet and their corresponding negative parity partners.A crucial aspect of this approach is the enabling of a chiral-invariant mass m 0 , which is possible because of the specific way the opposite parity baryon fields transform.Consequently, the chirally symmetric phase is not solely composed of massless baryons.Our model was calibrated to reproduce nuclear matter properties at  saturation density.By exploring the parameter space, we identified certain parameter sets that satisfy the mass and radius constraints imposed by recent astronomical observations without the need to include a quark matter phase.The early onset of the N(1535) resonance, which occurs immediately after the phase transition, inhibits strangeness formation and leads to hyperonization at higher densities compared to standard chiral restoration schemes for baryonic matter, which constitutes another solution to the "hyperon puzzle" problem.
We found a relatively subtle first-order transition into the chirally restored phase, with position of the critical density and energy density discontinuity mainly dictated by the chiral-invariant mass m 0 .The corresponding surface tension is one order of magnitude smaller than prior results using similar baryon-meson approaches.It exhibits a strong dependence on the value of m 0 , given its direct correlation to the alteration in effective baryonic mass across the phase transition.We employed the Maxwell construction, which remains applicable even for low surface tension values.Its validity stems from the energy associated with forming mixed phases, predominantly governed by the Coulomb contribution arising from the global charge neutrality condition.This condition favors the existence of a sharp interface between phases at the critical chemical baryon potential.
Investigating the impacts of parity doubling and hyperonization on the EoS led us to conclude that parity doubling significantly influences stellar radii of the most massive stars, generally causing their reduction.However, this effect is not big enough to be resolved by NICER measurements since the change in radii is smaller then the one-sigma ellipses.
By utilizing both the metastable and stable branches of the equation of state, we compared stars containing a metastable core to those composed exclusively of the most stable homogeneous solutions.Despite the apparently short lifespans for these metastable cores at finite temperatures due to their low surface tension values, the phase conversion and subsequent contraction of stars could potentially influence dynamical phenomena in neutron stars with masses around 1.3 solar masses or greater.In dynamical scenarios, where local regions could destabilize due to pressure or lepton fluctuations, the transition between chirally restored and chirally broken matter depends on weak decay timescales.For example, neutron star vibrations induce pressure changes with responses based on phase conversion timescales.Similarly, after compact star mergers and in core-collapse supernovae, lepton-rich zones stabilize chirally broken symmetry areas.As systems cool and densities drop, these stable regions may destabilize and shift to the restored phase in a time frame controlled by weak interactions.A detailed calculation of these timescales will be left for a future investigation.
One generalization of this model concerns the physical interpretation of the chiral-invariant mass m 0 .Possible dynamical ways of generating this term come from interactions with tetraquark condensates and the gluon condensate [12,58].Additionally, the inclusion of the ζ condensate as an independent parameter, or the consideration of different monotonic functions ζ(σ) = f (σ) instead of simple identification, is desirable.
Another relevant issue involves the inclusion of the vacuum contribution, commonly referred to as the "Dirac sea."This contribution is usually disregarded within relativistic mean field models applied to neutron stars and nucler physics, since vacuum effects were shown to have minimal impact on the equation of state [73].However, in a variety of models involving chiral symmetry restoration, it is known that a first-order chiral transition can be smoothed to a crossover when vacuum terms are included, bringing qualitative modifications [74][75][76][77], for a recent work on chiral density waves where the nucleonic vacuum contributions are included in a similar model see [78].
Our parity doubling framework can be extended to nonzero temperatures and scenarios beyond beta equilibrium.A crucial aspect in this context is the examination of the direct URCA process, which involves rapid neutrino emission in neutron stars.Neutrons transform into protons, releasing an electron and an electron antineutrino, while protons capture electrons, becoming neutrons and emitting electron neutrinos.This cooling mechanism carries away thermal energy via diffusion and radiation of neutrinos and antineutrinos.Its occurrence depends on specific conditions, with a critical proton fraction denoted as Y DU p .Notably, parity doubling modifies this critical value leading to an overall reduction of this critical proton threshold, as suggested in Ref. [79].
Expanding this model to finite temperatures also offers opportunities to explore the chiral phase transition in core-collapse supernovae, the evolution of protoneutron stars and neutron star mergers.In core-collapse supernovae, gravitational pressure cause the collapse of massive stars.The interplay between chiral symmetry restoration and exotic hadronic state formation significantly influences the explosion mechanism and subsequent nucleosynthesis.Neutron star mergers involve extreme conditions like high temperatures and rapid density changes.Understanding the state of nuclear matter formed in such processes, its phase transitions and their impact on gravitational wave emission and heavy element production through r-process nucleosynthesis is pivotal.This is especially timely given the ongoing operations of LIGO/Virgo/Kagra in the current O4 run [80], promising increased statistics and precision for gravitational wave signals from neutron star mergers.

Scalar mesons -baryon sector
In the nonlinear realization of SU(3) chiral symmetry, the various interaction terms of baryons with mesons have the same structure, except for the difference in Lorentz space [31,81].For a general meson field W they read where ΨOΨ F := Tr ( ΨOW Ψ − ΨOΨW ) and BOΨ D := Tr ( ΨOW Ψ + ΨOΨW ) − 2 3 Tr W .By a redefinition the parameters g W 8 , g W 1 and α W Eq. ( 38) is rewritten as For the interaction between baryons and scalar (and pseudoscalar) mesons, one takes W = X and O = 1.The baryon degrees of freedom are parameterized as usual The scalar and pseudoscalar meson nonets are encoded in the field X = Σ + iΠ = T a (σ a + iπ a ), where T a = λ a /2 for a = 0, . . ., 8, with the Gell-Mann matrices λ a for a = 1, . . ., 8 and λ 0 = 2/3 1.This is usually reparametrized as In the scalar sector, one can trade σ 0 and σ 8 for nonstrange and strange scalar fields by the transformation At mean field level (⟨π a ⟩ = 0, ⟨σ 0 ⟩ = σ 0 and ⟨σ 8 ⟩ = σ 8 ), the Lagrangian of Eq. ( 8) in terms of the transformed fields of Eq. ( 9) leads to the masses of Eqs.

VII. APPENDIX B
In this appendix, we provide an explanation of the parameter-fixing procedure used in this model.Most of the calculations are derived from Appendix A of Ref. [36].However, it is still valuable to demonstrate how those ideas are applied to the current model.Aside from the similarities, certain differences also become apparent.
At saturation density (where there are no hyperons), we have n B = n n + n p , n I = n n − n p , and we shall need µ B = (µ n + µ p )/2, µ I = (µ n − µ p )/2 (which is true in general).The stationarity equations ( 28 For isospin-symmetric matter, we have n I = 0 and thus also ρ = 0.In this case, we write the solutions for ω and ϕ as with and With lim x→0 f (x) = 1 we recover the case without quartic vector meson interactions, d = 0.The pressure of isospinsymmetric nuclear matter at saturation equals the vacuum pressure P vac = 0 and can be written as where µ * B = µ * n = µ * p = k 2 F + M 2 0 , and the Fermi momentum is given in terms of the saturation density via n 0 = 2k 3 F /(3π 2 ), while the stationarity equation (46a) at saturation for symmetric nuclear matter is The relation µ * B = µ 0 − g N ω ω 0 − g N ϕ ϕ 0 together with Eq. ( 47) can be used to write g N ω in terms of saturation properties and the coupling constant d, where µ 0 = 922.7 MeV is the value of µ B at saturation.The incompressibility is defined as where we assume isospin symmetry, n I = 0. We compute the derivative from For the first term, we insert Eq. ( 47) and take the derivative with respect to n B , while for the second term we can simply follow appendix A of Ref. [36].This yields .
Again, we can check that we recover the d = 0 result with f ′ (x) ∝ x for small x.
The symmetry energy is defined as where the derivative is taken at fixed n B and evaluated at n I = 0.This derivative is computed from We find from Eq. (46d) For the derivative of the second term in Eq. ( 57) we can again simply follow appendix A of Ref. [36].Consequently, With Eq. ( 47) we can use this result to write g N ρ as The slope parameter of the symmetry energy is defined as where the derivative is taken at fixed n I = 0, i.e., we can simply take the derivative of Eq. ( 59) with respect to n B .We compute We now use which follows from Eq. (46b), and we express ∂M ∂n B in terms of K with the help of the first line of Eq. ( 55).This yields at saturation The values of the hyperon potential depths which are given, explicitly U i = m i,+ (σ 0 ) − m i,+ (f π ) + g iω ω 0 + g iϕ ϕ 0 , (i = Λ, Σ, Ξ). (65) are fixed by the vacuum mass splitting between parity partners for N, Λ and Σ baryons (Eqs.(10a), (10b) and (10c), for ζ = σ = f π ).The mass of the negative parity partners are usually assumed to be m vac N,− = 1535 MeV, m vac Σ,− = 1750 MeV and m vac Λ,− = 1670 MeV.

Figure 1 :
Figure 1: From left to right the parameter space for L = 65, 70 and 80 MeV.In all panels the blue region corresponds to parameters that satisfy the astrophysical constraints, the red region corresponds to parameters that lead to crossovers.The gray region leads to a unstable vacuum potential with a4 < 0. The asterisks mark the sample sets of parameters of TableI.

Figure 2 :
Figure 2: Scalar condensate as a function of baryon chemical potential for the four parameter choices of TableI.Solid (Dashed) lines correspond to stable (metastable) solutions in each case.

2 ]Figure 3 :
Figure3: The surface tension as a function of m0 in the domain wall setup for both lepton chemical potentials at the phase transition.The region shaded in yellow corresponds to the values of the invariant mass that satisfy astrophysical constraints.On the other hand, the region shaded in blue indicates values of m0 for which nuclear matter exists solely as a metastable solution and the chirally broken phase used in the surface tension calculation is the vacuum phase.

Figure 4 :
Figure 4: Surface tension of bubbles in the spinodal region for m0 = 250, 350, 450 MeV.Red squares mark the corresponding domain wall solution at the phase transition.

Figure 5 :
Figure 5: Free energy densities of mixed phases (purple curve) and homogeneous phases (red curve) relative to the Gibbs constructed solutions (black horizontal line) for Σ = 0.32 MeV/fm 2 .The dotted line marks the critical baryon chemical potential.

Figure 6 :
Figure 6: Particle fractions for M0/MN = 0.65 and m0 = 460 MeV.On the left panel, parity partners are decoupled and matter consists solely of the baryon octet and leptons.On the right panel all parity partners are included.The blue band represents the discontinuity associated with the phase transition, while the red vertical lines indicate the maximum density reached within the interiors of the most massive stable stars built using the respective equations of state.

3 ]Figure 7 :Figure 8 :
Figure 7: Comparing parity doubling and hyperonization effects on the EoS.Left panel: EoS for nucleons and N(1535) states, revealing pressure reduction with parity doubling onset due to increased degrees of freedom post-transition.Parity doubling maintains consistent EoS slope through identical vector meson coupling.Right panel: Hyperonization-induced softening, stronger with rising strange baryon presence.Both positive and negative parity baryons were included.Weak hyperon-vector coupling eases repulsion, lowering pressure and speed of sound.All EoS were computed for the same parameters and couplings.

Figure
Figure9: Y-axis: gi,ωω + gi,ρρ + g i,ϕ + µ i,charge , with i = n, p, Σ.A vertical red line indicates the maximum density achieved in the most massive stable star.Additionally, a vertical dotted line represents the onset of the Σ − baryon, the sole strange baryon onset in our model.The plotted quantity serves to quantify the increased repulsion within the system as the number of the ith species is elevated for a given density.

Figure 10 :
Figure10: Left panel: mass-radius relations for the stable and metastable branches, distinguished by the orange and blue curves, respectively.The shaded orange band represents the mass range within which both configurations coexist.The gray ellipses correspond to constraints from NICER with black bars marking radii of 1.4 and 2.0 solar mass stars.Right panel: baryon density vs. star radius for stars with an identical baryonic mass of 1.3 solar masses.These stars are built using solutions from both branches of the equation of state, resulting in distinct density-radius profiles at the core.

Figure 11 :
Figure 11: Gravitation mass (left panel) and radius (right panel) for star configurations with the same baryonic mass.Both branches correspond to the metastable and stable EoS obtained for parameters M0/MN = 0.65 and m0 = 460 MeV.