Speed of sound and quark confinement inside neutron stars

Several observations of high-mass neutron stars (NSs), as well as the first historic detection of the binary neutron star merger GW170817, have delivered stringent constraints on the equation of state (EoS) of cold and dense matter. Recent studies suggest that, in order to simultaneously accommodate a 2M⊙ NS and the upper limit on the compactness, the pressure has to swiftly increase with density and the corresponding speed of sound likely exceeds the conformal limit. In this work, we employ a unified description of hadron-quark matter, the hybrid quark-meson-nucleon (QMN) model, to investigate the EoS under NS conditions. We show that the dynamical confining mechanism of the model plays an important role in explaining the observed properties of NSs.


Introduction
The equation of state (EoS) of strongly interacting matter is one of the key observables characterizing properties of matter under extreme conditions. It encodes information on the phase structure of Quantum Chromodynamics (QCD). One way of delineating the EoS is through the speed of sound, c 2 s ≡ dP/d . It has been conjectured that the speed of sound is bounded by c 2 s < 1/3 (conformal bound). This conjecture was in fact confirmed in ab initio calculations of lattice QCD (LQCD) at finite temperature and small net-baryon density [1][2][3][4][5][6]. However, the domain of low temperature and high net-baryon density is still inaccessible via the LQCD methods due to the infamous sign problem. This corner of the QCD phase diagram is of utmost importance for the understanding of the extraterrestrial observations, particularly for the study of neutron stars (NSs), their mergers [7] and supernovae [8]. Until recently, the progress in constraining the EoS at low temperature and high density was driven mostly by several discoveries of high-mass NSs [9][10][11][12]. These observations have set an important constraint on the maximum mass of a NS. More remarkably, the first ever detection of gravitational waves from the compact star merger GW170817 [13], the second detection from GW190425 [14], as well as the NICER observation of the millisecond pulsar PSR J0030+451 [15,16] delivered simultaneous measurement of masses and radii. Perturbative calculations of cold but dense QCD show that the speed of sound complies with the conformal bound [17], although they are reliable only at densities far beyond those realized in astrophysical objects. Contrary to these calculations, as well as the LQCD predictions, several recent analyses provided compelling reasons to expect that the conformal bound has to be violated at densities realized in the interior of NSs in order to support the observed NS properties [18,19].
The advancement of nuclear theory over the years has also tightened the constraints on the EoS over a wide range of densities. This has been achieved by systematic analyses of new astrophysical observations within simplistic approaches, such as the constant-speed-of-sound (CSS) model [20] or multipolytropic class of EoSs [21][22][23]. Recently, the interplay between the high-mass constraint (which requires high pressure) and the upper limit on the compactness from GW170817 event (which favors soft pressure) was used to derive a lower bound constraint on the maximal value of the speed of sound in the cold and dense matter EoS of a NS [24]. This new constraint strengthens previous expectations that the conformal bound is likely to be violated at densities realized inside NSs.
The lower bound of the speed of sound derived in [24] complies with constraints valid at various density regimes. This is particularly useful for determining a class of effective models in which the low-density and high-density regimes are not treated independently, but rather combined in a consistent unified framework. To this end, we employ the hybrid quark-meson-nucleon (QMN) model [25][26][27][28][29] to quantify the EoS of cold and dense matter under NS conditions. 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, before the quark deconfinement takes place. Quark degrees of freedom are included on top of hadrons, but their unphysical onset is prevented at low densities. This is achieved by an auxiliary scalar field which couples to both nucleons and quarks. This field serves as a momentum cutoff in the Fermi-Dirac distribution functions, thus it suppresses the unphysical thermal fluctuations of fermions, with the strength linked to the density. Our main focus is put on the role of the dynamical quark confinement in constraining the EoS of cold and dense matter under NS conditions. This paper is organized as follows. In Section 2, we introduce the hybrid quarkmeson-nucleon model. In Section 3, we discuss the obtained numerical results on the equation of state under neutron-star conditions and neutron-star relations, and we confront them with recent observations and constraints. Finally, Section 4 is devoted to the summary and conclusions.

