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

Abstract

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

1. Introduction

For the investigation of matter under extreme conditions and the structure of its phase diagram the equation of state (EoS) is the key target. In the region of finite temperature and vanishing baryon density ab initio calculations using Monte Carlo simulations of lattice QCD [1] provide a benchmark for developing phenomenological approaches that can be tested, e.g., in heavy-ion collision experiments. Up to now, however, the sign problem prevents the application of lattice QCD methods to the region at low temperature and high baryon density where the possible existence of a first-order phase transition with a critical endpoint has been conjectured. To elucidate the QCD phase structure in this domain inaccessible to terrestrial experiments and present techniques of lattice QCD simulations, valuable information comes from progress in observing the mass–radius (M-R) relationship of compact stars, due to its one-to-one correspondence with the EoS of compact star matter [2] via the solution of the Tolman–Oppenheimer–Volkoff (TOV) equations [3,4]. For the extraction of the compact star EoS via Bayesian analysis techniques using mass and radius measurements as priors, see Refs. [5,6,7]. In particular, in the era of multi-messenger astronomy, it shall soon become possible to constrain the sequence of stable compact star configurations in the mass–radius plane inasmuch that a benchmark for the EoS of cold and dense matter can be deduced from it.

Among the modern observatories for measuring masses and radii of compact stars, the gravitational wave interferometers of the LIGO-Virgo Collaboration (LVC) and the X-ray observatory NICER on-board the International Space Station provide new powerful constraints besides those from radio pulsar timing. In this work, we pay special attention to the state-of-the-art results from the recent measurement of the high mass ${2.17}_{-0.10}^{+0.11}{M}_{\odot }$ for PSR J0740+6620 by the NANOGrav Collaboration [8] and the compactness derived from the tidal deformability measurement for the binary compact star merger GW170817 [9] in its mass range ($1.16$–$1.60M\odot$ for the low-spin prior).

In the study of cold and dense QCD and its applications, commonly used are separate effective models for the nuclear and quark matter phases (two-phase approaches) with a priori assumed first-order phase transition, typically associated with simultaneous chiral and deconfinement transitions. Within this setting, for a constant-speed-of sound model of high-density (quark) matter, a systematic classification of hybrid compact star solutions has been given in [10], which gives a possibility to identify a strong first-order transition in the EoS by the fact that the hybrid star branch in the mass–radius diagram becomes disconnected from the branch of pure neutron stars. However, already before this occurs, a strong phase transition manifests itself by the appearance of an almost horizontal branch on which the hybrid star solutions lie, as opposed to the merely vertical branch of pure neutron stars. In the literature, this strong phase transition has been discussed as due to quark deconfinement [11,12,13]. This conclusion may however be premature since strong phase transitions with a large latent heat occur also within hadronic matter, for instance due to chiral symmetry restoration within the hadronic phase [14]. In the present work, we employ the hybrid quark–meson–nucleon (QMN) model [14,15,16] to explore the implications of dynamical sequential phase transitions at high baryon density on the phenomenology of neutron stars. To improve the description of nuclear matter properties at the saturation density, we extend the previous hybrid QMN model by including a six-point scalar interaction. Our main focus is on the role of the chiral symmetry restoration within the hadronic branch of the EoS.

This paper is organized as follows. In Section 2, we introduce the hybrid quark–meson–nucleon model. In Section 3, we discuss the obtained numerical results on the equation of state under neutron-star conditions. In Section 4, we discuss the obtained neutron-star relations. In Section 5, we present possible realizations of the low-temperature phase diagram. Finally, Section 6 is devoted to summary and conclusions.

2. Hybrid Quark–Meson–Nucleon Model

In this section, we briefly introduce the hybrid QMN model for the QCD transitions at finite temperature, density, and arbitrary isospin asymmetry for the application to the physics of neutron stars [14,15,16].

The hybrid QMN model is composed of the baryonic parity doublet [17,18,19] and mesons as in the Walecka model, as well as quark degrees of freedom as in the standard quark–meson model. The spontaneous chiral symmetry breaking yields the mass splitting between the two baryonic parity partners, while it generates an entire mass of a quark. In this work, we consider a system with $N_f=2$; hence, relevant for this study are the lowest nucleons and their chiral partners, as well as the up and down quarks. The hadronic degrees of freedom are coupled to the chiral fields $\left(\sigma ,\mathit{\pi }\right)$, isosinglet vector field (${\omega }_{\mu }$), and isovector vector field (${\mathit{\rho }}_{\mu }$). The quarks are coupled solely to the chiral fields. The important concept of statistical confinement is realized in the hybrid QMN model by considering a medium-dependent modification of the particle distribution functions.

