#  Implicit correlations within phenomenological parametric models of the
neutron star equation of state

Isaac Legred  Department of Physics, California Institute of Technology,
Pasadena, California 91125, USA  LIGO Laboratory, California Institute of
Technology, Pasadena, CA 91125, USA  Katerina Chatziioannou  Department of
Physics, California Institute of Technology, Pasadena, California 91125, USA
LIGO Laboratory, California Institute of Technology, Pasadena, CA 91125, USA
Reed Essick  Perimeter Institute for Theoretical Physics, Waterloo, Ontario
N2L 2Y5, Canada  Philippe Landry  Canadian Institute for Theoretical
Astrophysics, University of Toronto, Toronto, Ontario M5S 3H8, Canada

###### Abstract

The rapid increase in the number and precision of astrophysical probes of
neutron stars in recent years allows for the inference of their equation of
state. Observations target different macroscopic properties of neutron stars
which vary from star to star, such as mass and radius, but the equation of
state allows for a common description of all neutron stars. To connect these
observations and infer the properties of dense matter and neutron stars
simultaneously, models for the equation of state are introduced. _Parametric
models_ rely on carefully engineered functional forms that reproduce a large
array of realistic equations of state. Such models benefit from their
simplicity but are limited because any finite-parameter model cannot
accurately approximate all possible equations of state. _Nonparametric
methods_ overcome this by increasing model freedom at the cost of increased
complexity. In this study, we compare common parametric and nonparametric
models, quantify the limitations of the former, and study the impact of
modeling on our current understanding of high-density physics. We show that
parametric models impose strongly model-dependent, and sometimes opaque,
correlations between density scales. Such inter-density correlations result in
tighter constraints that are unsupported by data and can lead to biased
inference of the equation of state and of individual neutron star properties.

##  I Introduction

The equation of state (EoS) of the dense matter inside neutron stars (NSs) is
uncertain at densities near and beyond nuclear saturation,
$\rho_{\mathrm{nuc}}=2.8\times 10^{14}\mathrm{g/cm^{3}}$, because it cannot be
precisely constrained by theoretical calculations or terrestrial experiments
Lattimer and Prakash (2016); Özel and Freire (2016); Oertel  _et al._ (2017);
Baym  _et al._ (2018); Machleidt and Entem (2011); Han  _et al._ (2019);
Chatziioannou (2020). Astronomical observations Abbott  _et al._ (2017);
Cromartie  _et al._ (2019); Miller  _et al._ (2019); Riley  _et al._ (2019);
Miller  _et al._ (2021); Riley  _et al._ (2021); Antoniadis  _et al._ (2016);
Abbott  _et al._ (2020a); Raaijmakers  _et al._ (2019) target the macroscopic
properties of NSs, such as their mass $M$, radius $R$, and dimensionless tidal
deformability $\Lambda$, which in turn can be used to constrain the EoS at
densities greater than nuclear saturation Lattimer and Prakash (2001);
Lindblom (1992, 2010).

A set of observations of different systems can be used to constrain a shared
underlying property through a hierarchical inference scheme. The hierarchical
formalism is derived in the context of combining data from different sources
while faithfully incorporating their uncertainties and potential observational
selection effects Loredo (2004); see, e.g., Refs. Chatziioannou (2020); Landry
_et al._ (2020). In the context of NS structure, the main objective is to
obtain a posterior for the EoS as a shared variable among many astrophysical
observations. The prior corresponding to this posterior is not necessarily
straightforward to define because the space of potential EoSs (i.e., the space
of possible functions obeying basic physical constraints) that relate the
pressure $p$ and the baryon density $\rho$, $p=p(\rho)$, is infinite
dimensional.111We can equivalently use $p(\varepsilon)$, with $\varepsilon$
the internal energy density. In the zero-temperature limit
${d\varepsilon}/{d\rho}=\left(p(\rho)+\varepsilon(\rho)\right)/{\rho}$.

The simplest way to define such a prior is through a _parametrization_ of the
EoS, which is a functional form of $p(\rho)$ that typically depends on a few
parameters. Common phenomenological models such as piecewise-polytrope Read
_et al._ (2009a), spectral Lindblom (2010); Lindblom and Indik (2014) and
speed-of-sound Greif  _et al._ (2019) parametrizations have been used to
effectively sample candidate EoSs for use in inference. The simplicity of a
closed-form parametric expression comes at the cost, though, of being unable
to faithfully represent many of the possible degrees of freedom in the true
EoS. While many of these models can accurately represent most EoSs derived
from effective nuclear interactions Lindblom (2010); Read  _et al._ (2009a),
it is not always clear how to extend these parametrizations toward more
general behavior in the EoS that may arise from phase transitions or new
physics. This limitation of phenomenological parametric models has been
recognized from the outset Read  _et al._ (2009a). However, in this study we
investigate another way that they may artificially restrict the inferred EoS.

Parametric models use only a few parameters, which means that the values of
$p(\rho)$ at different densities are often _correlated_. These correlations
represent a source of _model dependence_ in the inference, which is
undesirable insofar as it does not reflect true prior knowledge of the EoS at
those densities. That is, the correlations induced by the choice of
parametrization can constitute strong, unintentional prior beliefs about the
EoS. This unwanted model dependence is a natural consequence of the
phenomenological nature of the parametric models.

An alternative method for constructing a prior on the space of EoSs, which we
call _nonparametric_ in what follows, targets more model flexibility by making
use of Gaussian processes (GPs) Landry and Essick (2019). This approach
produces a multivariate Gaussian distribution for the function
$\phi=\log\left((c/c_{s})^{2}-1\right)$, where $c_{s}$ is the speed of sound
and $c$ is the speed of light. By conditioning the prior only weakly on
existing nuclear-theory models, we generate a _model-agnostic_ prior process
for EoSs. The chosen correlations between $\phi(p_{i})$ and $\phi(p_{j})$, or
equivalently between the values of the EoS at different densities, are set by
a _kernel_ function, which is in turn described by a few parameters. Following
Essick  _et al._ (2020a), we consider a variety of possible kernel parameters
to probe a range of different correlations and thus maximize model freedom.
This approach allows us, in principle, to model any function $p(\rho)$, and
furthermore to probe a wide range of inter-density correlations and high-
density EoS behavior.

Of course, completely unrestricted freedom in the EoS is neither desirable nor
realistic, as certain physical constraints should be encoded into the EoS
prior. For example, an EoS must be causal,

$\frac{dp}{d\varepsilon}=c_{s}^{2}<c^{2},$ (1)

and thermodynamically stable,

$\frac{dp}{d\varepsilon}=c_{s}^{2}>0.$ (2)

Imposing these constraints in the prior is desirable as it excludes unphysical
models from the analysis.222Some analyses allow the EoS to be slightly acausal
at times; see Appendix B for more discussion.

In this paper, we examine common parametric and nonparametric EoS models to
determine the extent to which each prior’s assumptions impact inference of the
EoS and NS properties. We find that the three parametric models we study
(spectral, piecewise-polytrope, and speed-of-sound) build additional inter-
density correlations into the EoS beyond what can be attributed to causality
and stability. These correlations between densities typically lead to more
stringent constraints than are strictly supported by the data. On the other
hand, the nonparametric model demonstrates the largest degree of model
independence, restricted primarily only by causality and thermodynamic
stability. We demonstrate that these strong, model-dependent inter-density
correlations have already impacted inferred microscopic and macroscopic NS
properties. Such effects are expected to become more severe as statistical
uncertainties decrease with more data that probe different NS densities.

The remainder of the paper is organized as follows. In Sec. II, we describe
our inference methods and our approach to investigating model dependence. In
Sec. III, we examine the EoS and NS properties inferred with current data and
show that they are influenced by correlations in the EoS prior. In Sec. IV, we
illustrate the main limitations of parametric EoS inference with a toy model.
In Sec. V, we quantify the implicit EoS correlations and demonstrate that the
nonparametric model displays the largest degree of model independence. In Sec.
VI, we study the correlations’ potential impact on upcoming EoS inference
using mock astrophysical measurements. In Sec. VII, we demonstrate that the
limitations identified in parametric models cannot be resolved by making small
modifications to the prior distributions. Finally, in Sec. VIII we discuss our
conclusions.

##  II Methods and Models

The posterior for the EoS depends on two elements: (i) the prior EoS process,
and (ii) the data. Our goal in this study is to assess the effect of the prior
as generated from different parametric and nonparametric models for the EoS.
We therefore always employ the same data, which we briefly describe in Sec.
II.2. The hierarchical likelihood corresponding to this data is described in
Refs. Landry  _et al._ (2020); Legred  _et al._ (2021). The EoS priors are
described in detail in Sec. II.1, where we discuss the different EoS models
and parameter priors that generate each EoS prior process.

###  II.1 EoS Prior

We wish to establish a prior process over candidate EoSs. By this we mean a
probabilistic measure on the space of potential EoSs. To do this, we use
several _models_ of the EoS. We distinguish _parametric_ models, which provide
a functional form for the EoS, from _nonparametric models_ which do not impose
such a functional form. We use three different phenomenological parametric
models, a _piecewise-polytrope_ Read  _et al._ (2009a) parametrization, a
_spectral_ parametrization Lindblom (2010), and a direct parametrization of
the _speed of sound_ Greif  _et al._ (2019). The spectral and piecewise-
polytrope parametrizations use a polytropic form for the EoS, so that

$p(\rho)=K\rho^{\Gamma}.$

In the piecewise-polytrope case, the polytropic index $\Gamma$ is a piecewise-
constant function of the pressure while in the spectral case $\log(\Gamma)$ is
expanded as a polynomial in pressure. In the speed of sound parametrization,
the speed of sound is expressed as a constant plus a Gaussian and a logistic
curve which asymptotes to $c^{2}/3$. Following past practice Carney  _et al._
(2018), we slightly relax the causality threshold and consider EoSs with
$c_{s}^{2}<1.1c^{2}$ for all parametric models. See Appendix B for more
details about each model and its implementation.

To establish a prior process, we must additionally supply a joint prior
probability distribution on the parameters of each model from which a _draw_
is a realization of the parameters and therefore a candidate EoS. For example,
in the spectral model, the parameters are coefficients in the spectral
expansion. In the piecewise-polytrope, the parameters are the value of the
polytropic index itself. For our headline results, we use standard priors for
the parametric models Carney  _et al._ (2018); Wysocki  _et al._ (2020),
except for the speed-of-sound model, which we adapt to increase access to
astrophysically relevant EoS; again, see Appendix B for details.

