New Neutron Star Equation of State with Quark–Hadron Crossover

It is expected that the cold dense matter in neutron stars is nucleonic (hadronic) at lower densities and becomes a strongly interacting quark phase at high densities. The equation of state of the matter governs the structure of neutron stars. In this paper we construct a much improved version, called QHC19, of the equation of state of cold dense matter based on a smooth crossover from the hadronic to quark phase (Masuda et al. 2013; Kojo et al. 2015; Baym et al. 2018). Such crossover equations of state, labeled QHC (for Quark–Hadron Crossover), are based on a realistic nucleonic matter equation of state at low to medium densities, a strongly interacting quark matter equation of state at high density, and a thermodynamically consistent (and highly constrained) interpolation between the two regimes. See also Annala et al. (2019) for a closely related approach to the equation of state.

The new features we take into account are: first, a consistent equation of state for nucleonic matter (Togashi & Takano 2013; Togashi et al. 2017) from neutron drip in the crust to about 2n₀ in the hadronic interior, where n₀ is the nuclear saturation density, ∼0.16 fm⁻³, corresponding to 2.7 × 10¹⁴ g cm⁻³; and second we identify and implement strong constraints on the two phenomenological parameters used in the high density quark phase, the universal quark repulsion of coupling strength g_V, and the Bardeen–Cooper–Schrieffer (BCS), or diquark pairing strength, H. The resulting QHC19 equation of state agrees, as we show below, strikingly well with the inference of the equation of state by the LIGO/Virgo collaboration from the binary neutron star merger event GW170817 (see Figure 8 below), and is fully consistent with their radius and tidal deformability determinations (Abbott et al. 2018, 2019). It also allows for a maximum neutron star mass ≃2.35 M_⊙, consistent with the recent observation of a 2.14 ± 0.1 M_⊙ neutron star in the pulsar PSR J0740+6620 (Cromartie et al. 2019). In addition, the QHC19 equation of state can be generalized to finite temperature for use in dynamical simulations of neutron star and neutron star–black hole mergers, as well as in other high-energy astrophysical phenomena, including core-collapse supernovae and cooling of proto-neutron stars (Furusawa et al. 2017). The full details of QHC are posted on the CompOSE archive at https://compose.obspm.fr/eos/140/.

This paper is organized as follows. In the following section, we first give a brief overview of QHC equations of state, and then present the essence of the Togashi equation of state for hadronic matter in Section 2.1, the Nambu–Jona–Lasinio (NJL) model for quark matter at higher density in Section 2.2, and the interpolation region in Section 2.3. In Section 3, we discuss the physical ranges of the phenomenological parameters in the NJL model based on the constraints of thermodynamic stability and the sound velocity not exceeding the speed of light. In Section 4, we discuss the structure of neutron stars described by the QHC19 equation of state, and conclude in Section 5.

At very low baryon density the system is characterized by nonrelativistic nucleons interacting via nuclear forces, while at high densities one must take quarks seriously as the physical degrees of freedom. The equation of state in the low density regime can be found by the techniques of nuclear matter theory, while at high densities, owing to the inability to calculate via lattice quantum chromodynamics (QCD) simulations the equation of state at finite baryon density, one must adopt a phenomenological model of high density quark model; we use the NJL model in particular. Furthermore the strongly interacting intermediate regime, across which the system transitions from nucleonic to quark degrees of freedom, is still not well understood. While we interpolate between the two regimes, the interpolation turns out to be highly constrained in order that the matter be thermodynamically stable.

2.1. The Hadronic Regime