In the mean-field approximation, the thermodynamic potential of the hybrid QMN model reads [14]

$$\mathsf{\Omega }=\sum _x{\mathsf{\Omega }}_x+V_{\sigma }+V_{\omega }+V_b+V_{\rho }.$$

(1)

where the summation goes over the positive-parity nucleons, i.e., proton ($p_{+}$) and neutron ($n_{+}$), their negative-parity counterparts, denoted as $p_{-}$ and $n_{-}$, and up (u) and down (d) quarks. The positive-parity nucleons are identified as the positively charged and neutral $N\left(938\right)$ states. The negative-parity states are identified as $N\left(1535\right)$ [20]. The kinetic part of the thermodynamic potential in Equation (1), ${\mathsf{\Omega }}_x$, reads

$${\mathsf{\Omega }}_x={\gamma }_x\int \frac{{\mathrm{d}}^{3}p}{{\left(2\pi \right)}^3}T\left[\ln \left(1-n_x\right)+\ln \left(1-{\overline{n}}_x\right)\right].$$

(2)

The spin degeneracy of the nucleons is ${\gamma }_{\pm }=2$ for both positive- and negative-parity states, while the color-spin degeneracy factor for quarks is ${\gamma }_q=2\times 3=6$. The functions $n_x$ are the modified Fermi–Dirac distribution functions for the nucleons

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

(3)

and for the quarks, accordingly

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

(4)

where b is the expectation value of the b-field, and $\alpha$ is a dimensionless model parameter [15,16].

The hybrid QMN model employs confinement/deconfinement mechanism in a statistical sense. The approach used in this model is to introduce IR and UV momentum cutoffs to suppress quarks at low momenta and hadrons at high momenta. This notion has been widely used in effective theories and Dyson–Schwinger approaches [21,22]. In the current approach, the cutoff is replaced with a medium-dependent quantity, which is expected from asymptotic freedom. Such an intrinsic modification of the cutoff is determined self-consistently when the cutoff is regarded as a vacuum expectation value of a scalar field (see Equation (10d)). The role of the $\alpha b_0$ parameter can be understood twofold. First, its lower values trigger the chiral phase transition at lower densities. Second, the chiral phase transition is stronger and the equation of state becomes stiffer for lower values of the parameter. This can be seen in the equation of state and corresponding speed of sound squared as functions of net-baryon-number density (see Section 3).

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 \left(E_x-{\mu }_x\right)}},\hfill \\ \hfill {\overline{f}}_x& =\frac{1}{1+e^{\beta \left(E_x+{\mu }_x\right)}},\hfill \end{array}$$

(5)

with $\beta$ being the inverse temperature, the dispersion relation $E_x=\sqrt{{\mathit{p}}^2+m_x^2}$. The effective chemical potentials for $p_{\pm }$ and $n_{\pm }$ are defined as1

$$\begin{array}{cc}\hfill {\mu }_{p_{\pm }}& ={\mu }_B-g_{\omega }\omega -\frac{1}{2}{g}_{\rho }\rho +{\mu }_Q,\hfill \\ \hfill {\mu }_{n_{\pm }}& ={\mu }_B-g_{\omega }\omega +\frac{1}{2}{g}_{\rho }\rho .\hfill \end{array}$$

(6)

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

$$\begin{array}{cc}\hfill {\mu }_u& =\frac{1}{3}{\mu }_B+\frac{2}{3}{\mu }_Q,\hfill \\ \hfill {\mu }_d& =\frac{1}{3}{\mu }_B-\frac{1}{3}{\mu }_Q.\hfill \end{array}$$

(7)

In Equations (6) and (7), ${\mu }_B$ and ${\mu }_Q$ are the baryon and charge chemical potentials, respectively. The constants $g_{\omega }$ and $g_{\rho }$ couple the nucleons to the $\omega$ and $\rho$ fields, respectively. The strength of $g_{\omega }$ is fixed by the nuclear saturation properties, while the value of $g_{\rho }$ can be fixed by fitting the value of symmetry energy [23]. The properties are shown in

