Parity Doubling and the Dense Matter Phase Diagram under Constraints from Multi-Messenger Astronomy

We extend the recently developed hybrid quark-meson-nucleon model by augmenting a six-point scalar interaction and investigate the consequences for neutron-star sequences in the mass-radius diagram. The model has the characteristic feature that, at increasing baryon density, the chiral symmetry is restored within the hadronic phase by lifting the mass splitting between chiral partner states (parity doubling), before quark deconfinement takes place. At low temperature and finite baryon density, the model predicts a first-, second-order chiral phase transition, or a crossover, depending on the expectation value of the scalar field, and a first-order deconfinement phase transition. We discuss two sets of free parameters, which result in compact-star mass-radius relations that are at tension with the combined constraints for maximum-mass ($2~M_\odot$) and the compactness (GW170817). We find that the most preferable mass-radius relations result in isospin-symmetric phase diagram with rather low temperature for the critical point of the chiral phase transition.


Introduction
For the investigation of matter under extreme conditions and the structure of its phase diagram the equation of state (EoS) is the key target. In the region of finite temperature and vanishing baryon density ab-initio calculations using Monte Carlo simulations of lattice QCD [1] provide a benchmark for developing phenomenological approaches that can be tested, e.g., in heavy-ion collision experiments. Up to now, however, the sign problem prevents the application of lattice QCD methods to the region at low temperature and high baryon density where the possible existence of a first-order phase transition with a critical endpoint has been conjectured. In order to elucidate the QCD phase structure in this domain inaccessible to terrestrial experiments and present techniques of lattice QCD simulations valuable information comes from progress in observing the mass-radius (M-R) relationship of compact stars, due to its one-to-one correspondence with the EoS of compact star matter [2] via the solution of the Tolman-Oppenheimer-Volkoff (TOV) equations [3,4]. For the extraction of the compact star EoS via Bayesian analysis techniques using mass and radius measurements as priors see Refs. [5][6][7]. In particular, in the era of multi-messenger astronomy, it shall soon become possible to constrain the sequence of stable compact star configurations in the mass-radius plane inasmuch that a benchmark for the EoS of cold and dense matter can be deduced from it.
Among the modern observatories for measuring masses and radii of compact stars the gravitational wave interferometers of the LIGO-Virgo Collaboration (LVC) and the X-ray observatory NICER on board the International Space Station provide new powerful constraints besides those from radio pulsar timing.
In this work we pay special attention to the state-of-the-art results from the recent measurement of the high mass 2.17 +0. 11 −0.10 M for PSR J0740+6620 by the NANOGrav Collaboration [8], the compactness derived from the tidal deformability measurement for the binary compact star merger GW170817 [9] in its mass range (1.16 − 1.60 M for the low-spin prior) as well as the still preliminary result of a simultaneous mass and radius measurement by the NICER experiment for PSR J0030+451 [10].
In the study of cold and dense QCD and its applications, commonly used are separate effective models for the nuclear and quark matter phases (two-phase approaches) with a priori assumed first-order phase transition, typically associated with simultaneous chiral and deconfinement transitions. Within this setting, for a constant-speed-of sound model of high-density (quark) matter, a systematic classification of hybrid compact star solutions has been given in [11], which gives a possibility to identify a strong first-order transition in the EoS by the fact that the hybrid star branch in the mass-radius diagram becomes disconnected from the branch of pure neutron stars. However, already before this occurs, a strong phase transition manifests itself by the appearance of an almost horizontal branch on which the hybrid star solutions lie, as opposed to the merely vertical branch of pure neutron stars. In the literature, this strong phase transition has been discussed as due to quark deconfinement [12][13][14]. This conclusion may however be premature since strong phase transitions with a large latent heat occur also within hadronic matter, for instance due to chiral symmetry restoration within the hadronic phase [15]. In the present work, we employ the hybrid quark-meson-nucleon (QMN) model [15][16][17] to explore the implications of dynamical sequential phase transitions at high baryon density on the phenomenology of neutron stars. In order to improve the description of nuclear matter properties at the saturation density, we extend the previous hybrid QMN model by including a six-point scalar interaction. Our main focus is on the role of the chiral symmetry restoration within the hadronic branch of the EoS. This paper is organized as follows. In Sec. 2, we introduce the hybrid quark-meson-nucleon model. In Sec. 3, we discuss the obtained numerical results on the equation of state under neutron-star conditions. In Sec. 4, we discuss the obtained neutron-star relations. In Sec. 5, we present possible realizations of the low-temperature phase diagram. Finally, Sec. 6 is devoted to summary and conclusions.

