# Equation of state at neutron-star densities and beyond from perturbative QCD

Oleg Komoltsev<EMAIL_ADDRESS>Faculty of Science and Technology,
University of Stavanger, 4036 Stavanger, Norway  Rahul Somasundaram
<EMAIL_ADDRESS>Theoretical Division, Los Alamos National Laboratory, Los
Alamos, New Mexico 87545, USA  Department of Physics, Syracuse University,
Syracuse, NY 13244, USA  Tyler Gorda<EMAIL_ADDRESS>Institut für
Theoretische Physik, Goethe Universität, Max-von-Laue-Str. 1, 60438 Frankfurt
am Main, Germany  Technische Universität Darmstadt, Department of Physics,
64289 Darmstadt, Germany  ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum
für Schwerionenforschung GmbH, 64291 Darmstadt, Germany  Aleksi Kurkela
<EMAIL_ADDRESS>Faculty of Science and Technology, University of
Stavanger, 4036 Stavanger, Norway  Jérôme Margueron<EMAIL_ADDRESS>Institut de Physique des 2 infinis de Lyon, CNRS/IN2P3, Université de Lyon,
Université Claude Bernard Lyon 1, F-69622 Villeurbanne Cedex, France
International Research Laboratory on Nuclear Physics and Astrophysics,
Michigan State University and CNRS, East Lansing, MI 48824, USA  Ingo Tews
<EMAIL_ADDRESS>Theoretical Division, Los Alamos National Laboratory, Los
Alamos, New Mexico 87545, USA

###### Abstract

We explore the consequences of imposing robust thermodynamic constraints
arising from perturbative Quantum Chromodynamics (QCD) when inferring the
dense-matter equation-of-state (EOS). We find that the termination density, up
to which the EOS modeling is performed in an inference setup, strongly affects
the constraining power of the QCD input. This sensitivity in the constraining
power arises from EOSs that have a specific form, with drastic softening
immediately above the termination density followed by a strong stiffening. We
also perform explicit modeling of the EOS down from perturbative-QCD densities
to construct a new QCD likelihood function that incorporates additional
perturbative-QCD calculations of the sound speed and is insensitive to the
termination density, which we make publicly available.

††preprint: LA-UR-23-32625

##  I Introduction

Studying the impact of perturbative quantum chromodynamics (pQCD) calculations
on the inference of the equation of state (EOS) of neutron stars (NSs) has
recently become a promising area of research. The fundamental theory of the
strong interaction, quantum chromodynamics (QCD) is non-perturbative at low
energies but becomes perturbative at asymptotically high energies, which are
reached either at high temperatures or at high densities, at about
$20-40n_{\mathrm{sat}}$ in baryon number density Kurkela  _et al._ (2010);
Kurkela and Vuorinen (2016); Gorda  _et al._ (2018, 2021, 2023a). Here,
$n_{\mathrm{sat}}\approx 0.16$ fm-3 is the nuclear saturation density.

The cores of stable NSs explore densities that are most likely limited to
about $5-8n_{\mathrm{sat}}$, depending on the EOS which is realized in their
cores Schaffner-Bielich (2020). This implies that the theoretical description
of dense matter in stable NSs requires the solution of QCD in a non-
perturbative regime. At low densities, up to $1-2n_{\rm sat}$, effective-
field-theory (EFT) approaches, such as chiral EFT Epelbaum  _et al._ (2009);
Machleidt and Entem (2011), are commonly employed to study the NS EOS Hebeler
_et al._ (2013); Tews  _et al._ (2018a); Drischler  _et al._ (2019). Such
calculations are typically extended to higher densities using model-agnostic
extrapolations Read  _et al._ (2009); Lindblom (2010); Lindblom and Indik
(2012); Hebeler  _et al._ (2013); Lindblom and Indik (2014); Kurkela  _et al._
(2014); Annala  _et al._ (2018); Most  _et al._ (2018); Tews  _et al._
(2018b); Landry and Essick (2019); Greif  _et al._ (2019); Annala  _et al._
(2020); Raaijmakers  _et al._ (2020); Essick  _et al._ (2020); Annala  _et
al._ (2022); Altiparmak  _et al._ (2022); Somasundaram  _et al._ (2023a);
Jiang  _et al._ (2023); Gorda  _et al._ (2023b); Annala  _et al._ (2023),
conditioned with various astrophysical inputs. Some of these analyses have
also included input from the fundamental QCD Lagrangian in the perturbative
regime (henceforth referred to as “QCD input”) in the form of pQCD
calculations at high densities Kurkela  _et al._ (2014); Annala  _et al._
(2018); Most  _et al._ (2018); Annala  _et al._ (2020); Somasundaram  _et al._
(2023a); Annala  _et al._ (2022); Han  _et al._ (2023); Jiang  _et al._
(2023); Gorda  _et al._ (2023b); Annala  _et al._ (2023). These works have
reported results that differ from those that do not account for this QCD
input. In particular, the works including the QCD input have reported
softening of the EOS at energy densities above $750$ MeV/fm3, for which there
seems to be no evidence without this input.

The connection between the perturbative regime of QCD and NSs was first
investigated in Ref. Kurkela  _et al._ (2014), by extending model-agnostic
extrapolation functions of the EOS up to the densities where pQCD calculations
become reliable and imposing the pQCD values as boundary conditions for these
extrapolations. More recently, Komoltsev and Kurkela Komoltsev and Kurkela
(2022) have suggested a new method to link both density regimes, allowing them
to ‘integrate backwards’, i.e., to propagate the pQCD constraints to lower
densities in a completely general, analytical, and model-independent manner
using only the thermodynamic potential and the conditions of causality and
mechanical stability. In practice, the method provides a necessary condition
for an extrapolated EOS that has to be fulfilled at all densities below the
validity range of pQCD calculations. This condition asserts that all the
points of the EOS can be connected to the pQCD values by some causal and
stable interpolation. The QCD input was shown to constrain the EOS down to
densities of $\sim 2.2-2.5n_{\mathrm{sat}}$, if otherwise only causality and
mechanical stability are considered beyond chiral EFT. Within this
construction, the only uncertainty in the constraints arises from the
truncation error of the pQCD computation of the EOS at perturbative densities
(for a comprehensive study on the pQCD truncation errors, see Ref. Gorda  _et
al._ (2023c)). However, at these densities, astrophysical observations of NSs
are naturally more informative.

While the construction suggested in Ref. Komoltsev and Kurkela (2022) fully
bypasses the need for modeling NSs (as the input to these constraints derives
from high-energy particle physics, not from NSs), in order to interface the
QCD input with these astrophysical observations several assumptions must be
made. These include extrapolations involved in nuclear theory, assumptions
involved when extracting NS properties from observational data, the prior
assumptions of how EOS models are constructed, and up to what terminating
number density $n_{\mathrm{term}}$ such modeling is performed. These
assumptions and their interplay with the pQCD calculations are an important
subject in discussions about the constraining power of astrophysical and QCD
inputs.

The constraining power of pQCD calculations was studied in detail in two
recent publications Gorda  _et al._ (2023b); Somasundaram  _et al._ (2023b),
and was later also investigated by Refs. Ecker and Rezzolla (2022); Annala
_et al._ (2023); Essick  _et al._ (2023); Brandes  _et al._ (2023); Musolino
_et al._ (2023); Mroczek  _et al._ (2023); Fan  _et al._ (2023); Tang  _et
al._ (2023). These works have employed astrophysical input in the form of
gravitational-wave data constraining the tidal deformability of NSs in binary
mergers Abbott  _et al._ (2017, 2018, 2019), radio measurements of NS masses
Antoniadis  _et al._ (2013); Cromartie  _et al._ (2019); Fonseca  _et al._
(2021), as well as simultaneous mass-radius measurements using X-ray
observations Steiner  _et al._ (2018); Nättilä  _et al._ (2017); Shaw  _et
al._ (2018); Miller  _et al._ (2019a); Riley  _et al._ (2019); Miller  _et
al._ (2021); Riley  _et al._ (2021). They discussed whether the QCD input
reduces the uncertainties of the inferred EOS when imposed on top of
astrophysical inputs or whether the QCD and astrophysical inputs disfavor
similar behaviors. In the latter case, the two inputs would corroborate one
another, strengthening the conclusion that would have been obtained using only
one.

In this work, we compare the analyses conducted in Refs. Gorda  _et al._
(2023b) and Somasundaram  _et al._ (2023b). These studies arrived at different
conclusions but used the same pQCD calculations and the same way of
constraining the EOS Komoltsev and Kurkela (2022), see Fig. 1. While both of
these studies agreed that the QCD input impacted the EOS at the densities
reached in massive NSs, Ref. Gorda  _et al._ (2023b) found that QCD input
provides considerable additional information on top of other sources of data
whereas Ref. Somasundaram  _et al._ (2023b) found this additional impact to be
marginal. The conclusions of these two studies seem to be in contradiction,
and the aim of this work is to investigate the origin of this apparent
discrepancy.

Figure 1: The construction of the pQCD constraint of Ref. Komoltsev and
Kurkela (2022). The explicit form of the EOS is specified until
$n=n_{\mathrm{term}}$. At this point, the pQCD constraint considers all
possibilities for connecting the EOS from $n_{\mathrm{term}}$ to
$n_{\mathrm{pQCD}}$, where pQCD calculations become reliable. The blue and red
solid lines show EOSs that are compatible with pQCD at $n_{\mathrm{term}}$,
having a pQCD likelihood function of 1. However, for the blue EOS, the
parameter space for the EOS at higher densities is large (indicated by the
blue-dashed lines) while the EOS above $n_{\mathrm{term}}$ are extremely
constrained for the red EOS (indicated by red-dashed lines). Finally, the
black line corresponds to an EOS which cannot be connected to pQCD at
$n_{\mathrm{term}}$ without violating causality or thermodynamic stability.
For this EOS, the pQCD likelihood function is 0.

In addition to this direct comparison between Refs. Gorda  _et al._ (2023b);
Somasundaram  _et al._ (2023b), we also investigate the EOS behavior beyond
$n_{\mathrm{term}}$ that is required in order for the EOS to remain consistent
with the QCD input to higher densities. We do so by generating ensembles of
extensions of a given EOS beyond $n_{\mathrm{term}}$, which we condition with
the QCD input at $15n_{\rm sat}$. In this context, we also introduce a concept
of a pQCD tension index, which quantifies the degree to which a given EOS is
in tension with the high-density pQCD constraint within the construction of
Ref. Komoltsev and Kurkela (2022), before being fully excluded.