Table 1. Properties of the nuclear ground state at ${\mu }_B=923$MeV and the symmetry energy used in this work.

${\mathit{\rho }}_0$ [fm${}^{-3}$]	$\mathit{E}/\mathit{A}-{\mathit{m}}_{+}$ [MeV]	K [MeV]	${\mathit{E}}_{\mathbf{sym}}$ [MeV]
$0.16$	$-16$	240	31

.

The effective masses of the parity doublers $m_{p_{\pm }}=m_{n_{\pm }}\equiv m_{\pm }$ are given by

$$m_{\pm }=\frac{1}{2}\left[\sqrt{{\left(g_1+g_2\right)}^2{\sigma }^2+4m_0^2}\mp \left(g_1-g_2\right)\sigma \right],$$

(8)

and for quarks, $m_u=m_d\equiv m_q$,

$$m_q=g_q\sigma .$$

(9)

The parameters $g_1$, $g_2$, and $g_q$ are Yukawa-coupling constants, $m_0$ is the chirally invariant mass of the baryons and is treated as an external parameter (for more details, see [14,16]). The values of those couplings can be determined by fixing the fermion masses in the vacuum (see

Table 2. Physical vacuum inputs used in this work.

${\mathit{m}}_{+}$ [MeV]	${\mathit{m}}_{-}$ [MeV]	${\mathit{m}}_{\mathit{\pi }}$ [MeV]	${\mathit{f}}_{\mathit{\pi }}$ [MeV]	${\mathit{m}}_{\mathit{\omega }}$ [MeV]	${\mathit{m}}_{\mathit{\rho }}$ [MeV]
939	1500	140	93	783	775

). The quark mass is assumed to be $m_{+}=3m_q$ in the vacuum. When the chiral symmetry is restored, the masses of the baryonic parity partners become degenerate with a common finite mass $m_{\pm }\left(\sigma =0\right)=m_0$, which reflects the parity doubling structure of the low-lying baryons. This is in contrast to the quarks, which become massless as the chiral symmetry gets restored.

The potentials in Equation (1) are as in the SU(2) linear sigma model,

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

(10a)

$$\begin{array}{cc}\hfill V_{\omega }& =-\frac{1}{2}{m}_{\omega }^2{\omega }_{\mu }{\omega }^{\mu },\hfill \end{array}$$

(10b)

$$\begin{array}{cc}\hfill V_b& =-\frac{1}{2}{\kappa }_b^2{b}^2+\frac{1}{4}{\lambda }_b{b}^4,\hfill \end{array}$$

(10c)

$$\begin{array}{cc}\hfill V_{\rho }& =-\frac{1}{2}{m}_{\rho }^2{\mathit{\rho }}_{\mu }{\mathit{\rho }}^{\mu },\hfill \end{array}$$

(10d)

where ${\lambda }_2={\lambda }_4{f}_{\pi }^2-{\lambda }_6{f}_{\pi }^4-m_{\pi }^2$, and $ϵ=m_{\pi }^2{f}_{\pi }$. $m_{\pi }$, $m_{\omega }$, and $m_{\rho }$ are the $\pi$, $\omega$, and $\rho$ meson masses, respectively, The pion decay constant is denoted as $f_{\pi }$. Their values are shown in Table 2. The constants ${\kappa }_b$ and ${\lambda }_b$ are fixed following Ref. [15]. The parameters ${\lambda }_4$ and ${\lambda }_6$ are fixed by the properties of the nuclear ground state (see Table 1). We note that the introduction of the six-point scalar interaction term in Equation (10a) is essential in order to reproduce the experimental value of the compressibility $K=240\pm 20$MeV [24].

