© ESO 2020

1. Introduction

The equation of state (EoS) is one of the key observables characterizing properties of matter under extreme conditions. In the context of strongly interacting matter, the EoS encodes information on the phase structure and the phase diagram of quantum chromodynamics (QCD).

At a finite temperature and low net-baryon density, the EoS, obtained in ab initio calculations within lattice QCD (LQCD), exhibits a smooth crossover from hadronic matter to a quark-gluon plasma, which is linked to the restoration of the chiral symmetry and color deconfinement (Aoki et al. 2006; Bazavov et al. 2019). LQCD provides the state-of-the-art results describing properties of the low baryon-density QCD matter that are also successfully used for the interpretation of data obtained from the relativistic heavy ion collisions (HIC; Ejiri et al. 2006; Karsch & Redlich 2011; Borsanyi et al. 2012; Bazavov et al. 2012, 2014, 2020).

The chiral symmetry restoration is explicitly identified in LQCD calculations through the scaling properties of the chiral condensate or the behavior of observables that are sensitive to chiral criticality, such as fluctuations of conserved charges (Bazavov et al. 2020). Recently, the chiral symmetry restoration was explored via the temperature dependence of parity doublers when approaching the QCD phase boundary. The LQCD findings (Aarts et al. 2015, 2017, 2019) exhibit a clear manifestation of the parity doubling structure for the low-lying baryons around the chiral crossover. The masses of the positive-parity ground states are found to be rather temperature-independent, while the masses of negative-parity states drop substantially when approaching the chiral crossover temperature, T_c. The parity doublet states become almost degenerate with a finite mass in the vicinity of the chiral crossover. Despite the fact that these LQCD results are still not obtained in the physical limit, the observed behavior of parity partners is likely an imprint of the chiral symmetry restoration in the baryonic sector of QCD. It is to be expected that similar results can also appear in the cold and dense nuclear matter when approaching the chiral symmetry restored phase. Indeed, such properties of the chiral partners can be described in the framework of the parity doublet model (De Tar & Kunihiro 1989; Jido et al. 2000, 2001). The model has been applied to hot and dense hadronic matter, neutron stars, as well as the vacuum phenomenology of QCD (Zschiesche et al. 2007; Benic et al. 2015a; Marczenko & Sasaki 2018; Marczenko et al. 2018, 2019; Motornenko et al. 2019, 2020; Mukherjee et al. 2017a,b; Dexheimer et al. 2013; Steinheimer et al. 2011a,b; Weyrich et al. 2015; Sasaki & Mishustin 2010; Yamazaki & Harada 2019; Ishikawa et al. 2019).

Until now, LQCD calculations at nonzero baryon chemical potential suffer from the sign problem, which stymies their application to the region of low temperature and high net-baryon density. This domain of the QCD phase diagram and the corresponding EoS are of particular importance for the understanding of the extraterrestrial experiments, particularly for the study of compact stellar objects, especially neutron stars (NSs), their mergers (Bauswein et al. 2019a), and supernovae (Fischer et al. 2018).

The progress in constraining the EoS at a low temperature and high density was driven mostly due to discoveries of various two-solar mass NSs (Demorest et al. 2010; Antoniadis et al. 2013; Fonseca et al. 2016; Cromartie et al. 2019). The modern observatories for measuring masses and radii of compact objects, the gravitational wave interferometers of the LIGO/Virgo Collaboration (LVC), and the X-ray observatory Neutron star Interior Composition Explorer (NICER) provide new powerful constraints on the neutron-star mass-radius profile. The first ever detection of gravitational waves from the compact star merger GW170817 (Abbott et al. 2018), the second detection from GW190425 (Abbott et al. 2020), as well as the NICER observation of the millisecond pulsar PSR J0030+451 (Riley et al. 2019; Miller et al. 2019) delivered simultaneous measurements of masses and radii.

Recent advances in nuclear theory have 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 (Alford et al. 2013) or multipolytropic class of EoSs (Hebeler et al. 2013; Read et al. 2009; Alvarez-Castillo & Blaschke 2017). Such schemes are instructive but not microscopic approaches. Although they provide interesting heuristic guidance, they cannot replace a realistic model for the EoS.

One of the uncertainties when modeling the EoS of dense nuclear matter is NS composition. Due to presently limited knowledge, their internal structure remains unclear. Different studies indicate that a variety of exotic constituents may appear in the core of NSs (Blaschke et al. 2016), such as hyperons or meson condensates (see, e.g., Glendenning et al. 2017; Oertel et al. 2017). However, the 2 M_⊙ constraint is hard to reconcile when additional degrees of freedom are taken into account, since it is known that additional degrees of freedom usually implies softening of the EoS, and consequently, eventually obtained maximal masses are far below the 2 M_⊙. This leads, for instance, to the famous hyperon puzzle (Bombaci 2017). One of the most prominent solutions to this problem would be a strong phase transition to deconfined quark matter (Baldo et al. 2003; Shahrbaf et al. 2020).

Sufficient stiffness of the EoS at high densities is required in order to obtain massive solutions of the Tolman Oppenheimer Volkoff (TOV) equations (Tolman 1939; Oppenheimer & Volkoff 1939). In phenomenological models, it is typically provided by repulsive interactions mediated by the exchange of vector mesons. The inclusion of repulsive vector interactions in quark matter allows for hybrid stars that fulfill the 2 M_⊙ constraint, as demonstrated in Klähn et al. (2007, 2013), Benic et al. (2015b), and Klähn & Fischer (2015). They are also known to be important for modeling the QCD phase diagram (Kitazawa et al. 2002).

At moderate densities, however, the quark matter EoS should be soft enough such that the transition to the deconfined phase eventually takes place. Therefore, a natural question to ask is whether or not massive NSs, such as the recently discovered PSR J0740+6620 with mass of $2.{14}_{-0.18}^{+0.20}{M}_{\odot }$ (Cromartie et al. 2019), could contain deconfined quark matter in their core. This problem has attracted a lot of attention over the past two decades (see, e.g., Trümper et al. 2004; Özel 2006; Alford et al. 2007; Klähn et al. 2007; Masuda et al. 2013; Lenzi & Lugones 2012; Orsaria et al. 2013, 2014; Fortin et al. 2020), but no decisive observational evidence has yet been provided.

The mechanism of quark confinement and its relation to the chiral symmetry breaking are of major importance in probing the hadron-quark phase transition, although it is nontrivial to embed their interplay into a single effective theory. For this reason, despite its serious shortcomings, the conventional approach is to use 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 (Bastian & Blaschke 2016). Major attempts to go beyond the two-phase approach are, for example, the Polyakov loop-extended Nambu–Jona-Lasinio (PNJL) models (Fukushima 2004, 2008; Ratti et al. 2006; Sasaki et al. 2007a; Roessner et al. 2007), or the Polyakov loop-extended quark-meson (PQM; Schaefer et al. 2007; Skokov et al. 2010, 2011; Mizher et al. 2010; Herbst et al. 2011). While Polyakov loop-based models provide dynamical means to mimic confinement of quarks, their caveat is the lack of hadrons as degrees of freedom at low temperature and/or density. Moreover, while at finite temperature and low density the Polyakov loop expectation value serves as an approximated order parameter for deconfinement, it is highly questionable that it remains so at high density, which is relevant for astrophysical applications.

In this work, we provide the effective dynamical chiral model to quantify the high-density EoS in the presence of color deconfinement and the chiral symmetry restoration.

Instead of using the Polyakov-loop-based model, we employ the hybrid quark-meson-nucleon (QMN) model (Benic et al. 2015a; Marczenko & Sasaki 2018; Marczenko et al. 2018, 2019) to explore the implications of dynamical sequential phase transitions at high baryon density on the phenomenology of neutron stars. The model has the characteristic feature that, upon 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.

The previous implementation of the hybrid QMN model to obtain the nuclear matter EoS leads to the conclusion, that the hadronic branch of the mass-radius relation is at tension with the two-solar-mass constraint (Marczenko et al. 2018, 2019). In order to improve the description of quark matter properties, we extended the previous hybrid QMN model by including an isoscalar-vector and isovector-vector coupling to quarks. We systematically study the role of the repulsive quark-vector coupling on the deconfinement phase transition and phenomenology of compact stars.

The hybrid QMN model naturally embeds the interplay between the quark confinement and the chiral symmetry breaking into a single field-theoretical framework, which makes the phenomena inherently connected, in contrast to the two-phase approach. These key phenomenological ingredients of the model may have importance not only in the context of modeling of the interior of NSs, but also in dynamical simulations of their mergers, as well as supernova collapse and explosions. Previous studies have pointed out the importance of the QCD phase transitions, where the additional latent heat can trigger the local production of neutrinos (Klähn & Fischer 2015; Fischer et al. 2016, 2018, 2020; Klähn et al. 2017; Bauswein et al. 2019a,b; Pereira et al. 2020). This also motivates further study and possible applications of the proposed model in the context of HIC at finite temperature, where a novel signature of chiral symmetry restoration within the dense hadronic sector in dilepton production was recently proposed (Sasaki 2020).