Finally, we also introduce a method for modifying the QCD input itself by
marginalizing over a prior model of high density EOSs constructed by
extrapolating down from $n<n_{\mathrm{pQCD}}$. Such marginalization has so far
only been performed within parametric models Han  _et al._ (2023); Jiang  _et
al._ (2023); Annala  _et al._ (2023) (cf. the discussion in Ref. Essick  _et
al._ (2023)). Importantly, such a construction also permits us to condition
the model with further pQCD information not used in the original construction.
In particular, the speed of sound exhibits very good convergence down the
densities as low as $\approx 25n_{\mathrm{sat}}$. By conditioning the high-
density prior with the speed of sound gives us a novel, more theory-informed
QCD input at the price of increased model dependence.

The paper is organized as follows. In Sec. II, we review the construction of
Ref. Komoltsev and Kurkela (2022) and define the QCD input. We also introduce
the pQCD tension index. In Sec. III, we discuss our computational setup,
including the prior EOS models used below $n_{\mathrm{term}}$ and the
astrophysical constraints that we use within our Bayesian analysis. In Sec.
IV, we discuss our results. Finally, in Sec. V, we discuss our findings and
limitations of our work.

##  II Review of the QCD input

Following Refs. Gorda  _et al._ (2023b) and Komoltsev and Kurkela (2022), we
implement the QCD input as follows. Given the EOS at two different points, at
a low density $\beta_{L}=(n_{L},p_{L},\mu_{L})$ and at a high density
$\beta_{H}=(n_{H},p_{H},\mu_{H})$, we may ask whether these two points can be
connected by _any_ stable and causal EOS. Here $n$ is the baryon density, $p$
is the pressure, and $\mu$ is the baryon chemical potential. This question can
be answered by exploring the bounds for EOS behavior between $\beta_{L}$ and
$\beta_{H}$.

The range of possible pressure differences between $(\mu_{L},n_{L})$ and
$(\mu_{H},n_{H})$ are limited by the following two EOSs in this interval

$\displaystyle n_{\mathrm{pmin}}(\mu)$
$\displaystyle=n_{L}\frac{\mu}{\mu_{L}},\quad\mu_{L}\leq\mu<\mu_{H}$ (1)
$\displaystyle n_{\mathrm{pmax}}(\mu)$
$\displaystyle=n_{H}\frac{\mu}{\mu_{H}},\quad\mu_{L}<\mu\leq\mu_{H}$ (2)

namely those which have one segment of $c_{s}^{2}=1$ followed by a phase
transition or vice versa, respectively. These EOSs bound the range of possible
pressure differences to be

$\displaystyle\Delta p_{\min}$
$\displaystyle=\frac{n_{L}}{2}\left(\frac{\mu_{H}^{2}}{\mu_{L}}-\mu_{L}\right),$
(3) $\displaystyle\Delta p_{\max}$
$\displaystyle=\frac{n_{H}}{2}\left(\mu_{H}-\frac{\mu_{L}^{2}}{\mu_{H}}\right).$
(4)

Given values for $\beta_{L}$ and $\beta_{H}$, the pressure difference between
these two points is fixed. Therefore the simple check (for detailed derivation
see Refs. Komoltsev and Kurkela (2022); Gorda  _et al._ (2023b))

$\displaystyle\Delta p=p_{H}-p_{L}$ $\displaystyle\in[\Delta p_{\min},\Delta
p_{\max}]$ (5)

determines whether the two points can be connected by a causal and stable EOS.
This condition does not make use of any specific form of the EOS between the
low- and high-density points but rather implies the existence of at least one
EOS connecting the two limits. The condition on the pressure difference can be
rewritten as

$\displaystyle 0\leq I_{\mathrm{pQCD}}\equiv\frac{\Delta p-\Delta
p_{\min}}{\Delta p_{\max}-\Delta p_{\min}}\leq 1\,,$ (6)

where we introduce a new quantity $I_{\mathrm{pQCD}}$, the pQCD tension index.
If the index $I_{\mathrm{pQCD}}\notin[0,1]$ the EOS is excluded because there
is no possible causal and stable interpolation. If the index obtains the
values $I_{\mathrm{pQCD}}=0\,(I_{\mathrm{pQCD}}=1)$, then the EOS is forced to
follow $n_{\mathrm{pmin}}\,(n_{\mathrm{pmax}})$, respectively. In particular,
if $I_{\mathrm{pQCD}}=1$, the EOS must contain a phase transition of size

$\Delta n=n_{H}\frac{\mu_{L}}{\mu_{H}}-\mu_{L}\,,.$ (7)

For typical values of $n_{L}=n_{\mathrm{term}}=n_{\mathrm{TOV}}$, we obtain
that the strength of phase transition reaches $\Delta n\sim
20n_{\mathrm{sat}}$, followed by a $c^{2}_{s}$ = 1 segment of $\sim
10n_{\mathrm{sat}}$. Note that if $I_{\mathrm{pQCD}}$ is not exactly
$I_{\mathrm{pQCD}}=0$ or $I_{\mathrm{pQCD}}=1$, then both extreme EOSs of Eqs.
(1) and (2) are excluded. For $0<I_{\mathrm{pQCD}}<1$, the EOS is bounded in
the $n$–$\mu$–$p$ -space by a non-trivial envelope whose projection in the
plane of energy-density $\varepsilon$ and pressure $p$ is shown in Fig. 1 (for
an explicit expression, see Ref. Komoltsev and Kurkela (2022)).

Following Ref. Gorda  _et al._ (2023b), we adopt the scale-averaging
interpretation of the pQCD uncertainties for $\beta_{H}$. We use a log-uniform
distribution for the dimensionless renormalization scale
$X=3\bar{\Lambda}/2\mu_{\mathrm{pQCD}}$ in the range [1/2, 2], which can be
written as $w(\log X)={\mathbf{1}}_{[\log(1/2),\,\log(2)]}(\log X)$. Here
$\bar{\Lambda}$ is the renormalization scale in the modified minimal
subtraction scheme Kurkela  _et al._ (2010). The QCD likelihood function then
reads

$\displaystyle P({\rm QCD|EOS})=\int$ $\displaystyle d(\log X)w(\log X)\times$
$\displaystyle\mathbb{1}_{[0,1]}\bigl{(}I_{\mathrm{pQCD}}(X,\rm EOS)\bigr{)}$
(8)

where $\mathbb{1}_{[0,1]}(I_{\mathrm{pQCD}}(X,\rm EOS))$ is an indicator
function discarding EOSs with $I_{\mathrm{pQCD}}\notin[0,1]$ for a fixed value
of $X$. In general, the index $I_{\mathrm{pQCD}}$ depends on all six
quantities $(n_{L},p_{L},\mu_{L})$ and $(n_{H},p_{H},\mu_{H})$, which for
fixed values of $\mu_{H}=\mu_{\rm QCD}$ reduces to
$I_{\mathrm{pQCD}}(X,\mathrm{EOS})=I_{\mathrm{pQCD}}(X,n_{L},p_{L},\mu_{L})$.

We note that the convergence of the perturbative series depends on the
unphysical renormalization scale. In particular, the convergence is slowest
for the smaller values of $X\approx 1/2$. This slow convergence renders the
$X=1/2$ results, which differ quantitatively from the higher values of $X$,
less trustworthy; see Ref. Gorda  _et al._ (2023c) for a detailed analysis of
the convergence of the pQCD results.

In the following sections we apply the QCD likelihood function at different
termination densities, meaning that $(n_{L},p_{L},\mu_{L})$ corresponds to
$(n_{\mathrm{term}},p_{\mathrm{term}},\mu_{\mathrm{term}})$. For the high
density input, $\beta_{H}=\beta_{\mathrm{pQCD}}$, we use the partial next-to-
next-to-next-to-leading order pQCD computation of Ref. Gorda  _et al._ (2021).

Note that in terms of the sound speed, the QCD input constrains the density
$n_{\mathrm{term}}$, which is given by a (logarithmic) integral of
$c_{s}^{-2}(\mu)$ and the pressure $p_{\mathrm{term}}$, which is the integral
of this integral. Hence for a given NS EOS, the input is more sensitive to
individual features in the speed of sound at lower densities and is less
sensitive to those near $n_{\mathrm{term}}$.

##  III Computational Setup

In order to study the dependence of our results on the differences in the
analyses between Refs. Gorda  _et al._ (2023b); Somasundaram  _et al._
(2023b), we use two different models to generate our EOS sets: the parametric
sound-speed model of Ref. Somasundaram  _et al._ (2023b) and the Gaussian-
process (GP) model of the speed of sound of Ref. Gorda  _et al._ (2023b). Note
that similar models have also been employed by Refs. Somasundaram  _et al._
(2023a); Tews  _et al._ (2018a); Landry and Essick (2019); Mroczek  _et al._
(2023).

###  III.1 Prior Set 1: The sound-speed model

To generate the CSM set, we follow Ref. Somasundaram  _et al._ (2023b). At the
lowest densities, we fix the EOS to the crust model of Ref. Douchin and
Haensel (2001) and do not explore uncertainties in the crust modeling. The
crust EOS is then matched with the EOS of the core via a cubic spline in the
$c^{2}_{s}$ versus $n$ plane. For the core EOS, below a certain breakdown
density $n_{b}$, we assume that the EOS can be accounted for by nucleonic
degrees of freedom, which in this work is described by the meta-model of Refs.
Margueron  _et al._ (2018a, b). The meta-model is a density-functional
approach, similar to the Skyrme approach, that allows one to incorporate
nuclear-physics knowledge directly encoded in terms of the Nuclear Empirical
Parameters (NEPs). We vary the NEPs in a rather broad range in order to
explore a wide range of nuclear density functionals that describe matter up to
$n_{b}$. The parameter $n_{b}$ is varied uniformly in the range
$1-2n_{\mathrm{sat}}$. Furthermore, we discard all models that do not satisfy
constraints provided by the chiral EFT calculations of Ref. Lynn  _et al._
(2016) up to nuclear saturation density.