Following the previous studies of the parity-doublet-based models [14,15,16,25], as well as recent lattice QCD results [26,27], we choose rather large values, $m_0=700,800$ MeV. We note that the additional mass, $m_0$, is not associated with spontaneous chiral symmetry breaking. Thus, it has to originate through another mechanism. Although it is unknown how $m_0$ is expressed in terms of the QCD condensates, the constraint $m_0\le 800$MeV is transmuted into the nucleon mass such that at most $15\%$ of the entire mass is generated by the spontaneous chiral symmetry breaking. This is best seen in the chiral limit, where no dimensionful parameters are present in the QCD Lagrangian, but the appearance of the QCD scale breaks the scale invariance. Thus, one expects that both give rise to the emergence of dynamical hadronic scales at low energies [28,29,30]. Thus, the chirally invariant mass, $m_0$, can be identified with the gluon condensate $\langle G_{\mu \nu }{G}^{\mu \nu }\rangle$.

The physical inputs and the model parameters used in this work are summarized in Table 1, Table 2 and

Table 3. Sets of the model parameters used in this work. The values of ${\lambda }_4$, ${\lambda }_6$ and $g_{\omega }$ are fixed by the nuclear ground state properties, and $g_{\rho }$ by the symmetry energy (see the text). The remaining parameters, $g_q$, ${\kappa }_b$, and ${\lambda }_b$, do not depend on the choice of $m_0$, and their values are taken from Ref. [16].

${\mathit{m}}_0$ [MeV]	${\mathit{\lambda }}_4$	${\mathit{\lambda }}_6{\mathit{f}}_{\mathit{\pi }}^2$	${\mathit{g}}_{\mathit{\omega }}$	${\mathit{g}}_{\mathit{\rho }}$	${\mathit{g}}_1$	${\mathit{g}}_2$	${\mathit{g}}_{\mathit{q}}$	${\mathit{\kappa }}_{\mathit{b}}$ [MeV]	${\mathit{\lambda }}_{\mathit{b}}$
700	33.74	13.20	5.60	8.10	13.75	7.72	3.36	155	0.074
800	21.50	8.25	7.27	7.92	12.91	6.88

. In-medium profiles of the mean fields are obtained by extremizing the thermodynamic potential in Equation (1). The gap equations are obtained as follows

$$\begin{array}{cc}\hfill \frac{\partial \mathsf{\Omega }}{\partial \sigma }& =-{\lambda }_2\sigma +{\lambda }_4{\sigma }^3-{\lambda }_6{\sigma }^5-ϵ+\sum _{x=p_{\pm },n_{\pm },u,d}{s}_x\frac{\partial m_x}{\partial \sigma }=0,\hfill \end{array}$$

(11a)

$$\begin{array}{cc}\hfill \frac{\partial \mathsf{\Omega }}{\partial \omega }& =-m_{\omega }^2\omega +g_{\omega }\sum _{x=p_{\pm },n_{\pm }}{\rho }_x=0,\hfill \end{array}$$

(11b)

$$\begin{array}{cc}\hfill \frac{\partial \mathsf{\Omega }}{\partial b}& =-{\kappa }_b^{2}b+{\lambda }_b{b}^3+\alpha \sum _{x=p_{\pm },n_{\pm }}{\widehat{\omega }}_x-\sum _{x=u,d}{\widehat{\omega }}_x=0,\hfill \end{array}$$

(11c)

$$\begin{array}{cc}\hfill \frac{\partial \mathsf{\Omega }}{\partial \rho }& =-m_{\rho }^2\rho +\frac{1}{2}{g}_{\rho }\sum _{x=p_{\pm }}{\rho }_x-\frac{1}{2}{g}_{\rho }\sum _{x=n_{\pm }}{\rho }_x=0,\hfill \end{array}$$

(11d)

where the scalar and baryon densities are

$$s_x={\gamma }_x\int \frac{{\mathrm{d}}^{3}p}{{\left(2\pi \right)}^3}\frac{m_x}{E_x}\left(n_x+{\overline{n}}_x\right),$$

(12)

and

$${\rho }_x={\gamma }_x\int \frac{{\mathrm{d}}^{3}p}{{\left(2\pi \right)}^3}\left(n_x-{\overline{n}}_x\right),$$

(13)

respectively. The boundary terms in the gap Equation (11c) are given as

$${\widehat{\omega }}_{\pm }={\gamma }_{\pm }\frac{{\left(\alpha b\right)}^2}{2{\pi }^2}T{\left[\ln \left(1-f_{\pm }\right)+\ln \left(1-{\overline{f}}_{\pm }\right)\right]}_{{\mathit{p}}^2={\left(\alpha b\right)}^2},$$

(14)

and