This paper is organized as follows. In Sect. 2, we introduce the hybrid quark-meson-nucleon model. In Sect. 3, we discuss the obtained numerical results on the equation of state under neutron-star conditions. In Sect. 4, we discuss the obtained neutron-star relations and confront the results with recent observations. Finally, Sect. 5 is features our summary and conclusions.

2. Hybrid QMN model

In this section, we briefly introduce the hybrid QMN model, which is capable of describing the chiral symmetry restoration and deconfinement phase transitions (Benic et al. 2015a; Marczenko & Sasaki 2018; Marczenko et al. 2018, 2019). The hybrid QMN model is composed of the baryonic parity doublet (De Tar & Kunihiro 1989; Jido et al. 2000, 2001) and mesons as in the Walecka model (Walecka 1974), as well as quark degrees of freedom as in the standard linear sigma model (Scavenius et al. 2001). 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, which are 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 Lagrangian of the model reads

$$\begin{array}{c}\hfill \mathcal{L}={\mathcal{L}}_N+{\mathcal{L}}_M+{\mathcal{L}}_q,\end{array}$$(1)

with ℒ_N, ℒ_M, ℒ_q denoting the nucleon, meson, and quark parts, respectively. The nucleon part (ℒ_N) of the Lagrangian reads

$$\begin{array}{cc}\hfill {\mathcal{L}}_N& =i{\overline{\psi }}_1\partial {\psi }_1+i{\overline{\psi }}_2\partial {\psi }_2+m_0({\overline{\psi }}_1{\gamma }_5{\psi }_2-{\overline{\psi }}_2{\gamma }_5{\psi }_1)\hfill \\ \hfill & +g_1{\overline{\psi }}_1(\sigma +i{\gamma }_5\mathit{\tau }·\mathit{\pi }){\psi }_1+g_2{\overline{\psi }}_2(\sigma -i{\gamma }_5\mathit{\tau }·\mathit{\pi }){\psi }_2\hfill \\ \hfill & -g_{\omega }^N{\displaystyle \sum _{k=1,2}}{\overline{\psi }}_k/\omega {\psi }_k-\frac{1}{2}{g}_{\rho }^N{\displaystyle \sum _{k=1,2}}{\overline{\psi }}_k\mathit{\tau }·\mathbf{/}\mathit{\rho }{\psi }_k,\hfill \end{array}$$(2)

where ψ_k is the set of baryonic chiral fields. Parameters g₁, g₂, and $g_{\omega }^N$, $g_{\rho }^N$ are the baryon-to-meson coupling constants, and m₀ is a mass parameter.

The mesonic part, ℒ_M, of the Lagrangian is introduced as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{L}}_M& =\frac{1}{2}{\left({\partial }_{\mu }\sigma \right)}^2+\frac{1}{2}{\left({\partial }_{\mu }\mathit{\pi }\right)}^2-\frac{1}{4}{\left({\omega }_{\mu \nu }\right)}^2-\frac{1}{4}{\left({\mathit{\rho }}_{\mu \nu }\right)}^2\hfill \\ \hfill & -V_{\sigma }-V_{\omega }-V_{\rho },\hfill \end{array}\end{array}$$(3)

where ω_μν = ∂_μω_ν − ∂_νω_μ is the field-strength tensor of the vector-isoscalar field, ${\mathit{\rho }}_{\mu \nu }={\partial }_{\mu }{\mathit{\rho }}_{\nu }-{\partial }_{\nu }{\mathit{\rho }}_{\mu }-g_{\rho }^N{\mathit{\rho }}_{\mu }\times {\mathit{\rho }}_{\nu }$ is the field-strength tensor of the vector-isovector field, and the potentials read

$$\begin{array}{cc}\hfill V_{\sigma }& =-\frac{{\lambda }_2}{2}\mathrm{\Sigma }+\frac{{\lambda }_4}{4}{\mathrm{\Sigma }}^2-\frac{{\lambda }_6}{6}{\mathrm{\Sigma }}^3-ϵ\sigma ,\hfill \end{array}$$(4a)

$$\begin{array}{cc}\hfill V_{\omega }& =-\frac{m_{\omega }^2}{2}{\omega }_{\mu }{\omega }^{\mu },\hfill \end{array}$$(4b)

$$\begin{array}{cc}\hfill V_{\rho }& =-\frac{m_{\rho }^2}{2}{\mathit{\rho }}_{\mu }{\mathit{\rho }}^{\mu },\hfill \end{array}$$(4c)

where $\mathrm{\Sigma }={\sigma }^2+{\mathit{\pi }}^2$, ${\lambda }_2={\lambda }_4{f}_{\pi }^2-{\lambda }_6{f}_{\pi }^4-m_{\pi }^2$, and $ϵ=m_{\pi }^2{f}_{\pi }$. m_π, m_ω, and m_ρ are the π, ω, and ρ meson masses, respectively, and f_π is the pion decay constant.

The parameters λ₄ and λ₆ are fixed by the properties of the nuclear ground state.

We note that the six-point scalar interaction term in Eq. (4a) is essential in order to reproduce the experimental value of the compressibility K = 240 ± 20 MeV (Shlomo et al. 2006; Motohiro et al. 2015). Numerical values of model parameters are summarized in Tables 1–3.

In the physical basis, the effective masses of the chiral partners, m_p± = m_n± ≡ m_±, are given by

$$\begin{array}{c}\hfill m_{\pm }=\frac{1}{2}[\sqrt{{(g_1+g_2)}^2{\sigma }^2+4m_0^2}\mp (g_1-g_2)\sigma ].\end{array}$$(5)

The positive-parity nucleons are identified as the positively charged and neutral N(938) states: proton (p₊) and neutron (n₊). Their negative-parity counterparts, denoted as p₋ and n₋, are identified as N(1535) (Tanabashi et al. 2018). From Eq. (5), it is clear that the chiral symmetry breaking generates only the splitting between the two masses.

When the chiral symmetry is restored, the masses of the baryonic parity partners become degenerate with a common finite mass m_±(σ = 0) = m₀, which reflects the parity doubling structure of the low-lying baryons. Following the previous studies of the parity-doublet-based models (Zschiesche et al. 2007; Benic et al. 2015a; Marczenko & Sasaki 2018; Marczenko et al. 2018, 2019; Motornenko et al. 2020, 2019; Mukherjee et al. 2017a,b; Dexheimer et al. 2013; Steinheimer et al. 2011a,b; Weyrich et al. 2015; Sasaki & Mishustin 2010; Yamazaki & Harada 2019; Ishikawa et al. 2019), as well as recent LQCD results (Aarts et al. 2017, 2019), we choose a rather large value: m₀ = 700 MeV. The couplings g₁, g₂ in Eq. (5) can be determined by fixing the fermion masses in the vacuum. Their values used in this work are summarized in Table 1. Having fixed the chirally invariant mass m₀, we obtain the values of the effective masses at saturation density, which are denoted as $m_{\pm }^{\ast }$. We find that the effective mass of the positive-parity states, $m_{+}^{\ast }$, lies within the experimentally determined range (Kapusta & Gale 2011; Glendenning et al. 2017). We note, however, that their exact values depend on the choice of the parameter m₀ (Zschiesche et al. 2007). Thus, in principle, the experimental value of $m_{+}^{\ast }$ could be used as an input parameter to constrain the model parameters. In Table 2, we show the obtained values of the effective masses at saturation density for both parity partners.

The quark part (ℒ_q) of the Lagrangian in Eq. (1) reads

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathcal{L}}_q& ={\displaystyle \sum _{k=u,d}}i{\overline{\psi }}_k\partial {\psi }_k+g_q{\overline{\psi }}_k(\sigma +i{\gamma }_5\mathit{\tau }·\mathit{\pi }){\psi }_k\hfill \\ \hfill & -g_{\omega }^q{\displaystyle \sum _{k=u,d}}{\overline{\psi }}_k/\omega {\psi }_k-\frac{1}{2}{g}_{\rho }^q{\displaystyle \sum _{k=u,d}}{\overline{\psi }}_k\mathit{\tau }·\mathbf{/}\mathit{\rho }{\psi }_k,\hfill \end{array}\end{array}$$(6)

where ψ_k is a set of fields for up (u) and down (d) quarks. Parameters $g_{\sigma }^q$, $g_{\omega }^q$ and $g_{\rho }^q$ are the quark-to-meson coupling constants. The quark effective mass, m_u = m_d ≡ m_q, is linked to the sigma field as

$$\begin{array}{c}\hfill m_q=g_q\sigma .\end{array}$$(7)

We note that in contrast to the baryonic parity partners (cf. Eq. (5)), quarks become massless as the chiral symmetry gets restored. The value of the coupling g_q in Eq. (7) can be determined by assuming the quark mass to be m_q = 1/3 m₊ in the vacuum.