We compare the prior processes generated by the parametric models to a prior
process from a _nonparametric_ model Landry and Essick (2019). While our
nonparametric implementation does not assume a specific functional form for
the EoS, it does parametrize the correlations between the sound speed at
different densities. These correlations are described by a kernel function. In
practice, we choose a large set of points, $p_{i}$, and then the variable
$\phi(p_{i})$ is sampled from a multivariate Gaussian distribution. By
changing the kernel’s parameters and conditioning on different nuclear models,
we can generate a range of GPs. We choose a model-agnostic prior, which is to
say we average over multiple GPs with different correlations, each loosely
informed by nuclear-theory models Landry and Essick (2019). We do this to
maximize the freedom of the model. See Appendix A for more details.

Due to its construction, the GP itself has parameters which control
correlations. Such parameters have been termed _hyperparameters_ Landry and
Essick (2019), though we avoid this terminology here in order to avoid
potential confusion with the term’s use in hierarchical inference. In
addition, the parametric models also have parameters which control the prior
process; in general such details are unique to each model. We instead focus
primarily on the prior process induced by each EoS model with its chosen
prior, returning to the subject of parameter distributions briefly in Sec.
VII. For now, we simply note that our model implementations are typical of
those used in the literature Read  _et al._ (2009b); Lindblom (2010); Carney
_et al._ (2018); Lackey and Wade (2015); Greif  _et al._ (2019). Lastly, we
stress that our distinction between parametric and nonparametric models lies
not in the existence of parameters but in the specification of a functional
form for the EoS. In particular, we compare models with small, fixed numbers
of parameters that are commonly used in the literature. For these models, the
choice of functional form significantly impacts the range of EoSs that can be
represented.

###  II.2 Data and Likelihood

Unless otherwise stated, all analyses in this paper make use of the same
astronomical data as Ref. Legred  _et al._ (2021). Specifically, we include
two mass-tidal deformability measurements from gravitational wave (GW)
detections of merging NSs Abbott  _et al._ (2018, 2020b), one heavy pulsar
mass measurement with radio data Antoniadis  _et al._ (2013), and two X-ray
observations of NS masses and radii Miller  _et al._ (2019); Riley  _et al._
(2019); Miller  _et al._ (2021); Riley  _et al._ (2021). In the latter case,
we also use the up-to-date radio mass measurement of the pulsar J0740+6620
Fonseca  _et al._ (2021). Given this data $d$, the posterior probability
density of a particular EoS $\varepsilon$ is

$P(\varepsilon|d,\mathcal{I})=\frac{P(d|\varepsilon,\mathcal{I})}{P(d|\mathcal{I})}P(\varepsilon|\mathcal{I}),$
(3)

where $\mathcal{I}$ is any additional information we may have about the
system, e.g., knowledge that the data originates from a NS, as is the case for
the pulsar observations but not for the GWs. Here
$P(\varepsilon|d,\mathcal{I})$ is the _posterior_ probability of the EoS given
the data, $P(d|\varepsilon,\mathcal{I})$ is the _likelihood_ of the
astrophysical data given the EoS, $P(\varepsilon|\mathcal{I})$ is the prior
probability of the EoS, and $P(d|\mathcal{I})=\int
P(d|\varepsilon,\mathcal{I})P(\varepsilon|\mathcal{I})\mathcal{D}\varepsilon$
is the total probability of observing this data marginalized over all EoS in
the prior, often called the _evidence_. For general astrophysical data,
$P(d|\varepsilon,\mathcal{I})$ must be computed by marginalizing over the
astrophysical distribution of masses, spins, sky locations, and distances for
individual events, which remains poorly constrained Alsing  _et al._ (2018);
Farr and Chatziioannou (2020); Chatziioannou and Farr (2020); Fishbach  _et
al._ (2020); Abbott  _et al._ (2021); Landry and Read (2021); Farah  _et al._
(2021). For the full expression, see Refs. Landry  _et al._ (2020);
Chatziioannou (2020).

The different data sets we use primarily inform the EoS at different
densities. The heaviest pulsar mass measurements serve to downweight EoSs
which cannot support the observed NS masses; these constraints tend to most
significantly impact inference near $\sim(4$-$6)\rho_{\mathrm{nuc}}$ and
typically favor a stiffer EoS. The X-ray data provide constraints on the NS
radius, and constraints so far have given information about the EoS mainly in
the region $\sim(1$-$4)\rho_{\mathrm{nuc}}$ Miller  _et al._ (2021); Landry
_et al._ (2020); Legred  _et al._ (2021); Raaijmakers  _et al._ (2021). The GW
observations provide constraints on the tidal deformabilities of the binary
components, which are dominated by the loudest event observed so far, GW170817
Abbott  _et al._ (2017); Landry  _et al._ (2020). In terms of densities, the
relevant scale constrained by this measurement is
$\sim(1$-$3)\rho_{\mathrm{nuc}}$ Landry  _et al._ (2020); Legred  _et al._
(2021). Future constraints with GWs are likely to lie in this density range,
as the fractional uncertainty in $\Lambda$ will be smallest for lower-mass NSs
with less dense cores and larger tidal deformabilities. In principle, nuclear
experiments or calculations could also be included in such an analysis, and
would mainly constrain the EoS near or below $\rho_{\mathrm{nuc}}$ Essick  _et
al._ (2021a, b); Pang  _et al._ (2021); Biswas (2021); Raaijmakers  _et al._
(2021). However, we do not incorporate any in this work.

##  III Impact of EoS model on current EoS constraints

Figure 1:  Symmetric 90% credible region for the pressure $p$ at each density
$\rho$ in units of the nuclear saturation density using the nonparametric and
spectral prior processes. We show results including all astrophysical data
(labeled “astro”, solid lines) and restricting to the heavy pulsars only
(labeled “psr”, dashed lines). The latter choice ensures that prior choices on
the $M_{\rm max}$ supported by each model are irrelevant. Other parametric
models are shown in Fig. 9. In all cases, we find that the $p$-$\rho$
posterior depends on the EoS model even when identical data and inference
schemes are employed.

Figure 2:  Prior (dashed) and posterior (solid) for the radius of a
1.4$M_{\odot}$ NS, $R_{1.4}$, and maximum mass, $M_{\rm max}$, of a NS for two
choices of the marginal $M_{\rm max}$ prior: default (left) and flat (right).
We show results with the nonparametric and spectral EoS models and contours
denote 90% credible regions. The black line in the 2-dimensional plot
represents a maximally-stiff $M$-$R$ curve Rhoades and Ruffini (1974) stitched
to a fiducial low density EoS Kalogera and Baym (1996). Due to the low-density
stitching, this “causality” line should be interpreted as a fuzzy boundary and
not a sharp line. Both panels demonstrate that the spectral and nonparametric
EoS models produce different $M_{\rm max}$ posteriors and that these
differences cannot be attributed to the marginal priors. They are instead
caused by correlations between low and high densities (equivalently, between
$M_{\rm max}$ and $R_{1.4}$) imposed by the models. The correlations in the
nonparametric case are due to causality, while the spectral case exhibits
additional correlations and model dependence.

Following the above prescription, we analyze the existing data using the four
different EoS priors and plot the resulting marginal posteriors for $p(\rho)$
across a wide range of densities. Figure 1 compares the spectral and
nonparametric models; similar plots for the other parametric models can be
found in Appendix C. The posteriors differ in their predictions for the EoS.
For instance, the spectral posterior is stiffer on average than the
nonparametric one, especially above $4\rho_{\mathrm{nuc}}$. Similar
differences have also been pointed out in Refs. Greif  _et al._ (2019); Miller
_et al._ (2021); Raaijmakers  _et al._ (2021), where multiple EoS models were
employed under identical analysis settings. Our goal here is to understand the
origin of these discrepancies.

Since it is difficult to glean information about inter-density correlations
from envelope plots like Fig. 1, we turn our attention to two macroscopic NS
properties that roughly correspond to the EoS behavior at high and low
densities: the maximum mass, $M_{\rm max}$, and the radius of a 1.4$M_{\odot}$
NS, $R_{1.4}$. In Fig. 2 (left panel), we plot the one- and two-dimensional
marginal prior and posterior for $M_{\rm max}$ and $R_{1.4}$. As expected, the
marginal posteriors differ, but so do the marginal priors. Indeed, the $M_{\rm
max}$ plot shows that both the spectral prior and posterior seem to have more
support for $M_{\rm max}$ around $2.2-2.5M_{\odot}$ than the nonparametric
case. However, this trend is reversed above $2.6\,M_{\odot}$. This observation
suggests that the difference between the nonparametric and spectral posteriors
cannot be trivially assigned to different marginal priors. To further
demonstrate this, in the right panel we plot the same variables, but now
reweighted to a flat marginal $M_{\rm max}$ prior. As expected, the two
posteriors differ, with the spectral model producing a narrower posterior.

To understand this discrepancy, we revisit the possible reasons $M_{\rm max}$
is constrained on the high side. Though an upper limit on $M_{\rm max}$ has
been proposed based on the analysis of the counterpart of GW170817 Margalit
and Metzger (2017), our analysis does not make use of it. Excluding an origin
due to data, the upper limit on $M_{\rm max}$ must be the result of the EoS
prior. The two-dimensional $M_{\rm max}$-$R_{1.4}$ panel indeed shows that the
upper limit on $M_{\rm max}$ is related to the upper limit on $R_{1.4}$ Abbott
_et al._ (2017, 2018): the $M_{\rm max}$-$R_{1.4}$ prior does not cover the
entire available region for either model, with larger $M_{\rm max}$ requiring
stiffer EoSs and larger $R_{1.4}$.

The fact that larger $M_{\rm max}$ requires large values of $R_{1.4}$ is not
unexpected from causality considerations. Indeed, the causality condition and
the pressure at twice saturation, $p_{2.0}$, set an upper limit on the value
of pressure at $5$ times saturation, $p_{5.0}$. Since $p_{2.0}$ and $p_{5.0}$
correlate with $R_{1.4}$ and $M_{\rm max}$ respectively Lattimer and Prakash
(2001), any causal EoS model should limit $M_{\rm max}$ for certain low
$R_{1.4}$ configurations Rhoades and Ruffini (1974). To quantify this, we
overplot the limiting $M_{\rm max}$-$R_{1.4}$ relation given by Ref. Kalogera
and Baym (1996): each point on the line represents a soft low-density EoS
stitched to an EoS with $c_{s}^{2}=1$ at different densities. This curve
should be interpreted approximately, as the exact causality threshold depends
on the details of the low-density EoS Drischler  _et al._ (2021a).