Hybrid quark-meson-nucleon model
In this section, we briefly introduce the hybrid QMN model for the QCD transitions at finite temperature, density, and arbitrary isospin asymmetry for the application to the physics of neutron stars [15][16][17].
The hybrid QMN model is composed of the baryonic parity doublers [18][19][20] and mesons as in the Walecka model, as well as quark degrees of freedom as in the standard quark-meson model. The spontaneous chiral symmetry breaking yields the mass splitting between the two baryonic parity partners, while it generates an entire mass of a quark. In this work, we consider a system with N f = 2; hence, relevant for this study are the lowest nucleons and their chiral partners, as well as the up and down quarks. The hadronic degrees of freedom are coupled to the chiral fields (σ, π), and the isosinglet vector field ω µ , while the quarks are coupled solely to the former. Another important feature of the hybrid QMN model is that it realizes the concept of statistical confinement by taking into account a medium-dependent modification of the distribution functions, where an additional real scalar field b is introduced.
In the mean-field approximation, the thermodynamic potential of the hybrid QMN model reads [15] Ω = ∑ x where the summation in the first term goes over the up (u) and down (d) quarks, as well as the nucleonic states with positive and negative parity. The positive-parity nucleons correspond to the positively charged 939 1500 140 93 783 775 Table 2. Physical vacuum inputs used in this work.  Table 3. Sets of the model parameters used in this work. The values of λ 4 , λ 6 and g ω are fixed by the nuclear ground state properties, and g ρ by the symmetry energy (see the text). The remaining parameters, g q , κ b , and λ b , do not depend on the choice of m 0 , and their values are taken from Ref. [17]. and neutral N(938) states, i.e., proton (p + ) and neutron (n + ), respectively. The negative-parity nucleons are identified as their counterparts, N(1535) [21], and are denoted as p − and n − . The kinetic term, Ω x , reads The factor γ ± = 2 denotes the spin degeneracy for the nucleons with positive/negative parity, and γ q = 2 × 3 = 6 is the spin-color degeneracy factor for the quarks. The functions n x are the modified Fermi-Dirac distribution functions for the nucleons and for the quarks, accordingly where b is the expectation value of the b-field, and α is a dimensionless model parameter [16,17]. From the definition of n ± and n q , it is evident that, to mimic the statistical confinement, the expected behavior of the b field is to have a nontrivial vacuum expectation value, in order to favor the hadronic degrees of freedom over the quark ones at low densities. On the other hand, it is expected that it vanishes at higher densities in order to suppress the hadronic degrees of freedom and to allow for the population of quarks. This is achieved by allowing b to be generated from a potential V b (see Eq. (10d)). The functions f x andf x are the standard Fermi-Dirac distributions, with β being the inverse temperature, the dispersion relation E x = p 2 + m 2 x . The effective chemical potentials for p ± and n ± are defined as 1 The effective chemical potentials for up and down quarks are given by In Eqs. (6) and (7), µ B , µ Q are the baryon and charge chemical potentials, respectively. The constants g ω , g ρ couple the nucleons to the ω and ρ fields, respectively. The strength of g ω is fixed by the nuclear saturation properties, while the value of g ρ can be fixed by fitting the value of symmetry energy [22]. The properties are shown in Table 1. 1 In the mean-field approximation, the non-vanishing expectation value of the ω field is the time-like component; hence we simply denote it by ω 0 ≡ ω. Similarly, we denote the non-vanishing component of the ρ field, time-like and neutral, by ρ 03 ≡ ρ.  The effective masses of the parity doublers m p ± = m n ± ≡ m ± are given by and for quarks, m u = m d ≡ m q , m q = g q σ.
The parameters g 1 , g 2 , g q are Yukawa-coupling constants, m 0 is the chirally invariant mass of the baryons and is treated as an external parameter (for more details see Ref. [15,17]). The values of those couplings can be determined by fixing the fermion masses in the vacuum (see Table 2). The quark mass is assumed to be m + = 3m q in the vacuum. When the chiral symmetry is restored, the masses of the baryonic parity partners become degenerate with a common finite mass m ± (σ = 0) = m 0 , which reflects the parity doubling structure of the low-lying baryons. This is in contrast to the quarks, which become massless as the chiral symmetry gets restored. The potentials in Eq. (1) are as in the SU(2) linear sigma model, where λ 2 = λ 4 f 2 π − λ 6 f 4 π − m 2 π , and = m 2 π f π . m π , m ω , and m ρ are the π, ω, and ρ meson masses, respectively, The pion decay constant is denoted as f π . Their values are shown in Table 2. The constants κ b and λ b are fixed following Ref. [16]. The parameters λ 4 and λ 6 are fixed by the properties of the nuclear ground state (see Table 1). We note that the introduction of the six-point scalar interaction term in Eq. (10a) is essential in order to reproduce the experimental value of the compressibility K = 240 ± 20 MeV [23].
Following the previous studies of the parity-doublet-based models [15][16][17]24], as well as recent lattice QCD results [25,26], we choose rather large values, m 0 = 700, 800 MeV. The physical inputs and the model parameters used in this work are summarized in Tables 1-3. In-medium profiles of the mean fields are obtained by extremizing the thermodynamic potential in Eq. (1).
In the next section we discuss the the obtained equations of state in the hybrid QMN model and its impact on the chiral phase transition, under the neutron-star conditions of β equilibrium and charge neutrality.  [8]. The blue band is the 2.01 ± 0.04 M constraint [27]. The green and purple bands show 90% credibility regions obtained from the GW170817 event [9] for the low-and high-mass posteriors. The orange region shows the preliminary constraint on the radius from the observation of PSR J0030+0451 [10].