The strength of $g_{\omega }^N$ is fixed by the nuclear saturation properties, while the value of $g_{\rho }^N$ can be fixed by fitting the value of symmetry energy (Glendenning et al. 2017). The properties of the nuclear ground state and the symmetry energy are shown in Table 2. 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, we treat the couplings $g_{\omega }^q$ and $g_{\rho }^q$ as free parameters.

We parametrize the couplings of quarks to the ω and ρ mesons with a single dimensionless parameter χ as follows:

$$\begin{array}{cc}\hfill g_{\omega }^q& =\chi g_{\omega }^N,\hfill \end{array}$$(8a)

$$\begin{array}{cc}\hfill g_{\rho }^q& =\chi g_{\rho }^N.\hfill \end{array}$$(8b)

We stress that the original hybrid QMN model was designed to saturate the Stefan-Boltzmann limit at high temperature and density (Benic et al. 2015a; Marczenko & Sasaki 2018). However, the introduction of an additional quark-vector coupling spoils this limit. On the other hand, in this work we considered vector fields only up to quadratic order (cf. Eq. (4)), thus, the resulting terms in the effective mean-field potential are also up to quadratic order in densities and do not lead to the acausal behavior of the sound velocity, meaning $c_{\mathrm{s}}^2\equiv \mathrm{\partial }P/\mathrm{\partial }ϵ<1$ is fulfilled (Kaltenborn et al. 2017). We note that it is generally possible to sustain the 2 M_⊙ constraint and fulfill the conformal bound, $c_{\mathrm{s}}^2\le 1/3$. This can be obtained, for example, in a class of constant-speed-of-sound equations of state (Alford et al. 2015).

The hybrid QMN model realizes the concept of statistical confinement through a medium-dependent modification of the Fermi-Dirac distribution functions, where an auxiliary scalar field, b (bag field), is introduced. The distribution functions for the quarks and antiquarks are replaced with

$$\begin{array}{cc}\hfill n_q& =\theta ({\mathit{p}}^2-b^2)f_q,\hfill \end{array}$$(9a)

$$\begin{array}{cc}\hfill {\overline{n}}_q& =\theta ({\mathit{p}}^2-b^2){\overline{f}}_q,\hfill \end{array}$$(9b)

where b is the expectation value of the b-field. The distribution functions for the nucleons and antinuclones are replaced with

$$\begin{array}{cc}\hfill n_{\pm }& =\theta ({\alpha }^2{b}^2-{\mathit{p}}^2)f_{\pm },\hfill \end{array}$$(10a)

$$\begin{array}{cc}\hfill {\overline{n}}_{\pm }& =\theta ({\alpha }^2{b}^2-{\mathit{p}}^2){\overline{f}}_{\pm },\hfill \end{array}$$(10b)

where α is a dimensionless model parameter, and $f,\overline{f}$ are the standard Fermi-Dirac distribution functions for the particle and antiparticle, respectively. As we demonstrate in Sect. 3, the parameter α plays a crucial role in tuning the order of the chiral phase transition (Benic et al. 2015a; Marczenko & Sasaki 2018).

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, so as to suppress quark degrees of freedom at low densities in the confined phase and to allow for their population at high densities in the deconfined phase.

This is achieved by allowing b to be generated from the potential V_b. Following Benic et al. (2015a), we take the potential of the minimal form

$$\begin{array}{c}\hfill V_b=-\frac{1}{2}{\kappa }_b^2{b}^2+\frac{1}{4}{\lambda }_b{b}^4.\end{array}$$(11)

The potential V_b for a positive ${\kappa }_b^2$ develops a nontrivial vacuum expectation value at $b_0=\sqrt{{\kappa }_b^2/{\lambda }_b}$. From Eqs. (9) and (10), one finds that the nucleons favor large b, whereas the quarks favor small b. The potential (Eq. (11)) is chosen such that, at a certain T and μ_B, a transition sets in, causing the bag-field expectation value to abruptly drop. As a consequence, at low T and μ_B, the quark degrees of freedom are suppressed, while the nucleons get suppressed at high T and μ_B. This characteristic behavior is associated with the deconfinement transition, which is a crucial feature of the model (Benic et al. 2015a).

In Eqs. (9) and (10), the functions f_x and ${\overline{f}}_x$ are the standard Fermi-Dirac distributions,

$$\begin{array}{cc}\hfill f_x& =\frac{1}{1+e^{\beta (E_x-{\mu }_x)}},\hfill \end{array}$$(12a)

$$\begin{array}{cc}\hfill {\overline{f}}_x& =\frac{1}{1+e^{\beta (E_x+{\mu }_x)}},\hfill \end{array}$$(12b)

with β being the inverse temperature and the dispersion relation $E_x=\sqrt{{\mathit{p}}^2+m_x^2}$. The effective chemical potentials for p_± and n_± are defined as¹

$$\begin{array}{cc}\hfill {\mu }_{p_{\pm }}& ={\mu }_B-g_{\omega }^N\omega -\frac{1}{2}{g}_{\rho }^N\rho +{\mu }_Q,\hfill \end{array}$$(13a)

$$\begin{array}{cc}\hfill {\mu }_{n_{\pm }}& ={\mu }_B-g_{\omega }^N\omega +\frac{1}{2}{g}_{\rho }^N\rho .\hfill \end{array}$$(13b)

The effective chemical potentials for up and down quarks are given by

$$\begin{array}{cc}\hfill {\mu }_u& =\frac{1}{3}{\mu }_B-g_{\omega }^q\omega -\frac{1}{2}{g}_{\rho }^q\rho +\frac{2}{3}{\mu }_Q,\hfill \end{array}$$(14a)

$$\begin{array}{cc}\hfill {\mu }_d& =\frac{1}{3}{\mu }_B-g_{\omega }^q\omega +\frac{1}{2}{g}_{\rho }^q\rho -\frac{1}{3}{\mu }_Q.\hfill \end{array}$$(14b)

In Eqs. (13) and (14), μ_B, μ_Q are the baryon and charge chemical potentials, respectively.

The concept of momentum cutoff is naturally supported by the asymptotic freedom in QCD, which suggests that active degrees of freedom are different depending on their momenta, meaning hadrons (quarks) with low (high) momenta. Such notion has been widely used in effective theories and the Dyson-Schwinger approach (Roberts et al. 2010, 2011). In Benic et al. (2015a), this idea was pushed further by introducing a momentum cut-off as a medium-dependent quantity. Such an intrinsic modification of the cut-off is determined self-consistently when the cut-off is regarded as an expectation value of a scalar field, which is b in our model. As demonstrated in Benic et al. (2015a), the b field plays a role similar to that of the Polyakov loop at finite temperature and vanishing chemical potential, and it is responsible for the onset of quark degrees of freedom around the crossover temperature. Although the b field is so far not anchored to any QCD symmetry, its introduction is a minimal and practical way to realize the physical situation where quarks are not activated deep in the hadronic phase. Therefore, b is likely related to gluons, especially the chromoelectric component. So far it is not known how one is connected to another, since there is no rigorous order parameter of confinement in QCD with light fermions. Thus, the b field cannot be regarded as a true order parameter. In recent studies, it has been shown that such nontrivial order parameter of the center flavor symmetry exists in a QCD-like theory compactified on a circle (Cherman et al. 2017; Aitken et al. 2017).

The parameter α manifests itself only in combination with the expectation value of the b field. Therefore, α relates the scales for the momentum distributions of nucleons (αb) and quarks (b). The lower bound for the α parameter is set so that the ultraviolet cut-off αb₀ for nucleon distribution functions and the infrared cut-off b₀ for quark distribution functions do not spoil the properties of the nuclear ground state. This sets the lower bound to be roughly αb₀ ≥ 300 (Benic et al. 2015a; Marczenko & Sasaki 2018). On the other hand, the upper bound on the parameter can be derived by requiring that the system correctly saturates the Stefan-Boltzmann limit, which yields the maximal value for αb₀ ≤ 450 MeV (Benic et al. 2015a; Marczenko & Sasaki 2018). Following our previous works, we chose four representative values within that interval: αb₀ = 350,  370,  400, and 450 MeV, so as to systematically study the influence of the quark-vector couplings on the equation of state and the resulting mass-radius relations.

Intuitively, the meaning behind the momentum cut-off b and the α parameter can be understood within a simple bag model picture, where the nucleon is a bound state of three noninteracting massless valence quarks, which are confined in a bag of radius r. Due to the uncertainty principle, the quarks cannot have momenta lower than q₀ ≥ 1/r. With the mass of the ground state nucleon being m₊ = 3q₀ = 939 MeV, this yields the value of the radius r = 0.63 fm. When the nuclear matter is dense enough, the nucleon wave functions start overlapping and the system exhibits the Pauli blocking effect. In order to estimate the density at which the Pauli blocking manifests itself, one can use the van der Waals excluded volume of hard-sphere nucleons (Andronic et al. 2012):

$$\begin{array}{c}\hfill v_{\mathrm{ex}}=\frac{16}{3}\pi r^3,\end{array}$$(15)