Nonetheless, the right panel shows that in the nonparametric case the prior
fills more of the physically allowable $M_{\rm max}$-$R_{1.4}$ parameter space
compared to the spectral prior. This indicates that the nonparametric prior
has non-negligible support in the entire physically allowed region even if
specific marginal priors might downweight some regions (left panel). The same
is not true for the spectral model which cannot access certain regions of the
$M_{\rm max}$-$R_{1.4}$ plane.333Figure 2 shows 90% contours. If we plotted
99% contours instead, the nonparametric model accesses even more of the
allowed space, while the spectral model remains restricted.  This demonstrates
that the correlation between $M_{\rm max}$-$R_{1.4}$ that appears in the
spectral model is not entirely due to causality considerations, but it is also
affected by the specifics of the model. The nonparametric model, however, is
able to produce EoSs which fall near the causality limit, indicating physics
rather than modeling artifacts are the primary limitation to model freedom.
Figure 12 presents a qualitatively similar conclusion for the piecewise-
polytrope and speed-of-sound models. The piecewise-polytrope exhibits similar
behavior to the spectral model, while the speed-of-sound parametrization
exhibits the opposite problem: its prior does not include low $M_{\rm max}$
for large $R_{1.4}$ values.

Crucially, the correlations we see in these corner plots are only those that
are apparent from the two-dimensional marginalized posteriors. They do not
reveal the many hidden correlations within the parametric EoS models that are
not as easily detected. It is possible for implicit correlations within the
EoS prior to bias the inference in ways that are not obvious in low-
dimensional projections.

##  IV Impact of inter-density correlations: Toy Model

To better understand the effect of such implicit correlations, we first
consider a simple toy model that demonstrates several of the issues with
parametric models and introduce our techniques for diagnosing them.

We consider several simple linear parametrizations of the pressure as a
function of energy density $p(\varepsilon)$. This allows us to examine the
prior processes induced by the assumption of linearity from various
perspectives. We contrast this to a GP prior process in the same context,
finding particularly striking differences in the effect of a precise
measurement of the pressure at one density on our uncertainty in the pressure
at other densities.

We begin with the simple parametric model of a linear relationship between the
pressure and the energy density

$p(\varepsilon)=p_{a}+c_{s}^{2}(\varepsilon-\varepsilon_{a}),$ (4)

parametrized by the pressure at $p_{a}=p(\varepsilon_{a})$ and the slope
$c_{s}^{2}$. Ignoring causality constraints, we choose what appear to be
uninformative priors

$\displaystyle p_{a}\sim\mathcal{N}(\mu_{a},\sigma_{a}^{2}),\quad
c_{s}^{2}\sim\mathcal{N}(\mu_{c_{s}^{2}},\sigma_{c_{s}^{2}}^{2}),$ (5)

and refer to this as the _point+slope_ process. The top panel of Fig. 3 shows
the envelope plot for this prior process, i.e., the marginal distributions of
the pressure at each energy density.

Figure 3:  68% ($\pm$ 1-$\sigma$) marginal credible regions for the pressure
at each density under the _point+slope_ (top, green), _two-point_ (middle,
red), and GP (bottom, blue) prior processes. Shaded regions correspond to
marginal distributions induced by the prior process at each energy density
(light colors) and conditioned distributions for $p(\varepsilon)$ given a
precise observation of $p_{c}$ (dark colors). Only the GP process “fills the
prior volume” rapidly as one moves away from the observation point
$\varepsilon_{c}$. Compare to Fig. 5.

The envelope plot appears reasonable. That is, the prior process assigns
approximately equal uncertainty to each pressure. However, the envelope plot
only shows the marginal distributions at each density. Fig. 4 shows the
correlations between pressures at different densities. From this we see that
the prior process actually imposes strong correlations between all pressures.
We quantify this correlation between $p_{a}$ and the pressure at some other
density $p_{b}\equiv p(\epsilon_{b})$ with the _mutual information_ Cover and
Thomas (2005), defined as

$I(p_{a},p_{b})\equiv\int
dp_{a}dp_{b}\,P(p_{a},p_{b})\ln\left(\frac{P(p_{a},p_{b})}{P(p_{a})P(p_{b})}\right),$
(6)

where $P(p_{a})=\int dp_{b}\,P(p_{a},p_{b})$ is the marginal distribution. For
the _point+slope_ parametrization, we compute

$I(p_{a},p_{b})=\frac{1}{2}\ln\left(1+\frac{\sigma_{a}^{2}}{\sigma_{c_{s}^{2}}^{2}(\varepsilon_{a}-\varepsilon_{b})^{2}}\right).$
(7)

The mutual information between $p_{a}$ and $p_{b}$ can be made arbitrarily
small only in the limit
$\sigma_{a}\ll\sigma_{c_{s}^{2}}|\varepsilon_{b}-\varepsilon_{a}|$; however,
this limit corresponds to vanishingly small marginal uncertainty for $p_{a}$.
We conclude that the assumption of a linear functional form can produce what
seems to be a reasonable envelope plot in Fig. 3, but nevertheless induces
model-dependent correlations between the pressure at different densities in
Fig. 4.

Figure 4:  Joint and marginal distributions for the pressures at the three
reference densities called out in Fig. 3 for the _point+slope_ (green), _two-
point_ (red), and GP (blue) prior processes. Contours in the joint
distribution represent 90% credible regions. While certain choices for the
parametrization and priors can uncorrelate pairs of variables, only the GP
prior process induces minimal correlations between all variables
simultaneously.

In an attempt to remove the correlation between $p_{a}$ and $p_{b}$, we
consider the alternative parametrization

$p(\varepsilon)=p_{a}+\frac{p_{b}-p_{a}}{\varepsilon_{b}-\varepsilon_{a}}(\varepsilon-\varepsilon_{a}),$
(8)

described by the pressures at the two reference densities. We assume priors

$\displaystyle p_{a}\sim\mathcal{N}(\mu_{a},\sigma_{a}^{2}),\quad
p_{b}\sim\mathcal{N}(\mu_{b},\sigma_{b}^{2}),$ (9)

and refer to this as the _two-point_ process. Figures 3 and 4 show envelope
and marginal distributions, respectively.

By construction, $I(p_{a},p_{b})=0$ for the two-point prior process, as also
seen in Fig. 4. However, the envelope plot shows that the marginal prior
actually tightens for pressures between the reference densities. That is, we
are able to remove the correlation between two pressures only at the expense
of asserting greater prior knowledge about other pressures. Additionally, Fig.
4 shows that there are still correlations between ($p_{a}$, $p_{b}$) and other
pressures. This hints at the fact that, when one assumes a specific functional
form, it may be possible to remove the correlations between a small number of
statistics, but it is generally difficult to make all correlations vanish
simultaneously or to avoid making strong assumptions about specific values of
the function.

In order to consider this effect more quantitatively, we introduce a
generalization of the mutual information that considers three pressures Cover
and Thomas (2005)

$\displaystyle I(p_{a},p_{b},p_{c})$ $\displaystyle\equiv\int
dp_{a}dp_{b}dp_{c}\,P(p_{a},p_{b},p_{c})\ln\left(\frac{P(p_{a},p_{b},p_{c})}{P(p_{a})P(p_{b})P(p_{c})}\right)$
$\displaystyle=\int dp_{a}dp_{b}\,P(p_{a},p_{b})\int
dp_{c}\,P(p_{c}|p_{a},p_{b})\ln\left(\frac{P(p_{c}|p_{a},p_{b})}{P(p_{c})}\right)$
$\displaystyle\quad\quad+I(p_{a},p_{b}).$ (10)

Even if one can choose parametrizations and priors such that $I(p_{a},p_{b})$
vanishes, there is another term when considering mutual information for three
pressures. In fact, for both the _point+slope_ and _two-point_ prior
processes, the integral over $p_{c}$ diverges as $P(p_{c}|p_{a},p_{b})$ is a
delta function (determined by the closed-form parametrization), while
$P(p_{c})$ is a Gaussian with finite width. We conclude that the assumption of
a linear relationship between the pressure and the density implies an infinite
amount of information about the allowed relationships between variables. One
cannot undo all these correlations at the same time by a clever choice of
marginal prior distributions, although it may be possible to undo some of
them. The failure of this reparametrization scheme anticipates the results of
our investigation of alternative parametric priors in Sec. VII.

In general, the only way to undo all correlations simultaneously is to add
more model freedom into the prior process. For example, one may add more
reference densities to an existing model and generate a piece-wise linear
prior process. However, there will always be some densities between the
(finite number of) reference densities, regardless of how many reference
densities are chosen. In each of those regions, the piecewise-linear model is
equivalent to our _two-point_ prior process, and the strong correlations
remain. One is then left with the question of how to extend the
parametrization to remove all correlations in a scalable way. We show that a
GP is a natural solution.

A GP, defined in terms of a mean function and a covariance kernel, describes
our uncertainty in the infinitely-many degrees of freedom in a function. With
the assumption of Gaussianity, we can easily marginalize away uninteresting
degrees of freedom, in our case retaining only the pressures on a dense grid
of energy densities, and the GP reduces to a high-dimensional multivariate
Gaussian distribution. Specifically, we consider the joint distribution
induced over $p_{a}$, $p_{b}$, and $p_{c}$ by a GP

$\vec{p}\sim\mathcal{N}(\vec{\mu},\Sigma),$ (11)

with mean $\vec{\mu}$ and covariance $\Sigma$, with matrix elements defined by
a covariance kernel

$\Sigma_{ij}=\mathrm{Cov}(p_{i},p_{j})=K(\varepsilon_{i},\varepsilon_{j}).$
(12)

A common choice is the squared-exponential kernel

$K_{\mathrm{se}}(\varepsilon_{i},\varepsilon_{j})=\sigma^{2}\exp\left(-\frac{(\varepsilon_{i}-\varepsilon_{j})^{2}}{l^{2}}\right),$
(13)

although more complicated kernels are also used Essick  _et al._ (2020a).444It
is worth noting that linear regression is a special case of a GP. That is, a
GP can reproduce the linear model with an appropriate choice of covariance
kernel:
$K(\varepsilon_{i},\varepsilon_{j})\propto\varepsilon_{i}\varepsilon_{j}$.

Figures 3 and 4 show a GP assuming a squared-exponential kernel with
parameters chosen to match the marginal distribution of the _two-point_ prior
process and $l\ll|\varepsilon_{a}-\varepsilon_{b}|$. We also obtain