Equation of state under neutron-star conditions
The composition of neutron-star matter requires β equilibrium, as well as the charge neutrality condition. To this end, we include electrons and muons as gases of free relativistic particles.
In Fig. 1 The mixed phases of the chirally broken and restored phases are shown between circles. We stress that the chiral and deconfinement phase transitions are sequential in the current model setup (see Ref. [17]). The latter happen at higher densities and are not shown in the figure.
In all cases, the low-density behavior is similar. In general, for low values of αb 0 , the chiral transition is of first order, determined as a jump in the σ-field expectation value, causing the parity-doublet nucleons to become almost equally populated with mass m ± = m 0 . Higher values of the α parameter yield weaker transitions at higher densities. For αb 0 = 450 MeV, the transition turns into a smooth crossover, defined as a peak in ∂σ/∂µ B . This behavior stays in correspondence to the case of isospin-symmetric matter, where higher value of α weakens the first-order chiral phase transition, which goes through a critical point, and eventually turns into a crossover transition [17]. The values of the net-baryon density range for the mixed phase of the chirally broken and restored phases are shown in Table 4.

TOV solutions for compact star sequences
We use the equations of state introduced in the previous section (see Fig. 1) to solve the general-relativistic TOV equations [3,4] for spherically symmetric objects, dP(r) dr = − ( (r) + P(r)) M(r) + 4πr 3 P(r) r (r − 2M(r)) , with the boundary conditions P(r = R) = 0, M = M(r = R), where R and M are the radius and the mass of a neutron star, respectively. Once the initial conditions are specified based on a given equation of state, namely the central pressure P c and the central energy density c , the internal profile of a neutron star can be calculated. In general, there is one-to-one correspondence between an EoS and the mass-radius relation calculated with it. In Fig. 2, we show the relationship of mass versus radius, for the calculated sequences of compact stars, together with the state-of-the-art constraints on the maximum mass for the pulsar PSR J0348-0432 [27] and PSR J0740+6620 [8]. We point out that the chiral phase transition leads to a softening of the EoS so that it is accompanied by a rapid flattening of the mass-radius sequence. Notably, the chiral transition for all values of αb 0 occurs in the high-mass part of the sequence, but below the 2 M constraint, at around 1.8 M .
In both panels, the three curves for αb 0 = 350, 370, 400 MeV consist of three phases; the chirally broken phase in the low-mass part of the sequence, the chirally restored phase in the high-mass part, and the mixed phase between filled circles. Similarly to the equation of state, increasing the value of α softens the chiral transition, which eventually becomes a smooth crossover for αb 0 = 450 MeV and consists only of branches with chiral symmetry being broken and restored, separated by a circle.
In Fig. 2, we show mass-radius relations obtained for different values of the chirally invariant mass m 0 . What is evident is that increasing the value of m 0 strengthens the chiral phase transition. This is seen twofold, as a shrinking of the mixed phases, as well as more abrupt flattening of chirally restored branches. For a larger m 0 the transition becomes strong enough to produce disconnected branches (see, e.g., the red, solid line in the right panel of Fig. 2). These, in turn, cause the maximal mass of the mass-radius sequences to decrease with increasing value of m 0 . Eventually, the equations of state become not stiff enough to reach the 2 M constraint. We note that such small maximal masses are result of the additional six-point interaction term considered in the thermodynamic potential of the hybrid QMN model (see Eq. (10)). For m 0 = 700 MeV, the most favorable parametrizations are αb 0 = 350 − 370 MeV, while for m 0 = 800 MeV none of the EoSs is stiff enough.
The end points of the mass-radius relation correspond to the onset of quark d.o.f. in each parametrization. This leads to the conclusion that the hadronic matter is not stiff enough to fulfill the two-solar-mass constraint. In general, possible resolution to this problem could be another phase transition. This is the case in the hybrid QMN model, which features sequential chiral and deconfinement phase transitions. However, in the current model setup, the equation of state in the deconfined phase is not stiff enough to sustain the gravitational collapse and the branches become immediately unstable. This is because quarks are not coupled with the vector field leading to a repulsive force. On the other hand, it is known that repulsive interactions tend to stiffen the equation of state. Hence, an additional repulsive force in the quark sector could possibly make the branch stiff enough in order to reach the 2 M constraint, and an additional family of stable hybrid compact stars would appear, with the possibility for the high-mass twin scenario advocated by other effective models [28][29][30]. We note that the obtained mass-radius relations stay in good agreement with the low-mass constraints derived from the recent neutron-star merger GW170817 for the low-and high-mass posteriors [9]. In Fig. 2, they are shown as green and purple regions, respectively. Also, the preliminary constraint on the radius from the NICER experiment is shown as the orange rectangular [10].

