Treating quarks within neutron stars

TL;DR

This paper investigates the presence of deconfined quark cores in neutron stars using various phase transition models, finding that quarks appear at certain masses and can be consistent with observational constraints, but distinguishing hybrid stars remains challenging.

Contribution

It systematically explores different quark-hadron transition scenarios in neutron stars and their compatibility with recent observational data, highlighting the conditions for quark core formation.

Findings

Quarks appear in neutron stars with masses around 1.0-1.6 solar masses.

Supported quark matter requires a squared speed of sound greater than 0.4.

Certain equations of state are consistent with GW170817 tidal deformability bounds.

Abstract

Neutron star interiors provide the opportunity to probe properties of cold dense matter in the QCD phase diagram. Utilizing models of dense matter in accord with nuclear systematics at nuclear densities, we investigate the compatibility of deconfined quark cores with current observational constraints on the maximum mass and tidal deformability of neutron stars. We explore various methods of implementing the hadron-to-quark phase transition, specifically, first-order transitions with sharp (Maxwell construction) and soft (Gibbs construction) interfaces, and smooth crossover transitions. We find that within the models we apply, hadronic matter has to be stiff for a first-order phase transition and soft for a crossover transition. In both scenarios and for the equations of state we employed, quarks appear at the center of pre-merger neutron stars in the mass range $\approx 1.0-1.6\,{\rm M}_{\odot}$, with a squared speed of sound $c^2_{\rm QM}\gtrsim 0.4$ characteristic of strong repulsive interactions required to support the recently discovered neutron star masses $\geq 2\,{\rm M}_{\odot}$. We also identify equations of state and phase transition scenarios that are consistent with the bounds placed on tidal deformations of neutron stars in the recent binary merger event GW170817. We emphasize that distinguishing hybrid stars with quark cores from normal hadronic stars is very difficult from the knowledge of masses and radii alone, unless drastic sharp transitions induce distinctive disconnected hybrid branches in the mass-radius relation.

Tables

Table 1. Table 1: RMFT coupling strengths. Values of the meson masses used are m σ = 660 subscript 𝑚 𝜎 660 m_{\sigma}=660 MeV, m ω = 783 subscript 𝑚 𝜔 783 m_{\omega}=783 MeV and m ρ = 770 subscript 𝑚 𝜌 770 m_{\rho}=770 MeV.

Model	$g_{\sigma }$	$g_{\omega }$	$g_{\rho }$	$\kappa$	$\lambda$
MS-A	12.819	12.258	12.079	0.02544	-0.02179
MS-B	11.369	10.143	9.446	0.05098	-0.03396
MS-C	10.026	7.961	8.492	0.10841	-0.00365
Model	$\zeta$	$\xi$	${\Lambda }_{\sigma }$	${\Lambda }_{\omega }$
MS-A	0.0001	1.0	0.001	0.05
MS-B	0.0001	1.0	0.001	0.05
MS-C	0.0001	1.0	0.001	0.05

Table 2. Table 2: Properties at the nuclear equilibrium density n 0 subscript 𝑛 0 n_{0} for EoSs used in this work compared to that of the APR EoS Constantinou et al. ( 2014 ) . Entries in this table are the Landau effective mass m ∗ / M superscript 𝑚 𝑀 m^{*}/M , isospin symmetric matter compression modulus K 0 subscript 𝐾 0 K_{0} , kinetic and interaction parts S 2 ​ k subscript 𝑆 2 𝑘 S_{2k} and S 2 ​ d subscript 𝑆 2 𝑑 S_{2d} of the total symmetry energy S 2 subscript 𝑆 2 S_{2} , and the corresponding parts of the symmetry energy slope parameter L 𝐿 L . In MS models, m ∗ = E F ∗ = k F 2 + M ∗ 2 superscript 𝑚 superscript subscript 𝐸 𝐹 superscript subscript 𝑘 𝐹 2 superscript 𝑀 absent 2 m^{*}=E_{F}^{*}={\sqrt{k_{F}^{2}+M^{*2}}} .

Property	APR	MS-A	MS-B	MS-C	Units
$n_0$	0.16	0.16	0.16	0.16	${\mathrm{fm}}^{-3}$
$m^{\ast }/M$	0.698	0.662	0.763	0.847
$K_0$	266	230	230	230	MeV
$S_{2k}$	9.79	18.55	16.09	14.49	MeV
$S_{2d}$	22.80	11.45	13.91	15.51	MeV
$S_2$	32.58	30.0	30.0	30.0	MeV
$L_k$	12.69	61.74	44.35	34.52	MeV
$L_d$	45.78	-13.40	8.65	30.88	MeV
$L$	58.47	48.34	53.00	65.40	MeV

Table 3. Table 3: Structural properties of nucleonic neutron stars with M = 1.4 ​ M ⊙ 𝑀 1.4 subscript M direct-product M=1.4\,{\rm M}_{\odot} and M max subscript 𝑀 max M_{\rm max} for the indicated EoSs. For each mass, the compactness parameter β = ( G ​ M / c 2 ​ R ) ≃ ( 1.475 ​ R ) ​ ( M / M ⊙ ) 𝛽 𝐺 𝑀 superscript 𝑐 2 𝑅 similar-to-or-equals 1.475 𝑅 𝑀 subscript M direct-product \beta=(GM/c^{2}R)\simeq(1.475\,R)(M/{\rm M}_{\odot}) ; n c subscript 𝑛 𝑐 n_{c} , P c subscript 𝑃 𝑐 P_{c} and y c subscript 𝑦 𝑐 y_{c} are the central values of the density, pressure and proton fraction, respectively.

Property	APR	MS-A	MS-B	MS-C	Units
$R_{1.4}$	11.74	13.21	12.41	11.85	km
${\beta }_{1.4}$	0.176	0.157	0.167	0.174
$n_{c,1.4}/n_0$	3.35	2.05	2.80	3.72
$P_{c,1.4}$	89.33	41.78	64.43	94.24	$\mathrm{MeV}{\mathrm{fm}}^{-3}$
$y_{c,1.4}$	0.11	0.104	0.106	0.106
$R_{\max }$	10.26	12.44	10.91	9.94	km
$M_{\max }$	2.185	2.63	2.21	1.83	$M_{\odot }$
${\beta }_{\max }$	0.314	0.312	0.299	0.273
$n_{c,\max }/n_0$	6.97	4.71	6.38	8.30
$P_{c,\max }$	884.69	498.32	632.66	664.60	$\mathrm{MeV}{\mathrm{fm}}^{-3}$
$y_{c,\max }$	0.16	0.14	0.14	0.128

Table 4. Table 4: Effective mass dependence of the L 𝐿 L - R 1.4 subscript 𝑅 1.4 R_{1.4} relation. Entries are as in Tables  2 and 3 , but organized differently.

Model	$m^{\ast }/M$	${\Lambda }_{\sigma }$	${\Lambda }_{\omega }$	$L\left(\mathrm{MeV}\right)$	$R_{1.4}\left(\mathrm{km}\right)$	$n_{c,1.4}/n_0$
MS-A	0.662	0.001	0.05	48.34	13.21	2.05
	0.662	0.001	0.01	85.49	13.81	2.00
	0.662	0.0	0.0	96.1	14.07	1.93
MS-B	0.763	0.001	0.05	52.99	12.41	2.80
	0.763	0.001	0.01	78.71	12.93	2.68
	0.763	0.0	0.0	86.08	13.25	2.53
MS-C	0.847	0.001	0.05	65.40	11.85	3.72
	0.847	0.001	0.01	77.53	12.41	3.35
	0.847	0.0	0.0	81.05	12.67	3.12

Table 5. Table 5: Mapping onto the CSS phase transition parameters Alford et al. ( 2013 ) for the stiff MS-A+vMIT/vNJL hybrid EoSs with Maxwell construction; see also Fig.  2 , panels (a)-(c). Meanings of the various entries are explained in the associated text.

$(\frac{B^{1/4}}{\mathrm{MeV}},a)$	$\frac{n_{\mathrm{trans}}}{n_0}$	$\frac{P_{\mathrm{trans}}}{{\epsilon }_{\mathrm{trans}}}$	$\frac{\Delta \epsilon }{{\epsilon }_{\mathrm{trans}}}$	$c_{\mathrm{QM}}^2$	$\frac{\Delta {\epsilon }_{\mathrm{crit}}}{{\epsilon }_{\mathrm{trans}}}$	$M\left(R\right)$
(159, 0.2)	1.77	0.084	0.33	0.407	0.626	Connected
(162, 0.2)	2.10	0.136	0.33	0.416	0.704	Connected
(165, 0.2)	2.34	0.180	0.38	0.424	0.77	Connected
(180, 0.0)	2.04	0.127	1.13	0.326	0.691	Disconnected
(170, 0.1)	2.08	0.133	0.63	0.380	0.70	Both
(155, 0.3)	2.08	0.134	0.13	0.442	0.701	Connected
$(\frac{B_{\mathrm{eff}}^{1/4}}{\mathrm{MeV}},\frac{G_v}{G_s})$	$\frac{n_{\mathrm{trans}}}{n_0}$	$\frac{P_{\mathrm{trans}}}{{\epsilon }_{\mathrm{trans}}}$	$\frac{\Delta \epsilon }{{\epsilon }_{\mathrm{trans}}}$	$c_{\mathrm{QM}}^2$	$\frac{\Delta {\epsilon }_{\mathrm{crit}}}{{\epsilon }_{\mathrm{trans}}}$	$M\left(R\right)$
(218.3, 0.15)	2.91	0.284	0.594	0.236	0.926	Connected 111This connected branch is tiny $\left(M_{\max }-M_{\mathrm{trans}}\right)\lesssim {10}^{-3}{M}_{\odot }$ (invisible on the magnified $M\left(R\right)$ plot; see Fig. 4 (c)) and thus hybrid stars are undetectable.

Table 6. Table 6: Summary of the minimum density and minimum neutron star mass when quarks start to appear in various treatments of phase transitions explored in our work; see also the indicated figures for detailed information.

Treatment	$n_{\mathrm{trans}}/n_0$	$M_{\mathrm{trans}}$	Figure reference
Maxwell	1.77	$0.97M_{\odot }$	Fig. 2 (b) & (c)
Gibbs	$\lesssim 1.5$	$\lesssim 0.6M_{\odot }$	Fig. 2 (e) & (f)
Interpolation	2.0, 1.5	$0.81,0.48M_{\odot }$	Fig. 5 (b) & (c)
Quarkyonic	2.31	$0.97,1.21M_{\odot }$	Fig. 6 (b) & (c)

Table 7. Table 7: Summary of different treatments with the introduction of quarks in the dense matter EoS. For a sharp first-order transition with Maxwell construction, the most readily compatible scenario is stiff hadronic matter undergoes phase transition into stiff ( c QM 2 ≳ 0.4 greater-than-or-equivalent-to subscript superscript 𝑐 2 QM 0.4 c^{2}_{\rm QM}\gtrsim 0.4 ) quark matter. For a crossover transition, soft HM → → \to stiff QM is necessary.

HM $\to$ QM	First-order Transition (Maxwell) 222See text for details if the Gibbs construction is applied. Gibbs construction satisfies many observational constraints such as $R_{1.4}$ and $M_{\max }$ due to the earlier onset of quarks. However, distinguishability from purely-hadronic stars is lost.	Crossover Transition
stiff to soft	✗ vMIT: cannot support $M_{\max }\ge 2M_{\odot }$	✗ unphysical decreasing function of $P\left(n_B\right)$
	✗ vNJL: $M_{\mathrm{trans}}\gtrsim 1.7M_{\odot }$, $\tilde{\Lambda }\left(ℳ=1.186M_{\odot }\right)>720$
soft to stiff	✗ no intersection for $P\left(\mu \right)$ 333Limited by the specific quark models applied here; in a generic parametrization (e.g. CSS) the soft HM $\to$ stiff QM is possible.	✓interpolation: $M_{\mathrm{trans}}\lesssim 1.0M_{\odot }$,
		✓quarkyonic: $M_{\max }\ge 2M_{\odot }$; $R_{1.4}$ and $M_{\mathrm{trans}}$ vary
soft to soft	✗ cannot support $M_{\max }\ge 2M_{\odot }$	✗ cannot support $M_{\max }\ge 2M_{\odot }$
stiff to stiff	✓vMIT: $M_{\max }\ge 2M_{\odot }$; $R_{1.4}$ and $M_{\mathrm{trans}}$ vary	✗ $\tilde{\Lambda }\left(ℳ=1.186M_{\odot }\right)>720$, $R_{1.4}>13\mathrm{km}$
	✗ vNJL: onset for quarks too high; immediately destabilize