Above $n_{b}$, we use the CSM and create a non-uniform grid in density between
$n_{b}$ and $25n_{\textrm{sat}}$ with $9$ grid points, where the 9th point is
fixed at $25n_{\textrm{sat}}$ and the other 8 points are drawn from a uniform
probability distribution in density. At each grid point, the squared sound-
speed is uniformly sampled between $0$ and $1$ in units of the speed of light.
Then, the sound speeds at each grid point are connected by linear segments,
generating a sound-speed profile $c^{2}_{s}(n)$ for the EOS. The density-
dependent speed of sound $c^{2}_{s}(n)$ can then be integrated to give the
EOS, i.e., the pressure $p(n)$, the energy density $\varepsilon(n)$, and the
baryon chemical potential $\mu(n)$, see Ref. Somasundaram  _et al._ (2023a)
for more details. The CSM set used here consists of $\approx 16,000$ EOSs
generated this way.

###  III.2 Prior Set 2: Gaussian-Process model

In the GP model, for densities below $n=0.57n_{\mathrm{sat}}$ we use the crust
model by Baym, Pethick, and Sutherland Baym  _et al._ (1971). Above the crust-
core transition density, we perform a GP regression in an auxiliary variable
$\phi(n)~{}\equiv~{}-\ln(1/c_{s}^{2}(n)-1)$, choosing the prior on this
variable from the multivariate Gaussian distribution ,

$\phi(n)\sim\mathcal{N}\bigl{(}-\ln(1/\bar{c}_{s}^{2}-1),K(n,n^{\prime})\bigr{)}\,,$
(9)

with a Gaussian kernel
$K(n,n^{\prime})=\eta\exp\bigl(-(n-n^{\prime})^{2}/2\ell^{2}\bigr{missing})$
Gorda  _et al._ (2023b). In our hierarchical model the hyperparameters of the
GP $\eta$, $\ell$, and $\bar{c}_{s}^{2}$ are themselves drawn from normal
distributions,

$\begin{split}\eta&\sim\mathcal{N}(1.25,0.2^{2})\,,\\\
\ell&\sim\mathcal{N}(1.0n_{\mathrm{sat}},(0.25n_{\mathrm{sat}})^{2})\,,\\\
\bar{c}_{s}^{2}&\sim\mathcal{N}(0.5,0.25^{2})\,.\end{split}$ (10)

For $n\leq 1.1n_{\mathrm{sat}}$, this prior is then conditioned on the chiral
EFT results from Ref. Hebeler  _et al._ (2013), with the average and
difference between the “soft” and “stiff” EOSs in that work taken to be the
mean and 90% credible interval for the conditioning Gorda  _et al._ (2023b).
The final EOSs are then once again obtained by integrating $c_{s}^{2}(n)$. The
GP used in this work consists of 120,000 EOSs generated in this manner.

Figure 2: (Top) The QCD likelihood function for $\mu_{\rm{QCD}}=2.4$ GeV (left
panels) and $\mu_{\rm{QCD}}=2.6$ GeV (right panels) for EOS posteriors
conditioned on astrophysical data. The EOSs are ordered according to their QCD
likelihood. The blue lines correspond to the GP, whereas the red lines
correspond to the CSM. The likelihood function is displayed for two different
termination densities, $n_{\mathrm{term}}=n_{\mathrm{TOV}}$ and
${n_{\mathrm{term}}=10n_{\mathrm{sat}}}$, labeled in the figure. (Middle)
Fraction of the astrophysical evidence (see Eq. 11) that is cut by the QCD
input as a function of the termination density. The dashed blue distribution
shows the posterior for the TOV densities for the GP whereas the dashed red
lines show the same for the CSM. (Bottom) Fraction of the astrophysical
evidence (see Eq. 11) that is cut by the QCD input as a function of the
termination density in units of the TOV density.

###  III.3 Astrophysical input

The astrophysical inputs that we consider in our analysis are the following:

  1. 1.

Radio-astronomy constraints: We consider the mass-measurements of PSR
J0348+0432 with $M=2.01\pm 0.04M_{\odot}$ Antoniadis  _et al._ (2013) and PSR
J1624$-$2230 with $M=1.928\pm 0.017M_{\odot}$ Fonseca  _et al._ (2016). We
approximate these uncertainties as normal distributions, which is valid to
good accuracy.

  2. 2.

Gravitational-wave data from binary NS mergers: We consider the two-
dimensional PDF for the tidal deformability and mass ratio at fixed chirp mass
extracted from GW170817 by the LIGO-Virgo collaboration Abbott  _et al._
(2019), using the PhenomPNRT waveform model with low-spin priors (see Figure
12 of Abbott  _et al._ (2019)).

  3. 3.

NICER observations: We consider the simultaneous mass and radius measurement
of PSR J0740$+$6620 using the NICER + XMM-Newton data. We have employed here
the PDF from Ref. Miller  _et al._ (2021) (see their Fig. 1, right panel).

After conditioning our priors with the astrophysical input, we obtain a
posterior sample of 4,000 EOSs for CSM and 46,000 EOSs for the GP model. These
inputs are identical to those used in Ref. Gorda  _et al._ (2023b) (referred
to as “Pulsars + $\tilde{\Lambda}$”) with one exception. We choose not to use
the additional “Occam” factor related to the prior mass distribution of NSs,
${(m_{\mathrm{TOV}}-m_{\min})^{-4}}$ (one power per mass integral), which
penalizes high-mass EOSs that are completely consistent with astrophysical
data; see Ref. Miller  _et al._ (2019b); Landry  _et al._ (2020) for a
discussion of this point. These inputs are the same as used in Ref.
Somasundaram  _et al._ (2023b) with the exception of the NICER measurement of
PSR J0030+0451. We have not included the latter in the present analysis as we
found its impact on estimating the pQCD constraining power to be small.

Ref. Gorda  _et al._ (2023b) analyzed the astrophysical as well as the QCD
inputs using a Bayesian approach whereas Ref. Somasundaram  _et al._ (2023b)
performed a hard cut analysis. We have checked that, as far as the treatment
of the astrophysical data is concerned, this difference does not alter the
conclusions of our paper regarding the impact of the QCD input. The role
played by the different treatments of the QCD input is discussed in the next
section.

##  IV Results

Next, we focus on the impact of the differences in methodology between Refs.
Gorda  _et al._ (2023b) and Somasundaram  _et al._ (2023b). The three major
differences between the setups in these works can be summarized as follows:

  1. D1:

Ref. Gorda  _et al._ (2023b) used a fully Bayesian approach and the
predictions are analyzed in terms of posterior probability distributions
whereas Ref. Somasundaram  _et al._ (2023b) performed a hard cut analysis and
focused on the envelope associated with all EOS models for all studied
quantities.

  2. D2:

The functional basis of the EOS modeling in the two works was different: Ref.
Gorda  _et al._ (2023b) employed the non-parametric GP model (Prior set 2),
while Ref. Somasundaram  _et al._ (2023b) employed the parametric CSM (Prior
set 1).

  3. D3:

The two works made different choices for $n_{\mathrm{term}}$: Ref. Gorda  _et
al._ (2023b) chose $n_{\mathrm{term}}=10n_{\mathrm{sat}}$ whereas Ref.
Somasundaram  _et al._ (2023b) chose
$n_{\mathrm{term}}=n_{\mathrm{TOV}}\approx 5-8n_{\mathrm{sat}}$, which is
different for each EOS. Hence, the works inferred the EOS in different density
ranges.

To investigate point D3, we will discuss results for varying
$n_{\mathrm{term}}$ expressed either in terms of multiples of the saturation
density $n_{\mathrm{sat}}$ or in terms of the EOS-dependent central density of
the maximum-mass neutron star, $n_{\mathrm{TOV}}$.

###  IV.1 Impact of the QCD input on the EOS

We show the QCD likelihood of Eq. (II) for the two EOS models, GP and CSM, in
the upper panels of Fig. 2, where both models terminate at
$n_{\mathrm{term}}=n_{\mathrm{TOV}}$. In particular, we show the sorted QCD
likelihood for a large sample of EOSs drawn from the posterior distributions
of the CSM and GP models after they have been conditioned on the astrophysical
data. EOSs for which the likelihood is $1$ are accepted for all values of the
pQCD renormalization scale $X$ whereas EOSs that have a likelihood of $0$
cannot be connected to pQCD with a causal and stable interpolation for any
value of $X\in[1/2,2]$. EOSs with likelihoods between the two limits can be
connected to pQCD for some but not all values of $X$, see Fig. 3 for more
details. Figure 2 clearly shows that the choice of EOS model has a weak impact
on the QCD likelihood, illustrating that the difference D2 plays a minor role
in the following analysis.

We observe that for $n_{\mathrm{term}}=n_{\mathrm{TOV}}$ and $\mu_{H}=2.6$ GeV
($\mu_{H}=2.4$ GeV) roughly 20% ($\sim 40\%$) of the EOSs cannot be connected
to pQCD for all values of $X$ and are therefore affected by the QCD input. For
a very small fraction of EOS, the QCD likelihood drops below $20\%$ but no EOS
has a vanishing QCD likelihood for this choice of $n_{\mathrm{term}}$. In
contrast, the results are very different for
$n_{\mathrm{term}}=10n_{\mathrm{sat}}$, which was chosen in Ref. Gorda  _et
al._ (2023b). Now, roughly $50-60\%$ of EOS have a vanishing QCD likelihood.
This shows the strong impact of the choice of $n_{\mathrm{term}}$ on the pQCD
constraining power and explains a large part of the difference in the
conclusions of Refs. Gorda  _et al._ (2023b); Somasundaram  _et al._ (2023b).
The strong differences between the two cases arise from the behavior the EOS
must adopt beyond $n_{\mathrm{TOV}}$ to remain consistent with the QCD input
at $10n_{\mathrm{sat}}$. This is discussed further in Section IV.3.

We quantify the impact of the QCD likelihood function as a function of
$n_{\mathrm{term}}$ by computing the fraction of the marginalized likelihood
(evidence) excluded by the QCD input

$\displaystyle 1-\frac{\sum_{i}w^{\rm astro}_{i}\cdot w^{\rm
QCD}_{i}}{\sum_{i}w^{\rm astro}_{i}}\,,$ (11)

where the index $i$ denotes the EOS, and the summation is performed over all
EOSs in the ensemble. The symbols $w^{\rm astro}$ and $w^{\rm QCD}$ correspond
to the likelihood arising from the astrophysical and QCD inputs, respectively.
In Ref. Pang  _et al._ (2023) this quantity was called the average consistency
(see also Eq. (25) of Ref. Gorda  _et al._ (2023b)), and corresponds to one
minus the average QCD likelihood in the posterior ensemble conditioned with
only the astrophysical input, which is indicated by the hatched area.