Isospin-symmetric phase diagram
The observational neutron-star data provide useful constraints on the structure of strongly interacting matter. Furthermore, they may constrain the phase diagram of isospin-symmetric QCD matter, which is of major relevance for the heavy-ion physics. In Fig. 3, we show the low-temperature part of the isospin-symmetric phase diagram obtained in the model in the (T, ρ B )-plane for m 0 = 700 MeV (left panel) and m 0 = 800 MeV (right panel). In both panels, The first-order liquid-gas phase transition (green, dashed-doubly-dotted line) develops a critical point around T = 16 MeV, and turns into a crossover above it, similarly to the pure parity doublet model. By construction, in the hybrid QMN model, the liquid-gas phase transition is common for all values of the α parameter [16,17]. Similar phase structure is seen in the chiral phase transition for low values of α. For m 0 = 700 MeV, the critical points appear around T = 19, 9 MeV for αb 0 = 350, 370 MeV, respectively. On the other hand, for αb 0 = 400, 450 MeV, the chiral transition proceeds as a smooth crossover at all temperatures. For m 0 = 800 MeV, the critical points are developed at T = 36, 26, 15, 1 MeV for αb 0 = 350, 370, 400, 450 MeV, respectively. Higher values of the temperatures for the critical points are essentially a result of much stronger chiral phase transition at zero temperature. We note that the most favorable parametrizations, i.e., for smaller value of m 0 , yield rather low temperature for the critical point of the chiral phase transition.

Conclusions
In this work, we investigated the hybrid QMN model for the equation of state of dense matter under neutron-star conditions and the phenomenology of compact stars. In particular, we focused on the implications of including six-point scalar interaction and studied the consequences of the realization of the chiral symmetry restoration within the hadronic phase.
We have found that the apparent softening of the EoS results in mass-radius relations with maximal mass at tension with the 2 M constraint within the hadronic branch, especially if the new PSR J0740+6620 with 2.17 M [8] is considered. We have shown that parametrizations of the model which yield large maximal mass (i.e. for smaller value of m 0 ), suggest rather low value of the temperature for the critical end point of the first-order chiral phase transition in the phase diagram, which may also be absent. In view of this, if would interesting to establish a constraint on the chirally invariant mass m 0 . Since the hybrid QMN model features sequential chiral and deconfinement phase transitions, one possible resolution to this could be the onset of quark d.o.f. Such scenario would be even further supported in view of the recent formulation of the three-flavor parity doubling [31,32] and further lattice QCD studies [26], where it was found that to a large extent the phenomenon occurs also in the hyperon channels. In general, the inclusion of heavier flavors is known to soften the equation of state and additional repulsive forces are needed to comply with the 2 M constraint. Additional stiffness from the quark side would play a role, which is not included in the current study. Work in this direction is in progress and the results will be reported elsewhere.