$${\widehat{\omega }}_q={\gamma }_q\frac{b^2}{2{\pi }^2}T{\left[\ln \left(1-f_q\right)+\ln \left(1-{\overline{f}}_q\right)\right]}_{{\mathit{p}}^2=b^2},$$

(15)

for the nucleons and quarks, respectively. Note that the terms in Equations (14) and (15) come into the gap Equation (11c) 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=-\mathsf{\Omega }+{\mathsf{\Omega }}_0$, where ${\mathsf{\Omega }}_0$ is the value of the thermodynamic potential in the vacuum. The net-baryon number density for a species x is defined as

$${\rho }_B^x=-\frac{\partial {\mathsf{\Omega }}^x}{\partial {\mu }_B},$$

(16)

where ${\mathsf{\Omega }}^x$ is the kinetic term in Equation (2) for the species x. The total net-baryon number density reads

$${\rho }_B={\rho }_B^{n_{+}}+{\rho }_B^{n_{-}}+{\rho }_B^{p_{+}}+{\rho }_B^{p_{-}}+{\rho }_B^u+{\rho }_B^d.$$

(17)

In the next section, we discuss the obtained equations of state in the hybrid QMN model and their impact on the chiral phase transition, under the neutron-star conditions of $\beta$ equilibrium and charge neutrality.

3. Equation of State under Neutron-Star Conditions

The neutron-star conditions require additional constraints to be imposed on the EoS under investigation. To this end, electrons and muons are included as gases of free relativistic particles. The first constraint is the $\beta$-equilibrium. This condition is an equilibrium among protons, neutrons, and charged leptons. It assumes that the energy of the system is minimized, the system is electrically neutral, and the total net-baryon number is conserved. $\beta$-equilibrium condition can be expressed in terms of chemical potentials,

$${\mu }_{n_{+}}={\mu }_{p_{+}}+{\mu }_{e/\mu },$$

(18)

where ${\mu }_{n_{+}}$, ${\mu }_{p_{+}}$, ${\mu }_e$, and ${\mu }_{\mu }$ are the neutron, proton, electron, and muon chemical potentials, respectively. The electric-charge neutrality constraint dictates that the overall charge density in a neutron star has to be zero,

$${\rho }_Q^{p_{+}}+{\rho }_Q^{p_{-}}+\frac{2}{3}{\rho }_Q^u-\frac{1}{3}{\rho }_Q^d-{\rho }_Q^e-{\rho }_Q^{\mu }=0,$$

(19)

where ${\rho }_Q^x$ is the charge density of a species x.

In Figure 1, we show the calculated zero-temperature equations of state in the mean-field approximation with $m_0=700$MeV (Figure 1, left) and $m_0=800$MeV (Figure 1, right), for different values of the $\alpha$ parameter, namely $\alpha b_0=350$ MeV (red, solid line), $\alpha b_0=370$ MeV (purple, dashed line), $\alpha b_0=400$ MeV (blue, dotted line) and $\alpha b_0=450$ MeV (black, dash-dotted line). The value $b_0$ denotes the vacuum expectation value of the b-field. The coexistence phases of the chirally broken and restored phases are shown between circles. We stress that the chiral and deconfinement phase transitions are sequential in the current model setup (see [16]). The latter happen at higher densities and are not shown in the figure.

In all cases, the behavior at low densities is similar. In general, for low values of $\alpha b_0$ (except $\alpha b_0=450$MeV), the chiral transition is of first order, determined as a jump in the $\sigma$-field expectation value. The parity-doublet nucleons become degenerate with mass $m_{\pm }=m_0$. The chiral phase transition becomes weaker for higher values of the $\alpha$ parameter. For $\alpha b_0=450$ MeV, the transition turns into a smooth crossover, defined as a peak in $\partial \sigma /\partial {\mu }_B$. This behavior agrees with the case of isospin-symmetric matter, where higher value of $\alpha$ causes the first-order chiral phase transition to weaken and eventually go through a critical point, and turn into a crossover transition [16]. The values of the net-baryon density range for the coexistence phase of the chirally broken and restored phases are shown in

