New Neutron Star Equation of State with Quark-Hadron Crossover

TL;DR

This paper introduces an improved neutron star equation of state, QHC19, featuring a smooth hadron-quark crossover, aligning well with observational data and providing new insights into dense matter physics.

Contribution

The paper develops a novel neutron star equation of state with a continuous hadron-quark transition, integrating advanced models for both regimes and constraining key interaction parameters.

Findings

QHC19 matches GW170817 data on neutron star properties.

Maximum neutron star mass predicted is 2.35 solar masses.

The model provides a consistent description of matter at all densities.

Abstract

We present a much improved equation of state for neutron star matter, QHC19, with a smooth crossover from the hadronic regime at lower densities to the quark regime at higher densities. We now use the Togashi et al.~equation of state (Togashi:2017), a generalization of the Akmal-Pandharipande-Ravenhall equation of state of uniform nuclear matter, in the entire hadronic regime; the Togashi equation of state consistently describes non-uniform as well as uniform matter, and matter at beta equilibrium without the need for an interpolation between pure neutron and symmetric nuclear matter. We describe the quark matter regime at higher densities with the Nambu--Jona--Lasinio model, now identifying tight constraints on the phenomenological universal vector repulsion between quarks and the pairing interaction between quarks arising from the requirements of thermodynamic stability and causal propagation of sound. The resultant neutron star properties agree very well with the inferences of the LIGO/Virgo collaboration, from GW170817, of the pressure vs. baryon density, neutron star radii, and tidal deformabilities. The maximum neutron star mass allowed by QHC19 is 2.35 $M_\odot$, consistent with all neutron star mass determinations.

Figures

Equations

$$\frac{E}{A}$$

$$\tilde{\Lambda }=\frac{16}{13}\frac{(M_1+12M_2)M_1^4{\Lambda }_1+(12M_1+M_2)M_2^4{\Lambda }_2}{(M_1+M_2{)}^5},$$

Full text

††slugcomment: Submitted to ApJ

1. Introduction

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 2$n_{0}$ in the hadronic interior, where $n_{0}$ is the nuclear saturation density, $\sim$ 0.16 fm*-3*, corresponding to $2.7\times 10^{14}$g cm*-3*); 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 Fig. 8 below), and is fully consistent with their radius and tidal deformability determinations (Abbott et al. 2018a ; Abbott et al. 2018b ). It also allows for a maximum neutron star mass $\simeq$ 2.35 M⊙, consistent with the recent observation of a 2.14 $\pm$ 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 subsec, 2.1, the Nambu–Jona-Lasinio (NJL) model for quark matter at higher density in subsec, 2.2, and the interpolation region in subsec. 2.3. In Sec. 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 Sec. 4, we discuss the structure of neutron stars described by the QHC19 equation of state, and conclude in Sec. 5.

2. Physics of QHC19

At very low baryon density the system is characterized by non-relativistic 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 Nambu–Jona–Lasinio 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_{0}$. 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 $\alpha=1-2Y_{\rm p}$: Here $Y_{\rm p}=n_{{\rm p}}/n_{\rm B}$, with $n_{\rm B}=n_{\rm n}+n_{\rm 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_{\rm B}\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle>}}{{\sim}}>}5.5n_{0}, 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 non-uniform 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_{\rm p}$ without an interpolation between PNM and SNM. Specifically we use the Togashi equation of state from density111In QHC18, we did in fact use the Togashi equation of state up to density 0.26 $n_{0}$, to construct a smooth match to APR extrapolated to beta equilibrium at higher density, despite APR not describing the inhomogeneous matter in the crust. 4.74 $\times 10^{-10}n_{0}$ to 2 $n_{0}$.

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_{0})/3n_{0}$ and $\alpha=1-2Y_{{\rm p}}$, following the convention in the CompOSE repository:

[TABLE]