Figure 3:  The QCD likelihood function at $n_{\mathrm{term}}=n_{\mathrm{TOV}}$
for the GP EOS posteriors conditioned on astrophysical data. The EOSs are
ordered according to their QCD likelihood. The blue line corresponds to
$\mu_{\mathrm{pQCD}}=2.4$ GeV and a scale-averaged prescription for $X$ in the
range [1/2,2], whereas the red line corresponds to $\mu_{\mathrm{pQCD}}=2.6$
GeV. The colored ticks and dotted lines correspond to the likelihood function
for fixed values of $X$ ($X$=1/2, 1 and 2), taking the form of a step
function.  Figure 4: Posterior probability density of the energy density vs.
pressure with the QCD input enforced at three different fixed termination
densities. The average TOV density for the GP posterior after enforcing
astrophysical information is approximately 5.5$n_{\mathrm{sat}}$, while the
maximum TOV density reaches approximately 8$n_{\mathrm{sat}}$ for the softest
EOS. The areas correspond to 68% credible regions when considering only
astrophysical information (dotted gray) and when also including the QCD input
at $\mathrm{n_{term}}$ (solid blue). The QCD input is imposed at 2.6 GeV with
a scale-averaging prescription.

This quantity is shown as the function of the termination density for the two
EOS models in multiples of $n_{\mathrm{sat}}$ in the middle panels and
relative to $n_{\mathrm{TOV}}$ in the bottom panels of Fig. 2. As the models
are extended to higher $n_{\mathrm{term}}$, the fraction of the astrophysical
evidence cut by QCD increases rapidly; if the models are extended to
$n=1.2n_{\mathrm{TOV}}$ this fraction increases to 0.4 (0.5 for
$\mu_{\mathrm{QCD}}=2.4$ GeV). At a fixed density of
$n_{\mathrm{term}}=6.5n_{\mathrm{sat}}$, roughly corresponding to the
approximate 68%-credible upper limit for the $n_{\mathrm{TOV}}$, the average
consistency is 0.4. However, we stress that the stiffest EOSs, i.e., those
with the lowest QCD likelihoods, are in the lower percentile of the
$n_{\mathrm{TOV}}$ range, and hence, for these EOS
$n_{\mathrm{TOV}}<6.5n_{\mathrm{sat}}$. These EOS are more strongly disfavored
for larger $n_{\mathrm{term}}$, when the EOS branch beyond $n_{\mathrm{TOV}}$
is probed. For higher termination densities, the average consistency levels
out and takes values of $\sim 0.8$ at $n_{\mathrm{term}}=10n_{\mathrm{sat}}$.

We, hence, find that the constraining power of the QCD input depends largely
on the choice of the termination density $n_{\mathrm{term}}$ (D3). In fact, we
find that a value of $n_{\mathrm{term}}\approx 0.8n_{\mathrm{TOV}}$
($n_{\mathrm{term}}\approx 4n_{\mathrm{sat}}$ when applied at a fixed density)
constitutes a critical value, in the sense that it is at this point that the
QCD input begins to show any effect at all beyond the astrophysical input. We
note that the fact that this critical value of $n_{\mathrm{term}}$ is so close
to $n_{\mathrm{TOV}}$ is a coincidence, in the sense that this critical value
of $n_{\mathrm{term}}$ is set by a combination of many factors, including the
value of $n_{\mathrm{pQCD}}$ at which pQCD is well-converged and the
constraining power of the astrophysical input (the latter applying only below
$n_{\mathrm{TOV}}$). In particular, if $n_{\mathrm{pQCD}}$ were significantly
lower, then it is possible that this critical value of $n_{\mathrm{term}}$
would also be lower.

In order to investigate the impact of point D1 above, we show in Fig. 3 the
detailed impact of $X$ on the QCD likelihood. For different values of $X$, we
indicate where the likelihood function becomes 1, i.e., all EOS to the right
of these lines are accepted by the QCD input at
$n_{\mathrm{term}}=n_{\mathrm{TOV}}$, while those on the left side are ruled
out. Considering $\mu_{\mathrm{pQCD}}=2.6$ GeV as in Refs. Gorda  _et al._
(2023b); Somasundaram  _et al._ (2023b), we find that nearly all EOS are
accepted by $X=1/2$, recovering the results of Ref. Somasundaram  _et al._
(2023b). In contrast, roughly 25% of EOSs are rejected for $X=1$ and 60% for
$X=2$. If $\mu_{\mathrm{pQCD}}$ is decreased to $\mu_{\mathrm{pQCD}}=2.4$ GeV,
about 10% of EOS are rejected already at $X=1/2$. A Bayesian treatment of the
scale parameter $X$ as performed in Ref. Gorda  _et al._ (2023b) averages the
impact of $X=1/2$ leading to the findings of Ref. Gorda  _et al._ (2023b).
Hence we conclude that different statistical treatment of including pQCD
calculations (D1) plays a secondary role in the differing conclusions of the
two studies.

###  IV.2 Impact of the QCD input on the $\varepsilon$–$p$ values

Figure 5: Posterior probability density of the energy density vs. pressure at
$n_{\mathrm{TOV}}$ for different termination densities. Each contour
corresponds to 68% credible regions for the prior (black solid), considering
only astrophysical input (dotted gray) and when also including QCD input
(different colors corresponds to different termination densities). The left
plot shows results using the GP approach, while the right one shows results
for the CSM approach. The QCD input is imposed at 2.6 GeV with a scale-
averaging prescription.

We now analyze the effect of the QCD input on the allowed values of energy
density and pressure. Figure 4 shows posterior credible regions of energy
density $\varepsilon$ and pressure $p$ at various fixed densities in units of
saturation density, incorporating either the astrophysical input alone or
additionally also the QCD input. The EOS models are terminated at different
densities $n_{\mathrm{term}}\in\\{4.5,5.5,6.5\\}n_{\mathrm{sat}}$, which
correspond roughly to the range of $n_{\mathrm{TOV}}$ shown above in Fig. 2.
We note that taking the termination density to be the density we are studying
represents the maximally conservative approach to applying QCD constraint at a
fixed density, as the QCD input is completely independent of the prior beyond
$n_{\mathrm{term}}$. As discussed in Ref. Gorda  _et al._ (2023b), the main
effect of the QCD input is to disfavor the stiffest EOSs, visible in the
figure. The effect of the input depends on the density; consistent with Fig.
2, at densities $\sim 6.5n_{\mathrm{sat}}$ (corresponding to the 68%-credible
upper limit on $n_{\mathrm{TOV}}$) there is a significant reduction of allowed
$\varepsilon$–$p$ values, particularly at the highest $p$ and $\varepsilon$.
At lower densities the constraint is weaker as expected; at
5.5$n_{\mathrm{sat}}$, corresponding to the mean of $n_{\mathrm{TOV}}$, the
reduction is still significant, while taking
$n_{\mathrm{term}}=4.5n_{\mathrm{sat}}$ only shifts the posterior slightly at
the highest pressures. This is in line with the prediction of Ref. Komoltsev
and Kurkela (2022) stating that QCD input ceases to give any constraint at
$n<2.2-2.5n_{\mathrm{sat}}$. From Fig. 4, we can also clearly see a difference
between the two NS EOS models, namely that the GP model gives more weight to
EOSs with larger values of $p$ and $\varepsilon$ than the CSM, visible in the
posterior for the astrophysical input. This explains the slightly larger
effect of the QCD input within the GP for the same $n_{\mathrm{term}}$ (point
D2 above). However, after the QCD input is applied, the two models give quite
similar posteriors.

We now move from studying the EOS at fixed densities to varying the
termination density in terms of $n_{\mathrm{TOV}}$. In Fig. 5, we show the
posterior distribution of energy density $\varepsilon(n_{\mathrm{TOV}})$ and
pressure $p(n_{\mathrm{TOV}})$ for
$n_{\mathrm{term}}\in\\{1,1.2,1.4\\}n_{\mathrm{TOV}}$. The figure shows the
68%-credible regions for the prior distribution before introducing the
astrophysical constraints as well as the regions after the EOSs have been
conditioned on the astrophysical inputs but before the QCD input has been
imposed. For both EOS models, the upper (left) edge of the credible region
contour is nearly unchanged by the astrophysical input. This demonstrates that
there is a region in the
$\varepsilon(n_{\mathrm{TOV}})$–$p(n_{\mathrm{TOV}})$—plane that is only
marginally constrained by the present astrophysical input.

Similar to the $\varepsilon$–$p$—distribution at fixed density, the effect of
the QCD input is to soften the EOS, disfavoring the highest pressures. We
observe that for $n_{\mathrm{term}}\in\\{1.2,1.4\\}n_{\mathrm{TOV}}$ the QCD
input sets the upper left edge of the posterior distribution. This is also the
case for $n_{\mathrm{term}}=n_{\mathrm{TOV}}$ in the GP model, while for the
CSM at $n_{\mathrm{term}}=n_{\mathrm{TOV}}$ the impact is marginal. Consistent
with Fig 2, the quantitative impact is smallest for
$n_{\mathrm{term}}=n_{\mathrm{TOV}}$ but the impact grows rapidly when
$n_{\mathrm{term}}$ is increased and is qualitative similar to setting
$n_{\mathrm{term}}$ to fixed density $10n_{\mathrm{sat}}$, shown in Fig. 5.
From both of these two figures, we observe that the choice of EOS prior
(either GP or CSM) has only a minor impact on the effect of the QCD
constraint. Finally, we also show contours corresponding to a fixed
$n_{\mathrm{term}}\in\\{10,15,20\\}n_{\mathrm{sat}}$, where we show only CSM
contours for the latter two values. We note that applying the QCD input for
$n_{\mathrm{term}}\in\\{10,15,20\\}n_{\mathrm{sat}}$ yields very consistent
credible regions.

