Ab-initio QCD Calculations Impact the Inference of the Neutron-star-matter Equation of State

The past years have seen great interest in the properties of ultradense nuclear matter in the cores of neutron stars (NSs). Driven by the fast progress in astronomical observations, many works have attempted to infer the properties of the equation of state (EoS) of NS matter using different combinations of observational and theoretical inputs. From the observational side, the EoS is most constrained by the discovery of massive NSs from pulsar timing (Demorest et al. 2010; Antoniadis et al. 2013; Fonseca et al. 2016; Cromartie et al. 2019; Fonseca et al. 2021), combined mass and radius measurements from X-ray pulse profiling (Miller et al. 2019; Riley et al. 2019; Miller et al. 2021; Riley et al. 2021), as well as measurements of the tidal deformability from gravitational-wave observations of binary NS mergers (Abbott et al. 2017b, 2018, 2019, 2020) and their associated electromagnetic counterparts (Abbott et al. 2017c).

In addition to inputs from NS observations, the EoS is also informed by ab-initio theoretical calculations. On the one hand, calculations within the framework of the chiral effective theory (CET) constrain the EoS at low baryon densities n, i.e., at and around the nuclear saturation density n ≲ 1.1n_s , where n_s ≈ 0.16 fm⁻³ (Hebeler et al. 2013; Drischler et al. 2020). On the other hand, the EoS can also be directly computed from the fundamental theory of quantum chromodynamics (QCD). These computations become reliable at large densities n ≳ 40n_s , well above the densities reached in NSs (Freedman & McLerran 1977; Fraga & Romatschke 2005; Kurkela et al. 2010; Kurkela & Vuorinen 2016; Gorda et al. 2018, 2021a, 2021b; Gorda & Säppi 2022).

The inference of the EoS from these astrophysical observations and theoretical inputs often proceeds by forming large ensembles of model-agnostic EoSs that are constrained or conditioned using the different inputs (Hebeler et al. 2013; Kurkela et al. 2014). The majority of these EoS-inference setups have used different parametric (Hebeler et al. 2013; Tews et al. 2018; Capano et al. 2020; Dietrich et al. 2020; Miller et al. 2020; Raaijmakers et al. 2020; Al-Mamun et al. 2021; Huth et al. 2022; Lim & Holt 2022) or nonparametric extensions (Landry & Essick 2019; Essick et al. 2020; Landry et al. 2020; Essick et al. 2021; Miller et al. 2021) of the CET EoS to NS star densities n ∼ 5–10n_s . There are also fewer works that additionally include the fundamental, high-density QCD constraint and interpolate through the full density range from 1.1n_s to 40n_s (Kurkela et al. 2014; Annala et al. 2018; Most et al. 2018; Annala et al. 2020; Altiparmak et al. 2022; Annala et al. 2022; Jiang et al. 2023).

The results from analysis with and without the QCD input differ. Notably, the works with the QCD input report a strong softening of the EoS and lower speed of sounds at energy densities of the order ε ∼ 750 MeV/fm³, not discernible in works without the QCD input nor present in nuclear models. These features have been interpreted as the onset of a quark-matter phase in Annala et al. (2020). It is, however, an open question whether these features occurring at NS densities are genuine predictions of QCD or whether they arise from attempting to interpolate the EoS through 2 orders of magnitude in density with a too restrictive interpolation function. Whether the QCD computations offer a significant and robust constraint to the EoS at NS densities determines whether QCD should be included in any complete EoS-inference setup.

Here, we propose a new strategy to incorporate the QCD input that avoids the pitfalls of the interpolation. This advancement is made possible by a recent work (Komoltsev & Kurkela 2022) that demonstrated that the requirement of causality and stability of the EoS imposes global constraints that feed the high-density QCD information down to lower densities in a completely model-independent way. It was found that the EoS is impacted by QCD down to n ∼ 2.2n_s before including NS observations. What is not known, however, is whether QCD offers significant information beyond the constraints already provided by NS observations.

In this work we construct a simple Bayesian-inference setup, which we use to study the interaction of the QCD input with the NS observations. We extrapolate the CET EoS to n = 10n_s using a Gaussian-process (GP) regression and condition the prior using NS observations with or without QCD input. Our results are exemplified in Figure 1, which demonstrates the additional information gained by including the QCD input in conditioning the pressure p as a function of energy density ε especially at high densities, leading to the aforementioned softening seen in previous works.

Zoom In Zoom Out Reset image size

In this section, we describe our framework based on GP regression that forms our prior P(EoS). Governed by Bayes' theorem, the prior is conditioned with the data to form a posterior process

$$\begin{eqnarray}&&P(\mathrm{EoS}\,| \,\mathrm{data})=\displaystyle \frac{P(\mathrm{EoS})\,P(\mathrm{data}\,| \,\mathrm{EoS})}{P(\mathrm{data})},\end{eqnarray} \tag{ 1 }$$

where P(data ∣ EoS) is the likelihood of an EoS given the data.

Here we use data that are provided by astrophysical observations as well as from theoretical calculations in QCD. These inputs are independent of each other, and the likelihood factorizes into four terms:

$$\begin{eqnarray}\begin{array}{rcl}P(\mathrm{data}\,| \,\mathrm{EoS}) & = & P(\mathrm{QCD}\,| \,\mathrm{EoS})P(\mathrm{Mass}\,| \,\mathrm{EoS})\\ & & \times \,P(\mathrm{NICER}\,| \,\mathrm{EoS})P(\widetilde{{\rm{\Lambda }}},\mathrm{BH}\,| \,\mathrm{EoS}).\end{array}\end{eqnarray} \tag{ 2 }$$

These likelihoods correspond to the aforementioned QCD input, mass measurements of PSR J0348+0432 and PSR J1624−2230 using pulsar timing, the simultaneous X-ray measurement of the mass and radius of PSR J0740+6620 using the Neutron-star Interior Composition Explorer (NICER) telescope, the tidal deformability measurement of GW170817 by LIGO/Virgo as well as the hypothesis supported by the electromagnetic counterpart of GW170817 that the resulting binary merger product is a black hole (BH). In many of the figures we consider the effect of pulsar measurements together, and we denote

$$\begin{eqnarray}&&P(\mathrm{Pulsars}\,| \,\mathrm{EoS})=P(\mathrm{Mass}\,| \,\mathrm{EoS})P(\mathrm{NICER}\,| \,\mathrm{EoS}).\end{eqnarray} \tag{ 3 }$$

We now discuss the prior and the novel QCD likelihood function in detail and give a brief description of how we implement each of the astrophysical measurements.

2.1. Prior: P(EoS)

In this section we describe our prior, which consists of an ensemble of EoSs that extrapolate the low-density EoS to relevant NS densities, up to n_L ≡ 10n_s . We choose to extrapolate the EoS using a GP regression (Landry & Essick 2019; Essick et al. 2020), which is conditioned with an (β-equilibrated) EoS computed in CET up to n = 1.1n_s (Hebeler et al. 2013). For densities below n = 0.57n_s we merge the GP to the BPS crust EoS (Baym et al. 1971).

As discussed in Komoltsev & Kurkela (2022), in order to make full use of the ab-initio calculations in QCD, it is necessary to reconstruct the full form of the EoS with the pressure as a function of number density p(n), instead of using only its reduced form as a function of the energy density p(ε). For this reason, our GP regression differs from that of Landry & Essick (2019) and Essick et al. (2020). Namely, we choose to extrapolate an auxiliary variable related to the speed of sound c_s ,

$$\begin{eqnarray}&&\phi (n)=-\mathrm{ln}\left(1/{c}_{s}^{2}(n)-1\right)\end{eqnarray} \tag{ 4 }$$

as a function of the baryon number density n instead of as function of the energy density ε. The logarithm in the above expression maps the causal range of ${c}_{s}^{2}\in [0,1]$ to the range ϕ(n) ∈ [−∞ , ∞ ], suitable for a GP.

Before conditioning the GP, we choose ϕ(n) to be normally distributed

$$\begin{eqnarray}&&\phi (n)\sim { \mathcal N }\left(-\mathrm{ln}\left(1/{\bar{c}}_{s}^{2}-1\right),K(n,n^{\prime} )\right),\end{eqnarray} \tag{ 5 }$$

with a Gaussian kernel $K(n,{n}^{{\prime} })=\eta {e}^{-(n-{n}^{^{\prime} }{)}^{2}/2{l}^{2}}$. The three hyperparameters—the variance η, the correlation length l, and the mean speed of sound squared ${\bar{c}}_{s}^{2}$—determine the shape of the EoS where no data are available. To allow for a large variety of EoSs in our prior, we construct a hierarchical model by drawing the hyperparameters themselves from judiciously chosen random distributions

$$\begin{eqnarray}\begin{array}{rcl}l & \sim & { \mathcal N }\left(1.0{n}_{s},{\left(0.25{n}_{s}\right)}^{2}\right),\quad \eta \sim { \mathcal N }(1.25,{0.2}^{2}),\\ {\bar{c}}_{s}^{2} & \sim & { \mathcal N }(0.5,{0.25}^{2}).\end{array}\end{eqnarray} \tag{ 6 }$$

The GP is then conditioned with the CET EoS, leading to the ${c}_{s}^{2}(n)$ displayed in Figure 2. We have taken the average and the difference of the "soft" and "stiff" EoS from Hebeler et al. (2013) as the mean and 90% credible interval of the training data used to condition the EoS.

Zoom In Zoom Out Reset image size

After conditioning, we draw a sample of EoSs from the GP and solve for the speed of sound c_s , the baryon chemical potential μ, the energy density ε

$$\begin{eqnarray}&&{c}_{s}^{2}(n)=\displaystyle \frac{1}{{e}^{-\phi (n)}+1}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\mu (n)={\mu }_{0}\exp \left({\int }_{{n}_{0}}^{n}{dn}^{\prime} {c}_{s}^{2}(n^{\prime} )/n^{\prime} \right)\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\varepsilon (n)={\varepsilon }_{0}+{\int }_{{n}_{0}}^{n}{dn}^{\prime} \mu (n^{\prime} ),\end{eqnarray} \tag{ 9 }$$

and the pressure p(n) = −ε(n) + μ(n)n.