Hybrid quark-meson-nucleon model
In this section, we briefly introduce the hybrid QMN model for the chiral symmetry restoration and deconfinement phase transitions [25][26][27][28][29]. The hybrid QMN model is composed of the baryonic parity doublet [30][31][32] and mesons as in the Walecka model [33], as well as quark degrees of freedom as in the standard linear sigma model [34]. The spontaneous chiral symmetry breaking yields the mass splitting between the two baryonic parity partners, while it generates the entire mass of a constituent quark. In this work, we consider a system with N f = 2; hence, relevant for this study are the positive-parity nucleons, i.e., proton (p + ) and neutron (n + ), and their negative-parity partners, denoted as p − and n − , as well as the up (u) and down (d) quarks. The fermionic degrees of freedom are coupled to the chiral fields (σ, π), the isosinglet vector-isoscalar field (ω µ ), and the vector-isovector field (ρ µ ). The important concept of statistical confinement is realized in the hybrid QMN model by introducing a medium-dependent modification of the particle distribution functions.
The thermodynamic potential of the hybrid QMN model in the mean-field approximation reads [29] Ω = x=p±,n±,u,d Where the summation goes over the fermionic degrees of freedom. The spin degeneracy factor, γ x for nucleons is γ ± = 2 for both positive-and negative-parity states, while the spin-color degeneracy factor for up and down quarks is γ q = 2 × 3 = 6. The kinetic part, Ω x , reads where the functions n x andn x are the modified Fermi-Dirac distributions for nucleons and for quarks respectively. The model embeds the concept of statistical confinement through the modified Fermi-Dirac distribution functions, where b is the expectation value of an auxiliary scalar field b and α is a dimensionless model parameter. As demonstrated in references [25][26][27][28][29], the parameter α plays also a crucial role in tuning the order of the chiral phase transition. From the definition of n ± and n q , it is evident that, in order to mimic the statistical confinement, the b field should have a nontrivial vacuum expectation value, to suppress quark degrees of freedom at low densities in the confined and to allow for their population at high densities in deconfined phase. From equations (3) and (4), one finds that the nucleons favor large b, whereas the quarks small b. The functions f,f are the standard Fermi-Dirac distribution functions for particle and antiparticle, respectively. β is the inverse temperature, and the dispersion relation E x = p 2 + m 2 x . The effective chemical potentials for p ± and n ± are defined as The European Physical Journal Special Topics Table 1. Properties of the nuclear ground state at µB = 923 MeV and the symmetry energy used in this work.
0.16 −16 240 31 Table 2. Physical vacuum inputs used in this work. The effective chemical potentials for up and down quarks are given by In equations (6) and (7), µ B , µ Q are the baryon and charge chemical potentials, respectively. The strength of g N ω is fixed by the nuclear saturation properties, while the value of g N ρ can be fixed by fitting the value of symmetry energy [35]. The properties of the nuclear ground state and the symmetry energy are shown in Table 1. On the other hand, the nature of the repulsive interaction among quarks and their coupling to the ω and ρ mean fields are still far from consensus. To account for the uncertainty in the theoretical predictions, one may treat the couplings g q ω and g q ρ as free parameters. As demonstrated in reference [29], the repulsive quark-vector interaction has consequences for the phenomenological description of compact stellar objects of masses around 2M . In the current work, however, we are interested in NSs with masses of 1.4M . Thus, we neglect the repulsive quark-vector interactions and set g q ω = g q ρ = 0 for simplicity of discussion.
The effective masses of the chiral partners, m p± = m n± ≡ m ± , are given by The positive-parity nucleons are identified as the positively charged and neutral N (938) states, i.e., proton (p + ) and neutron (n + ). Their negative-parity counterparts, denoted as p − and n − are identified as N (1535) [36]. From equation (8), it is clear that the chiral symmetry breaking generates only the splitting between the two masses. When the chiral symmetry is restored, the masses become degenerate with a common finite mass m ± (σ = 0) = m 0 , which reflects the parity doubling structure of the low-lying baryons. Following the previous studies of the parity-doublet-based models [25][26][27][28][29][37][38][39][40][41][42][43][44][45][46][47], as well as recent lattice QCD results [48,49], we choose a rather large value, m 0 = 700 MeV. The couplings g 1 and g 2 in equation (8) can be determined by fixing the fermion masses in the vacuum. Their values used in this work are summarized in Table 2.
The quark effective mass, m u = m d ≡ m q , is linked to the sigma field as Table 3. Numerical values of the model parameters. The values of λ4, λ6 and g N ω are fixed by the nuclear ground state properties, g N ρ by the symmetry energy, and gq is fixed by the vacuum quark mass (see the text). The remaining parameters, κ b and λ b are fixed following reference [26]. We note that in contrast to the baryonic parity partners (cf. Eq. (8)), quarks become massless as the chiral symmetry gets restored. The value of the coupling g q in equation (9) can be determined by assuming the quark mass to be m q = 1/3 m + in the vacuum. The potentials in equation (1) read where λ 2 = λ 4 f 2 π − λ 6 f 4 π − m 2 π , and = m 2 π f π . m π , m ω , and m ρ are the π, ω, and ρ meson masses, respectively, and f π is the pion decay constant. The parameters λ 4 and λ 6 are fixed by the properties of the nuclear ground state. Numerical values of all model parameters are summarized in Table 3.
In-medium profiles of the mean fields are obtained by extremizing the thermodynamic potential in equation (1). In the grand canonical ensemble, the thermodynamic pressure is obtained from the thermodynamic potential as P = −Ω + Ω 0 , where Ω 0 is the value of the thermodynamic potential in the vacuum. The net-baryon number density for a species x is defined as where Ω x is the kinetic term in equation (2). The total net-baryon number density reads In the following section the above hybrid QMN model equation of state of strongly interacting matter will be applied to identify properties of compact stellar objects such as NSs. The allowed range for the α parameter is αb 0 = 300 − 450 MeV [25,26], where b 0 denotes the vacuum expectation value of the b-field. Following our previous works, we choose four representative values within that interval: αb 0 = 350, 370, 400, 450 MeV, to systematically study the phenomenology of compact stellar objects.

