Ab-initio QCD calculations impact the inference of the neutron-star-matter equation of state

We demonstrate that ab-initio calculations in QCD at high densities offer significant and nontrivial information about the equation of state of matter in the cores of neutron stars, going beyond that which is obtainable from current astrophysical observations. We do so by extrapolating the equation of state to neutron-star densities using a Gaussian process and conditioning it sequentially with astrophysical observations and QCD input. Using our recent work, imposing the latter does not require an extrapolation to asymptotically high density. We find the QCD input to be complementary to the astrophysical observations, offering strong additional constraints at the highest densities reached in the cores of neutron stars; with the QCD input, the equation of state is no longer prior dominated at any density. The QCD input reduces the pressure and speed of sound at high densities, and it predicts that binary collisions of equal-mass neutron stars will produce a black hole with greater than $95\%$ ($68\%$) credence for masses $M \geq 1.38 M_\odot$ ($M \geq 1.25 M_\odot$). We provide a Python implementation of the QCD likelihood function so that it can be conveniently used within other inference setups.


INTRODUCTION
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;Cromartie et al. 2019;Fonseca et al. 2016Fonseca et al. , 2021, combined mass and radius measurements from Xray pulse profiling Miller et al. 2019;Riley et al. 2021;Miller et al. 2021), as well as measurements of the tidal deformability from gravitational-wave observations of binary NS mergers (Abbott et al. 2017b(Abbott et al. , 2018(Abbott et al. , 2019(Abbott et al. , 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 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 −3 (Hebeler et al. 2013;Drischler et al. 2020). On the other hand, the EoS can also be computed directly 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. 2018Gorda et al. , 2021aGorda & 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;Dietrich et al. 2020;Capano et al. 2020;Raaijmakers et al. 2020;Al-Mamun et al. 2021;Miller et al. 2020;Huth et al. 2022;Lim & Holt 2022) or non-parametric extensions (Landry & Essick 2019;Landry et al. 2020;Essick et al. 2020;Miller et al. 2021;Essick 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 Altiparmak et al. 2022;Kurkela et al. 2014;Annala et al. 2018Annala et al. , 2020Annala et al. , 2022Jiang et al. 2022).
The results from analysis with and without the QCD input differ. Notably, the works with the QCD in- The bands correspond to 67%-credible intervals conditioned with different inputs. Pulsars+ Λ is conditioned with mass measurements of PSR J0348+0432 and PSR J1614−2230, with the combined mass and radius measurement of PSR J0740+6620, as well as the measurement of binary tidal deformability Λ from GW170817. Pulsars+ Λ+BH further assumes that the binary merger product of GW170817 was a black hole. Pulsars+ Λ+QCD additionally includes the QCD input, which excludes the largest pressures at high densities. The effect of the black-hole hypothesis is negligible if the QCD input is included. The grey band denotes a 67%credible interval for the maximal densities reached in stable non-rotating NSs.
put report a strong softening of the EoS and lower speed of sounds at energy densities of the order ε ∼ 750 MeV/fm 3 , not discernible in works without the QCD input nor present in nuclear models. These features have been interpreted as the onset of a quarkmatter 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 two 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 modelindependent 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 Fig. 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.

SETUP
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 where P (data | EoS) is the likelihood of an EoS given the data.
Here we use data that is provided by astrophysical observations as well as from theoretical calculations in QCD. These inputs are independent of each other and the likelihood factorises into 4 terms: (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 mass and radius of PSR J0740+6620 using the 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 We now discuss the prior and the novel QCDlikelihood function in detail, and give a brief description of how we implement each of the astrophysical measurements.
2.1. Prior: P (EoS) Figure 2. A sample of 10k EoSs drawn from the prior distribution conditioned with the low-density CET EoS. The coloring of the individual EoSs corresponds to the likelihood assigned according to Eq. (2), imposing the information arising from all NS observations and the QCD input (i.e., Pulsars+ Λ+BH+QCD).
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, pressure as a function of number density p(n), instead of using only its reduced form as a function of energy density p(ε). For this reason, our GP regression differs from that of Landry & Essick (2019); Essick et al. (2020). Namely, we choose to extrapolate an auxiliary variable related to the speed of sound c s , as a function of baryon number density n instead of as function of energy density ε. The logarithm in the above expression maps the causal range of c 2 s ∈ [0, 1] to the range ϕ(n) ∈ [−∞, ∞], suitable for a GP.
Before conditioning the GP, we choose ϕ(n) to be normally distributed ϕ(n) ∼ N − ln 1/c 2 s − 1 , K(n, n ′ ) , with a Gaussian kernel K(n, n ′ ) = ηe −(n−n ′ ) 2 /2l 2 . The three hyperparameters-the variance η, the correlation length l, and the mean speed of sound squaredc 2 sdetermine 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 The GP is then conditioned with the CET EoS, leading to the c 2 s (n) displayed in Fig. 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.
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 ε c 2 s (n) = 1 e −ϕ(n) + 1 (7) and the pressure p(n) = −ε(n) + µ(n)n.
We then construct sequences of non-rotating 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 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.

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 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 ⃗ β ≡ {p, n, µ}. 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,a).
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Λ, 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Λ ∝ µ, but the optimal coefficient is undetermined. Formally, the dependence on the renormalization scale is of 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 have 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 ⃗ β QCD (X) = {p QCD (µ H , X), n QCD (µ H , X), µ H }. This correlated set is evaluated at a fixed scale µ H , parameterized by the (dimensionless) renormalization scale X ≡ 3Λ/(2µ H ). The joint probability distribution for the set β at the scale µ H , defined as β H , is then obtained through a weighted average 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 with 1 S 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 (Schneider 2003;Rebhan & Romatschke 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 2 s,QCD (µ 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. , 2021a. 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 ⃗ β be able to reach the QCD  values at high densities using a stable and causal EoS imposes non-trivial constraints on the allowed values these quantities can take at NS densities. Here, we implement these constraints by excluding any GP extrapolations whose endpoint ⃗ β L = {p L , n L , µ L } 2 , 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 ⃗ β 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 ∆p min and ∆p max that can be reached with a thermodynamically stable EoS whose speed of sound squared does not exceed c 2 s,lim . For a given ⃗ β L the minimal and maximal pressure differences 2 Note that here the ⃗ β-triplet with subscript L indicates the last provided point of the EoS, while in Komoltsev & Kurkela (2022) the subscript L denoted the CET calculation. The reason is that we are effectively taking the last provided point of the EoS as the new low-density point.
are (Komoltsev & Kurkela 2022) ∆p min (c 2 s,lim ) = c 2 s,lim Note that the functions ∆p min and ∆p max implicitly depend on µ L and µ H , though we do not show this dependence in their arguments. Furthermore, when taking c 2 s,lim = 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 ⃗ β L leads to the condition ∆p / ∈ [∆p min , ∆p max ]. Accounting for the scale-variation error, the QCD likelihood function reads This integral is evaluated by substituting in Eq. (10) and integrating over ⃗ β H to resolve the δ-function. In practice we perform Monte-Carlo integration by randomly drawing X values from the distribution w(ln X) and counting the frequency of the last extrapolated point ⃗ β 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 2 s,lim < 1 where P 0 (c 2 s,lim ) is a prior distribution. Given that the QCD value for c 2 s is very close to its conformal value c 2 s = 1/3, we may average the weight over different values of c 2 s,lim with, e.g., a flat prior function for the maximal speed of sound reached between n L < n < n H P 0 (c 2 s,lim ) = 1 [1/3,1] (c 2 s,lim ).
In this work we will concentrate solely on the most conservative prediction given by Eq. (14) with c 2 s,lim = 1 and leave the exploration of the less-conservative likelihood functions for future studies.

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 with M min = 0.5M ⊙ . With this, our full prior is given by

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 nuisancemarginalized) 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 The final likelihood is then the product of two likelihoods for PSR J0348+0432 and PSR J1614−2230.

Radius measurements
In the case of PSR J0740+6620, we construct our likelihood function using the (nuisance-marginalized) joint

Tidal deformability
The LIGO/Virgo Collaborations give the accurate measurement of the chirp mass as well as the joint posterior density function P (d GW | Λ, q) (see Fig. 12 of Abbott et al. (2019), lowspin prior) for the binary-tidal-deformability parame-ter Λ and mass ratio of the binary merger components q = M 1 /M 2 , marginalized over all the other nuisance parameters (including the chirp mass). For this measurement, we require a joint mass prior for the binary P 0 (M 1 , M 2 | EoS), which we take as Here, the integral may be calculated by integrating from the q = 1 value of M 2 ≈ 1.362(1)M ⊙ .

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 modelling of NS mergers, consistent with the formation of a BH (Margalit & Metzger 2017;Rezzolla et al. 2018;Ruiz et al. 2018;Shibata et al. 2017Shibata 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 1 + N 2 = N remnant (q) + N ejecta (q), where N 1 and N 2 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 (non-rotating) 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 (cf. Eq. (22)) 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 1 and M chirp . This condition is however not independent from the constraint from the tidal deformability and both likelihoods depend on the same mass ratio q since it concerns the same event. Hence when using the BH hypothesis, we must impose both inputs simultaneously by defining the likelihood as

RESULTS
We have generated a prior sample of 120k EoSs 4 from the prior process, and for each of these EoS we have computed the corresponding likelihood functions appearing in Eq. (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, where here the likelihood includes only the selected inputs (these resampled ensembles are used only in Fig. 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); Raaijmakers et al. (2019). However, we find the comparison of posterior samples more intuitive. 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 Fig. 3a). 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 Λ, the combined pulsar measurements, and the QCD input all offer complementary information to one another, illustrated in Fig. 3b. Finally, we find that the input that most corroborates the QCD input is the BH hypothesis, see Fig. 3c. 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 Fig. 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).
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 Λ measurement from GW170817 also drives the ε-p process towards softer values but limits the pressure at signif-icantly smaller densities than the QCD input. This illustrates the complementary nature of these inputs seen earlier in Fig. 3b. 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 Fig. 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+disc and NS+disc systems. The effect of the QCD constraint becomes dominant at densities above ε ∼ 750 MeV/fm 3 (450 MeV/fm 3 , 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 Fig. 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 Fig. 6 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 Fig. 5. 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 poste- Pulsars+ +QCD Pulsars+ Prior Figure 8. The posterior probability that a binary-NS merger with chirp mass M chirp and mass ratio q leads to a BH. The purple area shows the prediction using information from pulsar and gravitational-wave observations but does not make an assumption about the binary merger product of GW170817. Including the QCD input (green area) leads to the prediction that a binary merger with M chirp = 1.186 leads to a BH with 90% credence.
rior 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 speed of sound c 2 s ≲ 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 Fig. 5. 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 2 s ≲ 1/3, is used in the conditioning.
Despite these changes on the ε-p and c 2 s processes, the QCD input has only a small effect in the M -R plane, as seen in Fig. 7a. 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 Fig. 7b.
As remarked above, after imposing the pulsar and Λ 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 Sec. 2.3.4 above for different potential binaries characterized by parameters M chirp and q, and show our results in Fig. 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 ⊙ ).
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] finding twin-star solutions. These works utilize more specific EoS parametrizations including explicit strong firstorder transitions, which are approximated by smooth but rapidly changing functions and have a small prior weight in our analysis.

CONCLUSIONS
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 3 .
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 non-trivial astrophysical modelling 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, non-overlapping 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 et al. 2023) 5 (Komoltsev 2022. Finally, after completing this study, another work appeared studying the impact of the QCD input on the NS EoS inference (Somasundaram et al. 2022a). 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 (Somasundaram et al. 2022b). Pushing the construction defined in Sec. 2.2 down to lower den-sity 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 phasetransition-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 ).