Figure 6:  (Left) The QCD constraint shown in Fig. 1 for 3 different EOS with
representative values of the pQCD tension index (see Eq. 6). Each EOS is
terminated at $n_{\mathrm{term}}=n_{\mathrm{TOV}}$ and accepted by QCD input
for a fixed $X$=1 and $\mu_{\rm QCD}=2.6$ GeV. The arrows illustrate the phase
transition and the subsequent $c_{s}^{2}=1$ segment, which an EOS is forced to
have when the tension index is close to 1. (Right) The distribution of the
pQCD tension index for EOSs that are accepted at $1.2n_{\mathrm{TOV}}$ and,
therefore, also at $n_{\mathrm{TOV}}$ (purple) as well as for EOS that are
accepted at $n_{\mathrm{TOV}}$, but rejected somewhere between
$n_{\mathrm{TOV}}$ and $1.2n_{\mathrm{TOV}}$ (blue). These distributions are
normalized to the number of EOSs in each distribution. Most of the EOSs from
the latter ensemble have a large pQCD tension index.  Figure 7: The behavior
an EOS with tension index $I_{\mathrm{pQCD}}=1$ at $n_{\rm TOV}$ (blue line)
has to follow above the termination density $n>n_{\mathrm{term}}$ to match
with pQCD for $X=1$ (black line). The purple band corresponds to pQCD
calculations of the speed of sound starting from $\mu_{\rm QCD}=2.2$ GeV, with
$X$ varying in the range [1/2,2].  Figure 8: Ensembles of possible extensions
for three representative EOSs that are allowed at
$n_{\mathrm{term}}=n_{\mathrm{TOV}}$. The EOSs correspond to different values
for the pQCD tension index $I_{\rm pQCD}(n_{\mathrm{TOV}})$. The extensions
are conditioned to fulfill the QCD constraint at 15$n_{s}$. In this figure,
the coloring of each EOS is given by the value of
$I_{\mathrm{pQCD}}(15n_{\mathrm{sat}})$. Figure 9: The distributions of the
averaged sound speed in the density interval $[8,15]n_{\rm sat}$ for the
possible EOS extensions of 100 different EOSs up to
$n_{\mathrm{term}}=n_{\mathrm{TOV}}$. For each EOS, we generate 1000
extensions. This distribution is shown for fixed values of the pQCD tension
index $I_{\mathrm{pQCD}}(n_{\mathrm{TOV}})=0.85,0.75,0.6,0.5$.

###  IV.3 Implications for the EOS beyond $n_{\mathrm{term}}$

So far, we have discussed the impact the QCD input has on the EOS for
$n<n_{\mathrm{term}}$. In this subsection, we discuss how the QCD constraint
affects the behavior of the EOS at densities $n>n_{\mathrm{term}}$. As shown
in Fig. 1, the QCD input affects the available parameter space for the EOS in
the $p$–$\varepsilon$ plane beyond $n_{\mathrm{term}}$,
$n_{\mathrm{term}}<n<n_{H}$. Some EOS have a large available parameter space
while others are tightly constrained. One possible measure to gauge this
freedom is the pQCD tension index $I_{\mathrm{pQCD}}$ defined before. The left
panel of Fig. 6 shows the allowed space in the $\varepsilon$–$p$ plane beyond
$n_{\mathrm{term}}$ for three EOSs with
$I_{\mathrm{pQCD}}(n_{\mathrm{TOV}})=0.75,0.85$, and $0.98$. The EOS with the
highest $I_{\mathrm{pQCD}}$ must undergo a rapid phase-transition-like
softening with $\Delta n\sim 20n_{\mathrm{sat}}$, followed by a segment of
$c_{s}^{2}\approx 1$ to remain consistent with pQCD calculations. For
$I_{\mathrm{pQCD}}=1$, this discontinuous behavior would be be completely
determined – see Fig. 7. Note that this phase transition is not responsible
for the onset of the unstable NS branch, which instead is determined by
general relativity within the density range modeled by the EOS. Instead, this
phase transition takes place right above $n>n_{\mathrm{TOV}}$. Therefore, the
onset of this phase transition coinciding with $n_{\mathrm{TOV}}$ is
accidental. If one wants to study the possibility of a phase transition
determining the onset of the unstable branch, phase transitions must be
included in the EOS models to a significant fraction, as in Ref. Gorda  _et
al._ (2023d). For $I_{\mathrm{pQCD}}=0.75$ a qualitatively similar softening
is required but there is more freedom in the EOSs for $n>n_{\mathrm{term}}$.
By contrast, an EOS with $I_{\mathrm{pQCD}}=0.5$ is largely unconstrained
above $n>n_{\mathrm{TOV}}$.

In order to visualize the behavior above $n_{\mathrm{term}}$ explicitly, we
have chosen representative EOSs that fulfill the QCD constraint at
$n_{\mathrm{term}}=n_{\mathrm{TOV}}$ but not necessarily above. For each of
these EOSs, we generate 5000 extensions between
$n_{\mathrm{TOV}}<n<15n_{\mathrm{sat}}$ using the GP approach. The extensions
are then conditioned with the QCD constraint at $15n_{\mathrm{sat}}$ so that
the behavior between $n_{\mathrm{TOV}}<n<15n_{\mathrm{sat}}$ remains
consistent with the QCD input. These extensions show how the EOS must evolve
beyond $n_{\mathrm{TOV}}$ given its behavior below $n_{\mathrm{TOV}}$. Figure
8 shows the extensions for representative EOSs with
$I_{\mathrm{pQCD}}=0.5,0.75$ and $0.85$ (for $I_{\mathrm{pQCD}}=0.98$, our
prior does not contain any valid samples). We have checked that these EOSs are
indeed representative of EOSs with these pQCD tension indices. We find a
difference in behavior above $n_{\mathrm{TOV}}$ based on the value of the pQCD
tension index at the TOV point $I_{\mathrm{pQCD}}(n_{\mathrm{TOV}})$. For the
highest index of $I_{\mathrm{pQCD}}(n_{\mathrm{TOV}})=0.85$ shown in Fig. 8,
we observe that the EOS must soften drastically above $n_{\mathrm{TOV}}$, with
$c_{s}^{2}(n)\leq 0.1$ all the way to $n=15n_{\mathrm{sat}}$. Furthermore, for
this EOS, the value of the index at $15n_{\mathrm{sat}}$ remains in the range
$I_{\mathrm{pQCD}}(15n_{\mathrm{sat}})\in[0.9,1]$, indicating that the EOSs
needs to remain very soft. For the smaller value of the pQCD tension index
$I_{\mathrm{pQCD}}(n_{\mathrm{TOV}})=0.75$, we similarly observe that the
index at $15n_{\mathrm{sat}}$ remains close to its value at $n_{\mathrm{TOV}}$
and the EOS must continue to remain soft. Finally, for the lowest
$I_{\mathrm{pQCD}}(n_{\mathrm{TOV}})=0.5$ shown, the EOS extensions can obtain
significantly different values of the pQCD tension index at
$15n_{\mathrm{sat}}$.

We quantify the degree of softening as a function of pQCD tension index in
Fig. 9, where for 100 EOSs in our posterior with fixed values of the pQCD
tension index $I_{\mathrm{pQCD}}(n_{\mathrm{TOV}})$, we show the distribution
of the average speed of sound of the extensions on the interval
$n\in[8,15]n_{\mathrm{sat}}$. From this figure we see that the average
$c_{s}^{2}$ for the extensions of EOSs with $I_{\mathrm{pQCD}}=0.85,0.75$ are
strongly skewed to small values, below $c_{s}^{2}\lesssim 0.03$ and
$c_{s}^{2}\lesssim 0.11$, respectively (95% credence). By contrast, the
extensions of EOSs with $I_{\mathrm{pQCD}}=0.5$ have a very large spread in
the average $c_{s}^{2}$, with a shape that follows the prior distribution,
indicating that there is no forced behavior in this case.

Figure 10: (Left) A sample of EOSs drawn from the hierarchical models of EOS
extensions with correlation length
$\ell\sim\mathcal{U}(1n_{\mathrm{sat}},20n_{\mathrm{sat}})$, the width
$\eta\sim\mathcal{N}(1.25,0.25^{2})$ and average
$\bar{c}_{s}^{2}=\mathcal{N}(0.3,0.3^{2})$. Negative values of $\eta$ and
$\bar{c}_{s}^{2}$ are removed. The dashed prior GP is anchored to the pQCD EOS
at $\mu_{\mathrm{pQCD}}=2.6$ GeV, while the solid magenta GP is conditioned
with the pQCD speed of sound on the interval $n\in[25,40]n_{\mathrm{sat}}$
(see main text). (Middle) The resulting QCD likelihood function at
$n_{\mathrm{term}}=5.5n_{\mathrm{sat}}$ for the GP conditioned on the pQCD
$c_{s}^{2}(n)$ for $n\in[25,40]n_{\mathrm{sat}}$. This termination density
corresponds approximately to the mean value of $n_{\mathrm{TOV}}$ for our
ensembles after applying the astrophysical input. Shown in blue and green are
the constraints arising from the chiral EFT calculation of Ref. Hebeler  _et
al._ (2013) and the $X=1$ pQCD calculations of Ref. Gorda  _et al._ (2021)
alone, using the construction of Ref. Komoltsev and Kurkela (2022). (Right)
The posterior $\varepsilon$–$p$ regions at $n_{\mathrm{TOV}}$ after
marginalizing over the EOS extensions between $n_{\mathrm{term}}$ and
40$n_{\mathrm{sat}}$.

In the previous sections, we have seen that the constraining power of the QCD
input increases as function of the termination density when going from
$n_{\mathrm{term}}=n_{\mathrm{TOV}}$ to
$n_{\mathrm{term}}=1.2n_{\mathrm{TOV}}$; see, e.g., Fig. 2. In Fig. 6 we show
the distribution of the pQCD tension indices
$I_{\mathrm{pQCD}}(n_{\mathrm{TOV}})$ for EOSs that are allowed at
$n_{\mathrm{term}}=1.2n_{\mathrm{TOV}}$ and for EOS that are allowed at
$n_{\mathrm{term}}=n_{\mathrm{TOV}}$ but rejected before
$n_{\mathrm{term}}=1.2n_{\mathrm{TOV}}$. The majority of the EOS that are
excluded at higher $n_{\mathrm{term}}$ are EOSs which have large indices
$I_{\mathrm{pQCD}}(n_{\mathrm{TOV}})$, forcing the EOS to undergo a rapid
softening beyond $n_{\mathrm{TOV}}$. This behavior cannot be excluded by the
rigorous physics principles behind our implementation of the QCD constraint,
but it would be interesting to understand what microphysical behavior could
lead to such behavior.

###  IV.4 Marginalization over EOS extensions beyond $n_{\mathrm{term}}$