We then construct sequences of nonrotating stars (M(n_central), R(n_central), Λ(n_central)) by solving the (perturbed) Tolman–Oppenheimer–Volkoff (TOV) equations (Hinderer 2008; Postnikov et al. 2010; Han & Steiner 2019) for the mass M, radius R, and tidal deformability Λ, with each EoS and for central densities up n_central = 10n_s . We identify the stable branches and find the maximal masses M_TOV supported by each EoS.

2.2. QCD Information: P(QCD∣EoS)

At high densities the asymptotic freedom of QCD allows one to access the EoS directly using perturbative calculations in QCD. These calculations give access to the full thermodynamic grand canonical potential as a function of the baryon number chemical potential Ω(μ), from which one can extract all other thermodynamic quantities. For our purposes, the most important physical quantities are the pressure p and the number density n evaluated at a given baryon chemical potential μ, which we group into a vector $\vec{\beta }\equiv \{p,n,\mu \}$. The current state-of-the-art result for the QCD grand canonical potential is the partial next-to-next-to-next-to-leading order (N3LO) perturbative-QCD computation of Gorda et al. (2021b, 2021a).

As is the case with any perturbative calculation, it is important to estimate the missing-higher-order (MHO) error arising from the truncation of the perturbative series at a finite order. This truncation renders the result dependent on a residual, unphysical renormalization scale $\bar{{\rm{\Lambda }}}$, which should be taken to be proportional to the physical scale in the problem to avoid large logarithms in the result. In the case of QCD at high density, one should take $\bar{{\rm{\Lambda }}}\propto \mu$, but the optimal coefficient is undetermined. Formally, the dependence on the renormalization scale is on the order of the first uncalculated term in series. This observation motivates the standard practice of estimating the missing-higher-order (MHO) terms by varying the renormalization scale around a fiducial scale by some fixed factor, though a Bayesian interpretation of the errors has been recently discussed in, e.g., Cacciari & Houdeau (2011), Duhr et al. (2021) ⁴ .

Here, we adopt the scale-averaging interpretation of Duhr et al. (2021) and interpret the perturbative result as a family of independent predictions for the correlated set ${\vec{\beta }}_{\mathrm{QCD}}(X)=\{{p}_{\mathrm{QCD}}({\mu }_{{\rm{H}}},X),{n}_{\mathrm{QCD}}({\mu }_{{\rm{H}}},X),{\mu }_{{\rm{H}}}\}$. This correlated set is evaluated at a fixed scale μ_H, parameterized by the (dimensionless) renormalization scale $X\equiv 3\bar{{\rm{\Lambda }}}/(2{\mu }_{{\rm{H}}})$. The joint probability distribution for the set β at the scale μ_H, defined as β_H, is then obtained through a weighted average

$$\begin{eqnarray}&&P({\vec{\beta }}_{{\rm{H}}})=\int d(\mathrm{ln}X)w(\mathrm{ln}X){\delta }^{(3)}({\vec{\beta }}_{{\rm{H}}}-{\vec{\beta }}_{\mathrm{QCD}}(X)).\end{eqnarray} \tag{ 10 }$$

In Cacciari & Houdeau (2011) it was argued that a physically motivated choice for the weight function w is a log-uniform distribution around a central scale

$$\begin{eqnarray}&&w(\mathrm{ln}X)={1}_{[\mathrm{ln}(1/2),\mathrm{ln}(2)]}(\mathrm{ln}X),\end{eqnarray} \tag{ 11 }$$

with 1_S being the indicator function on a set S. We adopt this distribution in the following. The range X ∈ [1/2, 2]—while in principle arbitrary—is suggested by phenomenological models (Rebhan & Romatschke 2003; Schneider 2003; Cassing 2007; Gardim & Steffens 2009) as well as the limit of QCD with a large number of flavors (Ipp & Rebhan 2003). Moreover, this range has been extensively used in previous studies and contains the minimally sensitive scales (defined through ∂_X f(X) = 0, for f(X)) of p_QCD(μ_H, X), n_QCD(μ_H, X) and ${c}_{s,\mathrm{QCD}}^{2}({\mu }_{{\rm{H}}},X)$ of the partial N3LO result used here.

As for the value of μ_H, we wish to use this perturbative information at a density that is as low as possible while staying in the regime where the MHO errors are under control. From these considerations, we adopt the scale μ_H = 2.6 GeV, which was chosen in Fraga et al. (2014), because the uncertainty estimation at this value is similar to that of CET at 1.1n_s (roughly ±24% variation around mean value). This has also been the standard choice in subsequent perturbative-QCD calculations (Gorda et al. 2018, 2021a, 2021b).

We now discuss how these QCD results can be used to construct a likelihood function that utilizes the perturbative calculation in a region where it is reliable but that can be imposed on an EoS at lower densities. The requirement that the triplet $\vec{\beta }$ be able to reach the QCD values at high densities using a stable and causal EoS imposes nontrivial constraints on the allowed values these quantities can take at NS densities. Here, we implement these constraints by excluding any GP extrapolations whose end point ${\vec{\beta }}_{L}=\{{p}_{L},{n}_{L},{\mu }_{L}\}$ ⁵ , with n_L ≡ 10n_s , cannot be connected to the corresponding QCD values with a stable and causal EoS. This requirement ensures that the EoS is consistent with the high-density limit also at all lower densities (assuming of course that the EoS is causal and stable up to ${\vec{\beta }}_{L}$). This very conservative requirement in effect allows the EoS to behave in a completely arbitrary way between the end point of the GP extrapolation and the density region where perturbative-QCD results are assumed to be reliable. In particular, it is allowed to have an arbitrary number of phase transition with arbitrary latent heats and may have arbitrarily large segments where the speed of sound is near the speed of light (without exceeding it).