Equations

$$\mathcal{L}$$

$$V(\sigma )$$

$$f(\sigma ,\omega )$$

$$S_2$$

$$m_{\rho }^{\ast 2}$$

$$L$$

$$L_k$$

$$L_d$$

$$\mathcal{L}$$

$$n_i$$

$${\epsilon }_Q$$

$$P_Q$$

$${\mu }_i$$

$${\epsilon }_Q$$

$$P_Q$$

$${\mu }_i$$

$$\mathcal{L}$$

$$\epsilon$$

$${\epsilon }_{\mathrm{int}}$$

$$m_i=m_{0,i}-4G_s⟨{\overline{q}}_i{q}_i⟩+2K⟨{\overline{q}}_j{q}_j⟩⟨{\overline{q}}_k{q}_k⟩,$$

$$⟨{\overline{q}}_i{q}_i⟩$$

$${\mu }_i={\frac{d\epsilon }{d{n}_i}|}_{n_j,n_k}=\sqrt{m_i^2+k_{Fi}^2}+4G_v{n}_i$$

$$P=\sum _i{n}_i{\mu }_i-\epsilon =P_{\mathrm{FG}}-B_{\mathrm{eff}}.$$

$${\epsilon }_Q$$

$$P_Q$$

$${\mu }_i$$

$$B_{\chi }^i$$

$$(1+x_e)+{\left[x_s^{2/3}+\frac{C_s-C_d}{r^{2/3}}\right]}^{3/2}+x_s$$

$${\left[(1+x_e{)}^{2/3}+\frac{C_u}{r^{2/3}}\right]}^{1/2}-{\left[x_s^{2/3}+\frac{C_s}{r^{2/3}}\right]}^{1/2}+(3x_e{)}^{1/3}$$

$$x_u=1+x_e\mathrm{and}{x}_d={\left[x_s^{2/3}+\frac{C_s-C_d}{r^{2/3}}\right]}^{3/2},$$

$$Q$$

$$n_{\mathrm{B}}$$

$$P(n)$$

$$f_{\pm }(n)$$

$$\epsilon (n)$$

$$\Delta \epsilon$$

$$g(n^{\prime })$$

$$k_{Fq}=\frac{(k_{FB}-\Delta )}{N_c}\Theta (k_{FB}-\Delta ),$$

$$\Delta =\frac{{\Lambda }_Q^3}{k_{FB}^2}+\kappa \frac{{\Lambda }_Q}{N_c^2}.$$

Full text

I introduction

The observation that the dense matter inside neutron stars might consist of weakly interacting quark matter owing to the asymptotic freedom of Quantum Chromodynamics (QCD) was first made by Collins and Perry Collins and Perry (1975). Since then, numerous explorative studies have been conducted to isolate neutron star observables that can establish the presence of quarks deconfined from hadrons. Starting from the QCD Lagrangian, lattice gauge simulations at finite temperature $T$ and net baryon number $n_{\rm B}=0$ naturally realize hadronic and quark degrees of freedom in a smooth crossover transition. However, lattice simulations for finite $n_{\rm B}$ at $T=0$, of relevance to neutron stars, have been thwarted due to the unsolved fermion sign problem and untenable imaginary probabilities. As a result, the possible phases of dense matter at $T=0$ have been generally explored by constructing equation of state (EoS) models of hadrons and quarks that are independent of each other although a few exceptions do exist.

Extensive studies of nucleonic matter in neutron stars for $n_{\rm B}\lesssim 0.5\,n_{0}$, where $n_{0}\simeq 0.16~{}{\rm fm}^{-3}$ is the isospin symmetric nuclear matter equilibrium density, have predicted the presence of a solid crust. Observations of the surface temperatures of accreting neutron stars in their quiescent periods have indeed confirmed the presence of a crust (see Ref. Meisel et al. (2018), and references therein). This region is characterized by a Coulomb lattice of neutron-rich nuclei surrounded by dripped neutrons with admixtures of light nuclei and a uniform background of electrons in chemical potential and pressure equilibrium in a charge-neutral state. Differences among different equations of state Baym et al. (1971); Negele and Vautherin (1973); Douchin and Haensel (2001); Sharma et al. (2015) are small and are of minor importance to the structure of stars more massive than $1\,{\rm M}_{\odot}$. In this work, we use the EoSs of Ref. Negele and Vautherin (1973) (for $0.001<n_{\rm B}<0.08~{}{\rm fm}^{-3}$) and Ref. Baym et al. (1971) (for $n_{\rm B}<0.001~{}{\rm fm}^{-3}$) to determine the structural properties of the star.

Models of the hadronic EoS for $n_{\rm B}>0.08~{}{\rm fm}^{-3}$ can be grouped into three broad categories: non-relativistic potential models, Dirac-Brueckner-Hartree-Fock models, and relativistic field-theoretical models. Microscopic many-body calculations in the first two of these categories (e.g., Brueckner-Hartree-Fock, variational, Greens’ function Monte Carlo, chiral effective field theory, as well as Dirac-Brueckner-Hartree-Fock) employ free-space two-nucleon interactions supplemented by three-nucleon interactions required to describe the properties of light nuclei as input. In contrast, coupling strengths of the two- and higher-body nucleon interactions mediated by meson exchanges are calibrated at $n_{0}$ in the relativistic field-theoretical models. Several schematic potential models based on zero- and finite-range forces also exist that take recourse in the Hohenberg-Kohn-Sham theorem Hohenberg and Kohn (1964); Kohn and Sham (1965) which assures that the ground state energy of a many-body system can be expressed in terms of local densities alone. Refinements in all of these approaches are guided by laboratory data on the bulk properties of isospin symmetric and asymmetric matter, such as the binding energy $\rm{BE}=-16\pm 1$ MeV Myers and Swiatecki (1966, 1996) at the saturation density $n_{0}=0.16\pm 0.01~{}{\rm fm}^{-3}$ Myers and Swiatecki (1966); Day (1978); Myers and Swiatecki (1996), compression modulus $K_{\rm nm}=240\pm 20$ MeV Garg (2004a); Colò et al. (2004); Shlomo et al. (2006), nucleon’s Landau effective mass $m^{*}/M=0.75\pm 0.1$ Bohigas et al. (1979); Krivine et al. (1980); Margueron et al. (2018), symmetry energy $S_{2}=28-35$ MeV Lattimer and Lim (2013); Tsang et al. (2012) and the symmetry energy slope parameter $L=60\pm 20$ MeV Lattimer and Lim (2013); Tsang et al. (2012) at saturation, etc. Low-to-intermediate energy (0.5-2 GeV) heavy-ion collisions have been used to determine the EoS for densities up to 2-3 $n_{0}$ through studies of matter, momentum, and energy flow of nucleons Gale et al. (1987); Prakash et al. (1988a); Welke et al. (1988); Gale et al. (1990); Danielewicz (2000); Danielewicz et al. (2002). The consensus has been that as long as momentum-dependent forces are employed in models that use Boltzmann-type kinetic equations, use of $K_{\rm nm}\sim 240\pm 20$ MeV, suggested by the analysis of the giant monopole resonance data Youngblood et al. (1999); Garg (2004b); Colo et al. (2004), fits the heavy-ion data as well Danielewicz et al. (2002).

The lack of Lorentz invariance in non-relativistic models leads to an acausal behavior at some high density particularly if contributions from three- and higher-body interactions to the energy are not screened in medium Bludman and Dover (1980); Prakash et al. (1988b). The general practice has been to enforce causality from thermodynamic considerations Nauenberg and Chapline (1973); Lattimer et al. (1990). In some cases, the reliability of non-relativistic models is severely restricted, sometimes only up to $2\,n_{0}$ as in the case of chiral effective field-theoretical (EFT) models owing to the perturbative scheme and the momentum cut-off procedure employed there Hebeler et al. (2010); Tews et al. (2013).

To explore consequences of the many predictions of these models at supra-nuclear densities, piecewise polytropic EoSs that are causal have also been extensively used to map out the range of pressure vs density relations (EoSs) that are consistent with neutron star phenomenology Lattimer and Prakash (2016); Steiner et al. (2016); Tews et al. (2018); Zhao and Lattimer (2018). The viability of these EoSs at supra-nuclear densities necessarily depends on the growing neutron star data to be detailed below.

The possibility of non-nucleonic degrees of freedom such as strangeness-bearing hyperons, pion and kaon condensates, and deconfined quarks above $n_{0}$ has also been examined in many of these models Lattimer and Prakash (2016); Chatterjee and Vidaña (2016); Vidaña (2018). At some $n_{\rm B}\gtrsim(2-4)\,n_{0}$, the presence of quark degrees of freedom has been invoked on the physical basis that the constituents of hadrons could be liberated as the compression in density progressively increases. First-principle calculations Baym and Chin (1976a, b); Freedman and McLerran (1978); Farhi and Jaffe (1984); Kurkela et al. (2010); Kurkela and Vuorinen (2016); Gorda et al. (2018); Alford et al. (2008) of the EoS of quark matter have thus far been limited to the perturbative region of QCD valid at asymptotically high baryon densities. The Nambu–Jona-Lasinio (NJL) model Nambu and Jona-Lasinio (1961), which shares many symmetries with QCD - but not confinement - has been used to mimic chiral restoration in quark matter Kunihiro (1989); Buballa and Oertel (1999); Buballa (2005). Also in common use are variations Klähn and Fischer (2015); Gomes et al. (2019) of the MIT bag model Baym and Chin (1976a).

Lacking knowledge about the nature of the phase transition, it has been common to posit a first-order phase transition in many recent studies Nandi and Char (2018); Paschalidis et al. (2018); Alvarez-Castillo et al. (2019); Wei et al. (2018); Gomes et al. (2019). Even in this case, the magnitude of the hadron-quark interface tension is uncertain Alford et al. (2001); Mintz et al. (2010); Lugones et al. (2013); Fraga et al. (2019). If the interface tension is regarded as being infinite, a Maxwell construction can be employed to determine the range of density for which chemical potential and pressure equality between the hadronic and quark phases exists Glendenning (1992). The other extreme case corresponds to a vanishing interface tension when a Gibbs construction is considered more appropriate. The Gibbs construction also corresponds to global charge neutrality instead of local charge neutrality, appropriate for matter with two conserved charges (baryon number and charge) Glendenning (2001).

Depending on the models used to calculate the EoSs of the hadron and quark phases, chemical potential and pressure equilibrium between the two phases may not be realized Baym et al. (2018). In such cases, several interpolatory procedures have been used to connect the two phases on the premise that at $n_{\rm B}>>n_{0}$, a purely hadronic phase is physically unjustifiable Masuda et al. (2013); Fukushima and Kojo (2016); Kojo et al. (2015); Baym et al. (2018). As a result, the hadron-quark transition becomes one of a smooth crossover with the proportion of each phase depending on the specific interpolation procedure used. This is in contrast to the Gibbs construction (which also renders the transition into a mixed phase to be smooth) in which the fraction of each phase is determined self-consistently.

Although differing in details, other examples of a smooth crossover transition are the chiral model of Ref. Dexheimer et al. (2015) and the quarkyonic model of Ref. McLerran and Reddy (2019). A quark phase with additional hadronic admixtures such as hyperons and meson condensates has also been explored Prakash et al. (1997). The precise manner in which the hadron-quark transition is treated influences the magnitudes of the mass and radius of the star. In addition, the behavior of the speed of sound with density affects the magnitude of tidal deformations. It is worth mentioning however, that stars with purely hadronic matter (HM) can sometimes masquerade as stars with quark matter (QM) Alford et al. (2005).

The objectives of this work are to seek answers to probing questions such as (a) What is the minimum neutron star (NS) mass consistent with the observational lower limit on the maximum mass ($M_{\rm max}$) that is likely to contain quarks? (b) What is the minimum physically reasonable density at which a hadron-quark transition of any sort can occur? (c) Which astronomical observations have the best potential to attest to the presence of quarks?