which corresponds to the closest packing density of ${\rho }_{\mathrm{c}}=1/{\mathit{v}}_{\mathrm{ex}}=0.238{\mathrm{fm}}^{-3}=1.49{\rho }_0$, where ρ₀ = 0.16 fm⁻³ is the saturation density. The density ρ_c naturally corresponds to the onset of the stiffening effect in the hybrid QMN model. At higher density, the wave functions of nucleons overlap and the quark momenta have to be rearranged, in order not to violate the Pauli exclusion principle. A simple argument how the Pauli principle may act in a dense gas of nucleons as three-quark bound states is that the higher quark momenta have to be populated. At zero temperature, the shell in the quark-momentum space, Δρ_q, which has to be filled, ranges from q₀ to q₀ + δq₀ and corresponds to the quark density

$$\begin{array}{c}\hfill \mathrm{\Delta }{\rho }_q=3\mathrm{\Delta }{\rho }_B=3{\gamma }_q{\displaystyle \underset{q_0}{\overset{q_0+\delta q_0}{\int }}}\mathrm{d}q\frac{q^2}{2{\pi }^2}=\frac{{\gamma }_q}{2{\pi }^2}{q}_0^2\delta q_0+\mathcal{O}(\delta q_0^2),\end{array}$$(16)

where γ_q is the spin-color-flavor degeneracy factor of quarks. From the above, the shift of the quark momentum, δq₀, can be related to the corresponding Pauli shift in the nucleon energy, ${\mathrm{\Delta }}^{\mathrm{P}\mathrm{a}\mathrm{u}\mathrm{l}\mathrm{i}}=3\delta q_0=18{\pi }^2/({\gamma }_q{q}_0^2)\mathrm{\Delta }{\rho }_B$. At densities exceeding the critical density, ρ_c, where the effect sets in, the Pauli shift corresponds to a density-dependent lowering of the chemical potential. The critical density can be translated to the critical Fermi momentum of nucleons,

$$\begin{array}{c}\hfill p_f=\sqrt[3]{\frac{3}{2}{\pi }^2{\rho }_c}=\sqrt[3]{\frac{9}{32}\pi }{q}_0,\end{array}$$(17)

and compared with the corresponding UV cut-off αb₀ in the hybrid QMN model. The value q₀ is identified with the IR cut-off in the quark momentum distribution due to confinement and the uncertainty relation. In the hybrid QMN model, it is given by b₀. Thus, one gets the UV cut-off in the nucleon distribution to be p_f = αb₀. From this, one gets $\alpha =\sqrt[3]{9\pi /32}\approx 0.959$. We note that the maximal value of the α parameter obtained in the hybrid QMN model is α_max = 2^−1/3 ≈ 0.794 (Benic et al. 2015a; Marczenko & Sasaki 2018). This simple analysis serves as a ballpark estimate for the parameter α and shows that it may be quantified by considering the Pauli exclusion principle on the quark level. The present estimate within an intuitive bag model picture and the van der Waals excluded volume leads to similar scaling as the one obtained in the hybrid QMN model, meaning smaller values of b₀ and q₀ lead to an earlier onset of the stiffening effect. A more realistic calculation of the quark Pauli blocking in nuclear matter can be found in Blaschke et al. (2020).

The thermodynamic potential in the mean-field approximation reads (Marczenko et al. 2018)

$$\begin{array}{c}\hfill \mathrm{\Omega }={\displaystyle \sum _{x=p_{\pm },n_{\pm },u,d}}{\mathrm{\Omega }}_x+V_{\sigma }+V_{\omega }+V_{\rho }+V_b,\end{array}$$(18)

where the summation goes over the fermionic degrees of freedom. The kinetic part, Ω_x, of the thermodynamic potential Eq. (18), reads

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_x={\gamma }_x{\displaystyle \int \frac{{\mathrm{d}}^{3}p}{{\left(2\pi \right)}^3}}T[ln(1-n_x)+ln(1-{\overline{n}}_x)],\end{array}$$(19)

where the functions n_x and ${\overline{n}}_x$ are the modified Fermi-Dirac distributions defined in Eqs. (9) and (10). The spin degeneracy factor, γ_x for nucleons is γ_± = 2 for both positive- and negative-parity states, while the spin-color degeneracy factor for quarks is γ_q = 2 × 3 = 6.

The physical inputs and the model parameters used in this work are summarized in Tables 1 and 2. In-medium profiles of the mean fields are obtained by extremizing the thermodynamic potential in Eq. (18), leading to the following gap equations:

$$\begin{array}{cc}\hfill \frac{\partial \mathrm{\Omega }}{\partial \sigma }& =\frac{\partial V_{\sigma }}{\partial \sigma }+{\displaystyle \sum _{x=p_{\pm },n_{\pm },u,d}}{s}_x\frac{\partial m_x}{\partial \sigma },\hfill \end{array}$$(20a)

$$\begin{array}{cc}\hfill \frac{\partial \mathrm{\Omega }}{\partial \omega }& =\frac{\partial V_{\omega }}{\partial \omega }+{\displaystyle \sum _{x=p_{\pm },n_{\pm },u,d}}{g}_{\omega }^x{\rho }_x,\hfill \end{array}$$(20b)

$$\begin{array}{cc}\hfill \frac{\partial \mathrm{\Omega }}{\partial \rho }& =\frac{\partial V_{\rho }}{\partial \rho }+{\displaystyle \sum _{x=p_{\pm },u}}\frac{g_{\rho }^x}{2}{\rho }_x-{\displaystyle \sum _{x=n_{\pm },d}}\frac{g_{\rho }^x}{2}{\rho }_x,\hfill \end{array}$$(20c)

$$\begin{array}{cc}\hfill \frac{\partial \mathrm{\Omega }}{\partial b}& =\frac{\partial V_b}{\partial b}+\alpha {\displaystyle \sum _{x=p_{\pm },n_{\pm }}}{\omega }_x-{\displaystyle \sum _{x=u,d}}{\omega }_x,\hfill \end{array}$$(20d)

where the scalar and baryon densities are

$$\begin{array}{c}\hfill s_x={\gamma }_x{\displaystyle \int \frac{{\mathrm{d}}^{3}p}{{(2\pi )}^3}}\frac{m_x}{E_x}(n_x+{\overline{n}}_x)\end{array}$$(21)

and

$$\begin{array}{c}\hfill {\rho }_x={\gamma }_x{\displaystyle \int \frac{{\mathrm{d}}^{3}p}{{(2\pi )}^3}}(n_x-{\overline{n}}_x),\end{array}$$(22)

respectively. The terms ω_x in the gap Eq. (20d) are given as

$$\begin{array}{c}\hfill {\omega }_{\pm }={\gamma }_{\pm }\frac{{(\alpha b)}^2}{2{\pi }^2}T{[ln(1-f_{\pm })+ln(1-{\overline{f}}_{\pm })]}_{{\mathit{p}}^2={(\alpha b)}^2}\end{array}$$(23)

and

$$\begin{array}{c}\hfill {\omega }_q={\gamma }_q\frac{b^2}{2{\pi }^2}T{[ln(1-f_q)+ln(1-{\overline{f}}_q)]}_{{\mathit{p}}^2=b^2}\end{array}$$(24)

for the nucleons and quarks, respectively. We note that the terms in Eqs. (23) and (24) come into the gap Eq. (20d) with opposite signs. This reflects the fact that nucleons and quarks favor different values of the bag field.

In the grand canonical ensemble, the thermodynamic pressure is obtained from the thermodynamic potential as P = −Ω + Ω₀, where Ω₀ is the value of the thermodynamic potential in the vacuum. The net-baryon number density for a species x is defined as

$$\begin{array}{c}\hfill {\rho }_B^x=-\frac{\partial {\mathrm{\Omega }}_x}{\partial {\mu }_B},\end{array}$$(25)

where Ω_x is the kinetic term in Eq. (19). The total net-baryon number density reads

$$\begin{array}{c}\hfill {\rho }_B={\rho }_B^{n_{+}}+{\rho }_B^{n_{-}}+{\rho }_B^{p_{+}}+{\rho }_B^{p_{-}}+{\rho }_B^u+{\rho }_B^d.\end{array}$$(26)

The particle-density fractions are defined as

$$\begin{array}{c}\hfill Y_x=\frac{{\rho }_B^x}{{\rho }_B}.\end{array}$$(27)

The symmetry energy and its slope at saturation density are given as

$$\begin{array}{c}\hfill E_{\mathrm{sym}}=\frac{1}{2}\frac{{\partial }^2(ϵ/{\rho }_B)}{\partial {\delta }^2}{|}_{\delta =0}\end{array}$$(28)

and

$$\begin{array}{c}\hfill L=3{\rho }_0\frac{\partial E_{\mathrm{sym}}}{\partial {\rho }_B}{|}_{{\rho }_B={\rho }_0},\end{array}$$(29)