The saturation energy $E_{0}$ of SNM is 16.1 MeV in Togashi versus. 16.0 MeV in APR, and the incompressibility $K_{0}$ = 245 vs. 267 MeV in APR. Near $n_{0}$, 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 vs. 34.0 MeV in APR, and the slope parameter, $L$ = 38.7 vs. 63.2 MeV in APR.222The APR numbers cited above are calculated from the fit function in the original APR paper (Akmal et al. (1998)). 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_{\odot}$ neutron star containing nuclear matter only (where in the crust in APR we use the Togashi equation of state for $n_{{\rm B}}\leq 0.26n_{0}$). The radii are $R$ = 11.6 km for Togashi vs. 11.5 km for APR, and the tidal deformability, $\Lambda$ = 360 vs. 268, In contrast, QHC19 (set D, Fig. 1) gives $R$ = 11.6 km and $\Lambda$ = 350 for a 1.4 $M_{\odot}$ neutron star.

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 $\simeq 0.625\ n_{0}$, found by comparing the energy per nucleon, for given average $n_{\mathrm{B}}$ and $Y_{\mathrm{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_{\mathrm{B}},Y_{\mathrm{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_{\rm B}/n_{0}$, from $4.7\times 10^{-10}$ to 38 $n_{0}$, $Y_{\rm p}$ from 0 to 0.65; and temperature $T$ from 0 to $\sim 400$ MeV.

2.2. Quark degrees of freedom

In cold dense matter above baryon densities some 5 $n_{0}$, quarks become the dominant degrees of freedom. The baryon chemical potential $\mu_{\rm B}$ is in the range $\sim$ 1.5-2 GeV, where the QCD strong coupling constant $\alpha_{s}$ is of order unity, too large to allow a perturbation theory calculation of the properties of quark matter. In this strong coupling regime of QCD, one must take into account non-perturbative effects including the generation of constituent quark masses owing to chiral symmetry breaking (strong quark–antiquark pairing; Hatsuda & Kunihiro (1994)), quark–quark pairing leading to color superconductivity (Alford et al. (2008)), and the mutual interplay of these effects (Hatsuda et al. (2006)). (Although QCD is an asymptotic-free gauge theory, the interactions between quarks and gluons become weak only at very high baryon densities (Collins et al. (1975)), well beyond those in neutron stars.)

As in QHC18, we describe the non-pertubative regime at baryon densities above $\sim 5n_{0}$ by the Nambu-Jona Lasinio (NJL) effective model (reviewed in Ref. Buballa (2005)), in which gluon degrees of freedom are taken into account only implicitly. The model contains quarks ($q$) of three flavors – up, down, and strange ($u,d,s$) – and of three colors, a flavor-dependent current mass $m_{q}$, as well as four-Fermi interactions in the scalar-pseudoscalar channel, the vector channel and the diquark channel, with coupling strengths $G$, $g_{V}$, and $H$, respectively. The Lagrangian describing these interactions is discussed in detail in Baym et al. (2018) where the vector repulsion, $g_{V}$, is needed for quark matter to be able to support heavy neutron stars. In perturbative QCD, $G$, $g_{V}$ and $H$ are all related to the single gluon exchange process between quarks, while in the non-perturbative regime, they are taken as independent parameters. The model also includes the (instanton-induced) six-quark interaction with coupling strengths $K$ and $K^{\prime}$ (Baym et al. (2018)). The new physics in QHC19 is that we now provide tight constraints on the parameters $g_{V}$ and $H$.

The calculation of the ground state energy in the NJL model, a sum of quark and lepton contributions, is reviewed in Baym et al. (2018). We take the parameters $G$, $K$ as well as the bare quark mass $m_{i}$ and the momentum cutoff $\Lambda_{\rm NJL}$ in the NJL model to have their standard vacuum values (the HK parameter set in Table I of Baym et al. (2018)). While in the standard description of the NJL model, the parameters $g_{V}$ and $H$ are determined by a Fierz transform of the one-gluon exchange interaction, we treat them as independent here as the equation of state and the subsequent neutron star structure are, as we see below, rather sensitive to their values. The effect of $K^{\prime}$ can be absorbed in the variation of $H$ and $g_{V}$ as far as the equation of state is concerned, so that we take $K^{\prime}=0$ throughout.

2.3. Quark–hadron crossover

We now discuss the transition from hadronic to quark degrees of freedom, a problem by no means satisfactorily resolved. The simplest picture is that of quark–hadron continuity or a quark–hadron crossover, discussed from the points of views of quark percolation (Baym (1979); Çelik et al. (1980); Satz (1998)), of QCD symmetry breaking (Schafer & Wilczek (1999); Hatsuda et al. (2006)), and of neutron star phenomenology (Masuda et al. (2013)). As the baryon density increases in this scenario, matter evolves from hadronic to strongly interacting quark matter with color superconductivity, without a discontinuous jump (or perhaps with a small jump) in the energy density. This picture is in strong contrast to the more conventional description of the transition from hadronic matter to weakly interacting quark matter via a first order transition (Baym & Chin (1976); Alford et al. (2013)). The quark–hadron crossover allows a relatively stiff equation of state capable of supporting neutron stars of masses greater than 2 $M_{\odot}$ with a substantial quark core (Masuda et al. (2013); Kojo et al. (2015)).

Since neither purely hadronic nor purely quark matter descriptions are reliable in the range 2-5 $n_{0}$, we construct, as in QHC18, the equation of state by a smooth interpolation between hadronic matter at $n_{\rm B}\leq 2n_{0}$ and quark matter at $n_{\rm B}\geq 5n_{0}$, expressing the thermodynamics in terms of the pressure, $P$, as a function of the baryon chemical potential, $\mu_{\rm B}$. (In the Appendix we examine interpolations with different choices of boundaries; with the boundary of the hadronic matter fixed to $2n_{0}$, and for the quark matter, to 5, 6, and 7 $n_{0}$.) The interpolation procedure,333We take a smooth interpolation as a baseline, but one can accommodate a first order phase transition by perturbing the pressure curve with a kink, although too large a first order transition tends to violate either the causality or two-solar mass constraints (Baym et al. (2018)). We emphasize that the interpolation uses the information on hadronic and quark matter only as boundary conditions, and does not extrapolate the equations of state of hadronic or quark matter into the interpolation region. Generally, extrapolated expressions for the thermodynamics, as well for example as transport coefficients, can behave unphysically and should be avoided where possible. reviewed in Baym et al. (2018), must satisfy the constraints obeyed by the pressure, that its first derivative, $\partial P/\partial\mu_{\rm B}=n_{\rm B}$, the baryon density, must be positive, and that it must be convex for all $\mu_{\rm B}$, i.e., $\partial^{2}P/\partial\mu_{\rm B}^{2}=\partial n_{\rm B}/\partial\mu_{\rm B}>0$; in addition the adiabatic sound velocity, $c_{s}=\sqrt{\partial P/\partial\varepsilon}=\sqrt{\partial\ln\mu_{\rm B}/\partial\ln n_{\rm B}}$, cannot exceed the speed of light, $c$. Given the hadronic pressure at $n_{\rm B}\leq 2n_{0}$, these constraints allow only a limited range of quark pressures at $n_{\rm B}\geq 5n_{0}$, and furthermore tightly constrain the strength of the universal quark repulsion, $g_{V}$, and the diquark pairing strength, $H$. Later we further impose the maximum neutron star mass constraint.

3. Constraints on Parameters

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}^{{\rm max}}/c)^{2}$ as a function of $g_{V}$ and $H$; the resolution in the figure is $\Delta g_{V}=\Delta H=0.01G$. The uncolored area violates either causality or thermodynamic stability (i.e., leads to an inflection point in $P(\mu_{\rm B})$ and is unphysical). As the figure indicates $g_{V}$ and $H$ are strongly correlated. The causality constraint bounds $g_{V}$ from above by $\simeq 1.3\,G$, implying that within this NJL description quark matter equations of state for $n_{\rm B}\geq 5n_{0}$ cannot be too stiff.

As we see in Fig. 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,444Parameters very close to the boundary of the colored region are unlikely in the sense that for a slight extension to finite temperatures or different lepton fractions the equation of state could easily violate the physical constraints. with $H$ chosen to give the least $c_{s}^{{\rm 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 Fig. 1.

In Fig. 2 we show the density dependence of $c_{s}^{2}$ 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_{\rm B}\gtrsim 5.5n_{0}$ is cured by this replacement. The bump in the QHC19 curves for $n_{\rm B}$ = 3-5 $n_{0}$ 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)).)

4. Neutron star structure

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_{{\rm max}}$, for given $g_{V}$ and $H$. The absolute maximum allowed mass is $\simeq$ 2.35 $M_{\odot}$ at $(g_{V},H)/G\simeq(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\gtrsim 0.60$-$0.75$, where the detailed value depends on $H$. The recently determined neutron star of mass 2.17$\pm$ 0.1 $M_{\odot}$ (Cromartie et al. (2019)) tightens the constraint to $g_{V}/G\gtrsim 0.9$. Beyond the minimal value of $H$, the maximum mass is reduced by increasing $H$ until the equations of state become acausal.

Figure 4 shows the central baryon density $n_{\rm B}^{{\rm c}}$ in units of $n_{0}$ for the maximum mass neutron star for each parameter set. The central density at $M_{{\rm max}}$ always exceeds $5.5n_{0}$. 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$.

The baryon density distribution of neutron stars for parameter set D is shown in Fig. 5. For $M>1.2M_{\odot}$ the core is beyond the pure nuclear regime. As $M$ increases to $\sim 2.0M_{\odot}$, the core density increases slowly to $\sim 4n_{0}$. Beyond this region the core density grows more rapidly, reflecting the peak structure of the sound velocity around $n_{\rm B}$ = 3-5 $n_{0}$ (Fig. 2); after passing the peak the equation of state softens.

The neutron star mass-radius relation is shown in Fig. 6 for the parameter sets A–D, together with the masses of PSR J1614–2230 ($M/M_{\odot}=1.928\pm 0.017$), PSR J0348–0432 ($M/M_{\odot}=2.01\pm 0.04$), and PSR J0740+6620 ($M/M_{\odot}=2.17\pm 0.10$). The curves for A–D largely overlap for $M\lesssim 1.4M_{\odot}$ where the corresponding density is $n_{\rm B}\lesssim 2.5$-$3n_{0}$. The radii at $M=1.4M_{\odot}$ are $\simeq 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\pm 1.4$ km (Abbott et al. 2018a ).

We look now at the dimensionless tidal deformability $\Lambda$ 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 $\Lambda$. The gravitational waveforms in a merger of two stars of masses $M_{1}$ and $M_{2}$ are sensitive to the combination,

[TABLE]

of the masses and $\Lambda_{i}$ of the individual neutron stars. The individual $\Lambda_{i}$ can be regarded as functions of the stellar $M_{i}$, and thus the independent variables in $\tilde{\Lambda}$ are $M_{1}$ and $M_{2}$.

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

Finally in Fig.8 we compare the QHC19 equation of state, $P$ vs. $m_{N}n_{\rm B}$, with that recently inferred by the LIGO/Virgo collaboration (Abbott et al. 2018a ). 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.

5. Summary and Conclusion

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 2 $n_{0}$, well describes neutron stars with masses up to 2.35 $M_{\odot}$, 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 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. (2001); Caplan et al. (2018)).

Hyperonic degrees of freedom play a role at higher hadronic densities. Ref. Togashi et al. (2016) extends the Togashi equation of state to hyperonic nuclear matter containing $\Lambda$ and $\Sigma^{-}$ hyperons at zero temperature. This study finds that the $\Lambda$ hyperon appears at 0.42 fm*-3* 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); Song et al. (2019).

In the quark regime, one needs to go beyond the phenomenological NJL model, 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 $\alpha_{s}(\mu_{\rm 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 non-perturbative 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\sim$ 400 MeV. On the other hand, thermal effects in the quark regime are highly suppressed by the pairing gaps up to $T\sim 100$ MeV. A finite temperature extension of the QHC19 equation of state in terms of the Helmholtz free energy, $F(n_{\mathrm{B}},Y_{\mathrm{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_{{\rm QM}}=5,6,$ and $7n_{0}$ while fixing the boundary of the Togashi equation of state at $n_{{\rm NM}}=2n_{0}$. The results are shown in Fig. 9. Overall, the extension of the interpolation range relaxes the constraint on the model parameters $g_{V}$ and $H$ extending their acceptable ranges by $\sim 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_{{\rm QM}}=5n_{0}$ but not for $7n_{0}$). 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_{0}$ to $7n_{0}$ does not significantly affect our conclusions in the main text. On the other hand, we find that reducing $n_{{\rm QM}}$ from $5n_{0}$ to $3$-$4n_{0}$, 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. We have looked at this variation to get a first impression of its effect. A more complete study is desirable but beyond the scope of this paper.