Toward providing answers to the above questions, we have undertaken a detailed study of the hadron-to-quark matter transition in neutron stars. Our focus is to study the sensitivity of outcomes on neutron star structure, principally mass-radius relations, in the different treatments of the phase transition. Results so obtained are then subjected to the constraints provided by precise measurements of heavy neutron stars Demorest et al. (2010); Antoniadis et al. (2013); Fonseca et al. (2016); Cromartie et al. (2019), bounds on the tidal deformability of neutron stars in the binary merger event GW170817 Abbott et al. (2017a, 2018); De et al. (2018); Abbott et al. (2019), and radius estimates of $1.4\,{\rm M}_{\odot}$ available from x-ray observations of neutron stars Steiner et al. (2016); Lattimer and Prakash (2016); Özel and Freire (2016).

Earlier studies in this regard have generally chosen one favored EoS in the hadronic sector and one approach to the quark matter EoS Steiner et al. (2000); Hanauske et al. (2001); Bhattacharyya et al. (2010); Klähn and Fischer (2015); Gomes et al. (2019); Wei et al. (2018). Contrasts between the Maxwell and Gibbs constructions have also been made in some of these works, but with the result that $R_{1.4}$ are typically larger than 14 km or more (characteristic of the use of mean-field theoretical (MFT) models) which is at odds with most of the available estimates. This work differs in that variations in the EoSs of both the hadronic and quark sectors are considered as well as a global view of the outcomes of different treatments of the transition is taken. By including terms involving scalar-vector and scalar-isovector interactions in MFT models, we show that values of $R_{1.4}$ more in consonance with data can be achieved. Additionally, we present an extension of the quarkyonic matter model of Ref. McLerran and Reddy (2019) to isospin asymmetric matter with the inclusion of interactions between quarks (not considered there) to enable calculations of beta-equilibrated neutron stars. This extension will be useful in applications involving compositional and thermal gradients in quarkyonic stars, such as their long-term cooling as well as quiescent cooling following accretion on them from a companion star and in investigating $f$-, $p$- and $g$- mode oscillations. Our in-depth study of the thermodynamics of quarkyonic matter sheds additional physical insight into the role that the nucleon shell plays in stiffening the EoS.

Our findings in this work reveal that several aspects of neutron star properties deduced from observations may have to be brought to bear in finding answers to the questions posed above. These properties include the masses $M$, radii $R$, periods $P$ and their time derivatives $\dot{P}$ and $\ddot{P}$, surface temperatures $T_{s}$ of isolated neutron stars and of those that undergo periodic accretion from companions, tidal deformations $\Lambda$ from the detection of gravitational waves during the inspiraling phase of neutron star mergers, etc. Currently, the accurately measured neutron star masses around and above $2\,{\rm M}_{\odot}$ Demorest et al. (2010); Antoniadis et al. (2013); Fonseca et al. (2016); Cromartie et al. (2019) pose stringent restrictions on the EoS. Even so, the EoS would be better restricted with knowledge of radii of stars for which the masses are also known, although this would not reveal the constituents of dense matter as the structure equations depend only on the pressure vs density relation $\varepsilon(P)$, and not on how it was obtained. In contrast, the surface temperatures of both isolated neutron stars and of quiescent cooling of accreting neutron stars are sensitive to the composition, but simultaneous knowledge of their masses and radii are yet unknown. The anomalous behavior of the braking indices $n=\Omega\,\ddot{\Omega}/\dot{\Omega}^{2}$, where $\Omega=2\pi/P$ is the spin rate, of several known pulsars Magalhaes et al. (2012); Hamil et al. (2015); Johnston and Karastergiou (2017) can also be put to good use in this connection.

The organization of this paper is as follows. In Sec. II, we present the models in the hadronic and quark sectors chosen for our study. The rationale for our choice and basic features of these models are highlighted here for orientation. We stress that our choices are representative, but not exhaustive. Results of neutron star properties for different treatments of the hadron-quark transition introduced in Sec. III are shown and discussed in Sec. IV. Our conclusions and outlook are contained in Sec. V. Appendix A contains details about the thermodynamics of nucleons in the shell of quarkyonic matter.

We use units in which $\hbar=c=1$.

II Equation of State Models

Nucleonic EoSs

To explore sensitivity to the hadronic part of the EoS, we use representative examples from both potential and relativistic mean field-theoretical (RMFT) models. In the former category, the EoS of Akmal, Pandharipande and Ravehall (APR) Akmal et al. (1998), which is a parametrization of the microscopic variational calculations of Akmal and Pandharipande Akmal and Pandharipande (1997), is chosen as its energy vs baryon density up to $2\,n_{0}$ closely matches those of modern EFT calculations of pure neutron matter and symmetric nuclear matter Hebeler et al. (2010); Tews et al. (2013). Moreover, it is compatible with current nuclear phenomenology from both structure (equilibrium density and energy, compression modulus, symmetry energy and its slope, etc.) and heavy-ion experiments Danielewicz et al. (2002) as well as with the latest constraints from astrophysical observations (largest known NS mass, upper limit on maximum NS mass, tidal deformability, NS radii, etc). Explicit expressions for the energy density $\varepsilon$, pressure $P$, compression modulus $K_{0}$, Landau effective mass $m^{*}/M$, symmetry energy $S_{2}$, and the symmetry energy slope parameter $L$ along with the coupling strengths of the various terms therein can be found in Ref. Constantinou et al. (2014). Recent fits of the APR calculations to the traditional Skyrme energy-density functional (EDF) can be found in Refs. Steiner et al. (2005); Schneider et al. (2019). The latter also details the calculation of a complete tabular EoS based on the original APR parametric form.

To provide contrast, we have constructed three EoSs, MS-A, MS-B and MS-C using the RMFT model of Müller and Serot Mueller and Serot (1996) employing terms that contain scalar-isovector and vector-isovector mixings as in Refs. Horowitz and Piekarewicz (2001a, b). The numerical results to be reported in this work are from these RMFT models; that is, we consider many-body effects at the Hartree level exclusive of quantum fluctuations in the meson fields. Fock (exchange) terms are beyond the scope of this paper. As demonstrated in Ref. Chin (1977), a simple re-parametrization of the couplings in the MFT Hartree models at $T=0$ yields very nearly the same $P$ vs $\varepsilon$ relations (and, hence masses, radii and tidal deformabilities of neutron stars) as models with the inclusion of Fock terms. Note that Fock terms and additional many-body contributions do influence thermal effects in a way that is not reproducible by re-parametrization; see e.g. Refs. Zhang and Prakash (2016); Constantinou et al. (2017, 2015). In this work however, we do not consider hot matter. Specifically, we have devised three new parametrizations for the coupling constants appearing in the MS Lagrangian such that consistency with contemporary experimental and observational data is achieved. Many other EoSs based on the MS model are currently in use; for an exhaustive list, see Ref. Oertel et al. (2017). Explicitly, the Lagrangian density for this model is

[TABLE]

with

[TABLE]

Expressions for the energy per particle $\varepsilon/n$, $P$, $K_{0}$, the Dirac effective mass $M^{*}$ and hence the sigma field $\sigma_{0}=(M-M^{*})/g_{\sigma}$ in the mean-field approximation can be found in Ref. Steiner et al. (2005). With input values of these quantities at $n_{0}$, the coupling strengths $g_{\sigma}$, $g_{\omega}$, $\kappa$ and $\lambda$ are straightforwardly determined by numerically solving the system of nonlinear equations containing these quantities. The strengths $\zeta$ and $\xi$, $\Lambda_{\sigma}$ and $\Lambda_{\omega}$ of the quartic $\omega$ and $\rho$ fields, remain as adjustable input parameters to control the high-density behavior. The density-dependent symmetry energy in this model is Horowitz and Piekarewicz (2001b)

[TABLE]

The first term on the right-hand side above contains effects of interaction through $\sigma$-meson exchange, whereas the second term includes those from the $\rho$-meson exchange along with $\rho$-$\sigma$ and $\rho$-$\omega$ mixing. The corresponding slope parameter at $n_{0}$ becomes

[TABLE]

Analogous expressions but without the term involving $\Lambda_{\sigma}$ can be found in Ref. Chen and Piekarewicz (2014). The strength $g_{\rho}$ may be fixed with a prescribed value of $S_{2}$ at $n_{0}$, which leaves one or a combination of $\Lambda_{\sigma}$ and $\Lambda_{\omega}$ to obtain a desired value of $L$. The values of the various couplings used in this work are listed in Table 1.

As noted in Refs. Horowitz and Piekarewicz (2001b); Steiner et al. (2005), the quartic and scalar-isovector and vector-isovector terms in Eq. (II) enable acceptable values Lattimer and Lim (2013) of the symmetry energy slope parameter $L$ at $n_{0}$ to be obtained. The reduction in $L$ from its generally large value found for RMFT models is made possible by the second term in $L_{d}$ of Eq. (II), the term in braces being positive definite. These density-dependent terms also influence the high-density behavior of these EoSs, leaving the near-nuclear-density behavior intact. Salient properties at $n_{0}$ for these nucleonic models are presented in Table 2. The values of $L$ in Table 2 are to be compared with those of the FSU models Chen and Piekarewicz (2014); Fattoyev et al. (2018) in the literature; see e.g. Fig. 2 and Table IV in Ref. Chen and Piekarewicz (2014): $L=60.5$ MeV for FSU (but it does not achieve $2\,{\rm M}_{\odot}$) and $L=112.8\pm 16.1$ MeV for FSU2 with $M_{\rm max}=2.07\pm 0.02\,{\rm M}_{\odot}$, $R_{\rm max}=12.2$ km and $R_{1.4}=14.42\pm 0.26$ km. In comparison to FSU2, the values of $L$ for the MS models of this work are significantly smaller, which result in smaller radii for the maximum mass and $1.4\,{\rm M}_{\odot}$ neutron stars (see Table 3 below).

Properties of nucleonic neutron stars

Structural properties of charge-neutral and beta-equilibrated neutron stars resulting from the chosen EoSs are listed in Table 3. Two of the three MS EoSs satisfy the requirement of supporting a star with mass $\geq 2\,{\rm M}_{\odot}$. The EoS of MS-C does not obey the $2\,{\rm M}_{\odot}$ constraint, but we have retained it in our analysis because, in conjunction with crossover transitions involving quark matter, masses well in excess of this observational limit can be obtained (see Secs. III and IV). Although the RMFT models employ terms that contain scalar-isovector and vector-isovector mixings as in Refs. Horowitz and Piekarewicz (2001a, b) to yield acceptable values of the symmetry energy slope parameter $L$ at $n_{0}$, the radii of neutron stars stemming from these models are somewhat larger than that of the APR model, but lie within the range of those extracted from data Lattimer and Lim (2013). The largest differences between the APR and RMFT models are in the central pressures of the maximum-mass stars. The proton fractions, $y_{c,{1.4}}$ and $y_{c,{\rm max}}$, are such that stars close to the maximum-mass stars allow the direct Urca processes with electrons and muons to occur Lattimer et al. (1991).

An examination of $L$ in Table 2 and $R_{1.4}$ and $R_{\rm max}$ in Table 3 would seem to imply an anti-correlation between these quantities for the MS models. That is, smaller values of $L$ appear to lead to larger values of $R_{1.4}$ and $R_{\rm max}$, which is a trend opposite to that observed for many EoS models. The reason for this reversal becomes clear when $L$’s corresponding to different $m^{*}$’s within the same model are compared, see Fig. 1 (b) and Table 4. In other words, the standard $L-R$ correlation holds within a class of MS models with the same effective mass, whereas there exists an anti-correlation between $m^{*}-R$ which, if not taken into account explicitly, manifests itself as a turnabout in $L-R$. Similar trends of correlation with $L$ and anti-correlation with $m^{*}$ are also seen in Ref. Hornick et al. (2018) which used the MS Lagrangian but without the term involving $\Lambda_{\sigma}$ in Eq. (2) as in Ref. Chen and Piekarewicz (2014). Fig. 7 of Ref. Hornick et al. (2018) suggests that, when both $m^{*}$ and $L$ are varied, $L$ and $R$ can appear correlated, anti-correlated, or uncorrelated. The latter two possibilities are due to the competing effects of $m^{*}$ and $L$ on neutron star radii. We have verified that nonrelativistic potential models also yield similar trends, which are not shown here for brevity.