$\displaystyle I(p_{a},p_{b})$
$\displaystyle=-\frac{1}{2}\ln\left(1-\frac{\sigma_{ab}^{4}}{\sigma_{aa}^{2}\sigma_{bb}^{2}}\right)$
$\displaystyle=-\frac{1}{2}\ln\left\\{1-\exp\left[-\frac{2(\varepsilon_{a}-\varepsilon_{b})^{2}}{l^{2}}\right]\right\\},$
(14)

which vanishes as $\exp[-2(\varepsilon_{a}-\varepsilon_{b})^{2}/l^{2}]$ in the
limit $|\varepsilon_{a}-\varepsilon_{b}|\gg l$. Furthermore,
$P(p_{c}|p_{a},p_{b})$ is a normal distribution, and the generalization of the
mutual information in Eq. (10) no longer diverges. If
$l\ll|\varepsilon_{b}-\varepsilon_{c}|,|\varepsilon_{a}-\varepsilon_{c}|$,
then $\Sigma_{ij}\rightarrow\sigma^{2}\delta_{ij}$ and
$P(p_{c}|p_{a},p_{b})\rightarrow P(p_{c})\ \forall\ (p_{a},p_{b})$. Therefore,
if $l\ll|\varepsilon_{b}-\varepsilon_{a}|$ as well,
$I(p_{a},p_{b},p_{c})\rightarrow 0$. We conclude, then, that the GP prior
process can be made to simultaneously produce reasonable envelope plots (broad
marginal distributions for all pressures) while retaining vanishingly small
correlations between (reasonably separated) pressures. This is in stark
contrast to the parametrized prior processes, where this is, in general, not
possible.

We demonstrate one more useful diagnostic in this toy model through the
conditioned distribution

$P(p_{i}|p_{j})=\frac{P(p_{i},p_{j})}{P(p_{j})}$ (15)

which shows how our knowledge of $p_{i}$ depends on $p_{j}$. Figure 3 shows
the envelope plots for the conditioned distributions corresponding to each of
our prior processes when we condition on $p_{c}$. We see that a constraint at
$\varepsilon_{c}$ is broadcast to nearby densities in all cases, but the
Gaussian process fills up the prior volume from the unconditioned marginal
distributions the fastest. This is a visual manifestation of the correlations
quantified by the mutual information. Indeed

$I(a,b)=\int daP(a)\int db\,P(b|a)\ln\frac{P(b|a)}{P(b)},$ (16)

is just the Kullback–Leibler divergence
$D_{\mathrm{KL}}\left(P(b|a)||P(b)\right)$ from the unconditioned marginal to
the conditioned marginal averaged over the possible $a\sim P(a)$.

In the case of realistic EoS inference, we have to consider even higher-
dimensional spaces. A natural measure of correlations, then, is a
generalization of the mutual information (sometimes called the total
correlation, multivariate constraint, or multiinformation Cover and Thomas
(2005))

$I(x_{1},\cdots,x_{N})\equiv-H(x_{1},\cdots,x_{N})+\sum_{i=1}^{N}H(x_{i})$
(17)

where $H(x)=-\int dxP(x)\ln P(x)$ is the entropy of the distribution $P(x)$.
Larger $H$ imply broader distributions. We will consider this statistic in the
context of real astrophysical constraints on the EoS in Sec. V. In general,
one can make $I$ small but still allow for very little model freedom (small
$H$). We therefore seek prior processes with both large
$H(x_{1},\cdots,x_{N})$ and small $I(x_{1},\cdots,x_{N})$. This is sometimes
captured in the variation of information, defined as $H-I$, but we find it
more useful to consider $H$ and $I$ separately.

##  V Impact of inter-density correlations: idealized measurement

Figure 5:  Similar to Fig. 1 but with a mock constraint injected directly into
$p(\rho=2\rho_{\operatorname{nuc}})$ for each EoS prior process. The posterior
after the simulated constraint is included (“astro+mock”) is overplotted on
the posterior with all current data (“astro”), and in a darker color. The
constraint at a single density affects the parametric posteriors over a much
wider range of density scales than the nonparametric one. The inset focuses
around $\rho=2\rho_{\operatorname{nuc}}$. The two black straight lines provide
an estimate of the constraints imposed by causality ($c_{s}^{2}<1$) and
thermodynamic stability $(c_{s}^{2}>0.1)$ around
$\rho=2\rho_{\operatorname{nuc}}$, subject to the heavy pulsar measurements.
The nonparametric posterior quickly “fills” more of the physically available
region after satisfying the mock constraint, while the parametric posteriors
do not.

Our toy model illustrates the potential impact of implicit correlations within
EoS models on the results of EoS inference. In order to quantify the
sensitivity of the parametric and nonparametric models to such model-dependent
correlations between density scales, we now consider simulated NS
observations. The macroscopic NS properties that astronomical observations
target (masses, radii, and tides) are determined by a range of NS densities,
it is therefore not straightforward to disentangle the effect of the data and
the model dependence in the constraint that a single astronomical observation
imposes on the EoS. Consequently, we begin with the same setup as Sec. IV: an
idealized direct measurement of $p(\rho)$ at a single density, while keeping
in mind that a realistic astronomical measurement would correspond to a
combination of many such constraints correlated across many densities.555An
example of how one may obtain direct constraints on the pressure from nuclear
experiments is demonstrated in Refs. Essick  _et al._ (2021a, b).

We consider a tight Gaussian constraint at $p_{2.0}$ with mean of $3.20\times
10^{34}\,\mathrm{dyn}/\mathrm{cm}^{2}$ based on a candidate EoS drawn from our
GP prior that is consistent with all current parametric posteriors near
$2\rho_{\operatorname{nuc}}$. We arbitrarily choose the standard deviation,
$2.61\times 10^{33}\,\mathrm{dyn}/\mathrm{cm}^{2}$ ($\sim 8\%$ relative
uncertainty). We then plot the corresponding envelope for $p(\rho)$ with this
mock constraint and all other real astronomical data for each model in Fig. 5.
In the nonparametric case, imposing this constraint pinches the $p$-$\rho$
envelope around $2\rho_{\operatorname{nuc}}$, but the uncertainty in the
$p(\rho)$ curve is unaffected beyond $\approx\pm 0.5\rho_{\mathrm{nuc}}$. All
the parametric models, though, change across several $\rho_{\mathrm{nuc}}$,
indicating that the EoS at many scales is informed significantly by the EoS
near $2\rho_{\mathrm{nuc}}$.

In each panel, the inset zooms in around the $\rho=2\rho_{\operatorname{nuc}}$
region; to guide the eye the two black lines provide a rough estimate of the
maximally causal ($c_{s}^{2}=1$) and minimally stable $(c_{s}^{2}=0.1)$ EoS
that can support the heavy pulsar observations (see Fig. 2 of Landry  _et al._
(2020)) around $\rho=2\rho_{\operatorname{nuc}}$. The two lines were obtained
by combining $dp/d\varepsilon=c_{s}^{2}$ and the first law of thermodynamics
with the approximation that $\varepsilon_{2.0}=c^{2}\rho_{2.0}$. The
nonparametric prior process contains EoS draws that approach this limiting
behavior near the constraint. Comparing the $4$ panels, the nonparametric
model fills more of the physically available space. The parametric models, on
the other hand, are clearly subject to additional correlations between
pressures besides causality and stability.

To quantify these correlations, we follow Sec. IV and compute the total
correlation between the pressures at several reference densities. Table 1
shows the total correlation ($I$) and joint entropy ($H$) between
$\ln\,p_{1.0}$, $\ln\,p_{1.5}$, $\ln\,p_{2.0}$, $\ln\,p_{3.0}$, and
$\ln\,p_{4.0}$ induced by the posterior process conditioned on the
astrophysical data as well as the astrophysical data and the mock constraint
on $p_{2.0}$.666We estimate the entropies via Monte Carlo sums over kernel
density estimates (KDEs) of the associated distributions. As such, the actual
correlations may be smoothed by the KDE, which may act as upper limits on the
estimates of the mutual information in some cases. We consider these pressures
as the central density of $M_{\rm max}$ stars may be as low as
$4\rho_{\mathrm{nuc}}$ Legred  _et al._ (2021), and therefore we focus on
pressures that are confidently relevant for NSs. Although the precise values
of $I$ and $H$ can be difficult to interpret, we notice some trends.

Table 1:  Total correlation ($I$) and entropy ($H$) of the joint distributions over $\ln\,p_{1.0}$, $\ln\,p_{1.5}$, $\ln\,p_{2.0}$, $\ln\,p_{3.0}$, and $\ln\,p_{4.0}$ induced by several processes as well as the entropy of the marginal distribution over only $\ln\,p_{2.0}$ ($H(\ln\,p_{2.0})$). The nonparametric processes consistently have smaller $I$ and (much) larger $H$ than any parametric process, implying much more model freedom. This is the case even though the entropy of the marginal distributions for $\ln\,p_{2.0}$ can be comparable.  | $I$ | $H$ | $H(\ln p_{2.0})$
---|---|---|---
| PSR | Astro | Astro+$p_{2.0}$ | PSR | Astro | Astro+$p_{2.0}$ | PSR | Astro | Astro+$p_{2.0}$
nonparametric | 3.7 | 3.1 | 2.9 | 0.7 | -1.0 | -2.5 | 1.0 | 0.5 | -1.1
spectral | 6.6 | 5.5 | 4.7 | -4.2 | -5.5 | -7.6 | 0.5 | 0.0 | -1.1
polytrope | 5.7 | 4.6 | 3.8 | -1.6 | -3.6 | -5.7 | 0.9 | 0.2 | -1.1
speed of sound | 5.0 | 4.7 | 4.3 | -2.6 | -4.3 | -7.1 | 1.0 | 0.6 | -1.1

Overall, the nonparametric process consistently has the largest joint entropy
and smallest total correlation, as desired. The parametric processes all have
approximately equal $I$, which are larger than the nonparametric process by
$\gtrsim 1$ nat. Additionally, the change in $I$ when we additionally
condition on mock constraint on $p_{2.0}$ is much smaller for the
nonparametric than for the parametric processes. This can be interpreted as
the constraint on $p_{2.0}$ removing some correlations from the parametric
processes by approximately fixing the value of $p_{2.0}$.

What is more, the parametric processes have much smaller joint entropies than
the nonparametric process in all cases. This is a manifestation of the reduced
model-freedom in the parametric processes as the nonparametric process
explores more combinations of pressures than any parametric process. Although
not exact, the exponential of the difference in entropies is an estimate of
the ratio of the effective number of pressure combinations supported in each
distribution: the nonparametric contains between 10 - 100 times as many
possible pressure combinations as the parametric processes.