respectively. In Eq. (28), ϵ is the energy density and the asymmetry parameter $\delta =({\rho }_B^{n^{+}}+{\rho }_B^{n^{-}}-{\rho }_B^{p^{+}}-{\rho }_B^{p^{-}}+{\rho }_B^d-{\rho }_B^u)/{\rho }_B$². The value of the slope of the symmetry energy at saturation density obtained in the model is L = 82 MeV. The value is rather large when compared to the recent chiral effective field theory estimate (Krüger et al. 2013). We note that similar values are found in other parity-doublet models at the mean-field-level thermodynamics (Motohiro et al. 2015), and it is consistent with the commonly considered range for the parameter (Oertel et al. 2017; Li & Han 2013).

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.

3. Equation of state

We started by generalizing the hybrid QMN EoS derived for the cold nuclear matter to account for the isospin effects, considering isospin-symmetric and highly asymmetric environments such as the interiors of NSs.

The composition of neutron-star matter requires β equilibrium with electrons (e) and muons (μ), included as free relativistic particles, as well as the charge neutrality condition. To this end, we solve the set of equations,

$$\begin{array}{cc}& {\mu }_e={\mu }_{\mu }={\mu }_{n_{+}}+{\mu }_{n_{-}}-{\mu }_{p_{+}}-{\mu }_{p_{-}}+{\mu }_d-{\mu }_u\hfill \end{array}$$(30a)

$$\begin{array}{cc}& {\rho }_{p_{+}}+{\rho }_{p_{-}}+\frac{2}{3}{\rho }_u-\frac{1}{3}{\rho }_d-{\rho }_e-{\rho }_{\mu }=0,\hfill \end{array}$$(30b)

at T = 0, where μ_x’s and ρ_x’s are the effective chemical potentials of given species (see Eqs. (13) and (14)) and their densities (see Eq. (22)), respectively.

In Fig. 1, we show the calculated zero-temperature equations of state in the mean-field approximation for selected parametrization, m₀ = 700 MeV and αb₀ = 370 MeV, for different strengths of the quark vector-interaction χ (see Eq. (8)), namely χ = 0 (red, solid line), χ = 0.1 (purple, dashed line), and χ = 0.2 (blue, dotted line), as functions of the baryon density in the units of saturation density. In the left panel of the figure, we show the isospin-symmetric matter EoSs. In each case, the behavior at low density is similar. Shown EoSs feature a common first-order chiral phase transition, determined as a jump in the σ-field expectation value, which causes the parity-doublet nucleons to become almost degenerate with the mass m_± = m₀. The chiral phase transition is triggered at roughly 3.25 ρ₀, with the mixed phase persisting up to 3.68 ρ₀. The mechanism of statistical confinement has a prominent twofold impact on the class of EoSs obtained in the model. First, the strength of the chiral phase transition depends on the choice of α. Namely, higher values yield weaker first-order transition, which turns into a second-order and eventually becomes a smooth crossover, defined as a peak in ∂σ/∂μ_B. Second, lower values of α result in a stiffening of the EoS, due to which, the chiral phase transition is triggered at lower densities (Marczenko & Sasaki 2018). We illustrate the stiffening of the EoS based on a simplified model in Appendix A.

Fig. 1. Thermodynamic pressure, P, for isospin-symmetric matter (left panel) and under the NS conditions of β-equilibrium and charge neutrality (right panel), as functions of the net-baryon number density, ρB, in units of the saturation density, ρ0 = 0.16 fm−3 for m0 = 700 MeV, αb0 = 370 MeV and different values of the parameter χ (see text for details). Left panel: the flow constraint (Danielewicz et al. 2002) is shown as the gray shaded region. In both panels, the inset plots show pressure for χ = 0 and different values of the parameter α, in the vicinity of the constraints obtained by Danielewicz et al. (2002) (left panel) and Hebeler et al. (2013) (right panel). We note that the inset figure in the right panel shows pressure as a function of energy density. In both panels, the first-order phase transitions are seen as plateaux of constant pressure.

We note that the stiffening effect of nuclear matter just above the saturation density is a feature that the present approach has in common with, for example, the relativistic density-functional models of nuclear matter with excluded nucleon volume (such as the one by Typel 2016). These are used, for example, in the effective multi-polytrope (Alvarez-Castillo & Blaschke 2017) and CSS (Alford & Sedrakian 2017) class of hybrid EoS. However, too extreme stiffening could represent a certain tension with the recent analysis of GW170817 (Abbott et al. 2018). In principle, this tension could be resolved, for example, by a strong phase transition occurring in the compact star mass range relevant for GW170817 (see Paschalidis et al. 2018).

In the left panel of Fig. 1, we also show the proton flow constraint (gray region) obtained in heavy ion collisions (Danielewicz et al. 2002), where the matter is assumed to be isospin-symmetric. The agreement with the EoSs obtained in the hybrid QMN model is fairly good. Shown EoSs go through the upper band of the constraint. Interestingly, the chiral phase transition happens within the baryon-density range of the flow constraint. The inset plot in the left panel of Fig. 1 shows the EoSs obtained for different values of the α parameter in the vicinity of the flow constraint. Equations of state for αb₀ = 350,  370, and 400 MeV feature the first-order chiral phase transition within this constraint. The softest EoS for αb₀ = 450 MeV features a smooth crossover transition slightly above the densities governed by the flow constraint, at roughly 5.4 ρ₀. The stiffening of the EoSs is evidently seen before the chiral phase transition takes place. This highlights the importance of the chiral phase transition for the density range accessible in heavy ion collisions. We note that the chiral symmetry restoration in the dense hadronic sector can possibly be identified in the production rate of dileptons, as proposed recently in Sasaki (2020).

As the density increases, the equation of state with no quark-vector coupling, that is, χ = 0, features another jump in baryon density, associated with the deconfinement of quarks, with a mixed phase from 9.65 − 16.04 ρ₀. When the inclusion of the quark-vector coupling is considered, meaning the χ = 0.1 and χ = 0.2 cases in Fig. 1, the onset of both quarks is systematically shifted toward higher densities, when compared to the case with vanishing coupling, with mixed phases between 11.19 − 17.28 ρ₀ and 19.18 − 24.83 ρ₀, respectively. Consequently, the hadronic part of the EoS is extended beyond the deconfinement in the case of vanishing χ. We note that all three equations of state remain the same up to the point where the quarks start to appear. We stress that the chiral and hadron-to-quark phase transitions are sequential.

In the right panel of Fig. 1, we show the corresponding EoSs under the neutron-star conditions for m₀ = 700 MeV and αb₀ = 370 MeV. Similarly to the isospin-symmetric case, the behavior at low densities is quite comparable. Shown EoSs feature a common first-order chiral phase transition. For comparison, the EoSs obtained for the remaining values of the parameter α are also shown in the vicinity of the chiral phase transitions. We note that the chiral phase transitions occur at systematically lower densities in comparison to the isospin-symmetric matter. However, as the density increases, the EoSs feature another two sequential jumps in baryon density. The first is associated with the onset of the down quark, and the second is associated with the onset of the up quark, after which the EoS is composed solely of quarks. The sequential appearance of quarks in the isospin-asymmetric matter stays in contrast to the case of symmetric case, where quarks deconfine simultaneously, owing to the isospin symmetry. This can be traced back to their different electric charge, as well as to the coupling to the ρ meson, which introduces different signs for both quarks in their effective chemical potentials (see Eq. (14)). Thus, the obtained EoSs under the neutron-star conditions feature sequential (i.e., flavor-dependent) deconfinement phase transition. We stress that the sequential deconfinement of quarks under the neutron-star conditions cannot be obtained in a class of pure quark models or even in a class of hadron-quark models, where the confined and deconfined phases are treated independently. This is due to the negative charge of the down quark, which cannot be neutralized with negatively-charged leptons. The effect of the finite quark-vector coupling is similar to the isospin-symmetric case. Namely, the onset of both quarks is systematically shifted toward higher densities when compared to the case with vanishing coupling. Consequent extension of the hadronic branch of the EoS is exhibited. In this case, the EoSs remain the same up to the point where the down quark appears. This behavior is also resembled by the density fractions shown in Fig. 2. The composition of the compact-star matter for m₀ = 700 MeV and αb₀ = 370 MeV and χ = 0 (top panel), χ = 0.1 (middle panel), χ = 0.2 (bottom panel), in terms of the particle-density fractions, Y_x, defined in Eq. (27), as a function of the baryon density in the units of saturation density is also illustrated in Fig. 2. Similarly to the EoS, the compositions at low density are the same, up to the point of the appearance of the down quark in the case with χ = 0 (green, dash-dotted line in the top panel). With a finite value of the parameter χ, the hadronic part of the EoS is extended, and consequently the deconfinement of quarks is shifted to higher densities. Because the down quark always appears before the up quark, due to lower Fermi momenta, in order to fulfill the charge neutrality constraint, its negative electric charge can only be neutralized by the positively-charged hadrons, meaning protons, p₊, and its chiral partners, p₋. Such a mixed phase, comprised of hadrons and the down quark, persists until the hadrons become suppressed and the up quark is populated, compensating for the positively-charged hadrons in order to meet the charge neutrality condition. We note that in the case of isospin-asymmetric matter, the phase transitions are triggered at lower densities. The baryon-density jumps associated with consecutive sequential transitions for the sets of model parameters used in this work are shown in Table 4.