Moreover, further examination of the $P$ versus $n_{\rm B}$ and $M$-$R$ relations for the MS models (Fig. 1) shows that the central densities of $1.4\,{\rm M}_{\odot}$ stars for the EoSs chosen are all $\geq 2\,n_{0}$ with that of the MS-C star being the largest. The symmetry energy slope parameter $L$ however, refers to that at $n_{0}$. The behaviors of the pressures (see panel (a) in this figure) at $n_{\rm B}\geq 2\,n_{0}$ for all of these EoSs are distinctly different from their corresponding behaviors at $n_{\rm B}\simeq n_{0}$. The $M$-$R$ curves in Fig. 1 (b) and Table 4 also clearly show how the value of $n_{c,1.4}$ differs in each of these cases. Evidently, the manner in which the size of a $1.4\,{\rm M}_{\odot}$ is built depends sensitively on the behavior of the EoS well above $n_{0}$. These features deliver the alert that the standard $L-R_{1.4}$ correlation involves more subtleties than generally thought.

Quark EoSs

For completeness, we briefly describe the quark matter EoSs considered in this work; details can be found in the references cited. Since the discovery of $2\,{\rm M}_{\odot}$ neutron stars Demorest et al. (2010); Antoniadis et al. (2013); Fonseca et al. (2016); Cromartie et al. (2019), the traditional MIT bag Baym and Chin (1976a) and NJL Kunihiro (1989) models have been supplemented with vector interactions Klähn and Fischer (2015) to achieve consistency with data. These models have been termed vMIT, vBag, vNJL, etc., and are outlined below. Common and different features of these models will be highlighted after a brief description of each model.

The bag model and its variations

The Lagrangian density of the MIT bag model is Baym and Chin (1976a)

[TABLE]

which describes quarks of mass $m_{i}$ confined within a bag as denoted by the $\Theta$ function. For three flavors $i=u,d,s$ and three colors $N_{c}=3$ of quarks, the number and baryon densities, energy density, pressure and chemical potentials in the bag model are Baym and Chin (1976a)

[TABLE]

The superscript $k_{Fi}$ in the integral signs is the Fermi wave number for each species $i$, at which the integration over $k$ is terminated at zero temperature. The first terms in $\varepsilon_{Q}$ and $P_{Q}$ are free Fermi gas contributions, $\varepsilon_{\rm FG}$ and $P_{\rm FG}$, respectively, the second terms are QCD perturbative corrections due to gluon exchange corresponding to ${\cal L}_{\rm int}$, and $B$ is the so-called bag constant which accounts for the cost in confining the quarks into a bag. The quark masses $m_{i}$ are generally taken to be current quark masses. Often, the $u$ and $d$ quark masses are set to zero (as at high density, $k_{Fi}$ in these cases far exceed $m_{i}$), whereas that of the $s$ quark is taken at its Particle Data Group (PDG) value. Refs. Baym and Chin (1976a, b); Freedman and McLerran (1978); Farhi and Jaffe (1984); Kurkela et al. (2010); Kurkela and Vuorinen (2016); Gorda et al. (2018) detail the QCD perturbative calculations of $\varepsilon_{\rm pert}$ and $P_{\rm pert}$, and the ensuing results for the structure of neutron stars containing quarks within the cores as well as self-bound strange quark stars. At leading order of QCD corrections, the results are qualitatively similar to what is obtained by just using the Fermi gas results with an appropriately chosen value of $B$ Prakash et al. (1990).

In recent years, variations of the bag model have been adopted Klähn and Fischer (2015); Gomes et al. (2019); Xia et al. (2019) to calculate the structure of neutron stars with quarks cores to account for $\geq 2\,{\rm M}_{\odot}$ maximum-mass stars. Termed as vMIT or vBag models, the QCD perturbative results are dropped and replaced by repulsive vector interactions between quarks in such works. We will provide some numerical examples of the vMIT model for contrast with other models as those of the vBag model turn out to be qualitatively similar.

The vMIT model

The form $\mathcal{L}_{\rm int}=-G_{v}\sum_{i}~{}\overline{\psi}\gamma_{\mu}V^{\mu}\psi+(m_{V}^{2}/2)V_{\mu}V^{\mu}$, where interactions among the quarks occur via the exchange of a vector-isoscalar meson $V^{\mu}$ of mass $m_{V}$, is chosen in Ref. Gomes et al. (2019). Here, the quark masses are chosen close to their current quark masses. Explicitly,

[TABLE]

where $n_{Q}=\sum_{i}n_{i}$, and the bag constant $B$ is chosen appropriately to enable a transition to matter containing quarks. Note that terms associated with the vector interaction above are similar to those in hadronic models. In the results reported below, we vary the model parameters in the range $B^{1/4}=(155-180)$ MeV and $a=(G_{v}/m_{V})^{2}=(0.1-0.3)~{}{\rm fm^{2}}$.

The vNJL model

In its commonly used form, the Lagrangian density for the vNJL model in the mean field approximation is

[TABLE]

Here, $q$ denotes a quark field with three flavors $u,d,s$ and three colors, $\hat{m_{0}}$ is the 3$\times$3 diagonal current quark mass matrix, $\lambda_{k}$ represents the 8 generators of SU(3), and $\lambda_{0}$ is proportional to the identity matrix. The four-fermion interactions are from the original formulation of this model Nambu and Jona-Lasinio (1961), whereas the flavor mixing, determinental interaction is added to break the $U_{A}(1)$ symmetry ’t Hooft (1986). The last term accounts for vector interactions Buballa (2005). As the constants $G_{s}$, $K$, and $G_{v}$ are dimensionful, the quantum theory is non-renormalizeable. Therefore, an ultraviolet cutoff $\Lambda$ is imposed, and results are meaningful only for quark Fermi momenta well below this cutoff.

The Lagrangian density in Eq. (8) leads to the energy density

[TABLE]

where the sums above run over $u,d,s$. The subscript $``0"$ denotes current quark masses and the superscript $\Lambda$ in the integral sign indicates that an ultraviolet cutoff $\Lambda$ is imposed on the integration over $k$. In both $\varepsilon_{\rm FG}$ [see Eq. (6)] and $\varepsilon_{\rm int}$, the quark masses $m_{i}$ are dynamically generated by requiring that $\varepsilon$ be stationary with respect to variations in the quark condensate $\langle\overline{q}_{i}q_{i}\rangle$:

[TABLE]

$(q_{i},q_{j},q_{k})$ representing any permutation of $(u,d,s)$. The quark condensate $\langle\overline{q}_{i}q_{i}\rangle$ is given by

[TABLE]

and the quark number density $n_{i}=\langle q_{i}^{\dagger}q_{i}\rangle$ is as in Eq. (6). Note that the integrals appearing in Eqs. (9)-(11) can all be evaluated analytically. Eqs. (10) and (11) render the dynamically generated masses $m_{i}$ density dependent, which tend to $m_{0,i}$ at high density mimicking the restoration of chiral symmetry in QCD.

To facilitate a comparison between the vMIT and vNJL models, Ref. Buballa and Oertel (1999) recommends a constant energy density $B_{0}=\varepsilon_{\rm int}|_{m_{u}=m_{d}=m_{s}=0}$ to be added to $\varepsilon_{\rm int}$ which makes the vacuum energy density zero. With this addition, the energy density takes the form $\varepsilon=\varepsilon_{\rm FG}+B_{\rm eff}$, with $B_{\rm eff}=B_{0}+\varepsilon_{\rm int}$.

The quark chemical potentials are

[TABLE]

using which the pressure is obtained from the thermodynamic identity

[TABLE]

To mimic confinement absent in the vNJL model, often a constant term $B_{\rm dc}$ is used with the replacement $B_{\rm eff}\rightarrow B_{\rm eff}-B_{\rm dc}$.

For numerical calculations, we use the HK parameter set Hatsuda and Kunihiro (1994): $\Lambda=631.4$ MeV, $G_{s}\Lambda^{2}=1.835$, $K\Lambda^{5}=0.29$, $m_{u,d}=5.5$ MeV, $m_{s}=135.7$ MeV and $B_{\rm dc}=0$.

The vBag model

In Ref. Klähn and Fischer (2015), vector interactions are used in the form of flavor-independent four-fermion interactions as in the NJL models (described below): $\mathcal{L}_{\rm int}=G_{v}\sum_{i}~{}(\overline{\psi}\gamma^{\mu}\psi)^{2}$. In this case Klähn and Fischer (2015),

[TABLE]

where the explicit forms of $\varepsilon_{\rm FG,i}$ and $P_{\rm FG,i}$ can be read off from Eq. (6). The effective bag constant $B_{\rm eff}$ in this model is composed of two parts: $B_{\rm eff}=\sum_{i}B_{\chi}^{i}-B_{\rm dc}$, where the flavor-dependent chiral bag constant

[TABLE]

where $m_{i}$ is the dynamically generated quark mass as in the NJL model, $m_{0,i}$ is the current quark mass, and $\Lambda$ is an ultraviolet cut-off on the integration over $k$. The quantity $B_{\rm dc}$ is tuned to control the onset of quark deconfinement.

Charge neutrality and beta-equilibrium conditions

Equilibrium with respect to weak-interaction processes $d\rightarrow ue^{-}\bar{\nu}_{e}$ and $s\rightarrow ue^{-}\bar{\nu}_{e}$ leads to the chemical potential equalities $\mu_{d}=\mu_{u}+\mu_{e}=\mu_{s}$ in neutrino-free matter. Charge neutrality requires that $2\,n_{u}-n_{d}-n_{s}-3\,n_{e}=0$. Together with the baryon number relation $n_{u}+n_{d}+n_{s}=3\,n_{\rm B}$, the simultaneous solution of the equations

[TABLE]

assures that quark matter with the three flavors $u,d,s$ is charge neutral and in beta equilibrium. In Eq. (II), $x_{i}=n_{i}/n_{\rm B}$ denote the particle concentrations, $r=n_{\rm B}/n_{0}$, and the factor $C_{i}={m_{i}^{2}}/({\pi^{4/3}n_{0}^{2/3}})$, $i=u,d,s$. Note that $C_{i}$ can depend on the density compression ratio $r$ through $m_{i}\equiv m_{i}(r)$ as in the vNJL model. The concentrations of the $u$ and $d$ quarks are given by

[TABLE]

respectively. Owing to the charges carried by the quarks, the electron concentration in quark matter is generally very small with increasing $r$.

Distinguishing features of the quark EoSs

The vMIT and vNJL models differ in important ways. Fashioned after the MIT bag for the nucleon, the vMIT model incorporates overall confinement of quarks within a giant bag Baym and Chin (1976a, b) through its density-independent (non-perturbative) bag constant $B$. Repulsive vector interactions in this model are $\propto n_{Q}^{2}$ with $n_{Q}=\sum n_{u,d,s}$ in the pressure and energy density. The kinetic energy is calculated with current quark masses, although use of constituent masses $m_{u,d,s}=M_{n}/3$ can also be found in the literature. Effects of interactions are included from perturbative QCD calculations, but often they are set to zero in favor of an altered value of $B$ to simulate the same effect.

The most important and distinguishing feature of the mean-field vNJL model is the chiral restoration of the quark masses present in the original QCD Lagrangian. Starting from the dynamically generated quark masses $m_{u}=m_{d}\simeq 350~{}{\rm MeV}$ and $m_{s}\simeq 525~{}{\rm MeV}$ in vacuum, the masses decrease steadily toward their current quark values with increasing density. In our numerical calculations, we have used the current quark masses $m_{0,u}=m_{0,d}=5$ MeV and $m_{0,s}=140$ MeV to conform to the values used in Ref. Buballa and Oertel (1999) (use of the current PDG values $m_{u}\simeq 2~{}{\rm MeV}$, $m_{d}\simeq 5~{}{\rm MeV}$, and $m_{s}\simeq 100~{}{\rm MeV}$ does not significantly affect the results). In addition to generating the quark masses, the scalar field energies involving the couplings $G_{s}$ and $K$ also enter the energy density (and hence the pressure). Vector interactions in the vNJL model are $\propto\sum_{i}n_{i}^{2}$, $i=u,d,s$; this differs from the vMIT model in that cross terms such as $n_{j}n_{k}$ are absent in the former case. The vNJL model lacks confinement, although a constant $B_{0}$ is added to the energy density so that $B_{\rm eff}=B_{0}+\varepsilon_{\rm int}$ to facilitate comparison with the $B$ of the vMIT bag model. $B_{\rm eff}$ is however density dependent, unlike the $B$ of the vMIT model.

In short, both the vMIT and vNJL models incorporate some aspects of the QCD Lagrangian, but only partially. Lacking a truly non-perturbative approach to QCD, we have explored both models as representative of the current status.