Results
The composition of NS 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 the left panel of Figure 1 Shown EoSs feature chiral phase transitions, defined as a jump in the σ-field expectation value, which causes the parity partners to become almost degenerate with mass m ± = m 0 . The mechanism of statistical confinement introduced in the previous section has a prominent impact on the stiffness of class of EoSs obtained in the model. Namely, higher values of the parameter α yield weaker first-order transitions, triggered at higher densities, which eventually becomes a smooth crossover. As the density increases, the EoSs feature another two sequential transitions. First, associated with the onset of the down quark, and second, associated with the onset of up quark, after which the matter is fully deconfined and comprised solely of quarks. We note that the sequential appearance stays in contrast to the isospin-symmetric case, where quarks deconfine simultaneously, owing to the isospin symmetry [29]. Such separation of the chirally broken and the deconfined phase might indicate the existence of a quarkyonic phase, where the quarks are partly confined to form a Fermi sphere, but the relevant degrees of freedom around the Fermi surface remain the nucleons with the restored chiral symmetry [50][51][52][53][54][55].
We note that the discussed EoSs are in good agreement with the maximum-mass constraint obtained by using a multi-polytrope ansatz for the EoS above the saturation density, shown as yellow-shaded region in the left panel of Figure 1. Interestingly, the chiral phase transitions and the deconfinement of down quark lie within the region. We note that, in principle, the inclusion of the quark-vector repulsive coupling has an impact only on the high density part of the EoS, when compared to the case with vanishing coupling. Namely, it extends the hadronic branch to higher densities, and Table 4. Baryon density ranges of the coexistence phases associated with the chiral restoration (top), onset of down (middle) and up (bottom) quark under the neutron-star conditions, in terms of saturation density units, ρ0, for m0 = 700 MeV and different values of αb0. In the cases were transitions proceed as smooth crossovers, a single value is given. simultaneously, shifts the appearance of the quarks (see [29] for details). The values of the density jumps associated with the chiral symmetry restoration, and consequent onset of down and up quarks featured in the class of EoS obtained in the model are shown in Table. 4. In the right panel of Figure 1, we show the calculated speeds of sound, c 2 s ≡ dP/d , in the units of the speed of light, as a function of the baryon number density, ρ B . For αb 0 = 350, 370, 400 MeV, the coexistence phases of chirally broken and restored phases due to first order phase transitions are seen as regions of vanishing speed of sound. For αb 0 = 450 MeV, the chiral phase transition proceeds as a smooth crossover and is seen as a dip around 3.5ρ 0 . We note that the sequential appearance of down and up quarks are triggered at higher densities and are not shown in the figure. The notable rapid increase of the speed of sound in each curve, before the chiral phase transition takes place, is a result of the stiffening mechanism that arises due to the statistical confinement implemented in the model (cf. Eq. (3)).
We use the obtained EoSs in the mean-field approximation to solve the generalrelativistic Tolman-Oppenheimer-Volkoff (TOV) equations for spherically symmetric objects at zero temperature, in order to obtain the mass-radius relations for NSs. In the left panel of Figure 2, we show the mass-radius relations calculated from the introduced EoSs, together with the state-of-the-art constraints: the high-mass measurement of the PSR J0740+6620 [12], two recent GW170817 [13] and GW190425 [14] events, and the mass-radius constraint obtained for J0030+0451 [16]. The agreement with all of the constraints is good. Notably, the chiral phase transitions (shown as circles) are featured within the mass and radius region accessible by GW190425, at around 1.8M . We note that the presented results are calculated for vanishing quark-vector interactions, i.e., χ = 0. Finite value of the parameter χ, leads to overall improvement of the phenomenological description of compact objects around 2M (see Ref. [29] for details).
In general, there is a tension between the constraints from high-mass measurements and gravitational-wave observations. On the one hand, it is long known that the existence of 2M NSs is supported by EoSs that are stiff enough at high densities above 2ρ 0 . On the other hand, the deformability constraint favors EoSs that are soft around ∼1 − 2ρ 0 . In [24], the interplay between constraints from these two types of measurements was used to derive a lower limit on the speed of sound in a 1.4M NS within a class of simplistic constant-speed-of-sound (CSS) EoSs. The constraint is shown in the right panel of Figure 2 as a function of radius, R 1.4 , of 1.4M NS (yellow-shaded region). The speed of sound monotonically decreases as R 1.4 increases. However, in principle, it is rather unlikely that the speed of sound is independent of density. Thus, if at some density the speed of sound is below the value needed to comply with the maximum-mass constraint, then it may have to be larger than the desired value of the constraint, at other densities. Therefore, the constraint places a lower estimate for the maximum of the speed of sound in dense NS matter. We note  The inner (outer) gray band shows the 68.3% (95.4%) credibility regions for the mass of PSR J0740+6620 [12]. The inner (outer) green and purple bands show 50% (90%) credibility regions obtained from the recent GW170817 [13] and GW190425 [14] events for the low-and high-mass posteriors. The inner (outer) black region corresponds to the mass and radius constraint at 68.2% (95.4%) obtained for PSR J0030+0451 by the group analyzing NICER X-ray data [16]. The circles mark the coexistence of the chirally broken and chirally restored phases. Right panel: The maximal value of the square of the speed of sound, c 2 s , within a 1.4M NS, as a function of its radius, R. The central values of the speed of sound (open symbols) and maximal values within 1.4M and corresponding radii are shown. The yellow-shaded region represents the lower bound for the maximum value of the speed of sound inside a 1.4M neutron star derived in reference [24]. The black horizontal line shows the conformal value c 2 s = 1/3. Results in both panels are obtained for m0 = 700 MeV, and four representative values of the parameter α.
that the constraint also implies that the maximal value of the speed of sound has to exceed the conformal bound, i.e., c 2 s = 1/3 (horizontal, black line). In the figure, we also show the maximal values of c 2 s (filled symbols) within 1.4M obtained for each parametrization in the hybrid QMN model, together with corresponding central values (open symbols). For αb 0 = 350, 370, 400 MeV, the values are not only above the conformal limit, but they also lie above the constraint. Notably, the maximal values of the speed of sound are obtained at densities where the stiffening of the EoS set in. This is seen in the right panel of Figure 1, where the maximal values of c 2 s in 1.4M are shown as filled symbols. On the other hand, the maximal value of the speed of sound for αb 0 = 450 MeV does not exceed the conformal value. In this case, c 2 s rises monotonically even beyond the central density of the 1.4M NS, and the stiffening sets in at higher densities (see right panel of Fig. 1). The obtained radii, central baryon densities, and the maximal speeds of sound of 1.4M neutron stars for all parametrizations are provided in Table 5. Seemingly, sufficient stiffening of the EoS at densities just above the saturation density, realizable in 1.4M NS, is a trademark that is required in order to comply with the constraint from reference [24]. In the hybrid QMN model, it is provided through the dynamical mechanism of confinement which strength in linked to the density. We note that this effect is also featured in other effective approaches, such as the relativistic density-functional models with excluded nucleon volume [56] and the class of CSS hybrid EoSs [57]. However, too extreme stiffening can become at certain tension with the analysis of GW170817 [13]. In principle, this tension could be resolved, e.g., by a strong phase transition occurring in the mass range relevant for GW170817 [58]. We note that because in the hybrid QMN model the chiral phase transition featured around 1.8M , the obtained 1.4M NSs are composed solely of nuclear matter with broken chiral symmetry. This is seen in the right panel of Figure 1, where the central values of the speed of sound (open symbols) of 1.4M NSs are obtained at baryon densities below the chiral phase transition region (plateaux of vanishing speed of sound). Therefore, the inclusion of the statistical confinement has important implications already at densities before the quarks deconfine. This is even more pronounced in a class of models with sequential chiral and deconfinement phase transitions, such as the hybrid QMN model. This may have important phenomenological implications for the study of multi-messenger astronomy and heavy ion collisions (HIC) [59].