As discussed in Komoltsev & Kurkela (2022), these constraints can be implemented in a robust way without specifying the form of the EoS in the range n_L < n < n_H. This is done by comparing the pressure difference between the last point of the extrapolated EoS and the perturbative-QCD result Δp = p_H −p_L to the minimal and maximal pressure differences ${\rm{\Delta }}{p}_{\min }$ and ${\rm{\Delta }}{p}_{\max }$ that can be reached with a thermodynamically stable EoS whose speed of sound squared does not exceed ${c}_{{\rm{s}},\ \mathrm{lim}}^{2}$. For a given ${\vec{\beta }}_{L}$ the minimal and maximal pressure differences are (Komoltsev & Kurkela 2022)

$$\begin{eqnarray}&&{\rm{\Delta }}{p}_{\min }({c}_{{\rm{s}},\mathrm{lim}}^{2})=\displaystyle \frac{{c}_{{\rm{s}},\ \mathrm{lim}}^{2}}{1+{c}_{{\rm{s}},\ \mathrm{lim}}^{2}}\left({\mu }_{{\rm{H}}}{\left(\displaystyle \frac{{\mu }_{{\rm{H}}}}{{\mu }_{L}}\right)}^{1/{c}_{{\rm{s}},\ \mathrm{lim}}^{2}}-{\mu }_{L}\right){n}_{L},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{p}_{\max }({c}_{{\rm{s}},\mathrm{lim}}^{2})=\displaystyle \frac{{c}_{{\rm{s}},\ \mathrm{lim}}^{2}}{1+{c}_{{\rm{s}},\ \mathrm{lim}}^{2}}\left({\mu }_{{\rm{H}}}-{\mu }_{L}{\left(\displaystyle \frac{{\mu }_{L}}{{\mu }_{{\rm{H}}}}\right)}^{1/{c}_{{\rm{s}},\ \mathrm{lim}}^{2}}\right){n}_{{\rm{H}}}.\end{eqnarray} \tag{ 13 }$$

Note that the functions ${\rm{\Delta }}{p}_{\min }$ and ${\rm{\Delta }}{p}_{\max }$ implicitly depend on μ_L and μ_H, though we do not show this dependence in their arguments. Furthermore, when taking ${c}_{{\rm{s}},\mathrm{lim}}^{2}=1$, we shall drop the argument for brevity. The most conservative QCD likelihood function is then constructed by assign a zero Bayesian weight to any EoS whose ${\vec{\beta }}_{L}$ leads to the condition ${\rm{\Delta }}p\notin [{\rm{\Delta }}{p}_{\min },{\rm{\Delta }}{p}_{\max }]$.

Accounting for the scale-variation error, the QCD likelihood function reads

$$\begin{eqnarray}&&P(\mathrm{QCD}\,| \,\mathrm{EoS})=\int d{\vec{\beta }}_{{\rm{H}}}P({\vec{\beta }}_{{\rm{H}}}){1}_{[{\rm{\Delta }}{p}_{\min },{\rm{\Delta }}{p}_{\max }]}({\rm{\Delta }}p).\end{eqnarray} \tag{ 14 }$$

This integral is evaluated by substituting in Equation (10) and integrating over ${\vec{\beta }}_{{\rm{H}}}$ to resolve the δ function. In practice we perform Monte Carlo integration by randomly drawing X values from the distribution $w(\mathrm{ln}X)$ and counting the frequency of the last extrapolated point ${\vec{\beta }}_{L}$ satisfying the condition for Δp.

The above conservative likelihood function does not in any way favor physically motivated EoSs in the density range n_L < n < n_H but merely asks if any causal interpolation exists. We may add further assumptions to the likelihood function by giving a higher weight to EoSs that can be connected to the QCD regime with an interpolation function that does not necessarily need to be extreme. One possible way to do so is to favor EoSs which can be connected to the QCD regime with EoSs that do not need to exceed some given limiting value of the speed of sound squared ${c}_{{\rm{s}},\ \mathrm{lim}}^{2}\lt 1$

$$\begin{eqnarray}\begin{array}{rcl}P(\mathrm{QCD}\,| \,\mathrm{EoS}) & = & \displaystyle \int {{dc}}_{{\rm{s}},\ \mathrm{lim}}^{2}{P}_{0}({c}_{{\rm{s}},\ \mathrm{lim}}^{2})\displaystyle \int d{\vec{\beta }}_{{\rm{H}}}P({\vec{\beta }}_{{\rm{H}}})\\ & & \times \,{1}_{[{\rm{\Delta }}{p}_{\min }({c}_{{\rm{s}},\mathrm{lim}}^{2}),{\rm{\Delta }}{p}_{\max }({c}_{{\rm{s}},\mathrm{lim}}^{2})]}({\rm{\Delta }}p),\end{array}\end{eqnarray} \tag{ 15 }$$