Note that in the vBag model, $B_{\chi}^{i}$ and $B_{\rm dc}$, and thus $B_{\rm eff}$, are independent of density. Unlike in the vNJL model in which all terms in the energy density and pressure are calculated with density-dependent dynamical masses $m_{i}$, the Fermi gas contributions in the vBag model are calculated with $m_{i}(k_{Fi}=0)$.

The striking similarity of the vMIT and vBag models is worthy of an explicit discussion. For the purpose of comparison, we can impose $a^{\rm vMIT}=G_{v}^{\rm vBag}$ (numerically), $B=B_{\rm eff}$ and set the quark masses the same. Then, the difference in the vector-interaction terms in $\varepsilon_{Q}$ and $P_{Q}$ in Eqs. (7) and (14) becomes apparent. Specifically, those of vMIT are $\propto n_{Q}^{2}$ with $n_{Q}=\sum_{i}n_{i}$ whereas those of vBag are $\propto\sum_{i}n_{i}^{2}$. These differences are caused by the associated terms in the respective Lagrangians. The corresponding terms in the chemical potentials will be $a^{\rm vMIT}n_{Q}$ for vMIT and $a^{\rm vMIT}n_{i}$ for vBag. When charge neutrality and beta equilibrium are imposed, even the Fermi gas parts in the two models will be different as the corresponding Fermi momenta will be different. Thus, although the two models look similar they are different because of the way interactions are treated.

Consequences of the vBag model on neutron star structure have been studied extensively in Refs. Klähn and Fischer (2015); Wei et al. (2018) (and in this work), and will not be repeated here.

III Treatment of Phase Transitions

First-order transitions

The manner in which the hadron-quark transition occurs is unknown. Even if the phase transition is assumed to be of first-order, description of the transition depends on the knowledge of the surface tension $\sigma_{s}$ between the two phases Alford et al. (2001); Mintz et al. (2010); Lugones et al. (2013); Fraga et al. (2019). In view of uncertainties in the magnitude of $\sigma_{s}$, two extreme cases have been studied in the literature.

Maxwell Construction

For very large values of $\sigma_{s}$, a Maxwell construction in which the pressure and chemical potential equalities, $P(H)=P(Q)$ and $\mu_{n}(H)=\mu_{n}(Q)$, are established between the two phases, hadronic (H) and quark (Q), has been deemed appropriate. In charge-neutral and beta-equilibrated matter, only one baryon chemical potential, often chosen to be $\mu_{n}$, is needed to conserve baryon number as local charge neutrality is implicit. The range of densities over which these equalities hold can be found using the methods detailed in Refs.  Constantinou et al. (2014); Lamb et al. (1983).

Gibbs Construction

For very low values of $\sigma_{s}$, a Gibbs construction in which a mixed phase of hadrons and quarks is present is more appropriate Glendenning (1992, 2001). The description of the mixed phase is achieved by satisfying Gibbs’ phase rules: $P(H)=P(Q)$ and $\mu_{n}=\mu_{u}+2\mu_{d}$. Further, the conditions of global charge neutrality and baryon number conservation are imposed through the relations

[TABLE]

where $f$ represents the fractional volume occupied by hadrons and is solved for at each $n_{\rm B}$. Unlike in the pure phases of the Maxwell construction, $Q(H)$ and $Q(Q)$ do not separately vanish in the Gibbs mixed phase. The total energy density is $\varepsilon=f\,\varepsilon(H)+(1-f)\,\varepsilon(Q)$. Relative to the Maxwell construction, the behavior of pressure vs density is smooth in the case of Gibbs construction. Discontinuities in its derivatives with respect to density, reflected in the squared speed of sound $c_{s}^{2}=dP/d\varepsilon$, will however be present at the densities where the mixed phase begins and ends.

The Maxwell and Gibbs constructions represent extreme cases of treating first-order phase transitions, and reality may lie in between these two cases. However, there are situations in which neither method can be applied as the required pressure and chemical potential equalities cannot be met for many hadronic and quark EoSs Baym et al. (2018). In such cases, an interpolatory method which makes the transition a smooth crossover has been used Masuda et al. (2013); Fukushima and Kojo (2016); Kojo et al. (2015); Baym et al. (2018); Li et al. (2018).

Crossover transitions

As it is not clear that a first-order phase hadron-to-quark transition at finite baryon density is demanded by fundamental considerations, crossover or second-order transitions have also been explored recently; see e.g. Refs. Baym et al. (2018); Dexheimer et al. (2015); McLerran and Reddy (2019). As details of results ensuing from the model of Ref. Dexheimer et al. (2015) have been recorded earlier in Refs. Dexheimer and Schramm (2010); Negreiros et al. (2010), we will only examine the cases of interpolated and quarkyonic models in what follows.

Interpolated EoS

We follow the simple recipe in Ref. Masuda et al. (2013) where the interpolated EoS in the hadron-quark crossover region is characterized by its central value $\bar{n}$ and width $2\,\Gamma$. Pure hadronic matter exists for $n\lesssim\bar{n}-\Gamma$, whereas a phase of pure quark matter is found for $n\gtrsim\bar{n}+\Gamma$. In the crossover region, $\bar{n}-\Gamma\lesssim n\lesssim\bar{n}+\Gamma$, strongly interacting hadrons and quarks coexist in prescribed proportions. The interpolation is performed for pressure vs baryon number density according to

[TABLE]

where $P_{H}$ and $P_{Q}$ are the pressure in pure hadronic and pure quark matter, respectively. The interpolated EoS for the crossover, Eq. (19), is different from that of the Gibbs construction within the conventional picture of a first-order phase transition in that the pressure equality between the two phases has been abandoned. Also, $f_{-}$ and $f_{+}$ are not solved for, but chosen externally. (Alternative forms of interpolation have also been suggested in Refs. Fukushima and Kojo (2016); Kojo et al. (2015), but do not qualitatively change the outcome.) The energy density $\varepsilon$ vs $n$ is obtained by integrating $P=n^{2}\partial(\varepsilon/n)/\partial n$:

[TABLE]

Quarkyonic matter

The transition to matter containing quarks in the model termed quarkyonic matter McLerran and Pisarski (2007); McLerran and Reddy (2019) is of second or higher order, depending on the behavior of the squared speed of sound $c_{s}^{2}=dP/d\varepsilon=d\ln\mu/d\ln n_{\rm B}=(n/\mu)(d^{2}P/d\mu^{2})^{-1}$ with $n_{\rm B}$. The order of the phase transition is not determined by the quarkyonic matter scenario a priori, but depends on its specific implementation. In the model proposed in Ref. McLerran and Reddy (2019), $c_{s}^{2}$ exhibits a kink at the onset of the transition, hence its derivative with respect to $n_{\rm B}$ is discontinuous. It is in this sense that the transition is of second order for Ref. McLerran and Reddy (2019) which may not be the case in other implementations of the quarkyonic matter scenario. This model is a departure from the first-order phase transition models insofar as once quarks appear, both nucleons and quarks are present until asymptotically large densities when the nucleon concentrations vanish. Keeping the structure of the quarkyonic matter model as in Refs. McLerran and Pisarski (2007); McLerran and Reddy (2019) in which isospin symmetric nuclear matter (SNM) and pure neutron matter (PNM) were considered, we present below its generalization to charge-neutral and beta-equilibrated neutron star matter. In quarkyonic matter, the appearance of quarks is subject to the threshold condition McLerran and Reddy (2019)

[TABLE]

where $k_{FB}$ is the baryon momentum, $N_{c}=3$ is the number of colors, and the momentum threshold $\Delta$ is chosen to be

[TABLE]

Above, $\Lambda_{Q}\sim\Lambda_{\rm QCD}\simeq 300-500$ MeV, and $\kappa\simeq 0.1-0.3$ is suitably chosen to preserve causality. In PNM, the transition density, $n_{\rm trans}$, for the appearance of quarks is $0.77\,(3.55)\,n_{0}$ for $\Delta=300\,(500)$ MeV and $\kappa=0.3$, where $n_{0}$ is the SNM equilibrium density. The corresponding values for $\kappa=0.1$ are $0.75\,n_{0}$ and $3.47\,n_{0}$, respectively, and show weak dependence of $n_{\rm trans}$ on $\kappa$. Unlike in the other approaches, the transition density at which quarks begin to appear in this model is independent of the EoSs used in the hadronic and quark sectors, being dependent entirely on $\Lambda_{\rm QCD}$ and large $N_{c}$ physics.

The total baryon density of quarkyonic matter is

[TABLE]

Notice that once quarks appear, the shell width $\Delta$ in which nucleons reside decreases with density as $n_{\rm B}^{-2/3}$, yielding the preponderance of quarks with increasing $n_{\rm B}$. Including leptonic (electron and muon) contributions $\varepsilon_{\ell}$, the total energy density is

[TABLE]

where $e_{k}$ is the single particle kinetic energy inclusive of the rest mass energy. The nucleonic part of the energy density for $n\gtrsim 0.5\,n_{0}$ can be taken from a suitable potential or field-theoretical model that is constrained by nuclear systematics near nuclear densities, and preserves causality at high densities. Below $0.5\,n_{0}$, the energy density is that of crustal matter as in e.g. Refs. Negele and Vautherin (1973); Baym et al. (1971). It is important to realize that the term $\varepsilon_{\rm int}(n_{n},n_{p})$ contributes in regions where $k_{FB}<\Delta$ as well as where $k_{FB}>\Delta$.

The chemical potentials and pressure are obtained from

[TABLE]

where the sum above runs over all fermions.

As with nucleons, an appropriate choice of the quark EoS is also indicated. Reference McLerran and Reddy (2019) set $\varepsilon_{\rm int}(q_{k},q_{\ell})=0$, and the quark masses $M_{q}$ were taken as $M_{n}/3$. The use of the nucleon constituent quark masses takes account of quark-gluon interactions to a certain degree as has been noted in the case of finite temperature QCD as well. This procedure however, omits density-dependent contributions from interactions between quarks. In our work, we will employ quark models (such as vMIT, vNJL) in which contributions from interacting quarks are included. Subtleties involved in the calculation of the kinetic part of the nucleon chemical potentials and in satisfying the thermodynamic identity are detailed in Appendix A.

This model has a distinct behavior for $c_{s}^{2}=dP/d\varepsilon$ vs $n_{\rm B}$ in that $c_{s}^{2}$ exhibits a maximum (its location controlled by $\Lambda_{Q}$, and the magnitude depending on both $\Lambda_{Q}$ and $\kappa$) before approaching the value of 1/3 characteristic of quarks at asymptotically large densities McLerran and Reddy (2019).

IV Results with Phase Transitions

The hadronic EoSs chosen in this study satisfy the available nuclear systematics near the nuclear equilibrium density (see Tables 1-3). Their supra-nuclear density behavior can however, be varied to yield a soft or stiff EoS by varying the parameters in the chosen model. Depending on the quark EoS examined such as vMIT, vNJL or of quarkyonic matter, the examination of a broad range of transitions into quark matter - soft-to-soft, soft-to-stiff, stiff-to soft and stiff-to-stiff - become possible. For both first-order and crossover transitions, we calculate the mass-radius curves and tidal deformabilities, and then discuss the results in view of the existing observational constraints. Of particular relevance to the zero-temperature EoS is the limit set by the data on the binary tidal deformability Flanagan and Hinderer (2008); Favata (2014)

[TABLE]

For each star, the dimensionless tidal deformability (or induced quadrupole polarizability) is given by Love (1909)

[TABLE]

where the second tidal Love number $k_{2}^{(1,2)}$ depends on the structure of the star, and therefore on the mass and the EoS. Here $G$ is the gravitational constant and $R_{1,2}$ are the radii. The computation of $k_{2}^{(1,2)}$ with input EoSs is described in Refs. Thorne and Campolattaro (1967); Hinderer (2008); Damour and Nagar (2009). For a wide class of neutron star EoSs, $k_{2}\simeq 0.05-0.15$ Hinderer (2008); Hinderer et al. (2010); Postnikov et al. (2010).

Combining the electromagnetic (EM) Abbott et al. (2017b) and gravitational wave (GW) information from the binary neutron star (BNS) merger GW170817, Ref. Margalit and Metzger (2017) provides constraints on the radius $R_{\rm ns}$ and maximum gravitational mass $M^{g}_{\rm max}$ of a neutron star:

[TABLE]