Conclusions
The progress in constraining the equation of state (EoS) of cold and dense matter under extreme conditions is driven mostly through modern multi-messengerastronomy observations. Several high-mass neutron star (NS) observations [9][10][11][12] and the first ever detection of the gravitational waves from the compact-star merger GW170817 [13] have delivered powerful constraints on the NS mass-radius profile independently. In fact, there is an apparent tension between the high-mass constraint (which requires high pressure) and the upper limit on the compactness from GW170817 event (which favors soft pressure). Recently, the interplay between them was used to derive a lower bound constraint on the maximal value of the speed of sound in the cold and dense matter EoS of a NS [24]. This new constraint strengthens previous expectations that the conformal bound is likely to be violated at densities realized inside NSs [18,19]. Because such constraint characterizes different density regimes, it is of particular use for a class of effective models in which the low-density and high-density regimes are not treated independently, but rather combined in a consistent unified framework, such as, e.g., [29,60,61].
In this study, we have utilized the hybrid quark-meson-nucleon (QMN) model to quantify the equation of state (EoS) of cold and dense matter. The model unifies the thermodynamics of quark and hadronic degrees of freedom. The interplay between the quark confinement and the chiral symmetry breaking is embedded in a dynamical way into a single unified framework. Within this approach, we have systematically investigated the EoS of cold and dense asymmetric matter under NS conditions. We have constructed the mass-radius relations based on solutions of the Tolman-Oppenheimer-Volkoff (TOV) equations.
We have shown that in order to comply with modern constraints from multimessenger astronomy, a rapid increase of pressure is required at densities inside a 1.4M NS. In the hybrid QMN model, such stiffening is naturally connected to the dynamical mechanism of confinement which strength in linked to the density. This result highlights the fact that the confinement plays a crucial role in the phenomenology of matter under extreme conditions, even at densities smaller than the density at which the system undergoes a hadron-to-quark phase transition. We note that in reference [26], it was shown that the confinement mechanism of the hybrid QMN model leads to strengthening of the chiral phase transition when compared to pure quark or hadronic models without confinement. This motivates further study and possible applications of the model, not only in the context of multi-messenger astronomy, but also in the context of heavy ion collisions (HIC) at finite temperature and density, where a novel signature of chiral symmetry restoration within the dense hadronic sector in dilepton production was recently proposed [59]. Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.