where ${P}_{0}({c}_{{\rm{s}},\ \mathrm{lim}}^{2})$ is a prior distribution. Given that the QCD value for c_s ² is very close to its conformal value ${c}_{s}^{2}=1/3$, we may average the weight over different values of ${c}_{{\rm{s}},\ \mathrm{lim}}^{2}$ with, e.g., a flat prior function for the maximal speed of sound reached between n_L < n < n_H

$$\begin{eqnarray}&&{P}_{0}({c}_{{\rm{s}},\mathrm{lim}}^{2})={1}_{[1/\mathrm{3,1}]}({c}_{{\rm{s}},\mathrm{lim}}^{2}).\end{eqnarray} \tag{ 16 }$$

In this work we will concentrate solely on the most conservative prediction given by Equation (14) with ${c}_{{\rm{s}},\mathrm{lim}}^{2}=1$ and leave the exploration of the less-conservative likelihood functions for future studies.

2.3. Astrophysical Observations

Here we briefly describe our implementations of the likelihoods for various NS measurements. We denote the various classes of observations d , and construct our likelihood functions by marginalizing over additional source-specific model parameters that we must introduce to characterize the properties of the stars. These in principle depend on the EoS and stellar populations, as well as on selection effects. In particular, in the following we will assume a flat prior distribution for the masses of stars, M_i , up to the EoS-dependent maximal mass M_TOV for each star i considered

$$\begin{eqnarray}&&{P}_{0}({M}_{i}\,| \,\mathrm{EoS})=\displaystyle \frac{{{\bf{1}}}_{[{M}_{\min },{M}_{\mathrm{TOV}}]}({M}_{i})}{{M}_{\mathrm{TOV}}-{M}_{\min }},\end{eqnarray} \tag{ 17 }$$

with ${M}_{\min }=0.5{M}_{\odot }$. With this, our full prior is given by P₀(M_i ∣ EoS)P(EoS). For a clear and more complete discussion of this point, we refer the reader to, e.g., Miller et al. (2020) and Landry et al. (2020) ⁶ .

2.3.1. Mass Measurements

We require that the maximal mass supported by the EoSs, M_TOV, is larger than the largest observed NS masses. We consider PSR J0348+0432 with M = 2.01 ± 0.04M_⊙ (Antoniadis et al. 2013) and PSR J1614−2230 with M = 1.928 ± 0.017M_⊙ (in order to avoid double counting with the NICER measurements we do not include PSR J0740+6620). For each of these stars i, we approximate the (per-source nuisance-marginalized) likelihood functions for the mass of the star M_i , P( d _M ∣ M_i ), given the measurements d _M as a normal distribution with central value M_i and standard deviation σ_i . We then give zero weight to an EoS if it cannot support a star of mass M_i , giving

$$\begin{eqnarray}&&\begin{array}{l}P({\mathrm{Mass}}_{i}\,| \,\mathrm{EoS})\propto \displaystyle \int {{dM}}_{i}P({d}_{{\boldsymbol{M}}}\,| \,{M}_{i}){P}_{0}({M}_{i}\,| \,\mathrm{EoS})\\ \ \ \approx \displaystyle \frac{1}{2({M}_{\mathrm{TOV}}-{M}_{\min })}\left(1+\mathrm{Erf}\left(\displaystyle \frac{{M}_{\mathrm{TOV}}-{M}_{i}}{\sqrt{2}{\sigma }_{i}}\right)\right).\end{array}\end{eqnarray} \tag{ 18 }$$

The final likelihood is then the product of two likelihoods for PSR J0348+0432 and PSR J1614−2230.

2.3.2. Radius Measurements

In the case of PSR J0740+6620, we construct our likelihood function using the (nuisance-marginalized) joint posterior distribution for mass and (equatorial circumferential) radius P( d _NICER ∣ M, R) given the observations d _NICER determined using combined measurements of NICER and XMM-Newton observatories (Figure 1 (right) of Miller et al. 2021)

$$\begin{eqnarray}&&\begin{array}{l}P(\mathrm{NICER}\,| \,\mathrm{EoS})\\ \ \ \propto \,\displaystyle \int {{dM}}_{i}P({{\boldsymbol{d}}}_{\mathrm{NICER}}\,| \,{M}_{i},R({M}_{i})){P}_{0}({M}_{i}\,| \,\mathrm{EoS}).\end{array}\end{eqnarray} \tag{ 19 }$$

2.3.3. Tidal Deformability

The LIGO/Virgo Collaborations give the accurate measurement of the chirp mass

$$\begin{eqnarray}&&{{ \mathcal M }}_{\mathrm{chirp}}\equiv \displaystyle \frac{{\left({M}_{1}{M}_{2}\right)}^{3/5}}{{\left({M}_{1}+{M}_{2}\right)}^{1/5}}=1.186(1){M}_{\odot },\end{eqnarray} \tag{ 20 }$$

as well as the joint posterior density function $P({{\boldsymbol{d}}}_{{\boldsymbol{GW}}}\,| \,\widetilde{{\rm{\Lambda }}},q)$ (see Figure 12 of Abbott et al. (2019); low-spin prior) for the binary-tidal-deformability parameter $\widetilde{{\rm{\Lambda }}}$ and mass ratio of the binary merger components q = M₁/M₂, marginalized over all the other nuisance parameters (including the chirp mass).

