Bayesian uncertainty quantification of perturbative QCD input to the neutron-star equation of state

The equation of state of neutron-star cores can be constrained by requiring a consistent connection to the perturbative Quantum Chromodynamics (QCD) calculations at high densities. The constraining power of the QCD input depends on uncertainties from missing higher-order terms, the choice of the unphysical renormalization scale, and the reference density where QCD calculations are performed. Within a Bayesian approach, we discuss the convergence of the perturbative QCD series, quantify its uncertainties at high densities, and present a framework to systematically propagate the uncertainties down to neutron-star densities. We find that the effect of the QCD input on the neutron-star inference is insensitive to the various unphysical choices made in the uncertainty estimation.


I. INTRODUCTION
The determination of the equation of state (EoS) of neutron-star (NS) cores is one of the grand questions of nuclear astrophysics [1,2]. The EoS determines many of the macroscopic properties of neutron stars and its features may give a unique inroad into determining the phase structure of Quantum Chromodynamics (QCD) at large baryon number densities [3][4][5][6].
At sufficiently high densities, well above the density range reached in stable NSs, perturbation theory gives a good approximation of the true EoS. It has furthermore been recently shown that these calculations-combined with the requirement that the EoS be mechanically stable, causal, and thermodynamically consistent (SCC) at all densities-give robust constraints to the EoS down to a few saturation densities n ∼ 2.3n s [44], with n s ≈ 0.16 fm −3 . The interaction between the astrophysical and the QCD constraints has also been studied, showing that the QCD input leads to a softening of the EoS at the highest densities reached inside the cores of stable NSs [3,4,33,[45][46][47][48][49] (cf. [50]). This feature has been inter-preted as a sign of loss of hadronic structure, and a phase change to quark matter [3][4][5]34].
The importance of the theoretical inputs in the EoS inference necessitates reliable and statistically interpretable uncertainty estimation of the calculations. In the low-density nuclear regime, theoretical uncertainty estimation including statistical uncertainty quantification has been an increasingly studied problem in recent years [1,[51][52][53][54][55][56][57][58][59][60][61]. Such quantified uncertainties are now routinely incorporated when inferring the NS-matter EoS. With this paper, we improve the treatment of corresponding uncertainties from the perturbative-QCD (pQCD) input used in the inference to produce statistically interpretable error estimates of the pQCD results and their impact at low densities.
The accuracy of the pQCD calculation is limited by the ignorance of the terms in the perturbative series beyond the last computed order, i.e., the missing-higher-order (MHO) terms. The current standard has been to estimate the MHO uncertainty of the result through its variation with respect to an unphysical renormalization scalē Λ. Explicitly, within theΛ MS renormalization scheme,Λ is varied by a factor of two around a central scale of 2µ/3, with µ the baryon chemical potential, to produce an uncertainty band for the result [37][38][39][40][41]. While this ad-hoc procedure is rooted in historical practice and experience in perturbative calculations, it lacks a well-defined statistical interpretation.
Recently, a Bayesian approach to estimate perturbative uncertainties has received attention in the high-energy community producing predictions for LHC physics [62][63][64][65][66]. These studies have differentiated between scale-variation uncertainties (i.e., those arising from setting the unphysical renormalization scale) and MHO corrections (i.e., those arising from the truncation of the perturbative series). For the MHO errors, various machine-learning-based models have been suggested to synthesize information from all the computed partial sums of the perturbative series (instead of solely from the highest order). Uncertainties due to scale-variation can also be folded in using different marginalization procedures.
In this work, we explore these different possibilities for quantifying the uncertainties of the high-density pQCD EoS at different chemical potentials µ using the MiHO code [66,67]. Additionally, we further develop the statistical framework to propagate the pQCD results to NS densities [33,44]. In particular, we discuss the marginalization of the chemical potential where the pQCD result is used. In this way we can combine the information of how constraining the pQCD results are with how convergent the pQCD series is at different chemical potentials. This extends the previous works that have so far considered only a fixed chemical potential for applying the pQCD results. The organization of the paper is outlined in the following Section.

II. OVERVIEW OF THE SETUP
In this section we present the overview of our setup of first estimating the perturbative uncertainties of the EoS at high densities using the Bayesian machine-learningbased framework MiHO [66,67] and then applying robust equation-of-state constraints [33,44] to translate the high-density information to NS densities.The different elements of the framework are introduced here while the details are discussed in the following sections.
Our goal is to determine the posterior probability of the (reduced) EoS p L (ϵ L ) at NS densities n L , given the first k + 1 terms of the perturbative series for the pressure 1 p (k) ≡ (p (0) , · · · , p (k) ) and number density n (k) ≡ (n (0) , . . . , n (k) ) at high densities where the perturbative description of QCD is reliable. Here, ϵ L , p L , and n L are the energy density, pressure, and baryon number density to which we wish to propagate the perturbative input; that is, at and around NS densities. The subscript L refers to low densities as opposed to the high densities, H, where the pQCD results converge and are directly valid. The first step is to convert the information about the perturbative series of the pressure p (k) and particlenumber density n (k) to a statistically interpretable probability distribution of pressure p H and density n H at a given large baryon number chemical potential µ H . To this end we use statistical machine-learning models that bound the higher-order terms {p (k+1) , p (k+2) , . . .} given the lowest k + 1 terms. We denote the resulting joint probability distribution quantifying the remaining MHO uncertainty as where X = 3Λ/(2µ) is related to the unphysical renormalization scaleΛ upon which the perturbative coefficients depend. We use the MiHO computer code to estimate these probabilities [67]. The current state-ofthe-art perturbative coefficients p (k) are reported in Section III and the details of the probability distribution of Eq. (2) is discussed in Section IV A. In order to remove the dependence on the unphysical renormalization scale X, one is forced to integrate over a sufficiently wide range of scales where P sa/sm (X|p (k) (µ H ), n (k) (µ H )) is an integration weight to be determined by a specific prescription. In this context we will discuss the scale-marginalization (sm) and scale-averaging (sa) prescriptions introduced in [65] and [66], respectively; for details, see Section IV B. This procedure is applied to the known non-perturbative results at high temperatures and zero chemical potential in Section IV C.
Given the EoS at high chemical potentials µ H -that is, the triplet of values p H , n H , and µ H -we can determine the region of allowed values ϵ L , p L at some lower density n L using the robust equation-of-state constraints introduced in [44]. These robust constraints are based on considering the most extreme extrapolations of the pQCD EoS that are allowed by mechanical stability, causality, and thermodynamic consistency to lower densities relevant to neutron stars. They can be expressed as a conditional probability 2 which takes a value 1 if the low density EoS at (n L , ϵ L , p L ) can be connected to the high density one at (µ H , n H , p H ) with any stable, causal, and consistent (SCC) extension and 0 otherwise. For details see Section V A. Combined with the joint probability distribution for the high-density EoS in Eq. (2), we can write down the distribution for ϵ L and p L at some given density n L and a matching scale µ H 2 We note that we use different thermodynamic variables at high and low densities. Namely, at low densities we characterized the EoS based on density n L whereas we use the chemical potential µ H at high densities. The reason for this choice is that the chemical potential is the natural variable in the pQCD calculation whereas the density is more directly relevant for the physics of NSs.
Finally, since the matching scale µ H is in principle arbitrary (as long as perturbative calculations remain reliable), we may marginalize over possible values of it similarly to what is done for the renormalization scale; the details are given in Section V B. The complete formula is then given by While optimally we would use the above expression to estimate the final posterior probability distribution, for practical reasons we have to make certain well-justified simplifications to it. Currently the MiHO computer code produces posterior distributions of one variable and is unable to compute joint probability distributions of two variables as required by Eq. (2). To estimate the joint probability distribution, we conservatively assume independent distributions for p H and n H such that where we have also used the fact that the series for n converges faster than the series for p in pQCD. Similarly, we estimate the integral weights arising from the scale marginalization procedure using only information about the perturbative series in pressure P sa/sm (µ H , X|p (k) , n (k) ) ≈ P sa/sm (µ H , X|p (k) ). (8) Both of these approximations are well justified given that the perturbative series for the number density n (k) converges faster than the corresponding series for the pressure p (k) . We leave improvements to the MiHO computer code and associated models discussed below in Section IV A to future work.
With these approximations the final result can be written as In the next three sections, we discuss the three types of input used in the master formula Eq. (9).

III. PQCD AT HIGH DENSITIES
In this section we discuss the current state-of-the-art pQCD results for the EoS at zero temperature T in β-equilibrium. We limit the discussion to 3 flavors of massless quarks, for which the β-equilibrium condition trivializes and the individual quark chemical potentials are given by µ u = µ d = µ s = µ/3 ≡ µ q , with µ being the baryon-number chemical potential. The approximation of zero quark masses is justified as we use the pQCD result only at very large chemical potentials, much larger than the strange quark mass µ q ≫ m s , such that the corrections from the strange quark masses are small [37,38,41]. At the same time, the chemical potential is below the charm threshold justifying restricting to 3 flavors.
In the following we consider only the results for the pressure because the pressure and density are linked to each other through the thermodynamic identity The perturbative expansion for the pressure reads where the different coefficients depend onΛ explicitly as well as through the strong coupling constant α s (Λ). The leading-order (LO) and next-to-leading-order (NLO) results are easily extracted as the T → 0 limit of results in [39] where C A = N c = 3, d A = N 2 c − 1 = 8 are the usual color factors and N f = 3 is the number of active fermion flavors.
The next-to-next-to-leading-order (N 2 LO) result consists of two types of terms, The first one is the hard contribution 3 , that arises from the 3-loop bubble diagrams. In these diagrams, the momenta are of order µ q . The second term in Eq. (14) is due to the soft divergence of QCD that is regulated by the in-medium screening of color charges, which leads to an enlarged set of diagrams with soft momenta ∼ α 1/2 s µ q . In this case, one must perform an allloop-orders resummation of a particular class of "ring diagrams" to correctly capture all of the physics from this soft momentum scale at this order in the coupling; for details, see [43]. The resummation is most conveniently performed using the hard-loop effective theory [68] leading to the result [39] p where m 2 E = 2αs(Λ) π N f µ 2 q is the in-medium screening mass scale.
Separately the two contributions in Eq. (14) are divergent at zero temperature, but the T dependence vanishes in the sum of the two contributions. Only the hard contribution has a physical UV divergence that is regulated by renormalization. Therefore, only the hard contribution explicitly depends on the unphysical renormalization scaleΛ. The dependence onΛ appears naturally in the logarithm together with the scale 2µ q . This factor of 2 can be traced to the coefficient of ϵ in performing Feynman integrals with integer powers of propagators in D = 4 − 2ϵ spacetime dimensions in theΛ MS -scheme. This is the motivation for choosing the central renormalization scaleΛ = 2µ q = 2 3 µ. The dimensionless quantity parameterizing the variation around the central scale is denoted by X, such thatΛ = 2µ q X, and the conventional scale variation by a factor of two is given by vary- At N 3 LO only the soft contribution to the pressure is known (this time at exactly zero T ) [42,43], We will not use this result unless explicitly stated, and when doing so we will denote this order with an asterisk -N 3 LO * -to remind the reader that not all contributions at this order are included.
In our results, we use the above expressions for the EoS, including the two-loop running of the strong coupling α s (Λ) in all results. We set the non-perturbative scale Λ QCD = 378 MeV, which fixes α s (2 GeV) = 0.2994. The renormalization-scale dependence of different perturbative orders of pressure and number density is depicted in Fig. 1 at different values of the baryon chemical potential. The fiducial range of scale variation is shown as vertical lines. We note that for all but the N 3 LO * results, there exists a smallest value of X for which the pressure remains positive. As µ is decreased, this smallest value of X increases, eventually approaching or moving into the fiducial range. We note that this diverging behavior at small X is the origin of sizeable scale variation errors in the conventional practice. The number density has a much smaller dependence on X and remains positive in the fiducial range even at µ = 2.2 GeV.

IV. ESTIMATING MHO AND RENORMALIZATION-SCALE UNCERTAINTIES
In this section, we discuss our handling of the uncertainties of the high-density QCD EoS. On the one hand, the truncated perturbative expansion Eq. (11) differs from the actual value by MHO terms. On the other hand, the truncated series depends on the choice of the renormalization scale, which is controlled by the parameter X. We use the recently introduced MiHO framework to estimate both uncertainties using Bayesian machinelearning techniques [66]. The main idea of the framework is to assume that the perturbative coefficients can be taken as independent draws from distributions arising from a statistical model of convergent series. Performing Bayesian inference on the available perturbative orders allows one to constrain the model parameters and construct the probability distribution for the next term in the series. The spread of the posterior distribution can then be used to quantify the uncertainty of the MHO terms. The scale uncertainty is incorporated by combining the probability distributions at different scales using a particular scale prescription, which will be specified below.

A. Posterior distribution at fixed scale
As explained in Section II [see Eq. (7)], we desire the probability distribution of the pressure and number den-sity given the first k + 1 terms in the perturbative expansion p (k) = (p (0) (µ, X), . . . p (k) (µ, X)) for fixed µ and X P MHO (p|p (k) (µ, X)).
For a convergent series the probability for the sum can be approximated by the posterior probability for the next partial sum in the series, which we will infer using two models: the geometric and abc models.
In practice we work with perturbative corrections normalized to the LO term. That is, we define a sequence of coefficients and δ k ≡ (δ 0 , . . . , δ k−1 , δ k ) with δ 0 = 1 by definition. The Bayesian model consists of a parameterized prior distribution for each δ k .
The geometric model assumes a flat distribution for with hidden model parameters a and c [65]. The normalized prior distribution for the δ k is That is, the sequence of perturbative coefficients is considered to be random variables, drawn from a uniform distribution δ k ∼ U [−ca k ,ca k ] whose width decreases geometrically with increasing order for 0 < a < 1.
The abc model introduced in [66] was proposed to allow for an asymmetric distribution with Note that in the abc model, the parameter a can take negative values ( −1 < a < 1), and therefore it can differentiate between alternating and non-alternating series. The model is then trained using the sequence of known terms, leading to posterior distributions for the model parameters a, (b), and c using Bayes's theorem P (a, c|δ k ) = P (δ k |a, c)P 0 (a)P 0 (c) where P (δ k |a, c) is the product of Eq. (21) [or Eq. (23)] for each term in δ k and P 0 (a) and P 0 (c) are judiciously chosen priors.
For the geometric model P 0 (a) and P 0 (c) are chosen to satisfy 0 < a < 1 and c ≥ 1: Here ϵ = 0.1 and ω = 1 are default constants defining the prior distributions. For the abc model the priors are chosen as The default constants for the abc model are (ϵ, ω, ξ, η) = (0.1, 1, 2, 0.1). 5 The marginalized likelihood (or evidence) for the geometric model is given by and similarly for the abc model. The marginalized distribution measures the (relative) confidence in the model's capabilities to reproduce the known terms in the series.
We use the trained model to find the posterior distribution of the next order in the series Then our desired posterior probability distribution for the pressure is approximated by [66] P MHO (p|p (k) (µ, X)) ≈ 1 In the top panels of Fig. 2 we show the posterior distributions for the pressure at fixed µ = 2.6 GeV and X = 1 for the geometric (left) and abc models (right). By construction, the LO distribution contains no useful information and is completely determined by priors. At NLO the distribution for the abc model is asymmetric, in contrast to the geometric model. Because the NLO correction is negative, the abc model infers an alternating series; hence, the distribution is skewed towards larger pressure values, i.e., it expects a positive N 2 LO correction. However, the actual N 2 LO correction is again negative. This causes the abc-model posteriors to become symmetric. In general both models predict similar posterior distributions with more input orders, as seen in the figure. In the bottom panels of Fig. 2 we show the 68% and 95% credible intervals (CIs) for the N 2 LO posterior distributions. 6 We see that the 68% CI fully incorporates the N 3 LO * corrections. The 95% CI for the abc model is noticeably wider than that of the geometric model. If the series were alternating as expected by the abc model, then the posterior CIs would be narrower than the ones of the geometric model.

B. Incorporating scale dependence
The truncated perturbative series depends on the unphysical renormalization scaleΛ = 2 3 µX, which is related to the characteristic hard scale of the process-the baryon number chemical potential µ-via the dimensionless scaling factor X. Consequently, the posterior probability of a Bayesian model is also sensitive to the choice of X. As there is no physical reason to prefer one scale choice over another, the scale dependence must be considered in estimating the observable's true value.
The scale-averaging prescription advocated in [66] treats all scales on equal footing and adds them coherently: Here the left-hand side should be understood as a scaleindependent probability distribution for the N (k+1) LO correction given terms up to N k LO. The weight function with F = 2, X 0 = 1 implements a log-uniform weight for the range 1/2 < X < 2. Note that the relative weight of the distributions does not depend on the convergence of the series. Alternatively, in the scale-marginalization prescription [65] the scale parameter X is treated as a hidden model parameter. The scale-independent probability is obtained by marginalization over X P sm (δ k+1 |δ k ) ≡ dXP (X|δ k (X))P (δ k+1 (X)|δ k (X)), (34) where the likelihood of X given δ k is P (X|δ k (X)) = P 0 (X)P (δ k (X)) dXP 0 (X)P (δ k (X)) .
Here P (δ k (X)) is the marginalized likelihood (evidence) given by Eq. (29) and P 0 (X) is now interpreted as a prior 6 Credible intervals are defined to contain the specified percentage of the probability distribution with the remaining probability split symmetry on either side of the CI.    for X. The marginalized likelihood is the largest for values of X for which the perturbative series is most convergent, i.e., favoring fastest apparent convergence [69,70].
In the case of pressure, the convergence is the fastest for the largest values of X, so that the weighting with the marginalized likelihood slightly favours higher values of X, see Fig. 2 for the marginalized likelihood as a function of X. In contrast, for the scale-averaging prescription, the integral Eq. (32) is dominated by the accumulation of probability distributions with slow X dependence, i.e., it accords with the principle of minimal sensitivity [70]. In Fig. 3 we display the scale-marginalized distributions for the pressure given in Eq. Notably, the perturbative series at low scale values does not affect the marginalized distribution. We also observe from this figure that the scale-marginalization prescription leads to distributions that are skewed to higher pressures than results of the naive scale-variation prescription, varying X ∈ [1/2, 2]. Scale-averaged distributions are also biased towards larger X, because there the scale dependence is slowing down. However, in this case convergence properties are not taken into account and the resulting distributions are wider.
In Fig. 4 we can see how the scale-marginalized distribution of the pressure (green line in Fig. 3) changes with the chemical potential µ. The green bands show the 68% and 95% CIs. The 68% CIs lies at the upper edge of the naive scale-variation estimate, shown in the figure as the hatched purple band. However the distributions have long power-law tails leading to 95% CIs with larger bands than scale variation, especially for large µ values.
C. Comparison with Lattice-QCD calculations at finite T and µ = 0 At high T and µ = 0, lattice computations of the EoS are available [71]. This enables the comparison between the precise nonperturbative value of the pressure p(T ) and the Bayesian inference from the perturbative calculations. In the high-T pQCD regime, due to the presence of thermally over-occupied, long-wavelength gluonic modes, the weak-coupling expansion for p is structured as a power series in g ≡ √ 4πα s , rather than in α s itself. Hence, in this context the terms p (k) denote terms of order O(g k ) = O(α k/2 s ), rather than O(α k s ). Note however, that the hard contributions (from momenta of order T ) are still structured as a power series in α s ∼ g 2 , since they arise from a finite number of individual bubble diagrams at the same naive order in the coupling. The resummed contributions on the other hand, arising from multiple Feynman diagrams featuring soft momenta ∼ α 1/2 s T appear at all orders starting at O(g 3 ). The p (k) are known FIG. 4. The pressure normalized by that of free Fermi gas of quarks as a function of chemical potential. The green bands correspond to the N 2 LO pQCD calculations whose uncertainties are estimated by the abc model using the scalemarginalization prescription for the renormalization scale X; the darker and the lighter bands represent 68% and 95%credible intervals, respectively. The relative confidence in the abc model is qualified by the marginalized likelihood P (µ) defined in Eq. (46) and illustrated by the black dashed line (and the fading of the green bands). The hatched purple band represents the standard error estimation of pQCD results obtained by renormalization scale variation by a factor of 2. Colored lines are the sample from the ensemble of NS EoSs used in [33] conditioned with astrophysical observations and QCD input for the scale-marginalization prescription for X in the range [1/2, 2] and µQCD in the range [2.2, 3] GeV. The coloring of individual EoSs corresponds to the posterior likelihood. The higher likelihood is associated with darker shades of red. up to k = 5, and can once again be extracted from [39] (see Appendix A).
As is well known [72] the perturbative expansion for the pressure at high-T is not well behaved without additional non-trivial resummations [73][74][75][76]. This can be readily seen from Fig. 5, where we show in different linestyles the different perturbative orders in the expansion. In particular, the result alternates above and below the free Stefan-Boltzmann (SB) value. In this case, due to the non-convergence of the series, the training of the statistical models is sensitive only to the g 5 term (i.e., if the model with given set of parameters can reproduce this term, it by construction can reproduce all the lower orders as well). Nevertheless, we show the CIs extracted with abc model with scale-marginalization (left panel) and scale-averaging prescriptions (right panel). Scalemarginalization leads to narrower CIs and only the 95% CI covers the non-perturbative results. Scale-averaging weights equally small X values which go downwards at small temperatures.
FIG. 5. Applying the statistically interpretable error estimation to QCD thermodynamics at zero chemical potential. In this case the pressure can also be accessed non-perturbatively using lattice-field-theory methods [71]. The green bands correspond to the 68% and 95%-CIs predicted by the abc model using (left) the scale-marginalization and (right) scale-averaging prescription for X in the range [1/2, 2]. The HTL-resummed perturbative result from [39] is shown with lines denoted by their order in the coupling constant g. The confidence in the model at various temperatures T is characterized by the marginalized likelihood P (T ) (displayed in arbitrary units).

V. PROPAGATING PQCD TO LOW DENSITIES
We turn now to a discussion of how pQCD results at high densities constrain the EoS at densities that are relevant for NSs. We use the recently introduced prescription from [44] to determine the allowed region for the low-density EoS given pQCD predictions at some fixed large value of µ H . The choice of this reference scale µ H is arbitrary as long as pQCD results are reliable. We will incorporate the uncertainty arising from the value of this scale into our predictions using a similar scale prescription as was done for the renormalization scaleΛ.

A. Constraints at low densities
Knowledge of the EoS at high densities n H and baryon number chemical potentials µ = µ H where the pQCD calculation is convergent imposes robust constraints to lower densities n L reached in cores in NSs [33,44]. The most conservative way of using the high-density information is to demand only that the EoS at lower densities can be connected to the EoS at µ H using an interpolation that is mechanically stable, causal, and thermodynamically consistent.
Given the triplet of values at high densities we can impose a condition on the corresponding triplet at lower densities ⃗ β L = p L , n L , µ L . Specifically, if the difference in pressure ∆p ≡ p H − p L does not lie in the with c 2 s,lim = 1, then there is no mechanically stable, causal and thermodynamically consistent way to connect the EoS at low densities to the EoS at high densities. Therefore any EoS with ∆p / ∈ [∆p min , ∆p max ] is excluded by pQCD.
These conditions restrict the area of allowed values the pressure p L and the energy density ϵ L can take at the given low density n L . In particular, to be able to match to the triplet β H , the energy density and pressure at the lower density n L , must be bounded by the curves p min (ϵ) < p < p max (ϵ) with ϵ being bounded by the intersections of the p min and p max curves denoted as ϵ min and ϵ max ϵ min ≡µ H n L − p H , See Fig. 6 for illustration. In this figure, we also show the corresponding region of p L and ϵ L values consistent with chiral effective field theory (cEFT) at low densities [77] (see Appendix B for the expressions for this corresponding boundary). Note that negative values of pressure are allowed by pQCD since it is possible to draw casual and stable EoS from these point to high-density limit (they might be, however, excluded for different reasons). As the information from a fixed high-density triplet β H is propagated down to a lower density n L , it gets spread over a wider and wider area in the ϵ, p plane. Specifically, the information becomes spread over the area bounded by the allowed curves in Eqs. (39) and (40), namely Now, the probability that a valid point (ϵ, p) lies within the bounded region at n L is unity, so that dϵ dp d 2 P (ϵ, p| ⃗ β H , n L ) dϵ dp , where d 2 P/(dϵ dp) is the differential probability of finding p and ϵ in the element of area. Taking the agnostic approach and not assigning higher likelihoods to less extreme EoSs, we may consider all allowed points equally likely, i.e. d 2 P (ϵ, p| ⃗ β H , n L ) dϵ dp where the constant is determined by the normalization condition of probability. Therefore the conditional probability is where 1 S is the indicator function on the set S. Note that for fixed µ H and n L the X dependence of the area is mild, since the dependence n H (X) is mild and the more X-sensitive p H (X) does not appear in the expression for the area.

B. Marginalization over µH
In [33] it was argued that one should use the perturbative information at the smallest possible density where MHO uncertainties are under control. However, the choice of µ H is arbitrary as long as perturbative results are reliable. Hence, here we shall incorporate the µ H dependence using the more sophisticated scalemarginalization prescription discussed above. Note that there are alternative ways to handle this µ H dependence such as the scale-averaging prescription.
In Eq. (44), we introduced the conditional probability for the SCC construction in [44] weighted with the inverse area A( ⃗ β H , n L ). The dependence of this area on µ H approximately scales as which follows from n (0) H ∼ µ 3 H . Therefore, one can see that applying the QCD input at larger values of µ H entails spreading the information over a larger allowed region of ϵ-p values. The area weight hence implies using smaller values of µ H in order to increase the constraining power. The scale-marginalization prescription introduces the opposite behavior (see the marginalized likelihood P (µ) in Fig. 4), i.e., increasing the value of µ H leads to more confidence in the perturbative calculation. This interplay allows a data-driven determination (based on p (k) ) of the relevant range of chemical potential that is at the same time as constraining as possible while being perturbatively reliable.
Note that we use a common marginalized likelihood P (δ k (X, µ H )) for both X and µ H . This is justifiable because in the procedure dictated by [44] one compares predictions at the same low-density point derived from different series labeled by µ H and X. Therefore Eq. (45) assigns weight according to the convergence properties, which depend both on µ H and X. We can also define the marginalized likelihood over X for µ H P (µ H ) ≡ dXP 0 (X)P (δ k (X, µ H )), (46) which is shown in Fig. 4 as the dashed black line. The interplay between scale-marginalization and the area weight can be seen in Fig. 7. The panel of 9 likelihoods displays P (ϵ L , p L |n L , p (k) ) calculated from Eq. (9) at fixed n L = 10n s using different prescriptions for the X and µ H variations. The range 1/2 ≤ X ≤ 2 is fixed for all likelihoods. Note that in these figures we have assumed also that the EoS can be causally connected to the low-density nuclear EoS [77] in a mechanically stable, causal, and thermodynamically consistent way (i.e., we have restricted to the purple shaded region in Fig. 6); this is purely to highlight the relevant region of the QCD likelihood function.
The first row reproduces the likelihood function used in [33], where a scale-averaging procedure for X was adopted, but µ H was kept fixed. The MHO uncertainty was not explicitly considered and only the highest available perturbative order (N3LO*) was used, corresponding to replacing the MHO distribution in Eq. (2) with δ-functions Going beyond [33] where only µ H = 2.6 GeV was used, we display in the first row results for three different chemical potentials: µ H = 2.2, 2.6, and 3 GeV corresponding to n H ≈ 23n s , 40n s , and 63n s . The smooth boundary of the likelihood function arises from the scale variation of the p max and p min lines as a function of X. For lower µ H values the allowed area is smaller as expected because P SCC is more constraining; however, the boundaries are less sharp because the X-variation of the pressure is stronger. The effect of the MHO and the scale-variation uncertainties for fixed µ H is demonstrated in the second row. Here, the MHO uncertainty is estimated using the abc model (using terms up to and including N 2 LO) and we have marginalized over X. We see (consistent with speculation in [33]) that the MHO and scale-variation uncertainties have only a minor effect to the final likelihood function; the main effect is to somewhat further blur the p max and p min boundaries compared to the first row.
Lastly, the third row constitutes a main result of this paper, in which the final QCD likelihood function using the master formula Eq. (9) is displayed. Compared to the second row, it includes the simultaneous marginalization over X and µ H from Eq. (45) with different µ H ranges. We observe that the interplay between the area weight and the marginalized likelihood makes the more aggressive range with lower values of µ H ∈ [2.2, 2.6] GeV (lower left panel) rather consistent with the more agnostic range µ H ∈ [2.2, 3] GeV (lower center panel). In particular it demonstrates the insensitivity to the upper limit of the µ H range. The most conservative range, µ H ∈ [2.6, 3] GeV (right lower panel), shows a somewhat larger allowed area in the ϵ-p plane; we note that this most conservative choice produces a likelihood rather close to the likelihood function used in [33], given by the µ H = 2.6 GeV panel in the top row.

VI. DISCUSSION AND CONCLUSIONS
In this paper, we have systematically discussed the effect of uncertainties of the perturbative-QCD input on the theoretically allowed ϵ-p region at neutron-star densities (n L = 10n s ) arising from the missing higherorder terms, the renormalization-scale variation, and the choice of reference density. To this end we employed Bayesian-inference techniques developed in high-energy physics [62,65,66]. The main results of this inference are displayed in the bottom row of Fig. 7, from which we conclude that pQCD constraints are robust with respect to different sources of uncertainties.
It remains a pertinent question how much the improved pQCD uncertainty estimation affects the global equationof-state (EoS) inference of neutron-star matter. To provide the first insights on this topic we interface the QCD likelihood function obtained here with the EoS inference setup of [33]. All the details of the EoS inference strictly follow those of [33] with the exception of the QCD likelihood function used. In Fig. 8 we show the effect of the pQCD input on the posterior EoS region that has already been conditioned using a set of astrophysical and cEFT constraints, shown in pink. The green hatched band shows the effect for a fixed value of µ H = 2.6 GeV with log-uniform weights for X ∈ [1/2, 2], which was denoted as "Pulsars+Λ+QCD" in [33]. The red dot-dashed band is the corresponding result when QCD is imposed at a larger chemical potential of µ H = 3 GeV, which is the most conservative choice discussed here. The blue dashed line corresponds to the case of scale marginalizing µ H in the range from 2.2 to 3 GeV and X in the range from 1/2 to 2 obtained using the abc model, which is the most comprehensive and agnostic uncertainty estimation we consider. We see that in all three of these cases, the pQCD input softens the EoS at high densities. Marginalizing over µ H only slightly widens the band at the high energy densities. Even though we allow µ H as high as 3.0 GeV (or n ≈ 64n s ) the band is only mildly affected.
In Fig. 4 in the introduction we have already shown EoSs in the form of pressure as a function of chemical potential, where the likelihood has been computed using the prescription of the blue dashed lines in Fig. 8 (with the same astrophysical constraints as above). We see from Fig. 4 that the QCD EoS constrained by astrophysical observations at low densities can be smoothly connected to our posterior distribution of the pQCD pressure at higher densities. Remarkably, even at the max-FIG. 7. The panel of the likelihoods P (ϵL, pL|nL, p (k) ) within the allowed region of Fig. 6 at fixed n = 10ns using different prescriptions. The first row corresponds to the fixed µH = 2.2, 2.6 and 3 GeV and the log-uniform weight for X, and it assumes pH and nH are given by the N 3 LO * pQCD values. The second row is obtained using the abc model with scale-marginalization prescription for X and fixed µH . The third row introduces the simultaneous marginalization over X and µH in the ranges [2.2, 2.6], [2.2, 3], [2.6, 3] GeV. The solid green (the dashed red) lines are the same as in Fig. 6 indicating the allowed region if EoS is connected to only pQCD (cEFT).
imal densities reached within stable NSs, the EoS is in rough agreement with the 68% CIs from the posterior distribution.
We conclude that the effect of the pQCD input on the NS-EoS inference is insensitive to the different prescriptions and choices used to estimate the uncertainties in the pQCD calculation. This increases the confidence in their use within neutron-star EoS inference. with ϵ being bounded by the intersections of the p χ min and p χ max curves denoted as ϵ χ min and ϵ χ max ϵ χ min ≡µ χ n L − p χ , The area bounded by these curves is given by (B5)