Table 4. Net-baryon density range of the coexistence phase of the chirally broken and restored phases in terms of saturation density units, ${\rho }_0$, for different values of $m_0$ and $\alpha b_0$ parameters. For the case of $\alpha b_0=450$MeV, the transitions are smooth crossovers for both values of $m_0$.

	$\mathit{\alpha }{\mathit{b}}_0$ [MeV]
${\mathit{m}}_{\mathbf{0}}$ [MeV]	350	370	400	450
700	$1.82$–$2.60$	$2.12$–$2.76$	$2.60$–$3.07$	$3.56$
800	$1.94$–$2.97$	$2.29$–$3.15$	$2.86$–$3.66$	$4.13$

.

In Figure 2, 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 are shown in between circles. As seen in the figure, the causality bound is preserved for all of the parameterizations. The apparent stiffening of the EoSs is a result of the modification of the Fermi–Dirac distributions (cf. Equation (3)) introduced in the hybrid QMN model. We note that it is in general possible to sustain the $2M_{\odot }$ constraint and fulfill the conformal bound, i.e., $c_s^2\le 1/3$. This can be obtained, e.g., in a class of constant-speed-of-sound equations of state [31].

4. TOV Solutions for Compact-Star Sequences

We use the equations of state introduced in the previous section (see Figure 1) to solve the general-relativistic TOV equations [3,4] for spherically symmetric objects,

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

(20a)

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

(20b)

with the boundary conditions $P(r=R)=0$ and $M=M(r=R)$, where R and M are the radius and the mass of a neutron star, respectively. Once the initial conditions are specified based on a given equation of state, namely the central pressure $P_c$ and the central energy density ${ϵ}_c$, the internal profile of a neutron star can be calculated.

In general, there is one-to-one correspondence between an EoS and the mass–radius relation calculated with it. In Figure 3 (left), we show the relationship of mass vs. central net-baryon number density, for the calculated sequences of compact stars, together with the state-of-the-art constraints on the maximum mass for the pulsar PSR J0348-0432 [32] and PSR J0740+6620 [8]. We point out that the chiral phase transition leads to a softening of the EoS so that it is accompanied by a rapid flattening of the sequence. Notably, the chiral transition for all values of $\alpha b_0$ occurs in the high-mass part of the sequence, but below the $2M_{\odot }$ constraint, at around $1.8M_{\odot }$.

In Figure 3 (left), the three curves for $\alpha b_0=350,370,400$MeV consist of three phases: the chirally broken phase in the low-mass part of the sequence, the chirally restored phase in the high-mass part, and the coexistence phase between filled circles. Similar to the equation of state, increasing the value of $\alpha$ softens the chiral transition, which eventually becomes a smooth crossover for $\alpha b_0=450$MeV and consists only of branches with chiral symmetry being broken and restored, separated by a circle.

In Figure 3, we show mass vs. central net-baryon density relations obtained for different values of the chirally invariant mass $m_0$. What is evident is that increasing the value of $m_0$ strengthens the chiral phase transition. This is seen twofold, as a shrinking of the coexistence phases and as more abrupt flattening of chirally restored branches. For a larger $m_0$, the transition becomes strong enough to produce disconnected branches (see, e.g., the red, solid line in the bottom right panel of Figure 3). These, in turn, cause the maximal mass of the sequences to decrease with increasing value of $m_0$. Eventually, the equations of state become not stiff enough to reach the $2M_{\odot }$ constraint. We note that such small maximal masses are result of the additional six-point interaction term considered in the thermodynamic potential of the hybrid QMN model (see Equation (10)). For $m_0=700$MeV (Figure 3, top), the most favorable parameterizations are $\alpha b_0=350-370$MeV, while for $m_0=800$MeV (Figure 3, bottom) none of the EoSs is stiff enough. In

Table 5. Maximal neutron-star masses in units of $M_{\odot }$ and corresponding radius in km (separated by comma) for different values of $m_0$ and $\alpha b_0$ parameters.

	$\mathit{\alpha }{\mathit{b}}_0$ [MeV]
${\mathit{m}}_{\mathbf{0}}$ [MeV]	350	370	400	450
700	$2.10,12.11$	$2.05,11.91$	$2.01,11.81$	$1.96,11.92$
800	$1.95,11.64$	$1.88,11.29$	$1.83,11.22$	$1.79,11.25$