For this measurement, we require a joint mass prior for the binary P₀(M₁, M₂ ∣ EoS), which we take as

$$\begin{eqnarray}&&{P}_{0}({M}_{1},{M}_{2}\,| \,\mathrm{EoS})=\displaystyle \frac{{1}_{[{M}_{2},{M}_{\mathrm{TOV}}]}({M}_{1}){1}_{[{M}_{\min },{M}_{\mathrm{TOV}}]}({M}_{2})}{1/2{\left({M}_{\mathrm{TOV}}-{M}_{\min }\right)}^{2}},\end{eqnarray} \tag{ 21 }$$

which is uniform on the set of possible binaries with ${M}_{1},{M}_{2}\in [{M}_{\min },{M}_{\mathrm{TOV}}]$ and M₁ ≥ M₂. Given the accurate measurement of ${{ \mathcal M }}_{\mathrm{chirp}}$, choosing a value for M₁ effectively enforces a δ function on M₂, determining the mass of the second merger component ${M}_{2}={M}_{2}({M}_{1},{{ \mathcal M }}_{\mathrm{chirp}})$. We then construct the likelihood function using the expression

$$\begin{eqnarray}\begin{array}{rcl}P(\widetilde{{\rm{\Lambda }}}\,| \,\mathrm{EoS}) & = & {\displaystyle \int }_{{M}_{1}\gt {M}_{2}}{{dM}}_{1}{P}_{0}({M}_{1},{M}_{2}\,| \,\mathrm{EoS})\\ & & \times \,P({{\boldsymbol{d}}}_{\mathrm{GW}}\,| \,\widetilde{{\rm{\Lambda }}},q({M}_{1},{M}_{2})).\end{array}\end{eqnarray} \tag{ 22 }$$

Here, the integral may be calculated by integrating from the q = 1 value of M₂ ≈ 1.362(1)M_⊙.

2.3.4. Black Hole Hypothesis

The final piece of observational input that we implement in our analysis is the hypothesis that the remnant in GW170817 collapsed to a BH. This conclusion is based on the short gamma-ray burst GRB170817A (Abbott et al. 2017a) associated with GW170817, whose properties are, in light of current numerical modeling of NS mergers, consistent with the formation of a BH (Margalit & Metzger 2017; Shibata et al. 2017; Rezzolla et al. 2018; Ruiz et al. 2018; Shibata et al. 2019). To implement this hypothesis, we proceed as follows.

Consider a possible initial binary with a mass ratio q for a fixed EoS. Throughout the merger process, the total baryon number of the system N(q) is conserved. Hence, N₁ + N₂ = N_remnant(q) + N_ejecta(q), where N₁ and N₂ are the baryon numbers of the two binary components. In order for the remnant to collapse to a BH for this value of q, it must be case that N_remnant(q) > N_TOV, the baryon number of the maximum stable (nonrotating) star. Following Annala et al. (2022), we ignore the small mass ejecta and consider a more conservative assumption that N(q) > N_TOV. To implement this condition within our Bayesian framework, we take as the likelihood (see Equation (22))

$$\begin{eqnarray}&&\begin{array}{l}P(\mathrm{BH}\,| \,\mathrm{EoS})={\displaystyle \int }_{{M}_{1}\gt {M}_{2}}{{dM}}_{1}{P}_{0}({M}_{1},{M}_{2}| \mathrm{EoS})\\ \ \ \times P({{\boldsymbol{d}}}_{{\boldsymbol{GW}}}\,| \,q({M}_{1},{M}_{2})){{\bf{1}}}_{[{N}_{\mathrm{TOV}},\infty ]}\left(N(q)\right),\end{array}\end{eqnarray} \tag{ 23 }$$

where P( d _GW ∣ q) is the 1D marginal posterior distribution for low-spin priors from Abbott et al. (2019). As above, here q is fixed by M₁ and ${{ \mathcal M }}_{\mathrm{chirp}}$.

This condition is however not independent from the constraint from the tidal deformability and both likelihoods depend on the same mass ratio q as it concerns the same event. Hence when using the BH hypothesis, we must impose both inputs simultaneously by defining the likelihood as

$$\begin{eqnarray}&&\begin{array}{l}P(\widetilde{{\rm{\Lambda }}},\mathrm{BH}\,| \,\mathrm{EoS})={\displaystyle \int }_{{M}_{1}\gt {M}_{2}}{{dM}}_{1}{P}_{0}({M}_{1},{M}_{2}\,| \,\mathrm{EoS})\\ \ \ \times \,P({{\boldsymbol{d}}}_{\mathrm{GW}}\,| \,\widetilde{{\rm{\Lambda }}},q({M}_{1},{M}_{2})){{\bf{1}}}_{[{N}_{\mathrm{TOV}},\infty ]}\left(N(q)\right).\end{array}\end{eqnarray} \tag{ 24 }$$

We have generated a prior sample of 120k EoSs ⁷ from the prior process, and for each of these EoS we have computed the corresponding likelihood functions appearing in Equation (2). Before turning to physical observables, we first discuss the extent to which the different constraints from the observations and QCD input either complement or corroborate each other. We do so by generating multiple posterior distributions, each conditioned with a subset of the inputs, and compare the overlap of these distributions.