Fig. 2. Composition of compact star matter as a function of the central baryon density in the units of saturation density for the present model for m0 = 700 MeV and αb0 = 370 MeV. Three different cases are considered for χ = 0 (top panel), χ = 0.1 (middle panel), χ = 0.2 (bottom panel). In each panel, the yellow shaded areas indicate the density regions of phase coexistence. The gray, solid vertical lines indicate the central densities of the associated maximum-mass neutron stars.

In Fig. 4, we show the speed of sound squared, $c_s^2=\mathrm{d}P/\mathrm{d}ϵ$, in units of the speed of light squared, as a function of the net-baryon number density. The coexistence phases due to first-order phase transitions are seen as plateaux of vanishing speed of sound. The notable swift increases of the speed of sound in each curve are results of the stiffening mechanism that arises due to the statistical confinement implemented in the model (cf. Eq. (10)).

We note that there is a theoretical uncertainty in the prediction for the ratio of the vector- and scalar-quark couplings in the class of Nambu–Jona-Lasinio (NJL) and the Polyakov loop-extended NJL models (Sasaki et al. 2007b; Fukushima 2008; Bratovic et al. 2013; Orsaria et al. 2013, 2014; Klähn et al. 2013). Interestingly, we find that the results are essentially independent of the exact value of the scalar coupling, g_q in the hybrid QMN model. This is because the value of the coupling is fixed to the properties in the vacuum, where quarks are not expected to be present as effective degrees of freedom. In the hybrid QMN model, this is controlled by the finite expectation value of the b-field, which suppresses the quark contribution to the gap equations at low density (see Sect. 2). On the other hand, the quarks are populated at high densities. There, however, the expectation value of the scalar mean field, σ, is of the order of a few MeV, and slowly vanishes as the density increases. Therefore, the quark contribution to the gap equation in the σ direction is negligible. Thus, contrary to other effective models, the laxity in the value of the coupling, g_q, is a direct consequence of the confinement mechanism implemented in the model, and seems natural from the phenomenological point of view. We note that at high temperatures it might not persist, owing to additional thermal fluctuations.

We stress that all EoSs discussed in this work fall into the region derived from the maximum mass constraint M_max ≥ 2.01(4) M_⊙ (Antoniadis et al. 2013) obtained by Hebeler et al. in Hebeler et al. (2013) using a multi-polytrope ansatz for the EoS above the saturation density. This is shown in Fig. 3 for the case of vanishing parameter χ. Interestingly, the chiral phase transition and the deconfinement of down quark lie within the region set by the constraint. This is in contrast to the symmetric-matter flow constraint, where only the chiral transition is featured within the constraint. The up quarks are populated only at higher values of the energy density. We note that for high enough value of the parameter χ, the deconfinement transition would shift to higher density and eventually fall out of the constraint.

Fig. 3. Thermodynamic pressure, P, under the NS conditions of β-equilibrium and charge neutrality, as a function of the energy density, ϵ, for m0 = 700 MeV and χ = 0, and four representative values of the parameter α (see text for details). The phase transitions are seen as plateaux of constant pressure. The gray shaded region marks the constraint obtained by Hebeler et al. (2013).

Fig. 4. Speed of sound squared as a function of the net-baryon number density, ρB, in units of the saturation density, ρ0 = 0.16 fm−3 for m0 = 700 MeV, αb0 = 370 MeV and different values of the parameter χ (see text for details). The first order phase transitions are seen as plateaux of vanishing speed of sound.

In general, the class of EoSs obtained in the hybrid QMN model for the case of isospin-asymmetric matter consist of four consecutive density regions, from low to high density: confined and chirally broken phase, confined and chirally restored phase, partially deconfined, and fully deconfined phase. The intermediate phases are for the novel states in which confined and deconfined matter may coexist. Such remarkable 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 surface, but the relevant degrees of freedom remain the nucleons with the restored chiral symmetry (Hidaka et al. 2008; McLerran et al. 2009; Andronic et al. 2010; McLerran & Reddy 2019; Jeong et al. 2020). Furthermore, in Blaschke et al. (2009), it was conjectured that a mixed phase, comprised of nucleons and down quarks, may have compelling nontrivial cooling and transport properties that are of interest for the thermal and rotational evolution of compact stars.

4. Properties of compact stars

In this section, we explore the physics of dynamical sequential phase transitions at high baryon density on the structure of neutron stars.

4.1. TOV solutions

We used the EoSs introduced in the previous section to solve the general-relativistic Tolman-Oppenheimer-Volkoff (TOV) equations (Tolman 1939; Oppenheimer & Volkoff 1939) for spherically symmetric objects,

$$\begin{array}{cc}\hfill \frac{\mathrm{d}P(r)}{\mathrm{d}r}& =-\frac{(ϵ(r)+P(r))(M(r)+4\pi r^{3}P(r))}{r(r-2M(r))},\hfill \end{array}$$(31a)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}M(r)}{\mathrm{d}r}& =4\pi r^2ϵ(r),\hfill \end{array}$$(31b)

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 the EoS and the mass-radius relation calculated from Eq. (31). In Fig. 5, we show the relationship of mass versus radius, for the calculated sequences of compact stars for m₀ = 700 MeV, for αb₀ = 350 MeV (upper left), αb₀ = 370 MeV (upper right), αb₀ = 400 MeV (lower left), and αb₀ = 450 MeV (lower left). In each panel, three calculated sequences for χ = 0 (red: solid line), χ = 0.1 (purple: dashed line), and χ = 0.2 (blue: dotted line) are shown. Also shown are the state-of-the-art constraints: the high-precision mass measurement of the high-mass pulsar PSR J0740+6620 (Cromartie et al. 2019), constraints from two recent GW170817 (Abbott et al. 2018) events, and the mass-radius constraint obtained for PSR J0030+0451 by the group analyzing NICER X-ray data (Miller et al. 2019). We note that the mass-radius relations obtained in the hybrid QMN model remain in good agreement with the aforementioned constraints. In each curve the onsets of down and up quarks are marked with diamonds and pentagons, respectively.

Fig. 5. Mass-radius sequences for compact stars as solutions of the TOV equations for m0 = 700 MeV, for αb0 = 350 MeV (upper left), αb0 = 350 MeV (upper right), αb0 = 400 MeV (lower left), αb0 = 450 MeV (lower right). The circles show the coexistence of the chirally broken and chirally restored phases. The onsets of up and down quarks are marked by pentagons and diamonds, respectively. We note that, in the upper-left figure, for χ = 0.2, the onset of up and down quarks is almost simultaneous. The inner (outer) gray band shows the 68.3% (95.4%) credibility regions for the mass of PSR J0740+6620 (Cromartie et al. 2019). The inner (outer) green and purple bands show 50% (90%) credibility regions obtained from the recent GW170817 (Abbott et al. 2018) event 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 (Miller et al. 2019).

Because the introduction of the quark-vector repulsive coupling shifts the onset of quarks to higher densities, when compared to the case with vanishing coupling (see Fig. 1), all three sequences in each of the panels share a chiral phase transition in the high-mass part (red circle), but still below the 2 M_⊙ constraint, at around 1.8 M_⊙. We note that increasing the value of αb₀ softens the chiral transition, which eventually becomes a smooth crossover for αb₀ = 450 MeV. We emphasize that an abrupt change in a mass-radius profile in the high-mass part of the sequence is, in general, a result of a phase transition. It leads to a softening of the EoS, meaning that it is accompanied by a rapid flattening of the mass-radius sequence. This is vividly seen in the case of the chiral phase transition on the mass-radius curves in all panels of Fig. 5. In a model with sequential transitions, a rapid flattening of the mass-radius curves is not necessarily associated with the deconfinement transition, and hence does not dictate the existence of quark matter in the core of a neutron star, in contrast to previous studies (Alvarez-Castillo et al. 2016). Interestingly, the chiral phase transition occurs at masses and radii accessible in the GW190425 merger event (Abbott et al. 2020). We note that this is similar to the proton flow constraint from HIC in the isospin-symmetric matter and the Hebeler constraint (cf. inset plots in Fig. 1).

The extension of the hadronic branch of an EoS due to a finite value of the quark vector-interaction χ is also reflected in the corresponding mass-radius sequence. Interestingly, the maximal mass is always reached within the hadronic branch of the sequence. For χ = 0, this happens just before the density jump associated with the onset of down quark is reached. This is the case, except for αb₀ = 370 and 400 MeV, where tiny fractions of down quarks appear just before the maximal mass is reached (see the top panel of Fig. 2). The appearance of the down quark makes the matter too soft to sustain from the gravitational collapse. For χ = 0.1 and 0.2, the hadronic branch extends beyond the density at which the maximal mass is reached and becomes gravitationally unstable. This is the case for all values of αb₀. Eventually, when down and up quark are sequentially populated, the matter is still not stiff enough to sustain from the collapse and turn into an additional family of stable hybrid compact stars. This is shown in Fig. 2, for m₀ = 700 MeV and αb₀ = 370 MeV, where for different values of χ, the densities at which the maximal mass is reached are marked by gray solid vertical lines. Clearly, the density for χ = 0.1 (middle panel) is greater than for χ = 0 (top panel), and the hadronic branch extends beyond it. For χ = 0.2, the hadronic branch extends even further, however the maximal mass stays the same. Thus, we conclude that a further increase of the quark-vector coupling does not support the maximal-mass constraint. In Table 5, we show the values of the maximal masses of a neutron star and corresponding radii obtained in each parametrization.