The termination of the explicit EOS models at $n_{\mathrm{term}}$ puts the EOS
below and above $n_{\mathrm{term}}$ in an asymmetric role; below the
termination density our prior contains a large number of non-extreme EOSs by
construction. Above $n_{\mathrm{term}}$, even the most extreme EOS extension
of Fig. 7 is realized without any penalty. This has been criticized in Ref.
Essick  _et al._ (2023), in which the authors point out that the construction
of Ref. Komoltsev and Kurkela (2022) is similar to performing a maximization,
rather than a marginalization over the likelihoods of different possible
extensions. As our final point in this work, we address the marginalization
over a set of extensions of the EOS beyond $n_{\mathrm{term}}$. We do so by
constructing a model of possible high-density EOSs down from the pQCD EOS at
$\mu_{\mathrm{pQCD}}=2.6$ GeV, distinct from the prior model of NS EOSs below
$n_{\mathrm{term}}$. Having an explicit model for EOSs also above
$n_{\mathrm{term}}$ allows us to incorporate further theoretical inputs. Since
the pQCD speed of sound is well converged at high densities, see Fig. 7, we
can further condition this model with the pQCD speed of sound. We may do so
either only at the pQCD matching point $n_{\mathrm{pQCD}}\approx
40n_{\mathrm{sat}}$(referred to below as “prior”), or over a larger range of
densities where the speed of sound is under good theoretical control,
$[25,40]n_{\mathrm{sat}}$ (“conditioned”).

For modeling the EOS above $n_{\mathrm{term}}$, we define another hierarchical
GP model anchored to the (scale-averaged) pQCD EOS defined with the following
hyperparameters:

$\begin{split}\ell&\sim\mathcal{U}(1n_{\mathrm{sat}},20n_{\mathrm{sat}})\,,\\\
\eta&\sim\mathcal{N}(1.25,0.25^{2})\,,\\\
\bar{c}_{s}^{2}&\sim\mathcal{N}(0.3,0.3^{2})\,,\end{split}$ (12)

see Fig. 10 (left). Because of the large interval in density between the
typical $n_{\mathrm{term}}$ and $n_{\mathrm{QCD}}\approx 40n_{\mathrm{sat}}$,
we choose a more general prior for the correlation length $\ell$ here than was
done in Ref. Gorda  _et al._ (2023b). We also choose the prior on
$\bar{c}_{s}^{2}$ to have a mean centered near the conformal value of $1/3$,
but allow it to vary broadly within the model. To be conservative with the
conditioned GP, we take the standard deviation of the GP at a given
$n\in[25,40]n_{\mathrm{sat}}$ to equal twice of that of the scale-averaged
pQCD calculation.

The procedure for applying this high-density prior as a model for EOS
extensions is as follows. First, we use either the low-density GP or CSM to
model the NS EOS up to a specified termination density $n_{\mathrm{term}}$.
Then, for this $n_{\mathrm{term}}$, we calculate a QCD likelihood function by
constructing the posterior distribution for different $\varepsilon$–$p$ values
that arise from marginalizing over the new high-density GP model above
$n_{\mathrm{term}}$. Concretely, we calculate the distribution by performing a
kernel-density estimation of the $\varepsilon$–$p$ values of the high-density
GP model at $n_{\mathrm{term}}$. This posterior distribution is then
reinterpreted as a new QCD likelihood function within the NS EOS inference at
$n_{\mathrm{term}}$ that marginalizes over the hierarchical model above, as
opposed to the construction from Ref. Komoltsev and Kurkela (2022), which
instead returns a weight of 1 for any $\varepsilon$–$p$ point within the
allowed region at $n_{\mathrm{term}}$ for a fixed $X$.

In Fig. 10 we show in the left panel a small sample of EOSs drawn from these
models of EOS extensions, showing a broad range of possible behaviors. The
conditioning of the magenta GP is clearly visible. In the middle panel we show
the QCD likelihood function resulting from the conditioned GP as applied at
$n_{\mathrm{term}}=5.5n_{\mathrm{sat}}$, compared to the constraints arising
from Ref. Komoltsev and Kurkela (2022) for the chiral EFT and pQCD
constraints. This value of $n_{\mathrm{term}}$ approximately corresponds to
the mean of the $n_{\mathrm{TOV}}$ for our ensembles using only the
astrophysical input (see Fig. 2). We see that the likelihood function is
mostly peaked in the interior of the $X=1$ pQCD region. As has been pointed
out already in Ref. Gorda  _et al._ (2023d), the impact of the QCD input
within the inference of the NS EOS occurs primarily along the upper boundary
of the pQCD region in Fig. 10, since the astrophysical constraints already
reduce the allowed $\varepsilon$–$p$ region in the neighborhood of the other
extreme boundaries (cf. also the discussion of the upper edge of the
$\varepsilon$–$p$ regions in Fig. 5).

Turning to the rightmost panel of Fig. 10, we see the effect of the
marginalization over the GP model of EOS extensions directly as the change in
the allowed $\varepsilon$–$p$ region at $n_{\mathrm{TOV}}$. As expected from
our observations of the central panel, we see a shift towards smaller $p$
values in this figure. The prior GP shows a similar shift to the effect of
using a termination density that is about 20-40% larger without the
marginalization (cf. the middle and right columns of Fig. 5), while the
conditioned GP shows a much more pronounced softening, more comparable to the
effect of the QCD input applied at $n_{\mathrm{term}}=10n_{\mathrm{sat}}$ (cf.
Fig. 5). In both cases, the quantitative effect is to shift the posterior
pressure distribution within the most massive NSs to smaller values, with the
conditioned GP shifted more strongly. That the upper limit for the conditioned
GP is significantly below that of the prior demonstrates that the additional
constraint is indeed due to the conditioning and not the assumed prior. This
decreased upper bound is natural, as the pQCD sound speed used in the
conditioning is inconsistent with the stiff high-density segment seen in Fig.
7, thus affecting the EOS at lower densities as well.

Lastly, we observe that the credible regions of $\varepsilon$–$p$ are almost
indistinguishable between $n_{\mathrm{term}}=n_{\mathrm{TOV}}$ and
$10n_{\mathrm{sat}}$, demonstrating independence of this parameter when using
the conditioned GP extension.

##  V Discussion and conclusions

Figure 11: The impact of QCD input at the densities reached in the cores of
maximally massive NS using different prescription for the EOS extension.
(Left) Pressure as a function of energy density. For those prescriptions with
$n_{\mathrm{term}}=n_{\mathrm{TOV}}$ the 68% credible intervals for the
pressure are extended only up to lower bound of the 95% credible interval of
the TOV energy densities. Above this point we cannot access the posterior
distribution $P(p|\varepsilon)$ but have access only to the conditional
distribution $P(p|\varepsilon,\varepsilon<\varepsilon_{\mathrm{TOV}})$.
(Right) The posterior $\varepsilon$–$p$ regions at $n_{\mathrm{TOV}}$.

In this paper, we have studied the impact of the high-density QCD input on the
modeling of the dense-matter EOS, and explained the different findings
presented in the literature. In particular, using the most general QCD input
we found that the constraining power of the QCD input depends sensitively on
the chosen termination density $n_{\mathrm{term}}$. For the conservative
choice of $n_{\mathrm{term}}=n_{\mathrm{TOV}}$, the impact of the QCD input is
largely removed, while the constraining power rises fast with
$n_{\mathrm{term}}$. This difference in constraining power just beyond
$n_{\mathrm{TOV}}$ arises from EOSs that undergo a change in behavior that is
qualitatively different from the assumed priors. In particular, we have shown
that most EOSs that are removed for $n_{\mathrm{term}}>n_{\mathrm{TOV}}$ have
a large pQCD tension index, leading to an abrupt softening above
$n_{\mathrm{TOV}}$ that is disfavored when the EOS priors employed here are
constructed to higher densities. Hence, to remain in agreement with the QCD
input, the EOS needs either to soften below $n_{\mathrm{TOV}}$ or it has to
exhibit a drastic softening immediately above $n_{\mathrm{TOV}}$ that persists
for an extended density range before stiffening again dramatically, see also
Ref. Annala  _et al._ (2020). This behavior is seen in the EOS extensions
depicted in Fig. 8. We reiterate that the strong softening coinciding with
$n_{\mathrm{TOV}}$ is accidental, as $n_{\mathrm{TOV}}$ is already set by the
behavior of the EOS at lower densities. While it cannot be excluded by the
generic thermodynamic considerations we employ, it would be interesting to see
what microphysical behavior would be necessary to realize such behavior within
dense matter. Note, however, that high speeds of sound above $n\gtrsim
25n_{\mathrm{sat}}$ are inconsistent with pQCD calculations.

We note that different choices for $n_{\mathrm{term}}$,
$n_{\mathrm{term}}=n_{\mathrm{TOV}}$ and $n_{\mathrm{term}}>n_{\mathrm{TOV}}$,
have their merits and answer slightly different questions: the choice of the
termination density eventually reflects the phenomenological question at hand.
On the one hand, the EOS is a fundamental prediction of QCD and exists
independently of NSs at densities above $n_{\mathrm{TOV}}$. The EOS above
$n_{\mathrm{TOV}}$ is required in modeling the post-merger phase of binary NS
mergers (see Refs. Ujevic  _et al._ (2023); Tootle  _et al._ (2022)111Note
that Ref. Ujevic  _et al._ (2023) found that multi-messenger signals from
binary NS merges are not effective at constraining the EOS above
$n_{\mathrm{TOV}}$, rendering the QCD input the only source of information in
this density range in the near future.) and sets a prediction for the quantity
that may be verified in the future using, e.g., lattice field theory222On this
point, see in particular two recent works in Refs. Moore and Gorda (2023);
Fujimoto and Reddy (2023).. In order to address these questions, the EOS must
be modeled beyond the densities reached in stable NSs, and
$n_{\mathrm{term}}>n_{\mathrm{TOV}}$ is required. On the other hand, one can
be more conservative and use $n_{\mathrm{term}}=n_{\mathrm{TOV}}$ if one
chooses to infer only the EOS along the stable NS branch. This choice avoids
prior-dependent modeling of the unstable branch, allowing for more robust
conclusions. Similarly, for an even more conservative modeling of NSs, one
could imagine terminating the EOS model at $n_{\mathrm{term}}=n(2.1M_{\odot})$
which is the mass of the heaviest observed neutron star. Whatever application
one has in mind, it is crucial that the QCD input in the form of Ref. Gorda
_et al._ (2023b) be applied at the highest density that will be studied. It is
inconsistent to apply the QCD input at $n_{\mathrm{term}}$ and then consider
the model for even higher densities $n>n_{\mathrm{term}}$.