In practice, for each selected subset of inputs, we resample the prior ensemble with an accept/reject step with probability proportional to the likelihood of the EoS,

$$\begin{eqnarray}&&\displaystyle \frac{P(\mathrm{data}\,| \,\mathrm{EoS})}{{\max }_{\mathrm{EoS}}[P(\mathrm{data}\,| \,\mathrm{EoS})]},\end{eqnarray} \tag{ 25 }$$

where here the likelihood includes only the selected inputs (these resampled ensembles are used only in Figure 3). The fraction of EoSs in such a posterior sample is then proportional to the marginal probability P(data) and illustrates how constraining the selected inputs are. We then examine the resulting posterior samples for different sets of selected inputs. If two posterior samples contain a large number of common EoSs, then the two observables corroborate each other, and the two inputs give similar information about the EoS; if the intersection is small, the two inputs offer complementary information about the EoS. Below, we report the fraction of EoSs remaining in a given posterior, normalized by the number of prior EoSs. An alternative way of characterizing this information would be, e.g., using the Kullback–Leibler divergence as done in Riley et al. (2019) and Raaijmakers et al. (2019). However, we find the comparison of posterior samples more intuitive.

Zoom In Zoom Out Reset image size

From our prior, 8% of the EoSs are consistent with combined NS observations (including the BH formation hypothesis). Similarly, 40% of the prior sample are consistent with QCD (see Figure 3(a)). Most importantly, only 4% of the EoSs are simultaneously allowed by both the NS observations and QCD input; that is, imposing the QCD input significantly constrains the EoS by removing half of the otherwise allowed EoSs in the posterior sample.

Further analyzing the overlap between different measurements, we find unsurprisingly that the pulsar mass measurements and the NICER measurement of PSR J0740+6620 give strongly overlapping information, with the latter being more constraining due to the information on the radius (not shown). We also find that the LIGO measurement of $\widetilde{{\rm{\Lambda }}}$, the combined pulsar measurements, and the QCD input all offer complementary information to one another, illustrated in Figure 3(b). Finally, we find that the input that most corroborates the QCD input is the BH hypothesis; see Figure 3(c). In fact, we see that after imposing the pulsar observations, all EoSs consistent with QCD are also consistent with the BH hypothesis in the resampled posterior distribution, indicating that the formation of a BH in GW170817 is a pre/postdiction of QCD. The effect of the BH hypothesis is negligible if the QCD input is included.

Moving on to physical observables, we show the impact of various constraints on the ε–p process in Figure 4, which displays samples drawn from our prior ensemble and colored by the different likelihood functions. We note in passing that considering only the QCD input (panel a) uses no astrophysical inputs and is thus dependent on neither general relativity nor the assumption that NSs are composed entirely of Standard Model matter. It therefore can be used as a starting point in studies testing these assumptions without using circular logic (Lope-Oter et al. 2019; Lope-Oter & Llanes-Estrada 2022).

Zoom In Zoom Out Reset image size

We see that the effect of the QCD input primarily decreases p(ε) for large ε, softening the EoS, which is qualitatively different from the effect from the pulsar input (panel d), which favors EoSs that are sufficiently stiff and limits the pressure from below. The $\widetilde{{\rm{\Lambda }}}$ measurement from GW170817 also drives the ε–p process toward softer values but limits the pressure at significantly smaller densities than the QCD input. This illustrates the complementary nature of these inputs seen earlier in Figure 3(b). As expected from the above discussion, the BH hypothesis corroborates the QCD input in the qualitative conclusion that the EoS must not be too stiff at high densities.

The combined effect of the different inputs is shown in Figure 1 above, where the prior is conditioned sequentially with different inputs. The QCD input offers a significant reduction of the uncertainty of the EoS, complementary to the astrophysical observations. This is particularly prominent if neglecting the BH hypothesis, which relies on models of jet formation from BH+disk and NS+disk systems. The effect of the QCD constraint becomes dominant at densities above ε ∼ 750 MeV fm⁻³ (450 MeV fm⁻³, if not including the BH hypothesis). Our result clearly demonstrates that the origin of the softening seen by previous studies utilizing the QCD constraint directly is indeed a robust prediction of QCD.

The gray band in Figure 1 corresponds to the density range reached in the cores of M_TOV stars. From the figure it is clear that the QCD input affects densities that can be reached within NSs. This is even more clearly visible in Figure 5 displaying a kernel density estimate of the pressure at four fixed energy densities. At the density n = 3n_s the QCD input mildly constrains the EoS, but starting already from n = 5n_s the effect of the QCD input is striking.

The softening due to QCD is also reflected in the posterior speed of sound, displayed in Figure 6. The NS observations primarily exclude EoSs that are too soft in the density range n ∼ 2–4n_s . At higher densities, the observations offer little information, and the posterior closely follows the prior. Imposing the QCD input however dramatically reduces the variation at densities n ≳ 3n_s forcing the equation of state to have speeds of sound ${c}_{s}^{2}\lesssim 0.6$. This is perhaps most dramatically seen in the distribution of the speed of sound in the centers of M_TOV stars, shown in Figure 6. Here, the QCD input strongly disfavors large speeds of sounds at the highest densities reached in NSs. Note that no information about the speed of sound of high-density QCD, ${c}_{s}^{2}\lesssim 1/3$, is used in the conditioning.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Despite these changes on the ε–p and c_s ² processes, the QCD input has only a small effect in the M–R plane, as seen in Figure 7(a). This renders it difficult to impose similar constraints on the EoS using M–R measurements only. The main effect of the QCD input in the M–R plane is to reduce the maximal masses reached in stable NSs; see Figure 7(b).