Additionally, the constraint on $p_{2.0}$ removes more entropy from each
parametric processes than from the nonparametric process. While we expect the
joint entropy to be smaller in all cases after measuring $p_{2.0}$ precisely,
the additional entropy lost in the parametric processes is associated with the
correlations between pressures. That is, knowledge of $p_{2.0}$ decreases our
uncertainty in other pressures within the parametric processes, something that
does not happen as strongly in the nonparametric process. This is apparent in
Fig. 5 as well.

As a final note, Table 1 also reports the entropy of the marginal
distributions over $\ln p_{2.0}$. We see smaller differences between these 1D
distributions, reinforcing the conclusion that the differences between the
nonparametric and parametric processes arise mainly from correlations between
multiple pressures.

##  VI Impact of inter-density correlations: Mock astrophysical observations

Different astronomical probes provide information about different density
scales, and therefore inter-density correlations are likely to matter even
more for realistic EoS inference than in the idealized case considered above.
Implicit correlations in EoS models could artificially give the appearance of
tension between observations of NSs or nuclear matter made via different
channels. This has already been shown to be relevant in comparisons of nuclear
experiments with astrophysical observations Essick  _et al._ (2021a, b). To
investigate this possibility, we now repeat the previous sections’ analysis
for a simulated set of astronomical observations.

We choose a candidate EoS with $M_{\rm max}=2.54\,\mathrm{M}_{\odot}$ and
$R_{1.4}=12.0\,\mathrm{km}$. This EoS is relatively soft at low densities and
stiff at high densities, but is consistent with the $90\%$ 1D marginal
pressure constraints for all of our models at all densities, except the speed-
of-sound parametrization above $4\rho_{\mathrm{nuc}}$. The combination of
macroscopic parameters lies outside the spectral 90% credible region in Fig.
2, motivating its use in studying how tension appears in an analysis when such
a mismatch arises.777Due to the broad prior of the nonparametric model,
finding a physically valid EoS with no support in the nonparametric
macroscopic or microscopic priors is much more challenging. We simulate 3
measurements of pulsar masses and radii with comparable uncertainty to the
recent measurement for J0740+6620 Miller  _et al._ (2021); Riley  _et al._
(2021). This observation incorporated radio data to constrain the pulsar mass
Fonseca  _et al._ (2021) and constrained the radius with x-ray data. We also
simulate 20 GW detections of binary NS mergers at A+ detector sensitivity
Abbott  _et al._ (2020c). Note, however, that we do not impose prior knowledge
of the NS nature of the components in our inference. The simulated pulsars are
drawn from a uniform-in-central-density distribution, while the simulated
binary NSs come from a uniform-in-mass distribution, under the condition that
the NS masses lie below $M_{\rm max}$.

Figure 6:  Inferred posterior for $R_{1.4}$ and $M_{\rm max}$ using the
nonparametric model and the spectral model and mock X-ray-radio (blue and
orange dashed-line, respectively) and GW (blue and orange solid-line,
respectively) observations. All posteriors also include all current
astrophysical data. The vertical and horizontal red lines shows the injected
value of $R_{1.4}$ and $M_{\rm max}$. The $M_{\rm max}$ posterior is only
weakly informed by the GW data as they typically cannot lead to a definitive
identification of a $>2M_{\odot}$ object as a NS, and it is thus similar to
that of Fig. 2.

We analyze each data set separately, folding the simulated measurements of
each type onto all current astrophysical data. We plot the inferred posteriors
in Fig. 6. We find that the nonparametric posterior for $R_{1.4}$ is centered
on the correct value ($12\,\mathrm{km}$) with either X-ray or GW data. In the
spectral case, though, while the GW measurements are consistent with the
correct value, the X-ray posterior is in tension at $90\%$ credibility.
Moreover, the GW and X-ray posteriors are less consistent with each other, an
observation that could lead to the erroneous conclusion of tension between
different EoS probes.

We can understand this as follows. We expect GW measurements of high-mass NSs
($\geq 1.7\mathrm{M}_{\odot}$) to be less informative than lower mass NSs, as
the absolute impact of tidal parameters on the signal is weaker for more
compact stars. In general, high mass NSs will most likely be indistinguishable
from black holes until the advent of next-generation detectors Chen  _et al._
(2020); Legred  _et al._ (2021). As such, high-mass systems offer little
information for either $M_{\rm max}$ or $R_{1.4}$, and thus the GW data
primarily probe only the low-mass/low-density part of the EoS. Indeed, Fig. 6
shows that each mock-GW $M_{\rm max}$ posteriors similar to the respective
posterior of Fig. 2, indicating that additional GW observations inform $M_{\rm
max}$ only weakly.

On the other hand, X-ray measurements have already proven capable of bounding
the radius of high-mass NSs Miller  _et al._ (2021); Riley  _et al._ (2021).
Additionally, X-ray detection of pulsations in a compact object proves it is a
NS, and thus its mass offers information about $M_{\rm max}$. Depending on the
mass distribution of observed events, X-ray probes could thus probe the EoS at
both low and high densities. Our mock X-ray data set contains one such NS with
mass $2.50\,M_{\odot}$. Figure 6, then, shows that when we use a parametric
model to fit all the X-ray data, biases can arise as no EoS in the prior
process can simultaneously reproduce the correct values for both $M_{\rm max}$
and $R_{1.4}$. The bias is smaller in the GW-based results as the data there
probe a narrower density range, resulting in the appearance of mild tension
between the two data sets. By extension, a newly observed GW signal for a
1.4$M_{\odot}$ NS would be in tension with the X-ray-based results, despite no
real astrophysical inconsistency. Recent concerns of tensions between PREX-II
Reed  _et al._ (2021) and astrophysical predictions may be influenced by a
similar mechanism, as noted in Refs. Essick  _et al._ (2021a, b).

##  VII Impact of parametric prior choices

These investigations show that the parametric EoS prior processes include
model-dependent inter-density correlations that influence the resulting
inference. Such prior processes are constructed based on two ingredients: (i)
a functional form for $p(\rho)$ (or an equivalent quantity) and (ii) a prior
for the parameters of the function. The former may be carefully engineered,
while the latter can be changed more easily. As in Sec. IV, it is therefore
reasonable to wonder if we can change the nature of the inter-density
correlations by a trivial change in the parameter prior, or whether the
correlations are inherent to the functional form. Below we argue for the
latter, as also demonstrated in Sec. IV.

Firstly, adding parameters does not necessarily always increase model freedom.
Adding a parameter to any model we have shown so far will require choosing a
distribution for that parameter, and a reasonable range will strongly depend
on the functional form. For the piecewise-polytrope model, this process is
somewhat easier, as additional adiabatic indices for new segments have clear
physical meaning. Therefore reasonable ranges can be chosen. For the spectral
model, with the addition of a new spectral component, there is no obvious
mapping of parameter values to physics and so tuning parameter ranges is much
harder.

Secondly, of the existing parameters in the models, we typically find that
only a few are meaningfully constrained. For example, the speed-of-sound model
has only a single gaussian bump, and thus current astrophysical data tightly
constrain this bump to be at densities low enough to produce pulsars
consistent with, e.g., Refs. Miller  _et al._ (2019); Riley  _et al._ (2019);
Landry  _et al._ (2020). This severely limits the flexibility of the model, as
the logistic term is, by itself, not strong enough to support realistic NSs.
In practice $a_{1}$ and $a_{2}$ are overconstrained in this model and $a_{4}$,
and $a_{5}$ are underconstrained. As a result, the speed-of-sound model is the
least flexible (and leads to the most stringent constraints) even though it
has the most parameters. We find similar behavior in the spectral and
piecewise-polytrope models, suggesting that the effective number of parameters
in the models is fewer than what is nominally stated.

Thirdly, due to the fine-tuning of the parametric models, attempting to
redefine priors on parameters is generally not an efficient way to expand
model freedom. As an example, we consider the spectral EoS where we find that
EoS candidates have a priori strong correlations between parameters in order
to satisfy causality and stability. These correlations were noted in Ref.
Wysocki  _et al._ (2020) and are also shown in Fig. 7. We find that the
$\gamma_{i}$ are alternately strongly correlated or anti-correlated with each
other. It is possible that other distributions, in particular distributions
that upweight EoSs further from the line of strongest correlation, reduce the
strength of inter-density correlations.

Figure 7:  Marginal 1- and 2-dimensional prior and posterior distributions for
the parameters of the spectral model, $\gamma_{i}$, as well as the maximum NS
mass, $M_{\rm max}$, and the radius, $R_{1.4}$. We show the default prior as
well as reweighted results that upweight more extreme values of $\gamma_{i}$.
In both cases, $M_{\rm max}$ and $R_{1.4}$ have very similar posteriors,
showing that the reweighting does not efficiently extend the coverage of the
prior toward the causality threshold in the $M_{\rm max}$-$R_{1.4}$ plane.
Additionally, extending the prior ranges for $\gamma_{i}$ is unlikely to
change the results as the posteriors are not limited by the prior ranges.

We test this by upweighting more extreme $\gamma_{i}$ values. The reweighting
procedure does indeed change the posterior distribution of $\gamma_{i}$
parameters, although the inferred distribution in $M_{\rm max}$-$R_{1.4}$ is
effectively unchanged. Moreover, the $M_{\rm max}$-$R_{1.4}$ prior does not
significantly extend into the previously excluded region closer to the
stability threshold. We find similar results with a different reweighting of
$\gamma_{i}$ that instead favors for central values. Figure 7 also shows that
extending the prior range on $\gamma_{i}$ will not extend the reach of the
spectral model as the parameter posteriors are not significantly affected by
the prior cutoffs. We reach similar conclusions with the piecewise-polytrope.

Overall, while it was possible to remove correlations between only a subset of
pressures within our toy models in Sec. IV, it is generally difficult to do
even that with real parametrized EoS models.

##  VIII Discussion and Conclusions

All models of the dense-matter EoS should contain some correlations between
density scales due to causality and thermodynamic stability requirements.
However, in this study we show that phenomenological parametric models such as
the spectral, piecewise-polytrope, and speed-of-sound models impose even
stronger correlations a priori. As a result, NS properties are constrained
more tightly in parametric models than in nonparametric ones in ways that are
not supported by the data. Regardless of whether these tighter constraints end
up being compatible with the true EoS, their emergence is attributable to what
are effectively model-dependent prior assumptions dictated by the
phenomenological nature of the parametrizations. Viewed in this way, they
deserve the same scrutiny as other prior choices imposed by the analyst.