Since the pQCD speed of sound is well converged to even lower densities than
the pressure Gorda  _et al._ (2021, 2023a), we have additionally constructed a
QCD likelihood function that incorporates this additional information. We have
done so by constructing a further GP model down from the pQCD EOS calculation
and conditioning it with the pQCD calculation of $c_{s}^{2}$. This also
addresses the points raised in Ref. Essick  _et al._ (2023), as this new
likelihood function involves marginalizing over the new GP above
$n_{\mathrm{term}}$. We find that the conditioned GP gives results that are
independent of $n_{\mathrm{term}}\geq n_{\mathrm{TOV}}$ (see Fig. 10) and
agrees with results obtained for $n_{\mathrm{term}}=10n_{\mathrm{sat}}$, see
Fig. 11. The conditioned likelihood function is made publicly available in
Ref. Komoltsev  _et al._ (2023). We also see in Fig. 11 that the unconditioned
GP leads to similar results as choosing
$n_{\mathrm{term}}=1.2n_{\mathrm{TOV}}$ or by restricting the pQCD tension
index $I_{\mathrm{pQCD}}(n_{\mathrm{TOV}})<0.75$.

Finally, we discuss the limitations of our study. An important finding of this
work is that the QCD input removes EOS that are very stiff around
$n_{\mathrm{TOV}}$. However, we allow for extreme stiffening of the EOS at low
densities, around $1-2n_{\mathrm{sat}}$, where EFT approaches can otherwise
provide constraints on the EOS and typically suggest the EOS to be soft Tews
_et al._ (2018a); Drischler  _et al._ (2020). In this work, we have not varied
the maximal density up to which we consider EFT input, which might bias our
EOS sets to be stiffer than if we marginalized over the maximal EFT density
Capano  _et al._ (2019); Essick  _et al._ (2020). We also have not
marginalized over the scale $\mu_{\mathrm{pQCD}}$ where we evaluate the pQCD
values; marginalizing over lower values than used here would provide a more
restrictive QCD input Gorda  _et al._ (2023c). Furthermore, we only considered
two EOS priors that model the EOS in terms of the speed of sound and are,
hence, rather similar. We have also not considered astrophysical data from
other X-ray sources such as quiescent low-mass X-ray binaries or photospheric-
radius–expansion burst, which decrease the upper bound on the radii with
respect to the NICER X-ray data Al-Mamun  _et al._ (2021). Finally, for the
purpose of marginalizing over EOS extensions above $n_{\mathrm{term}}$, it
would be interesting to investigate other prior choices for the EOS model,
such as GP models with a density dependent correlation length.

###### Acknowledgements.

We thank Reed Essick, Joonas Hirvonen, Débora Mroczek, Jorge Noronha, Jaki
Noronha-Hostler, Krishna Rajagopal, and Rachel Steinhorst for useful
discussions. This work has been supported in part by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) project-ID
279384907–SFB 1245, by the State of Hesse within the Research Cluster ELEMENTS
(projectID 500/10.006), and by the ERC Advanced Grant “JETSET: Launching,
propagation and emission of relativistic jets from binary mergers and across
mass scales” (Grant No. 884631) (T.G.). R.S. acknowledges support from the
Nuclear Physics from Multi-Messenger Mergers (NP3M) Focused Research Hub which
is funded by the National Science Foundation under Grant Number 21-16686, and
by the Laboratory Directed Research and Development program of Los Alamos
National Laboratory under project number 20220541ECR. J.M. is supported by the
CNRS/IN2P3 NewMAC project, and are also grateful to PHAROS COST Action MP16214
and to the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) of the
Université de Lyon for its financial support within the program
Investissements d’Avenir (ANR-11-IDEX-0007) of the French government operated
by the National Research Agency (ANR). The work of I.T. was supported by the
U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under
contract No. DE-AC52-06NA25396, by the Laboratory Directed Research and
Development program of Los Alamos National Laboratory under project number
20230315ER, and by the U.S. Department of Energy, Office of Science, Office of
Advanced Scientific Computing Research, Scientific Discovery through Advanced
Computing (SciDAC) NUCLEI program. Author contributions: O.K. performed most
of the numerical analysis, with contributions from R.S. T.G. wrote the scripts
for the marginalization and conditioning of the high-density GP. All authors
contributed to the preparation and revision of the manuscript.