We stress that the above structure persists to large extent when the chiral invariant mass m₀ changes. In order for the pure quark matter to be present in the cores of neutron stars, the repulsive interactions among them should be strong. On the other hand, the self-consistent treatment of the interactions in the hybrid QMN model also shifts the onset of the quark matter in the stellar sequence to higher densities. This eventually prevents the neutron stars with quark matter cores from existence in the gravitationally stable branch of the sequence. This leads to a general conclusion that the existence of hybrid stars is rather excluded in the current hybrid QMN model.

4.2. Tidal deformability

The dimensionless tidal deformability parameter Λ can be computed through its relation to the Love number k₂ (Hinderer 2008; Damour & Nagar 2009; Binnington & Poisson 2009; Yagi & Yunes 2013; Hinderer et al. 2010),

$$\begin{array}{c}\hfill \mathrm{\Lambda }=\frac{2}{3}{k}_2{C}^{-5},\end{array}$$(32)

where C = M/R is the star compactness parameter, with M and R being the total mass and radius of a star. The Love number k₂ reads

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill k_2& =\frac{8C^5}{5}{(1-2C)}^2[2+2C(y-1)-y]\hfill \\ \hfill & \times (2C[6-3y+3C(5\mathit{y}-8)]\hfill \\ \hfill & +4C^3[13-11\mathit{y}+C(3y-2)+2C^2(1+y)]\hfill \\ \hfill & +3{(1-2C)}^2[2-y+2C(y-1)ln(1-2C)]{)}^{-1},\hfill \end{array}\end{array}$$(33)

where y = Rβ(R)/H(R). The functions H(r) and β(r) are given by the following set of differential equations:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\beta }{\mathrm{d}r}& =2{(1-2\frac{M(r)}{r})}^{-1}\hfill \\ \hfill & H\{-2\pi [5ϵ(r)+9P(r)+\frac{\mathrm{d}ϵ}{\mathrm{d}P}(ϵ(r)+P(r))]\hfill \\ \hfill & +\frac{3}{r^2}+2{(1-2\frac{M(r)}{r})}^{-1}{(\frac{M(r)}{r^2}+4\pi rP(r))}^2\}\hfill \\ \hfill & +\frac{2\beta }{r}{(1-2\frac{M(r)}{r})}^{-1}\hfill \\ \hfill & \{\frac{M(r)}{r}+2\pi r^2(ϵ(r)-P(r))-1\},\hfill \end{array}$$(34)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}H}{\mathrm{d}r}& =\beta .\hfill \end{array}$$(35)

The above equations have to be solved along with the TOV Eqs. (31). The initial conditions are H(r → 0)=c₀r² and β(r → 0)=2c₀r, where c₀ is a constant, which is irrelevant in the expression for the Love number k₂.

In Fig. 6, we plot the tidal deformability parameters Λ₁ vs. λ₂ of the high- and low-mass members of the binary merger together with the 50% and 90% fidelity regions obtained by the LVC analysis of the GW170817 event (Abbott et al. 2018). Because the extracted most probable masses of the members are below the mass where the chiral phase transition occurs (cf. Fig. 5), the Λ₁ vs. Λ₂ diagrams are the same for given m₀ and αb₀, irrespective of the value of the parameter χ. We note that the tidal deformability parameter requires sufficiently soft EoS around the saturation density. This is seen in the figure, where the smallest tidal deformability and thus the best agreement with the constraint is obtained for the highest value of αb₀ = 450 MeV, which corresponds to the softest EoS at low density. We remark that this is in similar fashion to the proton flow constraint. On the other hand, the 2 M_⊙ requires a sufficiently stiff EoS at higher densities. In Table 5, we show the maximal masses obtained for each set of parameters. Inversely to the tidal deformability, the most massive stars are obtained for the stiffest EoS, namely for αb₀ = 350 MeV. In general, the two constraints are of exclusive character. The interplay between them could be further used to fix the allowed range of external model parameters. This would be of particular use in a class of effective models in which the low- and high-density regimes are not treated independently, but rather combined in a consistent unified framework.

Fig. 6. Relation between tidal deformabilities Λ1 and Λ2 of two compact stars that merged in the GW170817 event (Abbott et al. 2018). For comparison, also shown are the 50% and 90% fidelity regions from the analysis of the GW signal by the LIGO-VIRGO collaboration (Abbott et al. 2018). The gray shading corresponds to the unphysical region Λ2 <  Λ1.

5. Conclusions

Despite the success of LQCD at a finite temperature and vanishing density, the nature of the chiral and deconfinement phase transitions at a low temperature and high density remains unknown. At present, the investigation of matter under such extreme conditions is almost exclusively considered within either purely hadronic or purely quark effective models, with only few approaches to bridge the gap between them. They are typically based on excluded volume approach or a Polyakov loop representing a phenomenological description of the transition from the confined to the deconfined phase (Steinheimer et al. 2011a; Berges et al. 2003; Meyer et al. 2000; Lawley et al. 2006; Dexheimer & Schramm 2010).

In this work, we utilized a superior approach to any previous modeling, the hybrid QMN model (Benic et al. 2015a), to quantify the EoS of cold and dense nuclear matter. This model goes beyond the common two-phase approach to the EoS and embeds the interplay between the quark confinement and the chiral symmetry breaking in a dynamical way into a single unified framework. We have extended the previous version of the model by including repulsive interactions among quarks mediated by the exchange of vector mesons. Within this approach, we systematically investigated the EoS of cold and dense isospin-symmetric matter, as well as asymmetric matter under NS conditions. We constructed the mass-radius relations based on solutions of the TOV equations.

We find that the class of EoSs obtained in the considered model feature sequential chiral and deconfinement phase transitions. We demonstrate that, as a consequence of the realization of the chiral symmetry restoration by parity doubling within the hadronic phase, a strong first-order phase transition invalidates the implication that a flattening, eventually even occurrence of a mass-twin phenomenon, of the mass-radius relation for compact stars at 2 M_⊙, could inevitably signal the deconfinement phase transition in compact stars (Alvarez-Castillo et al. 2016). We verified that such flattening appears due to the chiral phase transition which also lies within the region of the proton flow constraint in isospin-symmetric matter (Danielewicz et al. 2002), as well as within the constraint derived in Hebeler et al. (2013) for the matter under NS conditions. Moreover, the chiral transition is featured in the mass region accessible by the recent GW190425 merger event (Abbott et al. 2020), that is, 1.46 − 1.87 M_⊙, which is of particular interest from the observational perspectives.

The EoS obtained in the hybrid QMN model, under the NS conditions, is comprised of four consecutive density regions, from low to high density: confined and chirally broken phase, confined and chirally restored phase, partially deconfined, and deconfined phase. In the intermediate phases, the confined and deconfined matter may coexist as what could be a possible manifestation of the conjectured quarkyonic phase (Hidaka et al. 2008; McLerran et al. 2009; Andronic et al. 2010; McLerran & Reddy 2019; Jeong et al. 2020). We stress that such a sequential structure of an EoS cannot be obtained in a class of pure quark models or even in a class of hadron-quark models, where the confined and deconfined phases are treated independently. Recent studies in this direction provide an interesting scenario for a rapid increase in the pressure within densities achieved in NSs (McLerran & Reddy 2019; Jeong et al. 2020). It would also be interesting to establish the cooling and transport properties of the hybrid QMN model EoS, especially the mixed phase comprised of nucleons and down quarks, as well as to establish their influence on the thermal and rotational evolution of compact objects.

We embedded the repulsive quark-vector interactions in the QMN model and studied their consequences for the phenomenological description of compact stellar objects. In particular, we analyzed a possible occurrence of a quark matter in the NS core. The inclusion of the quark-vector repulsion implies a shift of the onset of quarks to higher densities. This results in the mass-radius diagram with the maximal mass around 2 M_⊙, which is within the chirally restored hadronic configuration.

In the present studies, the transition to pure quark matter is likely to appear in the part of the stellar sequence that is already gravitationally unstable. We note, however, that the presence of the attractive diquark interaction can induce scalar diquark condensation in the quark phase and reduce the density of the onset of the quark matter (Klähn et al. 2007). This comes with a softening of the EoS which, however, can be compensated by a stronger vector meson coupling (Klähn et al. 2013; Alvarez-Castillo et al. 2019). Moreover, in recent LQCD studies, it was found that parity doubling also occurs in the hyperon channels (Aarts et al. 2019). The inclusion of heavier flavors is known to soften the equation of state. Therefore, it would be of further interest to consider the QMN equations of state that include the attractive diquark interaction and the hyperon degrees of freedom in order to test a model’s compliance with the 2 M_⊙ constraint. Work in this direction is already in progress.