, we show the values of the maximal masses of neutron star and corresponding radii obtained in each parameterization. In Figure 4, we show the radial profiles of energy density (Figure 4, top) and pressure (Figure 4, bottom), for a $2.01M_{\odot }$ neutron star, calculated for $m_0=700$MeV and $\alpha b_0=370$MeV. The chiral phase transition happens at roughly $7.4$km from the center of the star and is reflected in a jump in energy density.

The end points of the mass–radius relation correspond to the onset of quark d.o.f. in each parameterization. This leads to the conclusion that the hadronic matter is not stiff enough to fulfill the two-solar-mass constraint. In general, a possible resolution to this problem could be another phase transition. This is the case in the hybrid QMN model, which features sequential chiral and deconfinement phase transitions. However, in the current model setup, the equation of state in the deconfined phase is not stiff enough to sustain the gravitational collapse and the branches become immediately unstable. This is because quarks are not coupled with the vector field leading to a repulsive force. On the other hand, it is known that repulsive interactions tend to stiffen the equation of state. Hence, an additional repulsive force in the quark sector could possibly make the branch stiff enough in order to reach the $2M_{\odot }$ constraint, and an additional family of stable hybrid compact stars would appear, with the possibility for the high-mass twin scenario advocated by other effective models [33,34,35].

We note that the obtained mass–radius relations stay in good agreement with the low-mass constraints derived from the recent neutron-star merger GW170817 for the low- and high-mass posteriors [9]. In Figure 3 (right), they are shown as green and purple regions, respectively.

5. Isospin-Symmetric Phase Diagram

The observational neutron-star data provide useful constraints on the structure of strongly interacting matter. Furthermore, they may constrain the phase diagram of isospin-symmetric QCD matter, which is of major relevance for the heavy-ion physics. In Figure 5, we show the low-temperature part of the isospin-symmetric phase diagram obtained in the model in the $(T,{\rho }_B)$-plane for $m_0=700$MeV (Figure 5, left) and $m_0=800$MeV (Figure 5, right). The liquid–gas phase transition (green, dashed-doubly-dotted line) is common for both values of $m_0$ by construction of the hybrid QMN model [15,16]. Its critical point shows up at around $T=16$MeV, above which it turns into crossover. A similar phase structure is developed for the chiral phase transition for low values of $\alpha b_0$. For $m_0=700$MeV, the critical points appear around $T=19,9$MeV for $\alpha b_0=350,370$MeV, respectively. On the other hand, for $\alpha b_0=400,450$MeV, the chiral transition proceeds as a smooth crossover at all temperatures. For $m_0=800$MeV, the critical points are developed at $T=36,26,15,1$MeV for $\alpha b_0=350,370,400,450$MeV, respectively. Higher values of the temperatures for the critical points are essentially a result of much stronger chiral phase transition at zero temperature. We note that the most favorable parameterizations, i.e., for smaller values of $m_0$, yield rather low temperature for the critical point of the chiral phase transition.

6. Conclusions

In this work, we investigated the hybrid QMN model for the equation of state of dense matter under neutron-star conditions and the phenomenology of compact stars. In particular, we focused on the implications of including six-point scalar interaction and studied the consequences of the realization of the chiral symmetry restoration within the hadronic phase.

We found that the apparent softening of the EoS results in mass–radius relations with maximal mass at tension with the $2M_{\odot }$ constraint within the hadronic branch, especially if the new PSR J0740+6620 with $2.17M_{\odot }$ [8] is considered. We have shown that parameterizations of the model which yield large maximal mass (i.e., for smaller value of $m_0$) suggest rather low value of the temperature for the critical end point of the first-order chiral phase transition in the phase diagram, which may also be absent. In view of this, if would interesting to establish a constraint on the chirally invariant mass $m_0$. Since the hybrid QMN model features sequential chiral and deconfinement phase transitions, one possible resolution to this could be the onset of quark d.o.f. Such a scenario would be even further supported in view of the recent formulation of the three-flavor parity doubling [36,37] and further lattice QCD studies [27], where it was found that, to a large extent, the phenomenon occurs also in the hyperon channels. In general, the inclusion of heavier flavors is known to soften the equation of state and additional repulsive forces are needed to comply with the $2M_{\odot }$ constraint. Additional stiffness from the quark side would play a role, which is not included in the current study. Work in this direction is in progress and the results will be reported elsewhere.