The concerns about implicit correlations are alleviated by GP-based
nonparametric models that enjoy extensive model freedom, restricted only by
causality and thermodynamic stability. They allow us to generate, with no
additional modeling effort, candidate EoSs with complex phenomenology that
could be associated with, e.g., a transition to quark matter in the cores of
NSs. For example, Refs. Tan  _et al._ (2020, 2021) study EoSs with complex
speed-of-sound phenomenology, while Refs. Han and Steiner (2019);
Chatziioannou and Han (2020); Han and Prakash (2020); Drischler  _et al._
(2021b); Li  _et al._ (2021); Drischler  _et al._ (2021a) consider strong
first-order phase transitions that result in a discontinuity in the speed-of-
sound and multiple stable branches. The GP prior process is able to recreate
such behaviors generically.

The parametric models are relatively easier to implement. However this might
come at the cost of fine tuning which makes it harder to sample from the prior
as many draws are unphysical. Extension to more complex phenomenology, such as
phase transitions, is less straightforward and might need tailored parametric
models Alford  _et al._ (2013); Alford and Han (2016); Han  _et al._ (2019)
unless the parametrization supports such behavior inherently. The piecewise-
polytrope specifically, as implemented here, can lead to priors and posteriors
with “kinks” Carney  _et al._ (2018); Lackey and Wade (2015), while Refs.
Raaijmakers  _et al._ (2018); Riley  _et al._ (2018) discuss its behavior in
cases where the observed NSs do not reach high enough densities to probe all
polytropic segments. More complicated parametric models exist (see Appendix
B), but as we argue in Sec. VII, improving parametric models by adding
parameters or extending the priors ranges is not always straightforward or
efficient. However, extreme extensions to these models (for example a
$\mathcal{O}(1000)$ parameter extension to the piecewise-polytrope), could
exhibit behavior that is closer to the nonparametric results than the few-
parameter models they generalize.

In conclusion, commonly-used parametric models of the EoS are hampered by
built-in and often opaque correlations between density scales. These
correlations already affect inferences based on these models, and these
effects will only become more severe with additional astrophysical data. The
impact of the EoS model on inference acts as an additional systematic error
that must be addressed to achieve highly informative EoS constraints Dudi  _et
al._ (2018); Gamba  _et al._ (2021); Chatziioannou (2021); Pratten  _et al._
(2021); Kunert  _et al._ (2021); Essick (2021). Our work shows that the
nonparametric GP-based model addresses this EoS model systematic and restores
model freedom by forgoing the use of specific functional forms for the EoS
itself and instead parametrizing a wide range of possible correlations
directly.

###### Acknowledgements.