## References

  * Kurkela  _et al._ (2010) A. Kurkela, P. Romatschke, and A. Vuorinen, Phys. Rev. D 81, 105021 (2010), arXiv:0912.1856 [hep-ph] .
  * Kurkela and Vuorinen (2016) A. Kurkela and A. Vuorinen, Phys. Rev. Lett. 117, 042501 (2016), arXiv:1603.00750 [hep-ph] .
  * Gorda  _et al._ (2018) T. Gorda, A. Kurkela, P. Romatschke, S. Säppi, and A. Vuorinen, Phys. Rev. Lett. 121, 202701 (2018), arXiv:1807.04120 [hep-ph] .
  * Gorda  _et al._ (2021) T. Gorda, A. Kurkela, R. Paatelainen, S. Säppi, and A. Vuorinen, Phys. Rev. Lett. 127, 162003 (2021), arXiv:2103.05658 [hep-ph] .
  * Gorda  _et al._ (2023a) T. Gorda, R. Paatelainen, S. Säppi, and K. Seppänen, Phys. Rev. Lett. 131, 181902 (2023a), arXiv:2307.08734 [hep-ph] .
  * Schaffner-Bielich (2020) J. Schaffner-Bielich,  _Compact Star Physics_ (Cambridge University Press, 2020).
  * Epelbaum  _et al._ (2009) E. Epelbaum, H.-W. Hammer, and U.-G. Meissner, Rev. Mod. Phys. 81, 1773 (2009), arXiv:0811.1338 [nucl-th] .
  * Machleidt and Entem (2011) R. Machleidt and D. R. Entem, Phys. Rept. 503, 1 (2011), arXiv:1105.2919 [nucl-th] .
  * Hebeler  _et al._ (2013) K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Astrophys. J. 773, 11 (2013), arXiv:1303.4662 [astro-ph.SR] .
  * Tews  _et al._ (2018a) I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astrophys. J. 860, 149 (2018a), arXiv:1801.01923 [nucl-th] .
  * Drischler  _et al._ (2019) C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev. Lett. 122, 042501 (2019), arXiv:1710.08220 [nucl-th] .
  * Read  _et al._ (2009) J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Phys. Rev. D79, 124032 (2009), arXiv:0812.2163 [astro-ph] .
  * Lindblom (2010) L. Lindblom, Phys. Rev. D82, 103011 (2010), arXiv:1009.0738 [astro-ph.HE] .
  * Lindblom and Indik (2012) L. Lindblom and N. M. Indik, Phys. Rev. D 86, 084003 (2012), arXiv:1207.3744 [astro-ph.HE] .
  * Lindblom and Indik (2014) L. Lindblom and N. M. Indik, Phys. Rev. D 89, 064003 (2014), [Erratum: Phys.Rev.D 93, 129903 (2016)], arXiv:1310.0803 [astro-ph.HE] .
  * Kurkela  _et al._ (2014) A. Kurkela, E. S. Fraga, J. Schaffner-Bielich, and A. Vuorinen, Astrophys. J. 789, 127 (2014), arXiv:1402.6618 [astro-ph.HE] .
  * Annala  _et al._ (2018) E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018), arXiv:1711.02644 [astro-ph.HE] .
  * Most  _et al._ (2018) E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich, Phys. Rev. Lett. 120, 261103 (2018), arXiv:1803.00549 [gr-qc] .
  * Tews  _et al._ (2018b) I. Tews, J. Margueron, and S. Reddy, Phys. Rev. C 98, 045804 (2018b), arXiv:1804.02783 [nucl-th] .
  * Landry and Essick (2019) P. Landry and R. Essick, Phys. Rev. D 99, 084049 (2019), arXiv:1811.12529 [gr-qc] .
  * Greif  _et al._ (2019) S. K. Greif, G. Raaijmakers, K. Hebeler, A. Schwenk, and A. L. Watts, Mon. Not. Roy. Astron. Soc. 485, 5363 (2019), arXiv:1812.08188 [astro-ph.HE] .
  * Annala  _et al._ (2020) E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Nature Phys. 16, 907 (2020), arXiv:1903.09121 [astro-ph.HE] .
  * Raaijmakers  _et al._ (2020) G. Raaijmakers _et al._ , Astrophys. J. Lett. 893, L21 (2020), arXiv:1912.11031 [astro-ph.HE] .
  * Essick  _et al._ (2020) R. Essick, I. Tews, P. Landry, S. Reddy, and D. E. Holz, Phys. Rev. C 102, 055803 (2020), arXiv:2004.07744 [astro-ph.HE] .
  * Annala  _et al._ (2022) E. Annala, T. Gorda, E. Katerini, A. Kurkela, J. Nättilä, V. Paschalidis, and A. Vuorinen, Phys. Rev. X 12, 011058 (2022), arXiv:2105.05132 [astro-ph.HE] .
  * Altiparmak  _et al._ (2022) S. Altiparmak, C. Ecker, and L. Rezzolla, Astrophys. J. Lett. 939, L34 (2022), arXiv:2203.14974 [astro-ph.HE] .
  * Somasundaram  _et al._ (2023a) R. Somasundaram, I. Tews, and J. Margueron, Phys. Rev. C 107, 025801 (2023a), arXiv:2112.08157 [nucl-th] .
  * Jiang  _et al._ (2023) J.-L. Jiang, C. Ecker, and L. Rezzolla, Astrophys. J. 949, 11 (2023), arXiv:2211.00018 [gr-qc] .
  * Gorda  _et al._ (2023b) T. Gorda, O. Komoltsev, and A. Kurkela, Astrophys. J. 950, 107 (2023b), arXiv:2204.11877 [nucl-th] .
  * Annala  _et al._ (2023) E. Annala, T. Gorda, J. Hirvonen, O. Komoltsev, A. Kurkela, J. Nättilä, and A. Vuorinen, (2023), arXiv:2303.11356 [astro-ph.HE] .
  * Han  _et al._ (2023) M.-Z. Han, Y.-J. Huang, S.-P. Tang, and Y.-Z. Fan, Sci. Bull. 68, 913 (2023), arXiv:2207.13613 [astro-ph.HE] .
  * Komoltsev and Kurkela (2022) O. Komoltsev and A. Kurkela, Phys. Rev. Lett. 128, 202701 (2022), arXiv:2111.05350 [nucl-th] .
  * Gorda  _et al._ (2023c) T. Gorda, O. Komoltsev, A. Kurkela, and A. Mazeliauskas, JHEP 06, 002 (2023c), arXiv:2303.02175 [hep-ph] .
  * Somasundaram  _et al._ (2023b) R. Somasundaram, I. Tews, and J. Margueron, Phys. Rev. C 107, L052801 (2023b), arXiv:2204.14039 [nucl-th] .
  * Ecker and Rezzolla (2022) C. Ecker and L. Rezzolla, Mon. Not. Roy. Astron. Soc. 519, 2615 (2022), arXiv:2209.08101 [astro-ph.HE] .
  * Essick  _et al._ (2023) R. Essick, I. Legred, K. Chatziioannou, S. Han, and P. Landry, Phys. Rev. D 108, 043013 (2023), arXiv:2305.07411 [astro-ph.HE] .
  * Brandes  _et al._ (2023) L. Brandes, W. Weise, and N. Kaiser, Phys. Rev. D 108, 094014 (2023), arXiv:2306.06218 [nucl-th] .
  * Musolino  _et al._ (2023) C. Musolino, C. Ecker, and L. Rezzolla, (2023), arXiv:2307.03225 [gr-qc] .
  * Mroczek  _et al._ (2023) D. Mroczek, M. C. Miller, J. Noronha-Hostler, and N. Yunes, (2023), arXiv:2309.02345 [astro-ph.HE] .
  * Fan  _et al._ (2023) Y.-Z. Fan, M.-Z. Han, J.-L. Jiang, D.-S. Shao, and S.-P. Tang, (2023), arXiv:2309.12644 [astro-ph.HE] .
  * Tang  _et al._ (2023) S.-P. Tang, M.-Z. Han, Y.-J. Huang, Y.-Z. Fan, and D.-M. Wei, (2023), arXiv:2311.13805 [astro-ph.HE] .
  * Abbott  _et al._ (2017) B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017).
  * Abbott  _et al._ (2018) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, _et al._ (The LIGO Scientific Collaboration and the Virgo Collaboration), Phys. Rev. Lett. 121, 161101 (2018).
  * Abbott  _et al._ (2019) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. X 9, 011001 (2019).
  * Antoniadis  _et al._ (2013) J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, M. H. van Kerkwijk, M. Kramer, C. Bassa, V. S. Dhillon, T. Driebe, J. W. T. Hessels, V. M. Kaspi, V. I. Kondratiev, N. Langer, T. R. Marsh, M. A. McLaughlin, T. T. Pennucci, S. M. Ransom, I. H. Stairs, J. van Leeuwen, J. P. W. Verbiest, and D. G. Whelan, Science 340, 1233232 (2013), https://www.science.org/doi/pdf/10.1126/science.1233232 .
  * Cromartie  _et al._ (2019) H. T. Cromartie, E. Fonseca, S. M. Ransom, _et al._ (NANOGrav), Nature Astron. 4, 72 (2019), arXiv:1904.06759 [astro-ph.HE] .
  * Fonseca  _et al._ (2021) E. Fonseca, H. T. Cromartie, T. T. Pennucci, P. S. Ray, A. Y. Kirichenko, S. M. Ransom, P. B. Demorest, I. H. Stairs, Z. Arzoumanian, L. Guillemot, A. Parthasarathy, M. Kerr, I. Cognard, P. T. Baker, H. Blumer, P. R. Brook, M. DeCesar, T. Dolch, F. A. Dong, E. C. Ferrara, W. Fiore, N. Garver-Daniels, D. C. Good, R. Jennings, M. L. Jones, V. M. Kaspi, M. T. Lam, D. R. Lorimer, J. Luo, A. McEwen, J. W. McKee, M. A. McLaughlin, N. McMann, B. W. Meyers, A. Naidu, C. Ng, D. J. Nice, N. Pol, H. A. Radovan, B. Shapiro-Albert, C. M. Tan, S. P. Tendulkar, J. K. Swiggum, H. M. Wahl, and W. W. Zhu, Astrophys. J. Lett. 915, L12 (2021).
  * Steiner  _et al._ (2018) A. W. Steiner, C. O. Heinke, S. Bogdanov, C. Li, W. C. G. Ho, A. Bahramian, and S. Han, Mon. Not. Roy. Astron. Soc. 476, 421 (2018), arXiv:1709.05013 [astro-ph.HE] .
  * Nättilä  _et al._ (2017) J. Nättilä, M. C. Miller, A. W. Steiner, J. J. E. Kajava, V. F. Suleimanov, and J. Poutanen, Astron. Astrophys. 608, A31 (2017), arXiv:1709.09120 [astro-ph.HE] .
  * Shaw  _et al._ (2018) A. W. Shaw, C. O. Heinke, A. W. Steiner, S. Campana, H. N. Cohn, W. C. G. Ho, P. M. Lugger, and M. Servillat, Mon. Not. R. Astron. Soc. 476, 4713 (2018), arXiv:1803.00029 [astro-ph.HE] .
  * Miller  _et al._ (2019a) M. C. Miller, F. K. Lamb, A. J. Dittmann, S. Bogdanov, Z. Arzoumanian, K. C. Gendreau, S. Guillot, A. K. Harding, W. C. G. Ho, J. M. Lattimer, R. M. Ludlam, S. Mahmoodifar, S. M. Morsink, P. S. Ray, T. E. Strohmayer, K. S. Wood, T. Enoto, R. Foster, T. Okajima, G. Prigozhin, and Y. Soong, Astrophys. J. Lett. 887, L24 (2019a).
  * Riley  _et al._ (2019) T. E. Riley, A. L. Watts, S. Bogdanov, P. S. Ray, R. M. Ludlam, S. Guillot, Z. Arzoumanian, C. L. Baker, A. V. Bilous, D. Chakrabarty, K. C. Gendreau, A. K. Harding, W. C. G. Ho, J. M. Lattimer, S. M. Morsink, and T. E. Strohmayer, Astrophys. J. Lett. 887, L21 (2019).
  * Miller  _et al._ (2021) M. C. Miller _et al._ , Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Riley  _et al._ (2021) T. E. Riley, A. L. Watts, P. S. Ray, S. Bogdanov, S. Guillot, S. M. Morsink, A. V. Bilous, Z. Arzoumanian, D. Choudhury, J. S. Deneva, K. C. Gendreau, A. K. Harding, W. C. G. Ho, J. M. Lattimer, M. Loewenstein, R. M. Ludlam, C. B. Markwardt, T. Okajima, C. Prescod-Weinstein, R. A. Remillard, M. T. Wolff, E. Fonseca, H. T. Cromartie, M. Kerr, T. T. Pennucci, A. Parthasarathy, S. Ransom, I. Stairs, L. Guillemot, and I. Cognard, Astrophys. J. Lett. 918, L27 (2021).
  * Douchin and Haensel (2001) F. Douchin and P. Haensel, Astron. Astrophys. 380, 151 (2001), arXiv:astro-ph/0111092 .
  * Margueron  _et al._ (2018a) J. Margueron, R. Hoffmann Casali, and F. Gulminelli, Phys. Rev. C 97, 025805 (2018a), arXiv:1708.06894 [nucl-th] .
  * Margueron  _et al._ (2018b) J. Margueron, R. Hoffmann Casali, and F. Gulminelli, Phys. Rev. C 97, 025806 (2018b), arXiv:1708.06895 [nucl-th] .
  * Lynn  _et al._ (2016) J. E. Lynn, I. Tews, J. Carlson, S. Gandolfi, A. Gezerlis, K. E. Schmidt, and A. Schwenk, Phys. Rev. Lett. 116, 062501 (2016), arXiv:1509.03470 [nucl-th] .
  * Baym  _et al._ (1971) G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. 170, 299 (1971).
  * Fonseca  _et al._ (2016) E. Fonseca _et al._ , Astrophys. J. 832, 167 (2016), arXiv:1603.00545 [astro-ph.HE] .
  * Miller  _et al._ (2019b) M. C. Miller, C. Chirenti, and F. K. Lamb,  (2019b), 10.3847/1538-4357/ab4ef9, arXiv:1904.08907 [astro-ph.HE] .
  * Landry  _et al._ (2020) P. Landry, R. Essick, and K. Chatziioannou, Phys. Rev. D 101, 123007 (2020), arXiv:2003.04880 [astro-ph.HE] .
  * Pang  _et al._ (2023) P. T. H. Pang, L. Sivertsen, R. Somasundaram, T. Dietrich, S. Sen, I. Tews, M. Coughlin, and C. Van Den Broeck, (2023), arXiv:2308.15067 [nucl-th] .
  * Gorda  _et al._ (2023d) T. Gorda, K. Hebeler, A. Kurkela, A. Schwenk, and A. Vuorinen, Astrophys. J. 955, 100 (2023d), arXiv:2212.10576 [astro-ph.HE] .
  * Ujevic  _et al._ (2023) M. Ujevic, R. Somasundaram, T. Dietrich, J. Margueron, and I. Tews, (2023), arXiv:2311.04809 [astro-ph.HE] .
  * Tootle  _et al._ (2022) S. Tootle, C. Ecker, K. Topolski, T. Demircik, M. Järvinen, and L. Rezzolla, SciPost Phys. 13, 109 (2022), arXiv:2205.05691 [astro-ph.HE] .
  * Moore and Gorda (2023) G. D. Moore and T. Gorda, (2023), arXiv:2309.15149 [nucl-th] .
  * Fujimoto and Reddy (2023) Y. Fujimoto and S. Reddy, (2023), arXiv:2310.09427 [nucl-th] .
  * Komoltsev  _et al._ (2023) O. Komoltsev, T. Gorda, and A. Kurkela, “Marginalized QCD likelihood function,”  (2023).
  * Drischler  _et al._ (2020) C. Drischler, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, Phys. Rev. Lett. 125, 202702 (2020), arXiv:2004.07232 [nucl-th] .
  * Capano  _et al._ (2019) C. D. Capano, I. Tews, S. M. Brown, B. Margalit, S. De, S. Kumar, D. A. Brown, B. Krishnan, and S. Reddy, Nature Astron. (2019), 10.1038/s41550-020-1014-6, arXiv:1908.10352 [astro-ph.HE] .
  * Al-Mamun  _et al._ (2021) M. Al-Mamun, A. W. Steiner, J. Nättilä, J. Lange, R. O’Shaughnessy, I. Tews, S. Gandolfi, C. Heinke, and S. Han, Phys. Rev. Lett. 126, 061101 (2021), arXiv:2008.12817 [astro-ph.HE] .