We show how modern astrophysical constraints on the maximum-mass (Cromartie et al. 2019), the tidal deformability from the binary merger GW170817 (Abbott et al. 2018), and recent simultaneous mass-radius constraints from the NICER experiment (Riley et al. 2019; Miller et al. 2019), compiled together, allow for a consistent determination of parameters characterizing the EoS of dense nuclear matter. We demonstrate that not only too soft (excluded by the maximum mass constraint) but also too stiff (the tidal deformability constraint) equations of state may be ruled out in the current approach. With the optimal choice of parameters in the QMN model, we introduced a class of EoSs that fulfill the above astrophysical constraints, as well as phenomenological constraints imposed by flow data from heavy ion collisions. This makes the present hybrid QMN model EoSs superior to previously proposed class of EoSs obtained in pure-hadronic or pure-quark models, as well as in hybrid models based on a two-phase approach. It would be of interest to apply a Bayesian analysis method for selecting the most probable equation of state under a set of modern constraints from compact star physics, and from data obtained in heavy ion collisions. For applications of the QMN model EoS to multi-messenger astronomy in simulations of merger or supernova events as well as heavy ion collision, it is necessary to perform a finite-temperature extension. This is a straightforward task in the presented field-theoretical model as first results for corresponding phase diagrams of the unified quark-hadron matter description convincingly demonstrate (Marczenko et al. 2018, 2019).

We provided a first step toward a unified quark-hadron matter EoS for multi-messenger astronomy that allows conclusions for the composition of neutron star interiors and the phase structure of matter under extreme conditions, in accordance with mass-radius measurements of compact stars and heavy ion collision phenomenology. The next steps include be the generalization of the QMN model to include color superconducting quark matter phases that may allow the existence of hybrid quark-hadron stars due to the lowering of the onset mass for deconfinement while fulfilling the maximum mass constraint for a sufficiently large coupling to vector meson fields.

Acknowledgments

This work was partly supported by the Polish National Science Center (NCN), under OPUS Grant No. 2018/31/B/ST2/01663 (K.R. and C.S.), Opus Grant No. UMO-2019/33/B/ST9/03059 (D.B.), and Preludium Grant No. UMO-2017/27/N/ST2/01973 (M.M.). D.B. received support from the National Research Nuclear University (MEPhI) within the Russian Academic Excellence Project under contract No. 02.a03.21.0005. K.R. also acknowledges stimulating discussions with Larry McLerran and the support of the Polish Ministry of Science and Higher Education.

References

Appendix A: Impact of statistical confinement on EoS

In this appendix, we provide an example of the stiffening of the EoS due to the statistical confinement in the hybrid QMN model. For illustration, let us consider two EoSs. For the first EoS, we take a noninteracting gas of massless fermions at zero temperature. In this case, the pressure, P₁(μ), can be evaluated analytically as a function of chemical potential μ, as

$$\begin{array}{c}\hfill P_1(\mu )=-{\displaystyle \underset{0}{\overset{\mu }{\int }}}\frac{\mathrm{d}p}{2{\pi }^2}{p}^2(p-\mu )=\frac{{\mu }^4}{24{\pi }^2},\end{array}$$(A.1)

where the upper integration limit is due to the Fermi-Dirac distribution function for a massless particle in the zero-temperature limit, θ(μ−p). For the second EoS, let us consider the modification of the distribution function as in Eq. (10). For simplicity, we also set α = 1. The pressure, P₂(μ), is then obtained as

$$\begin{array}{c}\hfill P_2(\mu )=-{\displaystyle \underset{0}{\overset{\mu }{\int }}}\frac{\mathrm{d}p}{2{\pi }^2}{p}^2(p-\mu )\theta (b-p)=\{\begin{array}{c}\frac{{\mu }^4}{24{\pi }^2}\text{, if}b\ge \mu \hfill \\ \frac{b^4}{24{\pi }^2}\text{, if}b<\mu \hfill \end{array}.\end{array}$$(A.2)

Thus, if b ≥ μ, the pressure P₂(μ) evaluates to P₁(μ). However, if b <  μ, the mechanism of statistical confinement sets in and the integration region is limited up to b, and pressure P₂(μ)< P₁(μ). Let us assume that in the region where b <  μ, the functional form of the b field is b = γμ, where γ is a positive constant, such that γ ∈ (0, 1). Then we get the following:

$$\begin{array}{c}\hfill P_2(\mu )=\frac{b^4}{24{\pi }^2}={\gamma }^4\frac{{\mu }^4}{24{\pi }^2}={\gamma }^4{P}_1(\mu )<P_1(\mu ).\end{array}$$(A.3)

Let us consider two values of chemical potential, μ₁ <  μ₂, such that P₁(μ₁)=P₂(μ₂). From this, it follows that ${\mu }_1^4={\gamma }^4{\mu }_2^4$. Because the density is defined as the derivative of pressure with regard to chemical potential, ρ = ∂P/∂μ, we get

$$\begin{array}{c}\hfill {\rho }_2({\mu }_2)={\gamma }^4\frac{{\mu }_2^3}{6{\pi }^2}=\frac{{\mu }_1}{{\mu }_2}\frac{{\mu }_1^3}{6{\pi }^2}=\frac{{\mu }_1}{{\mu }_2}{\rho }_1({\mu }_1)<{\rho }_1({\mu }_1).\end{array}$$(A.4)

The last inequality holds by the assumption that μ₁ <  μ₂. Thus, the pressure P₂ increases faster than P₁ as a function of density ρ. This is the stiffening effect observed in the class of EoSs obtained in the hybrid QMN model. We note that, conversely to the case discussed in this appendix, the appearance of quarks always provides an extra contribution to the pressure, thus softening the corresponding EoS.

All Tables

All Figures

Fig. 1. Thermodynamic pressure, P, for isospin-symmetric matter (left panel) and under the NS conditions of β-equilibrium and charge neutrality (right panel), as functions of the net-baryon number density, ρB, in units of the saturation density, ρ0 = 0.16 fm−3 for m0 = 700 MeV, αb0 = 370 MeV and different values of the parameter χ (see text for details). Left panel: the flow constraint (Danielewicz et al. 2002) is shown as the gray shaded region. In both panels, the inset plots show pressure for χ = 0 and different values of the parameter α, in the vicinity of the constraints obtained by Danielewicz et al. (2002) (left panel) and Hebeler et al. (2013) (right panel). We note that the inset figure in the right panel shows pressure as a function of energy density. In both panels, the first-order phase transitions are seen as plateaux of constant pressure.

In the text

Fig. 2. Composition of compact star matter as a function of the central baryon density in the units of saturation density for the present model for m0 = 700 MeV and αb0 = 370 MeV. Three different cases are considered for χ = 0 (top panel), χ = 0.1 (middle panel), χ = 0.2 (bottom panel). In each panel, the yellow shaded areas indicate the density regions of phase coexistence. The gray, solid vertical lines indicate the central densities of the associated maximum-mass neutron stars.

In the text

	Fig. 3. Thermodynamic pressure, P, under the NS conditions of β-equilibrium and charge neutrality, as a function of the energy density, ϵ, for m0 = 700 MeV and χ = 0, and four representative values of the parameter α (see text for details). The phase transitions are seen as plateaux of constant pressure. The gray shaded region marks the constraint obtained by Hebeler et al. (2013).
In the text

	Fig. 4. Speed of sound squared as a function of the net-baryon number density, ρB, in units of the saturation density, ρ0 = 0.16 fm−3 for m0 = 700 MeV, αb0 = 370 MeV and different values of the parameter χ (see text for details). The first order phase transitions are seen as plateaux of vanishing speed of sound.
In the text

Fig. 5. Mass-radius sequences for compact stars as solutions of the TOV equations for m0 = 700 MeV, for αb0 = 350 MeV (upper left), αb0 = 350 MeV (upper right), αb0 = 400 MeV (lower left), αb0 = 450 MeV (lower right). The circles show the coexistence of the chirally broken and chirally restored phases. The onsets of up and down quarks are marked by pentagons and diamonds, respectively. We note that, in the upper-left figure, for χ = 0.2, the onset of up and down quarks is almost simultaneous. The inner (outer) gray band shows the 68.3% (95.4%) credibility regions for the mass of PSR J0740+6620 (Cromartie et al. 2019). The inner (outer) green and purple bands show 50% (90%) credibility regions obtained from the recent GW170817 (Abbott et al. 2018) event 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 (Miller et al. 2019).

In the text

	Fig. 6. Relation between tidal deformabilities Λ1 and Λ2 of two compact stars that merged in the GW170817 event (Abbott et al. 2018). For comparison, also shown are the 50% and 90% fidelity regions from the analysis of the GW signal by the LIGO-VIRGO collaboration (Abbott et al. 2018). The gray shading corresponds to the unphysical region Λ2 <  Λ1.
In the text