where $R_{1.3}$ is the radius of a $1.3\,{\rm M}_{\odot}$ neutron star and its numerical value above corresponds to $M^{g}_{\rm max}=2.17\,{\rm M}_{\odot}$. These estimates have been revisited in a recent analysis of Ref. Shibata et al. (2019) where a weaker constraint on the upper limit of the maximum mass $M^{g}_{\rm max}\lesssim 2.3\,{\rm M}_{\odot}$ has been reported. Combining the total mass measurement of $2.74^{+0.04}_{-0.01}\,{\rm M}_{\odot}$ from GW170817 with (empirical) universal relations between the baryonic and the maximum rotating and non-rotating masses of neutron stars, Ref. Rezzolla et al. (2018) constrains the maximum non-rotating neutrons star mass in the range $2.01^{+0.04}_{-0.01}\,{\rm M}_{\odot}\leq M_{\rm max}^{\rm nonrot}\leq 2.16^{+0.17}_{-0.15}\,{\rm M}_{\odot}$.

First-order transitions: Maxwell vs. Gibbs

We first survey the allowed parameter space for valid first-order phase transitions, namely, a critical pressure exists above which quark matter is energetically favored. We then proceed with both Maxwell and Gibbs constructions, calculating quantities to be compared with observational constraints. Our results are summarized in Figs. 2-4. Where possible, we also characterize the behavior of the hadron-to-quark transition with quantities introduced in the “Constant-Sound-Speed (CSS)” approach in Ref. Alford et al. (2013). This approach can be viewed as the lowest-order Taylor expansion of the high-density EoS about the transition pressure $P_{\rm trans}$, by specifying the discontinuity in energy density $\Delta\varepsilon$ at the transition, and the density-independent squared sound speed $c^{2}_{\rm QM}$ in quark matter. This generic parametrization has been widely used in recent studies on the manifestation of a first-order phase transition with Maxwell construction in neutron star phenomenology, see e.g. Paschalidis et al. (2018); Burgio et al. (2018); Christian et al. (2019); Montana et al. (2019). Despite different choices of the baseline hadronic EoS, comparison between separate works is afforded by mapping onto the CSS parameter space. To facilitate such a comparison, we list in Table 5 the corresponding CSS parameter values for calculations from our physically based models.

MS-A + vMIT (stiff $\to$ soft/stiff)

Fixing the hadronic EoS to be the stiff model MS-A, we choose in the vMIT model six parameter sets of ($B^{1/4},a$), where $B$ is the bag constant and $a=(G_{v}/m_{v})^{2}$ measures the strength of vector interactions between quarks. The bag constant is adjusted so that the transition to quark matter occurs at $n_{\rm trans}=1.5-2.4\,n_{0}$, and the finite vector coupling $a$ stiffens the quark matter EoS. Soft hadronic EoSs are not applied, as they either (with softer quark EoSs) violate the $M_{\rm max}\geq 2\,{\rm M}_{\odot}$ constraint or (with stiffer quark EoSs) cannot establish a valid first-order phase transition, i.e., there is no intersection between the two phases in the $P$-$\mu$ plane. We note that this limitation (hadronic matter being stiff) does not necessarily hold if a generic parameterization such as CSS is utilized instead of specific quark models to perform first-order transitions.

In the vMIT model, the sound speed varies little even with the inclusion of vector repulsive interactions within the star (see Fig. 2 (b) and (e)) and can be approximated as being density independent. The mass-radius topology with the Maxwell construction is determined by the three parameters ($P_{\rm trans}/\varepsilon_{\rm trans}$, $\Delta\varepsilon/\varepsilon_{\rm trans}$, $c^{2}_{\rm QM}$) in CSS, giving rise to either connected, disconnected (i.e. twin stars or third-family stars) or both branches of stable hybrid stars; $P_{\rm trans}$ and $\varepsilon_{\rm trans}$ are the pressure and energy density in hadronic matter at the transition, respectively, $\Delta\varepsilon$ is the discontinuity in energy density at $P_{\rm trans}$, and $c^{2}_{\rm QM}$ is the squared speed of sound in quark matter just above $P_{\rm trans}$. The threshold value $\Delta\varepsilon_{\rm crit}$ below which there is always a stable hybrid branch connected to the purely-hadronic branch is given by $\Delta\varepsilon_{\rm crit}/\varepsilon_{\rm trans}=\frac{1}{2}+\frac{3}{2}P_{\rm trans}/\varepsilon_{\rm trans}$ Seidov (1971); Schaeffer et al. (1983); Lindblom (1998). The relevant quantities for the mapping between the stiff MS-A+vMIT model (Maxwell) and the CSS parametrization are listed in Table 5.

After extensively varying all parameters and calculating the corresponding mass-radius relations, we find that $a=0.18$ is most likely the smallest value (corresponding to $c^{2}_{\rm QM}\approx 0.4$) that barely ensures $M_{\rm max}\simeq 2\,{\rm M}_{\odot}$. When $a$ is increased from zero, the energy density discontinuity becomes progressively smaller ($\Delta\varepsilon/\varepsilon_{\rm trans}\lesssim 0.5$) and eventually the twin-star solutions disappear, roughly at $a\geq 0.15$. Within the range $a=0.18-0.3$, the $M(R)$ curves of stable hybrid stars obtained are continuous, and quarks can appear at $1.0\leq M_{\rm trans}\leq 1.8\,{\rm M}_{\odot}$, pertinent to the range of component masses in BNS mergers. For too large vector interaction couplings e.g. $a=0.5$, the onset for quarks is beyond the central density of the maximum-mass hadronic star, and thus no stable quark cores would be present even though QM is sufficiently stiff.

Fig. 2 (c) shows that requiring $M_{\rm max}\geq 2\,{\rm M}_{\odot}$ excludes certain twin-star solutions obtained from EoSs with zero (gray dash-dot-dotted) and small (orange dotted) repulsive vector interactions between quarks, mainly due to the insufficient stiffness of the quark matter EoS. By invoking very stiff EoSs with $c_{s}^{2}\rightarrow 1$ in the quark sector and using the CSS parametrization coupled with hadronic EoSs at low density, recent works have reported twin stars compatible with the constraint $M_{\rm max}\geq 2\,{\rm M}_{\odot}$ and bounds on $\tilde{\Lambda}$ from GW170817 (see e.g. Refs. Christian et al. (2019); Sieniawska et al. (2019); Burgio et al. (2018); Han and Steiner (2019); Montana et al. (2019) and references therein). Moreover, the typical neutron star radius $R_{1.4}$ can be observationally constrained by radius estimates from x-ray emission and/or tidal deformability ($\Lambda$) measurements in pre-merger gravitational-wave detections. For hybrid EoSs with a sharp phase transition, the value of $R_{1.4}$ or $\Lambda_{1.4}$ is sensitive to the onset density $n_{\rm trans}$, above which $M(R)$ and $\Lambda(M)$ deviate from normal hadronic EoSs without a sharp transition. We demonstrate this effect in Fig. 3 by confronting calculated tidal deformations with inferred bounds from the first BNS event GW170817 Abbott et al. (2017a, 2018); De et al. (2018).