Zoom In Zoom Out Reset image size

As remarked above, after imposing the pulsar and $\widetilde{{\rm{\Lambda }}}$ inputs, the QCD input implies that the binary merger product in GW170817 was a BH. Similar predictions can be made for other potential binary NS events with different mass configurations. We conduct the same analysis as in Section 2.3.4 above for different potential binaries characterized by parameters ${{ \mathcal M }}_{\mathrm{chirp}}$ and q, and show our results in Figure 8. We find that the QCD input enhances the posterior probability of forming a BH in any given event, e.g., predicting that binary collisions of equal-mass NSs will produce a BH with greater than 95% (68%) credence for masses M ≥ 1.38M_⊙ (M ≥ 1.25M_⊙).

Zoom In Zoom Out Reset image size

Finally, we note that, while our prior ensemble contains a large number of EoSs that admit twin stars with multiple stable branches, the mass measurements give a vanishing weight to EoSs with multiple stable branches. We note however that this finding is not in contradiction with other works (e.g., Benic et al. 2015; Christian & Schaffner-Bielich 2020; Christian & Schaffner-Bielich 2022) finding twin-star solutions. These works utilize more specific EoS parameterizations including explicit strong first-order transitions, which are approximated by smooth but rapidly changing functions and have a small prior weight in our analysis.

In this work, we have characterized the impact from ab-initio calculations of the part of the Standard Model of particle physics describing the strong nuclear force, quantum chromodynamics, on the inference of the equation of state of ultradense matter in neutron-star cores. We have demonstrated that the QCD input significantly constrains the equation of state even after all the current astrophysical observations are taken into account, reducing uncertainties in the determination of the equation of state particularly at the highest densities reached in neutron stars. We found that imposing the QCD constraint requires the equation of state to significantly soften at energy densities ε ≳ 450 MeV fm⁻³.

Our work offers an explanation for the discrepancy between works that have observed the softening and those which have not as a result of either including the QCD input or not. By using the most conservative possible construction to impose the QCD input, we have demonstrated that this softening is a robust prediction of the fundamental theory of the strong interactions.

Of the observational inputs we have discussed, the black hole formation hypothesis has most overlap with QCD, with QCD giving a much stronger constraint of the two. The hypothesis requires nontrivial astrophysical modeling of the formation of the jet whose reliability may be difficult to quantify. Therefore is it is particularly fortunate that the QCD input offers an independent constraint supporting the conclusion arising from the black hole hypothesis. Similarly, future multimessenger observations of binary mergers with different masses offers a possible way to empirically test the QCD prediction.

We emphasize the importance of including all possible inputs when inferring the equation of state. All the measurements and theoretical calculations come with different, nonoverlapping systematics and uncertainties. Therefore, a convincing global understanding of the matter in the most extreme conditions is only reached through combining multiple, mutually consistent inputs.

Lastly, we note that a strategy that encompasses all the possible inputs is the one most likely to find a conflict between them. Assuming that all the systematics are under control, such a discrepancy would signal a failure of the underlying assumptions that NSs are described by the Standard Model and general relativity. This underlines the importance of including all possible reliable inputs in order fully exploit these extreme objects as a tool for fundamental discovery.

In order to facilitate an easy implementation of the QCD input to other, possibly more sophisticated, inference setups, we provide the QCD likelihood function used in this work in an easy-to-use Python implementation located on Zenodo (Komoltsev 2022; Komoltsev et al. 2023 ⁸ ).

Finally, after completing this study, another work appeared studying the impact of the QCD input on the NS EoS inference (Somasundaram et al. 2023). We are able to reproduce their results within our analysis by applying the constraints from Komoltsev & Kurkela (2022) at the maximum central density reached along a given EoS, rather than 10n_s , and they are able to reproduce our results within their analysis when applying the constraints at 10n_s (R. Somasundaram et al., private communication). Pushing the construction defined in Section 2.2 down to lower density obviously decreases the constraining power of the QCD input. We note that while this can be done, it requires the EoS to have significantly different behavior between the maximal central densities reached in NSs and 10n_s, compared to the prior used in this work and others. In particular, it requires a very strong phase-transition-like behavior and significant density ranges where the speed of sound is close to the speed of light. While this is interesting in its own right, we note that the effect of such phase-transition behavior is also not included for the low densities (for the first systematic studies see Gorda et al. 2022).

The authors thank Eemeli Annala, Bjorn Auestad, Aleksas Mazeliauskas, Joonas Hirvonen, Alex Nielsen, Joonas Nättilä Sanjay Reddy, Achim Schwenk, and Aleksi Vuorinen for useful discussions. The authors especially thank Tore Kleppe for discussions about Gaussian processes. This work is supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project ID 279384907—SFB 1245 and by the State of Hesse within the Research Cluster ELEMENTS (Project ID 500/10.006). The authors are listed in alphabetical order by last name.