In the previous version of our equation of state, QHC18 (https://compose.obspm.fr/eos/139/), we adopted the Akmal, Pandharipande, and Ravenhall (APR) equation of state (Akmal et al. 1998) to describe the liquid nuclear matter interior of neutron stars up to a density 2n₀. APR is the result of a thorough variational calculation, including realistic two- and three-body nuclear forces, of the ground state energy of pure neutron matter (PNM) and symmetric nuclear matter (SNM) containing equal densities of neutrons and protons. However, APR has three drawbacks. First, to describe matter in beta equilibrium containing a finite proton fraction, one generally makes a quadratic interpolate between the PNM and SNM limits in terms of α = 1 − 2Y_p; here Y_p = n_p/n_B, with n_B = n_n + n_p the total number density of baryons. The second is that descriptions of nonuniform nuclear matter in the crust do not use the same physics as APR uses in the liquid nuclear matter interior, and thus joining different equations of state in the crust to the liquid interior is inherently inconsistent. Finally APR predicts that the sound velocity exceeds the speed of light at n_B ≳ 5.5n₀, densities accessible in the core of high mass neutron stars (this drawback also shared with the Togashi equation of state).

Togashi et al. (Togashi & Takano 2013; Togashi et al. 2017) have recently extended the APR approach, to construct a thermodynamically consistent equation of state that overcomes these drawbacks of APR. The Togashi equation of state, which we now employ in QHC19, treats the nonuniform matter in the crust and the uniform matter at medium density within the same framework, and at the same time describes all values of the proton fraction Y_p without an interpolation between PNM and SNM. Specifically we use the Togashi equation of state from density⁸ 4.74 × 10⁻¹⁰n₀ to 2n₀.

We briefly summarize the key features of the Togashi equation of state in comparison to APR. Both are based on a Hamiltonian including the two-body Argonne V18 potential extracted by fitting two-nucleon experimental scattering data, as well as the more empirical three-body Urbana IX potential fit in part to light nuclei. Both APR and Togashi et al. solve the many-body problem by choosing a variational wave function with parameters determined by minimizing the total energy. The general form of the wave function is that of a free Fermi gas multiplied by factors that build in two-particle correlations dependent on the relative spin, isospin, and orbital angular momentum of the two particles. Despite differences in detail both the Togashi and APR calculations of the energy of PNM and of SNM are in good agreement (Kanzawa et al. 2007; Togashi & Takano 2013).

An illustrative comparison of the Togashi and APR equations of state near the saturation point in SNM can be made by expanding the energy per baryon in powers of x = (n − n₀)/3n₀ and α = 1 − 2Y_p, following the convention in the CompOSE repository:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{E}{A} & = & m-{E}_{0}+\displaystyle \frac{{x}^{2}}{2}{K}_{0}+{ \mathcal O }({x}^{3})\\ & & +\,{\alpha }^{2}\left(J+{xL}+{ \mathcal O }({x}^{2})\right)+{ \mathcal O }({\alpha }^{4}).\end{array}\end{eqnarray} \tag{ 1 }$$

The saturation energy E₀ of SNM is 16.1 MeV in Togashi versus 16.0 MeV in APR, and the incompressibility K₀ = 245 versus 267 MeV in APR. Near n₀, SNM in the Togashi equation of state is slightly softer than in APR, as indicated by the smaller values of the symmetry energy J = 29.1 versus 34.0 MeV in APR, and the slope parameter, L = 38.7 versus 63.2 MeV in APR.⁹ In neutron stars, the Togashi equation of state is slightly stiffer at higher densities than APR, as can be seen in the calculation of a 1.4 M_⊙ neutron star containing nuclear matter only (where in the crust in APR we use the Togashi equation of state for n_B ≤ 0.26n₀). The radii are R = 11.6 km for Togashi versus 11.5 km for APR, and the tidal deformability, Λ = 360 versus 268. In contrast, QHC19 (set D, Figure 1) gives R = 11.6 km and Λ = 350 for a 1.4 M_⊙ neutron star.

Zoom In Zoom Out Reset image size

In the neutron star crust, Togashi et al. assume a single species of heavy spherical nuclei forming a body-centered cubic (BCC) lattice surrounded by a gas of nucleons (Oyamatsu 1993; Shen et al. 2011), and construct the equation of state by the Thomas–Fermi method. The crust extends to density ≃0.625n₀, found by comparing the energy per nucleon, for given average n_B and Y_p, for inhomogeneous matter with that for homogeneous matter. As reported in Kanzawa et al. (2009), the properties of neutron star crusts calculated with the Togashi equation of state are consistent with those in previous studies (Baym et al. 1971; Negele & Vautherin 1973; Douchin & Haensel 2001).

The Togashi equation of state outlined above can be generalized to finite T in terms of the Helmholtz free energy, F(n_B, Y_p, T), (Togashi et al. 2017) thus providing a unified framework to treat both homogenous and inhomogeneous matter with thermodynamic consistency. The Togashi equation of state is available at http://www.np.phys.waseda.ac.jp/EOS/ and in the CompOSE archive at https://compose.obspm.fr/eos/105/ as a table of thermodynamic quantities over a wide range of baryon density, n_B/n₀, from 4.7 × 10⁻¹⁰ to 38n₀, Y_p from 0 to 0.65; and temperature T from 0 to ∼400 MeV.

2.2. Quark Degrees of Freedom

2.3. Quark–Hadron Crossover

The quark model parameters that have a particularly important impact on the pressure and hence neutron star structure are g_V and H. We first determine the ranges of these parameters that are consistent with the constraint that the sound velocity, c_s, never exceeds the speed of light. Figure 1 shows ${({c}_{s}^{\max }/c)}^{2}$ as a function of g_V and H; the resolution in the figure is Δg_V = ΔH = 0.01G. The uncolored area violates either causality or thermodynamic stability (i.e., leads to an inflection point in P(μ_B) and is unphysical). As the figure indicates g_V and H are strongly correlated. The causality constraint bounds g_V from above by ≃1.3G, implying that within this NJL description quark matter equations of state for n_B ≥ 5n₀ cannot be too stiff.

As we see in Figure 1, the maximum sound velocity for given g_V is least in the middle of the allowed region of H. We focus here on four parameter sets,¹¹ with H chosen to give the least ${c}_{s}^{\max }$ for given g_V: (g_V, H)/G = (0.6, 1.43) (set A), (0.8, 1.49) (set B), (1.0, 1.55) (set C) and (1.2, 1.61) (set D), indicated in Figure 1.

In Figure 2 we show the density dependence of c_s² for the sets A–D together with that in APR and Togashi. Replacing a hadronic with a quark equation of state at high density typically reduces the speed of sound. For the case of APR and Togashi, the violation of causality (c_s > c) for n_B ≳ 5.5n₀ is cured by this replacement. The bump in the QHC19 curves for n_B = 3–5n₀ is a consequence of connecting a soft low density equation of state with a stiffer one at high density, without violating causality. (By contrast, the sound velocity at finite temperature and zero baryon density dips between the hadronic and quark regimes; Asakawa & Hatsuda 1997.)

Zoom In Zoom Out Reset image size

We now indicate the astrophysical consequences of the QHC19 equation of state for the physical range of g_V and H. Figure 3 shows the maximum neutron star mass, M_max, for given g_V and H. The absolute maximum allowed mass is ≃2.35 M_⊙ at (g_V, H)/G ≃ (1.30, 1.65), which is at the causal boundary. The constraint that the equation of state must be able to support stars of at least two solar masses sets the lower bound, g_V/G ≳ 0.60–0.75, where the detailed value depends on H. The recently determined neutron star of mass 2.17 ± 0.1 M_⊙ (Cromartie et al. 2019) tightens the constraint to g_V/G ≳ 0.9. Beyond the minimal value of H, the maximum mass is reduced by increasing H until the equations of state become acausal.

Zoom In Zoom Out Reset image size

Figure 4 shows the central baryon density ${n}_{{\rm{B}}}^{{\rm{c}}}$ in units of n₀ for the maximum mass neutron star for each parameter set. The central density at M_max always exceeds 5.5n₀. The central density falls with g_V because the repulsion disfavors large density, while increasing H tends to increase the central density; although larger H stiffens the equation of state of quark matter, it softens it in the interpolated region, with the net effect that the central density is larger for larger H.

Zoom In Zoom Out Reset image size

The baryon density distribution of neutron stars for parameter set D is shown in Figure 5. For M > 1.2M_⊙ the core is beyond the pure nuclear regime. As M increases to ∼2.0M_⊙, the core density increases slowly to ∼4n₀. Beyond this region the core density grows more rapidly, reflecting the peak structure of the sound velocity around n_B = 3–5n₀ (Figure 2); after passing the peak the equation of state softens.

Zoom In Zoom Out Reset image size

The neutron star mass–radius relation is shown in Figure 6 for the parameter sets A–D, together with the masses of PSR J1614–2230 (M/M_⊙ = 1.928 ± 0.017), PSR J0348–0432 (M/M_⊙ = 2.01 ± 0.04), and PSR J0740+6620 (M/M_⊙ = 2.17 ± 0.10). The curves for A–D largely overlap for M ≲ 1.4M_⊙ where the corresponding density is n_B ≲ 2.5–3n₀. The radii at M = 1.4M_⊙ are ≃11.6 km, largely due to the relative softness of the Togashi equation of state. The radii are consistent with the LIGO/Virgo inference 11.9 ± 1.4 km (Abbott et al. 2018).

Zoom In Zoom Out Reset image size

We look now at the dimensionless tidal deformability Λ of neutron stars (Hinderer 2008) calculated with the QHC19 equation of state. The deformability has a strong impact on the phase evolution of gravitational waves from neutron star mergers, especially at the stage where two neutron stars are not touching but are close enough for the tidal fields to deform each other. The deformability is strongly correlated with the neutron star compactness, M/R; less compact stars deform more easily, and thus have larger Λ. The gravitational waveforms in a merger of two stars of masses M₁ and M₂ are sensitive to the combination,

$$\begin{eqnarray}&&\tilde{{\rm{\Lambda }}}=\displaystyle \frac{16}{13}\displaystyle \frac{({M}_{1}+12{M}_{2}){M}_{1}^{4}{{\rm{\Lambda }}}_{1}+(12{M}_{1}+{M}_{2}){M}_{2}^{4}{{\rm{\Lambda }}}_{2}}{{\left({M}_{1}+{M}_{2}\right)}^{5}},\end{eqnarray} \tag{ 2 }$$

of the masses and Λ_i of the individual neutron stars. The individual Λ_i can be regarded as functions of the stellar M_i, and thus the independent variables in $\tilde{{\rm{\Lambda }}}$ are M₁ and M₂.

In Figure 7 we show $\tilde{{\rm{\Lambda }}}$ calculated with QHC19 for three representative chirp masses, M_chirp = (M₁M₂)^3/5 (M₁ + M₂)^−1/5, equal to 1.188 M_⊙ (corresponding to GW170817), 1.045M_⊙, and 1.306M_⊙, as a function of the mass ratio q = M₁/M₂. For equal masses (q = 1.0) these chirp masses correspond to M₁ = M₂ =1.37M_⊙, 1.2M_⊙, and 1.5M_⊙. For the mass ratio q = 0.7–1.0 the parameter sets A–D give very similar curves for $\tilde{{\rm{\Lambda }}};$ we display only that for set D. The calculated $\tilde{{\rm{\Lambda }}}$ for ${{ \mathcal M }}_{\mathrm{chirp}}=1.188{M}_{\odot }$ is consistent with the LIGO/Virgo constraints: for low spin $\tilde{{\rm{\Lambda }}}={300}_{-190}^{+500}$ (symmetric interval) and $\tilde{{\rm{\Lambda }}}={300}_{-230}^{+420}$ (high posterior density interval), and for high spin, $\tilde{{\rm{\Lambda }}}=0\mbox{--}630$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, in Figure 8 we compare the QHC19 equation of state, P versus m_Nn_B, with that recently inferred by the LIGO/Virgo collaboration (Abbott et al. 2018). The latter takes into account the two solar mass constraint (which should now be larger), causality, the inferred tidal deformability, and assumes the nuclear equation of state at low density to be SLy4. The QHC19 equation of state is remarkably consistent, especially for parameter sets C and D, with the LIGO/Virgo inference, except at low density due to the difference between the LIGO/Virgo-assumed SLy4 and the Togashi nuclear equation of state.

The new QHC19 equation of state presented here, which includes the physics of strongly interacting quarks at high densities and the Togashi equation of state in the hadronic regime below 2n₀, well describes neutron stars with masses up to 2.35 M_⊙, and is in good agreement with inferences by LIGO/Virgo, from the GW170817 neutron star merger, of the equation of state as well as radii and tidal deformabilities. The requirements of thermodynamic stability and causal sound propagation tightly constrain the strength of the universal repulsive interaction and the pairing interaction between quarks. We await further data on binary neutron star mergers from the present LIGO/Virgo O3 run, as well as from the NICER X-ray timing observatory on the the masses and radii of neutron stars, to further constrain the parameters of the QHC19 equation of state.

There is still room for fundamental improvement, however, of the equation of state, from the crust through to dense quark matter. The present description of the nuclei in the crust neglects quantum effects such as pairing, as well as the shell structure of the neutron rich nuclei present. How the usual shell closures at neutron number 50 and 82 are modified in very neutron rich nuclei, as a consequence of the nuclear tensor force, remains an open question. In addition one must take into account nuclear "pasta" phases in the innermost crust (Ravenhall et al. 1983; Hashimoto et al. 1984; Watanabe et al. 2000; Caplan et al. 2018).

Hyperonic degrees of freedom plays a role at higher hadronic densities. Togashi et al. (2016) extends the Togashi equation of state to hyperonic nuclear matter containing Λ and Σ⁻ hyperons at zero temperature. This study finds that the Λ hyperon appears at 0.42 fm⁻³ for relatively attractive hyperon interactions. This density is well into the interpolation region in this work, where confidence in the validity of hadronic descriptions begins to break down. Improving our understanding of hadronic interactions, which should be based on quark dynamics at short distance, would allow us to narrow the interpolation domain and to directly discuss the impact of strangeness. In this respect the large uncertainty in hyperon interactions makes it desirable to investigate further effects of hyperons in dense matter with state-of-the-art hyperon–nucleon and hyperon–hyperon interactions from laboratory experiments and from lattice QCD simulations (Hatsuda 2018).

A major problem in QCD is the lack of a theory of the crossover from hadronic to quark matter, containing both hadronic and quark degrees of freedom. Efforts in this direction include, e.g., McLerran & Reddy (2019) and Song et al. (2019).

In the quark regime, one needs to go beyond the phenomenological NJL model with a momentum cutoff, which does not contain explicit gluons. One needs to understand microscopically how the parameters of QHC19—the vector coupling g_V and the pairing strength H—arise from QCD in terms of the strong interaction running coupling α_s(μ_B), and how they depend on baryon density. A step in this direction has been to relate g_V to the QCD quark exchange energy at moderate densities, including effects of nonperturbative chiral and diquark condensates as well as a gluon effective mass (Song et al. 2019).

Finally, an equation of state at non-zero temperature is needed for simulations of neutron star–neutron star and neutron star–black hole mergers. The primary effects of temperature are at lower densities where the Togashi equation of state is already available up to T ∼ 400 MeV. On the other hand, thermal effects in the quark regime are highly suppressed by the pairing gaps up to T ∼ 100 MeV. A finite temperature extension of the QHC19 equation of state in terms of the Helmholtz free energy, F(n_B, Y_p, T), is under investigation and will be reported elsewhere (Baym et al. 2019).

We are grateful to Vicky Kalogera and collaborators for very useful discussions at the Aspen Center for Physics on equations of state for neutron stars and gravitational wave observations. The research of author G.B. was supported in part by National Science Foundation grant No. PHY1714042; T.H. was partially supported by the RIKEN iTHEMS program and JSPS Grant-in-Aid for Scientific Research (S), No. 18H05236; T.K. was supported by NSFC grant 11650110435 and 11875144; S.F. was supported by JSPS KAKENHI (grant Nos. JP17H06365 and 19K14723); and H.T. by JSPS KAKENHI (Nos. 18K13551 and 18H04598). This work was carried at the Aspen Center for Physics, which is supported by National Science Foundation grant No. PHY-1607611.

Our interpolation method assumes fixed locations of the boundaries for hadronic and quark matter. As the onset density of quark matter is the major uncertainty, we study upper boundaries n_QM = 5, 6, and 7n₀ while fixing the boundary of the Togashi equation of state at n_NM = 2n₀. The results are shown in Figure 9. Overall, the extension of the interpolation range relaxes the constraint on the model parameters g_V and H extending their acceptable ranges by ∼10%. Meanwhile wider interpolation range also excludes part of the allowed domain for a narrower interpolation range (for example (g_V, H)/G = (1.3, 1.65) is allowed for n_QM = 5n₀ but not for 7n₀). This is simply an artifact due to the insufficiency of our interpolation function with sixth-order polynomials; in order to cover a wider range we need to include higher order polynomials in the interpolation. We have checked this statement by increasing the number of polynomials whose coefficients, not fixed by the boundary conditions, are varied for some range. These analyses show that the results presented in the main context are good representatives; the change of the quark matter boundary from 5n₀ to 7n₀ does not significantly affect our conclusions in the main text. On the other hand, we find that reducing n_QM from 5n₀ to 3–4n₀, allowing quark matter at more dilute densities, significantly reduces the acceptable domain for (g_V, H) that is compatible with the constraints from high mass neutron stars. The variations considered here give a first impression of the effect of changing the interpolation boundaries. A more complete study is desirable but beyond the scope of this paper.

Zoom In Zoom Out Reset image size