We thank Les Wade for useful discussions on the implementation of the spectal
model in LALSuite. R.E. thanks the Canadian Institute for Advanced Research
(CIFAR) for support. Research at Perimeter Institute is supported in part by
the Government of Canada through the Department of Innovation, Science and
Economic Development Canada and by the Province of Ontario through the
Ministry of Colleges and Universities. P.L. is supported by the Natural
Sciences and Engineering Research Council of Canada (NSERC). This research has
made use of data, software and/or web tools obtained from the Gravitational
Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO
Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration.
Virgo is funded by the French Centre National de Recherche Scientifique
(CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the
Dutch Nikhef, with contributions by Polish and Hungarian institutes. This
material is based upon work supported by NSF’s LIGO Laboratory which is a
major facility fully funded by the National Science Foundation. The authors
are grateful for computational resources provided by the LIGO Laboratory and
supported by National Science Foundation Grants PHY-0757058 and PHY-0823459.

##  Appendix A Description of the nonparametric EoS model

Figure 8:  Example EoS draws from the nonparametric prior plotted in terms of
the pressure $p$ versus the density $\rho$ (top), the mass $M$ versus the
radius $R$ (middle), and the speed of sound $c_{s}$ versus the density $\rho$
(bottom). We only draw EoS with nonnegligible contribution to the posterior.
For reference, we also plot the 90% symmetric credible intervals for the
posterior using only heavy pulsar observations and all astrophysical data.

Our GP is tailored to incorporate a variety of possible correlation lengths,
established by the form of the kernel function. Each GP draw is a realization
of a multivariate Gaussian distribution, which is loosely conditioned on
nuclear models. The GP from which EoS candidates are sampled has a covariance
which is governed by a kernel function through the parameters $\sigma$ and
$\ell$ that control the strength and length of the correlations respectively
Essick  _et al._ (2020a). The EoS prior process includes EoSs drawn from
multiple underlying GPs with different parameters: $\log\sigma\in U(1,10)$ and
$\ell\in U(0.1,0.9)$. However, our GPs’ kernels contain additional terms as
well. See Ref. Essick  _et al._ (2020a) for more details. In total we use
$\sim 2\times 10^{6}$ draws, of which $\sim 3\times 10^{5}$ contribute to the
prior nontrivially.

In Ref. Miller  _et al._ (2021), a single GP is used to generate EoS
realizations using the same method. This single GP is more tightly bound to
the mean realization and nuclear models than corresponding piecewise-polytrope
and spectral models due to the values of $\sigma=1$ and $\ell=1$ chosen. Such
values correspond to stronger correlations and over larger lengthscales than
any GP we employ here. This demonstrates that although the GP model is
flexible, it is not necessarily agnostic. This can be useful, for example, to
examine the validity of a set of related nuclear models given astrophysical
data Essick  _et al._ (2020b).

Figure 8 shows example draws from our prior process plotted on top of
posteriors for various parameters. The candidate EoSs exhibit a wide range of
behavior as is perhaps most evident in the bottom panel.

EoS prior process | Parameter | Prior
---|---|---
Spectral | $r_{0}$ | U(-4.37722, 4.91227)
$r_{1}$ | U(-1.82240, 2.06387)
$r_{2}$ | U(-0.32445, 0.36469)
$r_{3}$ | U(-0.09529, 0.11046)
Piecewise-polytrope | $\log p_{1}$ | U(33.6, 35.4)
$\Gamma_{1}$ | U(1.9, 4.5)
$\Gamma_{2}$ | U(1.1, 4.5)
$\Gamma_{3}$ | U(1.1, 4.5)
Speed-of-sound | $a_{1}$ | U(0.5, 1.5)
$a_{2}$ | U(1.3, 5)
$a_{3}$ | U(0.05, 3)
$a_{4}$ | U(1.5, 21)
$a_{5}$ | U(0.1, 1)
Table 2:  List of parameters and corresponding priors on which each parametric
EoS prior process depends.

##  Appendix B Description of the parametric EoS models

Figure 9:  Symmetric 90% credible region for the pressure $p$ at each density
$\rho$ in units of the nuclear saturation density (left) and the radius $R$ as
a function of the mass $M$. From top to bottom we show results with the
spectral, piecewise-polytrope, and speed-of-sound parametric models. At each
panel we overplot the corresponding nonparametric result for comparison. The
spectral $p$-$\rho$ panel is identical to Fig. 1, but we show it for
completeness. As in Fig. 1 we show results with all astrophysical data
(labeled “astro”, solid lines) and restricting to the heavy pulsars only
(labeled “psr”, dashed lines).

Figure 10:  The speed of sound squared as a function of baryon density in the
spectral (top), and piecewise-polytrope (middle) models, and speed-of-sound
(bottom) models, compared to our nonparametric model.

###  B.1 Piecewise-polytrope parametrization

In the piecewise-polytrope approach, consistent with Refs. Read  _et al._
(2009a); Carney  _et al._ (2018), the polytropic exponent is a piecewise
constant function, which changes value at 2 predetermined densities

$p(\rho)=\begin{cases}K_{1}\rho^{\Gamma_{1}}:\rho<\rho_{1}\\\
K_{2}\rho^{\Gamma_{2}}:\rho_{1}<\rho<\rho_{2}\\\
K_{3}\rho^{\Gamma_{3}}:\rho_{2}<\rho\end{cases}$ (18)

Here $\rho_{1}=10^{14.7}\mathrm{g}/\mathrm{cm}^{3}$ and
$\rho_{2}=10^{15}\mathrm{g}/\mathrm{cm}^{3}$ are fixed via an optimization for
a set of candidate EoSs following from nuclear models Read  _et al._ (2009a).
The parameter $K_{1}$ is chosen to give some value $p_{1}\equiv p(\rho_{1})$,
and $K_{2}$ and $K_{3}$ are then fixed by continuity. Therefore
$\\{\Gamma_{1},\Gamma_{2},\Gamma_{3},p_{1}\\}$ are the parameters in this
model. Their corresponding priors are given in Table 2. Extensions to this
model with more polytropic segments or allowing the transition densities to
vary are proposed in Refs. Steiner  _et al._ (2010, 2016); Raithel  _et al._
(2016); O’Boyle  _et al._ (2020).

When $\\{\Gamma_{1},\Gamma_{2},\Gamma_{3},p_{1}\\}$ are sampled from a uniform
distribution then the resulting total EoS will be neither necessarily causal
or stable. Therefore we have to enforce these constraints after the fact,
specifically we sample a set of parameters, compute the corresponding EoS, and
save it only if it obeys causality and stability.888In practice, we impose a
weaker causality constraint ($c_{s}\leq 1.1c$) for our parametric models.
Overall, we retain $\sim 1.6\times 10^{5}$ EoSs. We verified that this number
is enough to efficiently characterize the posterior by confirming that we get
consistent results with half as many draws. For the computation of the
$p(\rho)$, and $\varepsilon(p)$ relations, we used LALSimulation, a subsection
of LALSuite LIGO Scientific Collaboration (2018). For checks of NS properties
such as causality, we used LALInference LIGO Scientific Collaboration (2018);
Veitch  _et al._ (2015). Our priors are slightly more restrictive than those
used in Ref. Lackey and Wade (2015) due to computational problems that arise
for candidates with the highest $\Gamma_{2},\Gamma_{3}$, which tend to
represent acausal EoS candidates anyway.

###  B.2 Spectral parametrization

In the spectral approach, the polytropic exponent is expanded in a series of
basis functions. Following the conventions of Ref. Lindblom (2010) which
introduced the spectral parameterization, we take $x\equiv p/p_{0}$ where
$p_{0}$ is the smallest pressure where the spectral parametrization will be
used; the parametrization is matched to some other EoS at this density which
serves as the low-density crust Gamba  _et al._ (2020). Then we set

$p(\rho)=\rho^{\Gamma(x)}$ (19)

with

$\Gamma(x)=\sum_{i=0}^{n}\gamma_{i}\left(\log(x)\right)^{i}$ (20)

In most of the literature, and for our purposes $n$ is set to 3. Note that the
overall scaling of $p(\rho)$ is fixed by $\gamma_{0}$, and again we have 4
total parameters $\\{\gamma_{0},\gamma_{1},\gamma_{2},\gamma_{3}\\}$. In
practice sampling individual parameters is impractical because generic
combinations of parameters produce unphysical EoSs, even if the parameter
ranges are chosen carefully. Instead, following Ref. Wysocki  _et al._ (2020),
we sample in a different parameter space $r=(r_{0},r_{1},r_{2},r_{3})$ and
under an affine map construct samples in $\gamma$. The prior on $r$ is given
in Table 2. Our analysis uses a total of $\sim 1.9\times 10^{5}$ draws from
the spectral model. We again use the LALSuite components LALSimulation and
LALInference LIGO Scientific Collaboration (2018); Veitch  _et al._ (2015),
with particular spectral components implemented by Ref. Carney  _et al._
(2018). The spectral EoS is stitched to a model of the SLy EoS just below
$0.5\rho_{\mathrm{nuc}}$ Carney  _et al._ (2018) (see Fig. 1).

###  B.3 Speed of sound parametrization

In this approach, the speed of sound is parametrized as a function of energy
density. Taking $z\equiv\varepsilon/(\rho_{\operatorname{nuc}}c^{2})$, we
write

$\frac{c_{s}^{2}(z)}{c^{2}}=a_{1}e^{-\frac{1}{2}\left(z-a_{2}\right)^{2}/a_{3}^{2}}+a_{6}+\frac{\frac{1}{3}-a_{6}}{1+e^{-a_{5}\left(z-a_{4}\right)}}$
(21)

with $a_{1},a_{2},a_{3},a_{4},a_{5}$ real parameters, and $a_{6}$ fixed by
matching to a low-density crust. In Ref. Greif  _et al._ (2019), the matching
is done to a chiral effective field theory at $\sim\rho_{\mathrm{nuc}}$ with
limits based on Fermi liquid theory enforced up to a density of
$1.5\rho_{\mathrm{nuc}}$. Since we do not wish to use more nuclear theory
information for this model than others, we instead stitch to SLy at a density
of $0.6\rho_{\mathrm{nuc}}$, comparable to the stitching density of the
spectral model. Because of this, the parameter ranges in our implementation
must be adjusted to generate realistic EoS candidates. The prior on each
parameter is given in Table 2. Our analysis uses a total of $\sim 1.6\times
10^{5}$ draws from this model. A similar model based on the speed-of-sound is
presented in Ref. Tews  _et al._ (2018).

###  B.4 Causality in parametric models

Figure 11: Comparison between using strictly causal ($c_{s}^{2}<c^{2}$)
parametric EoSs with each model (gray) and the headline results allowing some
violation of the causal limit ($c_{s}^{2}<1.1c^{2}$). When restricting to only
causal EoSs, the issues of model dependence and insufficient coverage of the
physically allowed $M_{\rm max}$-$R_{1.4}$ space are more severe, especially
for the piecewise-polytrope.

Because of the relatively large uncertainties, we follow Ref. Carney  _et al._
(2018) in not excluding parametric EoSs until they have a large violation of
the speed of sound $c_{s}>1.1c$. This is the standard criteria used in
LALinference for the piecewise-polytrope and spectral models as part of
determining if an EoS is physical. For consistency we extend it to the speed-
of-sound parametrization as well. The primary motivation for this is to allow
a possibly acausal EoS to represent another, causal EoS which is not modeled
effectively by the prior on EoSs Carney  _et al._ (2018). In addition, the
LALinference implementation of the spectral and piecewise-polytrope models
enforces this criterion only up to the central density of the maximum mass NS.
In the speed-of-sound model, we require the EoS to be causal (or approximately
causal) everywhere. The nonparametric model obeys exact causality ($c_{s}\leq
c$) at all densities.

These choices were made for consistency with past work Carney  _et al._
(2018); Abbott  _et al._ (2018), but we still find that this extra model
freedom does not enable to spectral and piecewise-polytropic models to fill in
the physically available $M_{\rm max}$-$R_{1.4}$ space up to the causality
threshold. Figure 11 shows that excluding the acausal models minimally affects
the spectral and speed-of-sound results. However, the piecewise-polytrope
results are noticeably tighter and the $M_{\rm max}-R_{1.4}$ allowed parameter
space is covered significantly less. This also explains why the nominal
piecewise-polytrope prior supports larger pressures than the nonparametric
prior in Fig. 9.

##  Appendix C Further results with the parametric models

Figure 12:  Similar to Fig. 2 but for the piecewise-polytrope (top panels) and
the speed-of-sound (bottom panels) models. As with the spectral case, we find
that the parametric model lead to tighter posteriors for $M_{\rm max}$
compared to the nonparametric model. The 2-dimensional plots show that this is
again due to model-dependent correlations in $M_{\rm max}$-$R_{1.4}$.

Most results presented the main body of this paper were obtained using the
spectral EoS model. In this Appendix we present similar results with the
piecewise-polytropic and the speed-of-sound models. Figure 9 shows the
pressure-density and mass-radius posteriors for all EoS models. As expected,
we find that the posteriors differ due to the different models. Figure 10
shows the posteriors for the speed of sound as a function of the density where
again the nonparametric case results in the less constrained results as an
outcome of larger model flexibility.

Figure 12 shows the equivalent of Fig. 2 for the piecewise-polytrope and the
speed-of-sound models. We again find that both parametric models lead to
tighter constraints on $M_{\rm max}$ for high values. The 2-dimensional plots
show that model-dependent correlations between the maximum mass and the radius
(equivalently low and high densities) exclude certain regions of the $M_{\rm
max}$-$R_{1.4}$ space in the parametric marginal priors. As with the spectral
model, there is a “gap” between the piecewise-polytrope prior and the
approximate causality threshold. The corresponding plots for the speed-of-
sound parametrization show the opposite behavior: the prior reaches the
causality threshold, but it fails to produce EoSs with large $R_{1.4}$ and
small $M_{\rm max}$.

## References

  * Lattimer and Prakash (2016) J. M. Lattimer and M. Prakash, Phys. Rept. 621, 127 (2016), arXiv:1512.07820 [astro-ph.SR] .
  * Özel and Freire (2016) F. Özel and P. Freire, Ann. Rev. Astron. Astrophys. 54, 401 (2016), arXiv:1603.02698 [astro-ph.HE] .
  * Oertel  _et al._ (2017) M. Oertel, M. Hempel, T. Klahn, and S. Typel, Rev. Mod. Phys. 89, 015007 (2017), arXiv:1610.03361 [astro-ph.HE] .
  * Baym  _et al._ (2018) G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka, Rept. Prog. Phys. 81, 056902 (2018), arXiv:1707.04966 [astro-ph.HE] .
  * Machleidt and Entem (2011) R. Machleidt and D. R. Entem, Phys. Rept. 503, 1 (2011), arXiv:1105.2919 [nucl-th] .
  * Han  _et al._ (2019) S. Han, M. A. A. Mamun, S. Lalit, C. Constantinou, and M. Prakash, Phys. Rev. D 100, 103022 (2019), arXiv:1906.04095 [astro-ph.HE] .
  * Chatziioannou (2020) K. Chatziioannou, Gen. Rel. Grav. 52, 109 (2020), arXiv:2006.03168 [gr-qc] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (LIGO Scientific Collaboration, Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Cromartie  _et al._ (2019) H. T. Cromartie _et al._ , Nature Astron. 4, 72 (2019), arXiv:1904.06759 .
  * Miller  _et al._ (2019) M. C. Miller, C. Chirenti, and F. K. Lamb, (2019), arXiv:1904.08907 [astro-ph.HE] .
  * Riley  _et al._ (2019) T. E. Riley _et al._ , Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * 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 _et al._ , Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Antoniadis  _et al._ (2016) J. Antoniadis, T. M. Tauris, F. Ozel, E. Barr, D. J. Champion, and P. C. Freire, (2016), arXiv:1605.01665 [astro-ph.HE] .
  * Abbott  _et al._ (2020a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 892, L3 (2020a), arXiv:2001.01761 [astro-ph.HE] .
  * Raaijmakers  _et al._ (2019) G. Raaijmakers _et al._ , Astrophys. J. Lett. 887, L22 (2019), arXiv:1912.05703 [astro-ph.HE] .
  * Lattimer and Prakash (2001) J. M. Lattimer and M. Prakash, Astrophys. J. 550, 426 (2001), arXiv:astro-ph/0002232 .
  * Lindblom (1992) L. Lindblom, Astrophys. J. 398, 569 (1992).
  * Lindblom (2010) L. Lindblom, Phys. Rev. D 82, 103011 (2010), arXiv:1009.0738 [astro-ph.HE] .
  * Loredo (2004) T. J. Loredo, AIP Conf. Proc. 735, 195 (2004), arXiv:astro-ph/0409387 .
  * Landry  _et al._ (2020) P. Landry, R. Essick, and K. Chatziioannou, Phys. Rev. D 101, 123007 (2020), arXiv:2003.04880 [astro-ph.HE] .
  * Read  _et al._ (2009a) J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Phys. Rev. D 79, 124032 (2009a), arXiv:0812.2163 [astro-ph] .
  * Lindblom and Indik (2014) L. Lindblom and N. M. Indik, Phys. Rev. D89, 064003 (2014), [Erratum: Phys. Rev.D93,no.12,129903(2016)], arXiv:1310.0803 [astro-ph.HE] .
  * 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] .
  * Landry and Essick (2019) P. Landry and R. Essick, Phys. Rev. D 99, 084049 (2019), arXiv:1811.12529 [gr-qc] .
  * Essick  _et al._ (2020a) R. Essick, P. Landry, and D. E. Holz, Phys. Rev. D 101, 063007 (2020a), arXiv:1910.09740 [astro-ph.HE] .
  * Legred  _et al._ (2021) I. Legred, K. Chatziioannou, R. Essick, S. Han, and P. Landry, Phys. Rev. D 104, 063003 (2021).
  * Carney  _et al._ (2018) M. F. Carney, L. E. Wade, and B. S. Irwin, Phys. Rev. D98, 063004 (2018), arXiv:1805.11217 [gr-qc] .
  * Wysocki  _et al._ (2020) D. Wysocki, R. O’Shaughnessy, L. Wade, and J. Lange, (2020), arXiv:2001.01747 [gr-qc] .
  * Read  _et al._ (2009b) J. S. Read, C. Markakis, M. Shibata, K. Uryu, J. D. Creighton, _et al._ , Phys. Rev. D 79, 124033 (2009b), arXiv:0901.3258 [gr-qc] .
  * Lackey and Wade (2015) B. D. Lackey and L. Wade, Phys. Rev. D 91, 043002 (2015), arXiv:1410.8866 [gr-qc] .
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ (LIGO Scientific Collaboration, Virgo Collaboration), Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Abbott  _et al._ (2020b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 892, L3 (2020b), arXiv:2001.01761 [astro-ph.HE] .
  * Antoniadis  _et al._ (2013) J. Antoniadis, P. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, _et al._ , Science 340, 1233232 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Fonseca  _et al._ (2021) E. Fonseca _et al._ , Astrophys. J. Lett. 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Alsing  _et al._ (2018) J. Alsing, H. O. Silva, and E. Berti, Mon. Not. Roy. Astron. Soc. 478, 1377 (2018), arXiv:1709.07889 [astro-ph.HE] .
  * Farr and Chatziioannou (2020) W. M. Farr and K. Chatziioannou, Research Notes of the American Astronomical Society 4, 65 (2020), arXiv:2005.00032 [astro-ph.GA] .
  * Chatziioannou and Farr (2020) K. Chatziioannou and W. M. Farr, Phys. Rev. D 102, 064063 (2020), arXiv:2005.00482 [astro-ph.HE] .
  * Fishbach  _et al._ (2020) M. Fishbach, R. Essick, and D. E. Holz, Astrophys. J. Lett. 899, L8 (2020), arXiv:2006.13178 [astro-ph.HE] .
  * Abbott  _et al._ (2021) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021), arXiv:2111.03634 [astro-ph.HE] .
  * Landry and Read (2021) P. Landry and J. S. Read, Astrophys. J. Lett. 921, L25 (2021), arXiv:2107.04559 [astro-ph.HE] .
  * Farah  _et al._ (2021) A. M. Farah, M. Fishbach, R. Essick, D. E. Holz, and S. Galaudage, (2021), arXiv:2111.03498 [astro-ph.HE] .
  * Raaijmakers  _et al._ (2021) G. Raaijmakers, S. K. Greif, K. Hebeler, T. Hinderer, S. Nissanke, A. Schwenk, T. E. Riley, A. L. Watts, J. M. Lattimer, and W. C. G. Ho, (2021), arXiv:2105.06981 [astro-ph.HE] .
  * Essick  _et al._ (2021a) R. Essick, P. Landry, A. Schwenk, and I. Tews, Phys. Rev. C 104, 065804 (2021a).
  * Essick  _et al._ (2021b) R. Essick, I. Tews, P. Landry, and A. Schwenk, Phys. Rev. Lett. 127, 192701 (2021b), arXiv:2102.10074 [nucl-th] .
  * Pang  _et al._ (2021) P. T. H. Pang, I. Tews, M. W. Coughlin, M. Bulla, C. Van Den Broeck, and T. Dietrich, Astrophys. J. 922, 14 (2021), arXiv:2105.08688 [astro-ph.HE] .
  * Biswas (2021) B. Biswas, Astrophys. J. 921, 63 (2021), arXiv:2105.02886 [astro-ph.HE] .
  * Rhoades and Ruffini (1974) C. E. Rhoades, Jr. and R. Ruffini, Phys. Rev. Lett. 32, 324 (1974).
  * Kalogera and Baym (1996) V. Kalogera and G. Baym, Astrophys. J. Lett. 470, L61 (1996), arXiv:astro-ph/9608059 .
  * Margalit and Metzger (2017) B. Margalit and B. D. Metzger, Astrophys. J. 850, L19 (2017), arXiv:1710.05938 [astro-ph.HE] .
  * Drischler  _et al._ (2021a) C. Drischler, S. Han, and S. Reddy, (2021a), arXiv:2110.14896 [nucl-th] .
  * Cover and Thomas (2005) T. Cover and J. Thomas,  _Elements of information theory_ (Wiley, 2005).
  * Abbott  _et al._ (2020c) B. P. Abbott _et al._ , Living Reviews in Relativity 23, 3 (2020c).
  * Chen  _et al._ (2020) A. Chen, N. K. Johnson-McDaniel, T. Dietrich, and R. Dudi, Phys. Rev. D 101, 103008 (2020), arXiv:2001.11470 [astro-ph.HE] .
  * Reed  _et al._ (2021) B. T. Reed, F. J. Fattoyev, C. J. Horowitz, and J. Piekarewicz, Phys. Rev. Lett. 126, 172503 (2021), arXiv:2101.03193 [nucl-th] .
  * Tan  _et al._ (2020) H. Tan, J. Noronha-Hostler, and N. Yunes, Phys. Rev. Lett. 125, 261104 (2020), arXiv:2006.16296 [astro-ph.HE] .
  * Tan  _et al._ (2021) H. Tan, T. Dore, V. Dexheimer, J. Noronha-Hostler, and N. Yunes, (2021), arXiv:2106.03890 [astro-ph.HE] .
  * Han and Steiner (2019) S. Han and A. W. Steiner, Phys. Rev. D 99, 083014 (2019), arXiv:1810.10967 [nucl-th] .
  * Chatziioannou and Han (2020) K. Chatziioannou and S. Han, Phys. Rev. D 101, 044019 (2020), arXiv:1911.07091 [gr-qc] .
  * Han and Prakash (2020) S. Han and M. Prakash, Astrophys. J. 899, 164 (2020), arXiv:2006.02207 [astro-ph.HE] .
  * Drischler  _et al._ (2021b) C. Drischler, S. Han, J. M. Lattimer, M. Prakash, S. Reddy, and T. Zhao, Phys. Rev. C 103, 045808 (2021b), arXiv:2009.06441 [nucl-th] .
  * Li  _et al._ (2021) A. Li, Z. Miao, S. Han, and B. Zhang, Astrophys. J. 913, 27 (2021), arXiv:2103.15119 [astro-ph.HE] .
  * Alford  _et al._ (2013) M. G. Alford, S. Han, and M. Prakash, Phys. Rev. D 88, 083013 (2013), arXiv:1302.4732 [astro-ph.SR] .
  * Alford and Han (2016) M. G. Alford and S. Han, Eur. Phys. J. A 52, 62 (2016), arXiv:1508.01261 [nucl-th] .
  * Raaijmakers  _et al._ (2018) G. Raaijmakers, T. E. Riley, and A. L. Watts, Mon. Not. Roy. Astron. Soc. 478, 2177 (2018), arXiv:1804.09087 [astro-ph.HE] .
  * Riley  _et al._ (2018) T. E. Riley, G. Raaijmakers, and A. L. Watts, Mon. Not. Roy. Astron. Soc. 478, 1093 (2018), arXiv:1804.09085 [astro-ph.HE] .
  * Dudi  _et al._ (2018) R. Dudi, F. Pannarale, T. Dietrich, M. Hannam, S. Bernuzzi, F. Ohme, and B. Brügmann, Phys. Rev. D98, 084061 (2018), arXiv:1808.09749 [gr-qc] .
  * Gamba  _et al._ (2021) R. Gamba, M. Breschi, S. Bernuzzi, M. Agathos, and A. Nagar, Phys. Rev. D 103, 124015 (2021), arXiv:2009.08467 [gr-qc] .
  * Chatziioannou (2021) K. Chatziioannou, (2021), arXiv:2108.12368 [gr-qc] .
  * Pratten  _et al._ (2021) G. Pratten, P. Schmidt, and N. Williams, (2021), arXiv:2109.07566 [astro-ph.HE] .
  * Kunert  _et al._ (2021) N. Kunert, P. T. H. Pang, I. Tews, M. W. Coughlin, and T. Dietrich, (2021), arXiv:2110.11835 [astro-ph.HE] .
  * Essick (2021) R. Essick, (2021), arXiv:2111.04244 [astro-ph.HE] .
  * Essick  _et al._ (2020b) R. Essick, I. Tews, P. Landry, S. Reddy, and D. E. Holz, Phys. Rev. C 102, 055803 (2020b), arXiv:2004.07744 [astro-ph.HE] .
  * Steiner  _et al._ (2010) A. W. Steiner, J. M. Lattimer, and E. F. Brown, Astrophys. J. 722, 33 (2010), arXiv:1005.0811 [astro-ph.HE] .
  * Steiner  _et al._ (2016) A. W. Steiner, J. M. Lattimer, and E. F. Brown, Eur. Phys. J. A 52, 18 (2016), arXiv:1510.07515 [astro-ph.HE] .
  * Raithel  _et al._ (2016) C. A. Raithel, F. Ozel, and D. Psaltis, Astrophys. J. 831, 44 (2016), arXiv:1605.03591 [astro-ph.HE] .
  * O’Boyle  _et al._ (2020) M. F. O’Boyle, C. Markakis, N. Stergioulas, and J. S. Read, Phys. Rev. D 102, 083027 (2020), arXiv:2008.03342 [astro-ph.HE] .
  * LIGO Scientific Collaboration (2018) LIGO Scientific Collaboration, “LIGO Algorithm Library - LALSuite,” free software (GPL) (2018).
  * Veitch  _et al._ (2015) J. Veitch _et al._ , Phys. Rev. D 91, 042003 (2015), arXiv:1409.7215 [gr-qc] .
  * Gamba  _et al._ (2020) R. Gamba, J. S. Read, and L. E. Wade, Class. Quant. Grav. 37, 025008 (2020), arXiv:1902.04616 [gr-qc] .
  * Tews  _et al._ (2018) I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astrophys. J. 860, 149 (2018), arXiv:1801.01923 [nucl-th] .


* Wagner [2011] Wagner, W (2011), “Stochastic models in kinetic theory,” Phys. Fluids 23, 030602.
  * Wagner [2004] Wagner, Wolfgang (2004), “Stochastic models and monte carlo algorithms for boltzmann type equations,” in  _Monte Carlo and Quasi-Monte Carlo Methods 2002: Proceedings of a Conference held at the National University of Singapore, Republic of Singapore, November 25–28, 2002_ (Springer) pp. 129–153.
  * Wang  _et al._ [1997] Wang, Q, KD Squires, M Chen, and JB McLaughlin (1997), “On the role of the lift force in turbulence simulations of particle deposition,” International Journal of Multiphase Flow 23 (4), 749–763.
  * Wiegel [1986] Wiegel, Frederik W (1986),  _Introduction to path-integral methods in physics and polymer science_ (World Scientific).
  * Yeung [2002] Yeung, PK (2002), “Lagrangian investigations of turbulence,” Annual review of fluid mechanics 34 (1), 115–142.
  * Yeung and Pope [1989] Yeung, Pui-Kuen, and Stephen B Pope (1989), “Lagrangian statistics from direct numerical simulations of isotropic turbulence,” Journal of Fluid Mechanics 207, 531–586.