With high accuracy, the chirp mass ${\mathcal{M}}=(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$, where $m_{1,2}$ are the masses of the merging neutron stars, was determined to be ${\mathcal{M}}=1.186_{-0.001}^{+0.001}{\rm M}_{\odot}$ in GW170817 Abbott et al. (2019). This event also revealed information on the binary tidal deformability, $\tilde{\Lambda}({\mathcal{M}}=1.186_{-0.001}^{+0.001}{\rm M}_{\odot})=300_{-230}^{+420}$ for low-spin priors (using a 90% highest posterior density interval). Furthermore, by assuming a linear expansion of $\Lambda(M)$, which holds fairly well for normal hadronic stars without sharp transitions, limits on the dimensionless tidal deformability of a $1.4\,{\rm M}_{\odot}$ NS were derived Abbott et al. (2018): $70\leq\Lambda_{1.4}\leq 580$ for low spin priors (at $90\%$ confidence level). This single detection of GW170817 rules out purely-hadronic EoSs that are too stiff and correlated with large tidal deformabilities, as shown in Fig. 3 (b) and (c). The stiff MS-A model by itself is incompatible with the estimated ranges of $\Lambda_{1.4}$ and $\tilde{\Lambda}$. The only solution to rescue such a stiff hadronic EoS is to introduce a phase transition at not-too-high densities, e.g., a possible smaller $\Lambda$ can be achieved in a hybrid star that already exists in the pre-merger stage. For Maxwell constructions, one of the six parameter sets, $(B^{1/4},a)=(159,0.2)$ (blue dash-dotted) with $n_{\rm trans}/n_{0}=1.77$ (see Table 5) is successful to survive the LIGO constraint. Together with the maximum-mass constraint, the parameter space for sharp phase transitions is severely limited.

The panels (e)-(f) of Figs. 2-3 represent results for the stiff MS-A+vMIT model with Gibbs constructions, for which the model parameters remain the same as in their Maxwell counterparts (panels (a)-(c)). The smooth feature of the Gibbs construction advances the appearance of quarks in the mixed phase to lower densities, while it defers the region of the purely quark phase to higher densities. These features are also manifested in the corresponding $\varepsilon(P)$ relation and its finite speed-of-sound behavior (Fig. 2 (e) and (f)). Effectively, the softening due to (Gibbs) phase transition occurs earlier, smoothly decreasing the NS radii and tidal deformabilities for a broader range of masses, which gives rise to increased compatibility with observational constraints. Three more parameter sets of the stiff MS-A+vMIT model that satisfy $M_{\rm max}\geq 2\,{\rm M}_{\odot}$ are now consistent with the tidal deformability constraint (Fig. 3 (e) and (f)), in contrast to the only candidate that qualifies in Maxwell constructions. In this respect, applying Gibbs construction is advantageous to enlarging the quark model parameter space that suitably satisfies the current constraints from observation (and also revives previously-excluded stiff hadronic models). However, the clear-cut distinction between hybrid and purely-hadronic branches in terms of $M(R)$ and $\Lambda(M)$ diminishes: the drastic effect from a sharp hadron/quark transition is toned down, and thus distinguishability of quarks with regard to global observables becomes less feasible if they take the form of a mixture with hadrons. This feature accentuates the significance of dynamical properties such as NS cooling and spin-down, and the evolution of merger products.

MS-A + vNJL (stiff $\to$ soft)

In the vNJL model, pressures at $\lesssim 2\,n_{0}$ exhibit an unphysical behavior (of being negative and/or decreasing with density) which forbids attempts to shift $n_{\rm trans}$ to low densities. If a finite vector coupling $G_{v}$ is introduced, the onset of quarks is typically reached at $n_{\rm trans}\gtrsim 2.3\,n_{0}$ ($M_{\rm trans}\gtrsim 1.7\,{\rm M}_{\odot}$), leading to a short stable hybrid branch that obeys $M_{\rm max}\geq 2\,{\rm M}_{\odot}$ because of the stiff hadronic EoS. We display one such example in Fig. 4 for both Maxwell and Gibbs constructions. Note that the speed of sound in the quark phase remains small, restricted by the fact that a too large $G_{v}$ (correlated with stiffer QM) significantly delays the onset for quarks which yields no stable hybrid stars. Some relevant points to note are:

(i) $M_{\rm trans}\gtrsim 1.7\,{\rm M}_{\odot}$ indicates that most likely there will be no quarks in e.g., the component neutron stars of a binary before they coalesce. Thus, tidal properties are not shown in Fig. 4 due to the high onset density for quarks: i.e., in this case BNS observables are irrelevant;

(ii) A small $G_{v}$ has little effect on stiffening quark matter ($c^{2}_{\rm QM}\lesssim 1/3$), which is not desirable in terms of supporting $2\,{\rm M}_{\odot}$ mostly by quarks; and

(iii) Gibbs construction helps maintaining slightly more quark content than Maxwell in the most massive stars, but quarks are effectively “invisible” even if they exist.

Note that the tidal deformability constraint rules out a very stiff hadronic EoS, e.g., MS-A. This stiffness in the hadronic EoS is nevertheless a prerequisite for vNJL to construct a valid first-order transition; stable hybrid stars that are consistent with observation do not exist in this scenario. There is no solution other than an alternative treatment, such as a crossover transition to which we turn below. It is noteworthy that the only successful scenario we find for first-order phase transitions to be compatible with observations is a stiff HM $\to$ stiff QM transition (see summary in Table 7). This conclusion agrees fairly well with those from other previous studies in which specific models of quark matter were used, see e.g. Nandi and Char (2018); Paschalidis et al. (2018); Alvarez-Castillo et al. (2019); Gomes et al. (2019).

Crossover transitions: Interpolatory procedures and quarkyonic matter

In obtaining the results shown below in Figs. 5 and 6, we have followed the methods detailed in Sec. III for constructing crossover hadron-to-quark transitions. Although the generalization of the quarkyonic matter model to beta-equilibrated stars is presented in that section, results shown here for this case are for pure neutron matter only, in order to provide a direct comparison with the results of Ref. McLerran and Reddy (2019).

Interpolated EoSs

The results shown for this case correspond to a smooth interpolation in the window $(\bar{n},\Gamma)=(3\,n_{0},n_{0})$ between the soft hadronic EoS MS-B and stiff quark EoSs in the vNJL model with $G_{v}/G_{s}=1.5,2.0,2.5$. Outside this window in density, pure hadronic and pure quark phases are expected to exist. Due to the abrupt cutoff imposed in the boundary condition, there is a finite jump in $c_{s}^{2}$ at the lower-end of the crossover window $\bar{n}-\Gamma=2\,n_{0}\equiv n_{\rm trans}$ below which only pure hadronic phase is present. At the higher-end and above, we continue to use the interpolated form. This is different from Ref. Masuda et al. (2013), where the interpolated form extended to all densities. As we will see below, the cutoffs are important to typical radii and thus could be significant.

Effects of introducing quarks above $n_{\rm trans}=2\,n_{0}$ through smooth interpolations in the EoS are shown in Fig. 5 (a)-(c). The maximum mass is primarily determined by the stiffness above $4\,n_{0}$, hence the use of large vector-coupling strengths in vNJL. Consequently, one can derive a constraint on $G_{v}/G_{s}$ from $M_{\rm max}\geq 2\,{\rm M}_{\odot}$ if other parameters are fixed, e.g., $G_{v}/G_{s}=1.5$ is probably ruled out. On the other hand, typical radii for $1.0-1.6\,{\rm M}_{\odot}$ stars are sensitive to the stiffness in the hadronic phase at $n_{\rm B}\lesssim 2\,n_{0}$, as well as to the choice of the threshold density. For instance, we have found that for $n_{\rm trans}=1.5\,n_{0}$ instead of $2\,n_{0}$, $R_{1.4}$ decreases by about $0.3$ km. Note that the hyperbolic construction results in admixtures of the hadron and quark EoSs in the interpolated region. This feature causes a finite discontinuity in $c_{s}^{2}(n_{\rm B})$ at low density, which is an artifact of the scheme. Alternative forms of interpolation suggested in e.g. Refs. Fukushima and Kojo (2016); Kojo et al. (2015) do not allow for spillovers into the region of interpolation. Use of such forms, however, does not qualitatively change the outcome: while $2\,{\rm M}_{\odot}$ NS can still be produced, the constraints on $R_{1.4}$ and $\tilde{\Lambda}$ cannot be easily transformed into constraints on the parameters of interpolated EoSs. If, however, a stiff hadronic matter EoS such as MS-A in Sec. III is applied, the resulting radius and tidal deformability are apparently too large and violate the condition $\tilde{\Lambda}({\mathcal{M}}=1.186\,{\rm M}_{\odot})\leq 720$ Abbott et al. (2019).

As can be seen from Fig. 5 (d)-(f), the softer MS-B EoS is by itself compatible with the current constraint on the binary tidal deformability. Implementing the crossover region through interpolation further enhances the compatibility. Better measurements of $\tilde{\Lambda}({\mathcal{M}})$ from multiple merger detections in the future might help in limiting the relevant interpolation parameters. Recall that such “soft HM $\to$ stiff QM” combination is usually forbidden in a first-order transition, given the absence of an intersection in the $P$-$\mu$ plane between pure hadronic and pure quark phases (see summary in Table 7).

Quarkyonic matter

In this case, we present results obtained by using the hadronic EoSs MS-B/C for pure neutron matter, and two-flavor quark EoSs with and without interactions between quarks when they appear. The main reason for the rapid increase in pressure at supra-nuclear densities and the attendant behavior of $c_{s}^{2}$ vs $n_{\rm B}$ is also elucidated in more detail than was done in Ref. McLerran and Reddy (2019).

In the quarkyonic picture, both the maximum mass and typical radii are larger than those obtained by EoSs with neutrons only. In fact, some EoSs that are too soft to survive the $M_{\rm max}\geq 2\,{\rm M}_{\odot}$ constraint can be rescued by the boost in stiffness once quarkyonic matter appears; see e.g. MS-C (PNM) in Fig. 6 (a)-(c). However, for a stiff neutrons-only EoS, if a transition into quarkyonic matter takes place, compatibility with binary tidal deformability constraint from GW170817 becomes reduced, because of the tendency to also increase $R$ and therefore $\Lambda$. These increases put the model at more risk of breaking the upper limit on $\Lambda$. This is evident in Fig. 6 (d)-(f), where the MS-B (PNM) EoS is at the edge of exclusion, and with quarkyonic matter the situation is slightly worse.

An examination of the behavior of $c_{s}^{2}$ vs $n_{\rm B}$ with and without quarks offers insights into the role played the presence of the shell for $k_{Fn}>\Delta$ in the quarkyonic model. Fig. 7 shows results of $c_{s}^{2}$ for the cases in which there is no shell (i.e. PNM throughout the star), neutrons only below and above $\Delta$ (i.e. with a shell but no quarks), and with the inclusion of quarks for $k_{Fn}>\Delta$. The results in this figure correspond to the neutron matter EoS used in Ref. McLerran and Reddy (2019) and the MS-C+vNJL model of this work with $G_{v}/G_{s}=0.5$. For the former case, values of $\Lambda_{Q}=420$ MeV and $\kappa=1$ were used to calculate the shell width $\Delta$, and the onset of quarks occurs at $n_{\rm trans}=0.37~{}{\rm fm}^{-3}$. This is to be compared with $n_{\rm trans}=0.24~{}{\rm fm}^{-3}$ obtained with $\Lambda_{Q}=380$ MeV and $\kappa=0.3$ in the EoS of Ref. McLerran and Reddy (2019). For the two-flavor vNJL model used in this connection, values of the parameters used were $\Lambda=631.4$ MeV and $G_{s}\Lambda^{2}=1.835$ as in Ref. Hatsuda and Kunihiro (1994).

The main differences between the models in Ref. McLerran and Reddy (2019) and this work are:

(i) For pure neutron matter (no quarks), the EoS of Ref. McLerran and Reddy (2019) becomes acausal for $n_{\rm B}/n_{0}\simeq 6$ owing to the term proportional to $n_{n}^{3}$ in its interacting part. As the central density of the star is $\simeq 6.74\,n_{0}$, this feature may be of some concern. However, the MS-B/C+vNJL models - being relativistically covariant - remain causal for all densities, and

(ii) Interactions between quarks are not included in the EoS of Ref. McLerran and Reddy (2019) except in the kinetic energy term with the use of $M_{q}=M_{n}/3$, whereas the MS-B/C+vNJL model uses density-dependent dynamically generated $u,d$ quark masses that steadily decrease with increasing density from their vacuum values of $\simeq M_{n}/3$. In addition, repulsive vector interactions between quarks were used in the vNJL models.

The above differences notwithstanding, the inner workings of the quarkyonic model - particularly, the influence of quarks - are apparent from Fig. 7 (a) and (c). Without the presence of quarks in the shell, the EoSs in both models are very stiff even to the point of being substantially acausal. The presence of quarks in the shell abates this undesirable behavior by softening the overall EoS (dash-dotted blue curves) relative to the case when only nucleons are present (dotted gray curves). With progressively increasing density, the density of nucleons is depleted within the shell whereas that of the quarks becomes predominant. As $c_{s}^{2}\rightarrow 1/3$ for quarks at asymptotically high densities, it exhibits a maximum (as well as a minimum) at some intermediate density. Note however, that compared to the case of no shell, pure neutron matter everywhere (black solid curves), the overall EoS of the quarkyonic matter is still stiffer within the central densities of the corresponding stars.

Insofar as $c_{s}^{2}$ is a measure of the stiffness of the EoS, the $M$-$R$ curves shown in Fig. 7 (b) and (d) reflect the corresponding $c_{s}^{2}$ vs $n_{\rm B}$ behavior. The presence of quarkyonic matter (dash-dotted blue curve) causes an increase in the $M_{\rm max}$ for both models. If only the neutron content of quarkyonic matter is considered (dotted gray curve), then the increase in $M_{\rm max}$ is more substantial. Similarly, the radii of both the maximum mass and $1.4\,{\rm M}_{\odot}$ stars are significantly larger in quarkyonic matter. Quantitative differences between the two cases can be attributed to the presence of interactions between quarks in the MS-B+vNJL model.

The hadron-to-quark transition density $n_{\rm trans}$, the peak value of the squared speed of sound $c_{s,\,{\rm max}}^{2}$, and the maximum mass $M_{\rm max}$ all depend on the choice of $\Lambda_{Q}$ and $\kappa$ used to calculate the shell width $\Delta$. Fig. 8 shows the variation of these quantities as a function of $\Lambda_{Q}$ with $\kappa=0.1$, $0.6$ and 1 for the MS-C+vNJL model chosen here. Intermediate values of $\kappa$ lead to results that lie within the boundaries shown in this figure. Note that high values of both $\Lambda_{Q}$ and $\kappa$ are required to ensure that $n_{\rm trans}\gtrsim 1.5\,n_{0}$ and $c_{s}^{2}<1$. This requirement decreases $M_{\rm max}$, but masses above the current constraint of $\gtrsim 2\,{\rm M}_{\odot}$ can be still obtained. In the absence of interactions between quarks (as in Ref. McLerran and Reddy (2019)), the window of $\Lambda_{Q}$ and $\kappa$ values that are usable is very small. We stress however, that the optimum choice of these parameters is model dependent in that if a different hadronic or quark EoS is used, the values of $\Lambda_{Q}$ and $\kappa$ can change.

On a physical level, low values of $\Lambda_{Q}$ and $\kappa$ lead to a substantial quark content in the star, but at the expense of $n_{\rm trans}\rightarrow n_{0}$ – a disturbing trend. Although quarks soften the overall EoS, the presence of the shell and the redistribution of baryon number between nucleons and quarks causes a substantial stiffening of the overall EoS, which in turn leads to very high values of $M_{\rm max}$. Conversely, very high values of $\Lambda_{Q}$ and $\kappa$ decrease the quark content which makes the overall EoS to be nearly that without quarks. This feature is generic to the quarkyonic model, which enables it to achieve maximum values consistent with the observational mass limit even when the EoS with hadrons only fails to meet this constraint.

The low transition densities and the extreme stiffening of the EoS caused by the shell in quarkyonic matter bear further investigation. Although inspired by QCD and large $N_{c}$ physics, the width of the shell is independent of the EoSs in both hadronic and quark sectors, at least in the initial stage of the development of the model. The energy cost in creating such a shell in dense matter is another issue that warrants scrutiny.

V Conclusion and Outlook

In this work, we have performed a detailed comparison of first-order phase transition and crossover treatments of the hadron-to-quark transition in neutron stars. For first-order transitions, results of both Maxwell and Gibbs constructions were examined. Also studied were interpolatory schemes and the second-order phase transition in quarkyonic matter, which fall in the class of crossover transitions. In both cases, sensitivity of the structural properties of neutron stars to variations in the EoSs in the hadronic as well as in the quark sectors were explored. The ensuing results were then tested for compatibility with the strict constraints imposed by the precise mass measurements of $2\,{\rm M}_{\odot}$ neutron stars, the available limits on the tidal deformations of neutron stars in the binary merger GW170817, and the radius estimates of $1.4\,{\rm M}_{\odot}$ stars inferred from x-ray observations. These independent constraints from observations are significant in that the lower limit on the maximum mass reflects the behavior of the dense matter EoS for densities $\gtrsim 4-6\,n_{0}$, whereas bounds on binary tidal deformability $\tilde{\Lambda}$ and estimates of $R_{1.4}$ depend on the EoS for densities $\gtrsim 2-3\,n_{0}$, respectively.

Table 6 provides a summary of the transition density $n_{\rm trans}/n_{0}$ for the appearance of quarks and the associated neutron star mass, $M_{\rm trans}$, for the EoSs and different treatments of the phase transition considered in this work. The entries in this table allow us to answer the first two of the three questions posed in the introduction:

(a) What is the minimum neutron star (NS) mass consistent with the observational lower limit of the maximum mass ($M_{\rm max}$) that is likely to contain quarks?

The answer to this question depends on both the low-density hadronic and high-density quark EoSs, as well as the order and the method of implementing the phase transition. Barring rare cases, such as in the Gibbs construction and interpolation (see why below), the minimum mass for phase transition is $M_{\rm trans}\gtrsim 1\,{\rm M}_{\odot}$.

(b) What is the minimum physically reasonable density at which a hadron-quark transition of any sort can occur?

Our results indicate the minimum density $n_{\rm trans}$ to be $\gtrsim 2\,n_{0}$, again excluding the rare cases. The reasons for the exclusions are as follows. In the Gibbs construction, valid in the extreme case of the interface tension between the phases being zero, the onset density is generally lower than that of the corresponding Maxwell construction. Depending on the softness or stiffness of the hadronic and quark EoSs, $n_{\rm trans}$ can approach near-nuclear densities. Such cases should be discarded as being in conflict with nuclear data near saturation. In the case of interpolation, the onset density is chosen a priori. In this approach, input values of $n_{\rm trans}\gtrsim 2\,n_{0}$ yield the minimum mass for phase transition $M_{\rm trans}\gtrsim 1\,{\rm M}_{\odot}$.

The values of $n_{\rm trans}$ and $M_{\rm trans}$ quoted in Table 6 do not conflict with experimental data on nuclei that probe densities near and below $n_{0}$. A mild tension however, exists with theoretical interpretations of low-to-intermediate energy heavy-ion data Danielewicz et al. (2002) which probe densities up to $3-4\,n_{0}$. We wish to note, however, that analysis of such data using Boltzmann-type kinetic equations has not yet been performed with quark degrees of freedom and their subsequent hadronization as in RHIC and CERN experiments.

Table 7 provides a summary of the generic outcomes of our study. If the hadron-to-quark transition is strongly first-order, as is the case for standard quark models such as vMIT and vNJL that we used, then the hadronic part needs to be relatively stiff to guarantee a proper intersection in the $P$-$\mu$ plane. For a hadronic EoS as stiff as MS-A, this combination brings tension with $\Lambda_{1.4}$ or $R_{1.4}$ estimates. Concomitantly, a too-high transition density that yields $M_{\rm trans}\gtrsim 1.7\,{\rm M}_{\odot}$ results in either very small quark cores or completely unstable stars that are indistinguishable from those resulting from the (stiff) purely-hadronic EoS. Thus, such hybrid EoSs can easily be ruled out. This is typical for NJL-type models; see e.g. Fig. 4. In contrast, lower transition densities (that yield $M_{\rm trans}\lesssim 1.0-1.6\,{\rm M}_{\odot}$) are capable of decreasing radii, and if accompanied by a stiff quark matter EoS with sizable repulsive interactions, these hybrid EoSs produce $M_{\rm max}\geq 2\,{\rm M}_{\odot}$. Figs. 2 and 3 show examples in which the vMIT model was used with Maxwell/Gibbs constructions.

Our analysis indicates that use of the Gibbs construction is beneficial in satisfying the current constraints from observation for many stiff hadronic EoSs, as it enlarges the parameter space of quark models. As similar $M(R)$ and $\Lambda(M)$ relations for hybrid and purely-hadronic stars can be obtained, the distinction between the two is, however, lost. This feature underscores the significance of dynamical properties such as neutron star cooling and spin-down, and the evolution of merger products.

To sum up the part about first-order phase transitions, current observational constraints disfavor weakly-interacting quarks at the densities reached in neutron star cores. Should a first-order transition into strongly-interacting quark matter (as described by the vMIT bag model or vNJL-type models) take place, the onset density is likely of relevance also to canonical neutron star masses in the range $1.0-1.6\,{\rm M}_{\odot}$.

One should keep in mind, however, that perturbative approaches to the quark matter EoS are not expected to hold in the density range $\approx 2-4\,n_{0}$. This limitation brings the validity of first-order phase transitions caused by such EoSs into question. In this regard, model-independent parameterizations circumvent the issue and have the advantage of translating observational constraints more generically. For instance, specific QM models prohibit the transition into soft hadronic matter, but in the CSS parameterization this restriction disappears and a much larger parameter space can be explored including soft HM $\to$ stiff QM Han and Steiner (2019). However, such parameterizations lack a physical basis and beg for the invention of a non-perturbative approach.

If the hadron-to-quark transition is a smooth crossover, as in the case of interpolatory schemes and in quarkyonic matter, the pressure in the transition region is stiffened unlike the sudden softening of pressure caused by a first-order transition. This stiffening is also reflected in a local peak in the sound velocity before the pure quark phase is entered. This stiffening is responsible for supporting massive stars that are compatible with the current lower limit of 2 ${\rm M}_{\odot}$.

It is also common that the onset density for quarks is somewhat low ($M_{\rm trans}\lesssim 1.0-1.6\,{\rm M}_{\odot}$) in these crossover approaches. This feature implies that all neutron stars we observe should contain some quarks admixed with hadrons. We find that at low densities soft hadronic EoSs are necessary, but above the transition changes in radii rely heavily on the methods of implementing the crossover in both the interpolation approach and in quarkyonic matter. Consequently, it is difficult to obtain physical constraints on the crossover EoSs from a better determination of the radius, e.g. $R_{1.4}$, or improved tidal deformability measurements. It is promising, however, in limiting parameters, e.g. the vector coupling strength $G_{v}/G_{s}$ in vNJL or $\kappa$ and $\Lambda$ in the quarkyonic model, pertinent to the required stiffening to satisfy the limits imposed by mass measurements of heavy neutron stars.

Regardless of the phase transition being first-order or crossover, our results suggest that the presence of quarks in the pre-merger component neutron stars of GW170817 is a viable possibility. If quarks only appear after the merger (before the remnant collapsing into a black hole), there is a valid soft HM $\to$ stiff QM first-order transition that cannot be captured by the vMIT bag or vNJL models. There are exceptions when the onset occurs close to the $2\,{\rm M}_{\odot}$ limit, so that quarks are precluded in cold beta-equilibrated NSs due to immediate collapse. While we rejected these solutions by default, these cases can however be relevant for the dynamic products of mergers where quarks may emerge temporarily Most et al. (2019); Chesler et al. (2019). Numerical simulations that involve quarks Most et al. (2019); Bauswein et al. (2019) will assist in identifying such cases during post-merger gravitational-wave evolution. Better understanding and progress in theory, experiments, and observations are required to clarify the situation.

Although the presence of quarks in neutron stars is not ruled out by currently available constraints, it is nearly impossible to confirm it even with improved determinations of radii from x-ray observations and tidal deformabilities from gravitational wave detections. This conundrum arises because purely hadronic EoSs can also satisfy the current constraints of $M_{\rm max},~{}R_{1.4}$ and $\tilde{\Lambda}$; i.e., the “masquerade problem” Alford et al. (2005) persists. Similarly, it will be difficult to identify the nature of the phase transition on the basis of $M$ and $R$ observations only, unless there is a sufficiently strong first-order transition that gives rise to separate branches of twin stars with discontinuous $M$-$R$ and/or $\tilde{\Lambda}$-${\mathcal{M}}$ relations.

We now turn to the third question in the introduction: (c) Which astronomical observations have the best potential to attest to the presence of quarks? Dynamical observables such as supernovae neutrino emission, thermal/spin evolution, global oscillation modes, continuous gravitational waves, dynamic collapse etc. that are sensitive to transport properties would potentially provide more distinct signatures of exotic matter in neutron stars Fischer et al. (2018); Aloy et al. (2019); Alford et al. (2019). In future work, it is worthwhile to achieve consistency with dynamical observables, particularly for the crossover scenarios of the transition to quark matter.

Acknowledgements

We thank Larry McLerran and Sanjay Reddy for helpful discussions relating to quarkyonic matter. This work was performed with research support from the National Science Foundation Grant PHY-1630782 and the Heising-Simons Foundation Grant 2017-228 for S.H., the U.S. Department of Energy Grant No. DE-FG02-93ER-40756 for M.A.A.M., S.L. and M.P. The research of C.C. is supported by the National Science Foundation Grant PHY-1748621.

Appendix A Thermodynamics of nucleons in the shell of quarkyonic matter

In this Appendix, we provide some details of the evaluation of the kinetic parts of the energy density, chemical potential and energy density and pressure for nucleons in the shell. The expressions we obtain will then be used to establishing the thermodynamic identity (TI) in the presence of a shell. For the evaluation of these quantities the relation

[TABLE]

where $\alpha$ is a parameter in the functions $\phi_{1},\phi_{2}$ and $F$ will be useful.

Energy density

The kinetic energy density of nucleons, neutrons to be specific, in the shell is

[TABLE]

where $\Delta=\Lambda_{Q}^{3}/k_{Fn}^{2}+\kappa\Lambda_{Q}/N_{c}^{2}$. In analytical form,

[TABLE]

For $k_{Fn}<\Delta$, the upper limit $U=k_{Fn}$, the lower limit $L=0$, and $e_{k}={\sqrt{k^{2}+M_{n}^{2}}}$, which leads to the familiar expression for spin-$\frac{1}{2}$ relativistic particles of mass $M_{n}$. For neutrons in the shell with $k_{Fn}>\Delta$, however, $U=k_{Fn}$ and $L=(k_{Fn}-\Delta)$ with $e_{k}={\sqrt{(k-\Delta)^{2}+M_{n}^{2}}}$.

Chemical potential

The associated chemical potential ensues from

[TABLE]

For the neutrons in the shell,

[TABLE]

For evaluating $dF_{1}/dk_{Fn}$, use of the relations

[TABLE]

in Eq. (LABEL:eqn:Leibnitz) yields

[TABLE]

The evaluation of $dF_{2}/dk_{Fn}$ proceeds along similar lines, but with

[TABLE]

with the result

[TABLE]

Putting these results together, we obtain after some simplification

[TABLE]

Pressure

The kinetic theory expression for a single species of spin-$\frac{1}{2}$ fermions is

[TABLE]

where $\beta=1/T$, $e$ is the single particle spectrum and $\mu$ the chemical potential. A partial integration on the right-hand side yields

[TABLE]

At finite $T$, the first term vanishes when $U=\infty$ and $L=0$, leaving the second term as the kinetic pressure. For $T\rightarrow 0$, and finite $L$ and $U$, however, we have

[TABLE]

The expressions for $P_{n}^{{\rm(kin)}}$, thus take different forms in the regions $k_{Fn}<\Delta$ and $k_{Fn}>\Delta$. For $k_{Fn}<\Delta$, the limits $U=k_{Fn}$ and $L=0$ yield the familiar kinetic theory expression

[TABLE]

where $e_{k}=\sqrt{k^{2}+M_{n}^{2}}$. The last two terms in Eqs. (32) and (LABEL:Pnkin) cancel, and thus in this region, $\varepsilon_{n}^{\rm(kin)}+P_{n}^{{\rm(kin)}}=n_{n}\mu_{n}^{{\rm kin}}$ (the TI) with $n_{n}=k_{Fn}^{3}/(3\pi^{2})$ and $\mu_{n}^{\rm(kin)}=\sqrt{k_{Fn}^{2}+M_{n}^{2}}$.

In the region $k_{Fn}>\Delta$ with $U=k_{Fn}$ and $L=(k_{Fn}-\Delta)$, the kinetic theory pressure becomes

[TABLE]

The first term above gives the contribution from the boundaries of the shell. Inserting the appropriate limits for the shell, this term reads as

[TABLE]

Thermodynamic identity

Collecting the results in Eqs. (LABEL:PnTI1) and (A) and evaluating $\varepsilon_{n}^{\rm(kin)}$ in Eq. (32) for the shell, we obtain the TI (after many cancellations)

[TABLE]

the last relation giving the number density of neutrons in the shell. The neutron chemical $\mu_{n}^{\rm kin}$ here is independently calculated from Eq. (39). We have verified that, numerical calculations of $\mu_{n}^{\rm kin}=d\varepsilon_{n}^{\rm(kin)}/dn_{n}$ and $P_{n}^{\rm(kin)}=n_{n}^{2}d(\varepsilon_{n}^{\rm(kin)}/n_{n})/dn_{n}$ in the shell yield the same results as evaluations from the analytical formulas above, which provide additional physical insight concerning the role played by the shell.

All expressions in this Appendix apply also to protons so that a generalization to a 2-component system consisting of neutrons and protons is straightforward.