#  A Detailed Examination of Astrophysical Constraints on the Symmetry Energy
and the Neutron Skin of 208Pb with Minimal Modeling Assumptions

Reed Essick<EMAIL_ADDRESS>Perimeter Institute for Theoretical
Physics, 31 Caroline Street North, Waterloo, Ontario, Canada, N2L 2Y5  Kavli
Institute for Cosmological Physics, The University of Chicago, Chicago, IL
60637, USA  Philippe Landry<EMAIL_ADDRESS>Nicholas & Lee Begovich
Center for Gravitational-Wave Physics & Astronomy, California State
University, Fullerton, 800 N State College Blvd, Fullerton, CA 92831  Achim
Schwenk<EMAIL_ADDRESS>Technische Universität Darmstadt,
Department of Physics, 64289 Darmstadt, Germany  ExtreMe Matter Institute
EMMI, GSI Helmholtzzentrum für Schwerionenforschung GmbH, 64291 Darmstadt,
Germany  Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, 69117
Heidelberg, Germany  Ingo Tews<EMAIL_ADDRESS>Theoretical Division, Los
Alamos National Laboratory, Los Alamos, NM 87545, USA

###### Abstract

The symmetry energy and its density dependence are pivotal for many nuclear
physics and astrophysics applications, as they determine properties ranging
from the neutron-skin thickness of nuclei to the crust thickness and the
radius of neutron stars. Recently, PREX-II reported a value of $0.283\pm
0.071$ fm for the neutron-skin thickness of 208Pb, $R_{\rm
skin}^{{}^{208}\text{Pb}}$, implying a symmetry-energy slope parameter $L$ of
$106\pm 37$ MeV, larger than most ranges obtained from microscopic
calculations and other nuclear experiments. We use a nonparametric equation of
state representation based on Gaussian processes to constrain the symmetry
energy $S_{0}$, $L$, and $R_{\rm skin}^{{}^{208}\text{Pb}}$ directly from
observations of neutron stars with minimal modeling assumptions. The resulting
astrophysical constraints from heavy pulsar masses, LIGO/Virgo, and NICER
favor smaller values of the neutron skin and $L$, as well as negative symmetry
incompressibilities. Combining astrophysical data with chiral effective field
theory ($\chi\mathrm{EFT}$) and PREX-II constraints yields
$S_{0}=33.0^{+2.0}_{-1.8}$ MeV, $L=53^{+14}_{-15}$ MeV, and $R_{\rm
skin}^{{}^{208}\text{Pb}}=0.17^{+0.04}_{-0.04}$ fm. We also examine the
consistency of several individual $\chi\mathrm{EFT}$ calculations with
astrophysical observations and terrestrial experiments. We find that there is
only mild tension between $\chi\mathrm{EFT}$, astrophysical data, and PREX-
II’s $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ measurement ($p$-value
$=12.3\%$) and that there is excellent agreement between $\chi\mathrm{EFT}$,
astrophysical data, and other nuclear experiments.

††preprint: LA-UR-21-21363

##  I Introduction

Knowledge of the nuclear symmetry energy is vital for describing systems with
neutron-proton asymmetry, ranging from atomic nuclei to neutron stars Lattimer
and Lim (2013); Tsang  _et al._ (2012); Huth  _et al._ (2021). The symmetry
energy is defined as the difference between the nuclear energy per particle in
pure neutron matter (PNM) and symmetric nuclear matter (SNM),

$S(n)=\frac{E_{\rm PNM}}{A}(n)-\frac{E_{\rm SNM}}{A}(n)\,.$ (1)

Pure neutron matter consists only of neutrons and resembles neutron-star
matter closely, while SNM consists of equal parts of protons and neutrons and
can be probed through the bulk energy of atomic nuclei. The value of
$S_{0}=S(n_{0})$, typically defined at nuclear saturation density
$n_{0}\approx 0.16\,\mathrm{fm}^{-3}$, and the density dependence of $S(n)$,
described by its slope parameter $L$ and curvature $K_{\rm sym}$,

$\displaystyle L=3n\left.\frac{\partial S(n)}{\partial n}\right|_{n_{0}}\,,$
(2) $\displaystyle K_{\rm sym}(n)=9n^{2}\left.\frac{\partial^{2}S(n)}{\partial
n^{2}}\right|_{n_{0}}\,,$ (3)

can be correlated to several observables in nuclear physics and astrophysics,
e.g., to the neutron-skin thickness of nuclei ($R_{\mathrm{skin}}$ Typel and
Brown (2001); Viñas  _et al._ (2014); Reinhard and Nazarewicz (2016); Mondal
_et al._ (2016)), their electric dipole polarizability ($\alpha_{D}$ Tamii
_et al._ (2011); Piekarewicz  _et al._ (2012); Roca-Maza  _et al._ (2013,
2015)), the radius ($R$) of neutron stars (NSs) Lattimer and Prakash (2001);
Steiner and Gandolfi (2012), and properties of the NS crust Neill  _et al._
(2021). This is because $L$ is related to the pressure of PNM at $n_{0}$,
where $d(E_{\mathrm{SNM}}/A)/dn=0$. Typical values for $S_{0}$ and $L$ from
nuclear experiments Lattimer and Lim (2013); Tsang  _et al._ (2012); Tamii
_et al._ (2011); Roca-Maza  _et al._ (2015); Russotto  _et al._ (2016) and
theory Huth  _et al._ (2021); Tews  _et al._ (2013); Hebeler  _et al._ (2013);
Drischler  _et al._ (2016); Lonardoni  _et al._ (2020); Drischler  _et al._
(2020a) are $30$–$35$ MeV and $30$–$70$ MeV, respectively.

In particular, the neutron-skin thickness of 208Pb,
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$, is strongly correlated with $L$
Typel and Brown (2001); Viñas  _et al._ (2014); Reinhard and Nazarewicz
(2016); Mondal  _et al._ (2016). Recently, the PREX collaboration determined
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ by measuring the parity-violating
asymmetry ($\rm A_{\rm PV}$) in the elastic scattering of polarized electrons
off ${}^{208}\mathrm{Pb}$. Using data from two experimental runs, PREX-I and
PREX-II, the PREX collaboration reported
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}=0.283\pm 0.071\,\mathrm{fm}$ (mean
$\pm$ standard deviation) Adhikari  _et al._ (2021). Using a correlation
between $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ and $L$, Ref. Reed  _et al._
(2021) inferred $L=106\pm 37\,\mathrm{MeV}$ from this measurement. Note that
Ref. Reinhard  _et al._ (2021) has found lower values of
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ and $L$ when folding in information
from other nuclear observables.

In recent work Essick  _et al._ (2021), we examined astrophysical constraints
on the symmetry energy, its density dependence, and
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ using a nonparametric inference
framework for the equation of state (EOS) Landry and Essick (2019); Essick
_et al._ (2020a). This framework is based on Gaussian Processes (GPs) that
simultaneously represent the uncertainty in the (infinitely many) functional
degrees of freedom of the sound speed in $\beta$-equilibrium as a function of
pressure. This approach avoids the modeling assumptions implicit in
parametrized EOS representations—e.g., speed-of-sound Alford  _et al._ (2013);
Tews  _et al._ (2018); Greif  _et al._ (2019), polytropic Hebeler  _et al._
(2013); Raithel  _et al._ (2016), or spectral Lindblom (2010); Lindblom and
Indik (2012) extension schemes—which attempt to capture the variability in the
EOS in terms of a number of parameters. Hence, our extraction of the symmetry
energy and the neutron-skin thickness allows for increased model freedom
relative to astrophysical inferences using explicit parameterizations of the
EOS (e.g., Refs. Alam  _et al._ (2016); Carson  _et al._ (2019); Biswas  _et
al._ (2020); Yue  _et al._ (2021)). Indeed, our approach reduces systematic
uncertainties from a priori modeling assumptions, which can otherwise be
difficult to quantify, and provides constraints obtained directly from the
astrophysical data.

In this paper, we provide a more detailed description of our method and
present additional new results for symmetry-energy parameters, the neutron-
skin thickness and NS properties. In Ref. Essick  _et al._ (2021), we
marginalized over four nuclear-theory calculations of the EOS from chiral
effective field theory ($\chi$EFT). Here, we examine the results for the
individual calculations and discuss what we can learn about nuclear
interactions from comparisons with astrophysical data. In general, we find no
significant tension between the PREX-II data and astrophysical observations,
primarily because $L$ is less strongly correlated with NS observables than has
typically been claimed Lattimer and Prakash (2001); Lattimer and Lim (2013).
Given current measurement uncertainties, there is only mild tension between
PREX-II and the $\chi\mathrm{EFT}$ predictions, while the latter agree very
well with measurements of the dipole polarizability of 208Pb
($\alpha_{D}^{{}^{208}\mathrm{Pb}}$) Tamii  _et al._ (2011); Roca-Maza  _et
al._ (2013, 2015). Finally, we show that allowing for a nonparametric high-
density extension of the EOS leads to a significantly weaker correlation of
the $L$ parameter with NS radii, which must be taken into account when
discussing the impact of a precise $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$
measurement on NS radii.

This paper is structured as follows. In Sec. II, we introduce the
nonparametric EOS inference scheme. In Sec. III, we explain how we extract the
nuclear parameters from the nonparametric EOS realizations. We then present
the results of the inference of microscopic and macroscopic dense-matter
properties in Sec. IV. In particular, we address the consistency of various
$\chi$EFT predictions with astrophysical observations and experimental
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ and
$\alpha_{D}^{{}^{208}\mathrm{Pb}}$ measurements. In Sec. V, we discuss
possible future areas of improvement and their expected impact before
concluding in Sec. VI.

##  II Methodology

We briefly review our GP-based nonparametric EOS inference scheme in Sec. II.1
before summarizing the astrophysical data used in our inference in Sec. II.2.
Section II.3 describes the $\chi\mathrm{EFT}$ calculations employed in this
work, against which we contrast the constraints obtained without nuclear-
theory input at low densities.

###  II.1 Nonparametric EOS Inference

To extract dense matter information from astrophysical observations of NSs, we
need a model for the NS EOS, i.e., the relation between energy density and
pressure in the stellar interior. In this work, we use the nonparametric
representation of the EOS introduced in Refs. Landry and Essick (2019); Essick
_et al._ (2020a) based on GPs that model the uncertainty in the correlations
between the sound speed in $\beta$-equilibrium at different pressures. By
construction, the GPs generate EOS realizations that are causal,
thermodynamically stable, and matched to a NS crust model (BPS Baym  _et al._
(1971)) at very low densities, $n<0.3n_{0}$. Although GPs can be constructed
to closely emulate the behavior of specific theoretical models, we instead
construct GPs that explore as much functional behavior as possible (see the
discussion of _model-informed_ vs. _model-agnostic_ priors in Refs. Landry and
Essick (2019); Essick  _et al._ (2020a)). That is, our GPs are not strongly
informed by a specific description of the microphysics; they are designed to
be theory-agnostic.

Our GPs are conditioned on a training set of tabulated EOSs from the
literature. In particular, we follow Ref. Essick  _et al._ (2020a) and
construct priors from mixture models of GPs separately conditioned on
hadronic, hyperonic and quark EOSs. We condition 50 GPs with agnostic
hyperparameters for each composition, and then marginalize over the
compositions to obtain our final prior; see Ref. Essick  _et al._ (2020a) for
more details. In this way, our prior emulates the functional behavior of
established EOSs on average. However, each process’s uncertainties are very
large, so that the EOS realizations we generate span a much wider range of
behavior than the training set. This includes EOSs that are much stiffer or
much softer than EOSs from the literature, as well as many that exhibit sharp
features reminiscent of strong phase transitions that can give rise to
multiple stable branches in the mass-radius relation. By sampling many EOS
realizations from the GPs, one obtains a discrete prior process over the EOS.
We typically draw $10^{4}$–$10^{6}$ EOS realizations for each prior we
consider.

Given this large set of EOS realizations, our analysis proceeds through a
Monte-Carlo implementation of a hierarchical Bayesian inference. Every EOS
from the prior is assigned a marginal likelihood from each astrophysical
observation. In turn, the likelihood for each observation is modeled as an
optimized kernel density estimate (KDE), and we directly marginalize over
nuisance parameters (e.g., the masses $M$) with respect to a fixed prior (see
Ref. Landry  _et al._ (2020) for more details). This results in a
representation of the posterior EOS process as a set of discrete samples with
weights equal to the product of the marginal likelihoods. The posterior
probability for an EOS realization $\varepsilon_{\beta}$ is then

$P(\varepsilon_{\beta}|\\{d\\})\propto
P(\varepsilon_{\beta})\prod_{i}P(d_{i}|\varepsilon_{\beta})\,,$ (4)

where $\\{d\\}=\\{d_{1},d_{2},\dots\\}$ is the set of observations,
$P(d_{i}|\varepsilon_{\beta})$ are the corresponding marginal likelihoods, and
$P(\varepsilon_{\beta})$ is the EOS realization’s prior probability.

###  II.2 Astrophysical Data

The nonparametric inference scheme can incorporate different types of
astrophysical observations Landry  _et al._ (2020), including the existence of
massive pulsars Antoniadis  _et al._ (2013); Cromartie  _et al._ (2020),
simultaneous $M$-$\Lambda$ measurements from compact binary mergers with
gravitational waves (GWs) Abbott  _et al._ (2017a, 2019) observed by the
Advanced LIGO Aasi  _et al._ (2015) and Virgo Acernese  _et al._ (2015)
interferometers, and simultaneous $M$-$R$ measurements from X-ray pulse-
profile modeling of Neutron Star Interior Composition Explorer (NICER) Miller
_et al._ (2019); Riley  _et al._ (2019) observations. We use these
astrophysical observations to constrain the GPs described in the previous
section.

For the masses of the two heaviest known NSs, measured via pulsar timing, we
model the likelihoods $P(d|m)$ as Gaussian distributions. For PSR J0740+6620
Cromartie  _et al._ (2020); Fonseca  _et al._ (2021) (respectively, PSR
J0348+0432 Antoniadis  _et al._ (2013)) the mean and standard deviation are
$2.08\pm 0.07\,M_{\odot}$ ($2.01\pm 0.04\,M_{\odot}$). The likelihood of an
EOS realization $\varepsilon_{\beta}$, given this observation, is

$P(d|\varepsilon_{\beta})\propto\int P(d|M)P(M|\varepsilon_{\beta})dM\,.$ (5)

We take the mass prior $P(M|\varepsilon_{\beta})$ to be flat up to the maximum
mass supported by the EOS realization, and take care to include the proper
normalization. This ensures that EOSs that predict a maximum mass far below
the pulsar mass are assigned zero likelihood, while among EOSs that support
greater masses the models that least overestimate the maximum mass relative to
the observation are favored (see Appendix of Landry  _et al._ (2020) and
discussion in Ref. Miller  _et al._ (2019)). In practical terms, this is
because the non-observation of pulsars with masses significantly above
$2.1\,M_{\odot}$ is informative in itself.

For $M$-$\Lambda$ measurements from GW170817 Abbott  _et al._ (2017a, 2019),
we model the likelihood $P(d|M_{1},M_{2},\Lambda_{1},\Lambda_{2})$ with an
optimized Gaussian KDE as explained in Ref. Essick  _et al._ (2020a). The
corresponding likelihood of an EOS realization $\varepsilon_{\beta}$ given
this observation is

$P(d|\varepsilon_{\beta})\propto\int\Big{[}P(d|M_{1},M_{2},\Lambda_{1},\Lambda_{2})P(M_{1},M_{2})\\\
\times\delta(\Lambda_{1}-\Lambda(M_{1}))\delta(\Lambda_{2}-\Lambda(M_{2}))\Big{]}dM_{1}dM_{2}\,.$
(6)

The mass prior is taken to be uniform. We do not truncate it at the maximum
mass supported by the EOS because we do not exclude a priori the possibility
that one of the components of the binary was a BH. Our analysis does not
incorporate the binary NS observation GW190425, as it was not loud enough to
yield a measurable matter signature and hence inform inference of the EOS.
Furthermore, we do not include light-curve models of electromagnetic
counterparts associated with GW events because of the systematic uncertainties
involved in interpreting the kilonova physics and its connection to the EOS
(see, e.g., discussions in Refs. Abbott  _et al._ (2017b); Coughlin  _et al._
(2018); Kasen and Barnes (2019); Siegel (2019); Metzger (2020); Dietrich  _et
al._ (2020); Stachie  _et al._ (2021); Raaijmakers  _et al._ (2021)).

Finally, we consider X-ray pulse-profile measurements of PSR J0030+0451’s mass
and radius assuming a three-hotspot configuration Miller  _et al._ (2019) (see
also Ref. Riley  _et al._ (2019), which yields comparable results Landry  _et
al._ (2020)). The likelihood $P(d|M,R)$ for this observation is also modeled
with an optimized Gaussian KDE Essick  _et al._ (2020a). Weighing an EOS
realization $\varepsilon_{\beta}$ by this likelihood, we obtain

$P(d|\varepsilon_{\beta})\propto\int P(d|M,R)P(M|\varepsilon_{\beta})dM\,.$
(7)

The mass prior should, in principle, extend only up to the maximum mass for a
given EOS realization because, like for the pulsar mass measurements, we know
that PSR J0030+0451 is a NS. However, for convenience we instead assume a NS
population model that truncates the mass prior for X-ray sources well below
the maximum NS mass. As discussed in Ref. Landry  _et al._ (2020), these two
prescriptions are effectively equivalent in the case of PSR J0030+0451 because
its mass is clearly smaller than the maximum mass of any viable EOS.

Nonetheless, we would need to truncate $P(M|\varepsilon_{\beta})$ at
$M_{\mathrm{max}}$ if we were to include the recent NICER+XMM Newton
observations of J0740+6620 Miller  _et al._ (2021); Riley  _et al._ (2021);
Raaijmakers  _et al._ (2021). We do not consider this measurement in the
present work because the NICER results for J0740+6620 were published after
Ref. Essick  _et al._ (2021) and the properties of this high-mass NS do not
influence significantly the EOS inference at $n_{0}$ (see also Refs. Pang  _et
al._ (2021); Legred  _et al._ (2021)), especially within our nonparametric
framework (see, e.g., Fig. 12). However, the updated mass measurement for
J0740+6620 reported in Ref. Fonseca  _et al._ (2021) is incorporated as one of
the two pulsar mass observations described above.

###  II.3 Chiral EFT Calculations

The nonparametric EOS prior based on a crust EOS with GP extensions to higher
densities can also be conditioned on theoretical calculations of the EOS for
densities above the crust and up to around $1-2n_{0}$, where nuclear theory
calculations are well controlled. At higher densities, our EOS framework still
uses the model independence of the GP construction Essick  _et al._ (2020b).
Following our previous work Essick  _et al._ (2021), we separately condition
the EOS on the uncertainty band obtained from four different calculations
based on $\chi$EFT interactions and marginalize over all four bands.

First, we consider quantum Monte Carlo calculations (QMC) using local
$\chi$EFT interactions up to next-to-next-to-leading order (N2LO) Lynn  _et
al._ (2016). These results, labeled
$\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$, are based on a nonperturbative
many-body method that is proven to be accurate for strongly correlated
systems, but are presently limited to N2LO due to nonlocalities entering at
higher order in $\chi$EFT. As a result, the
$\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$ band has somewhat larger
uncertainties. In addition, we consider two calculations based on many-body
perturbation theory (MBPT) Tews  _et al._ (2013); Drischler  _et al._ (2019)
using nonlocal $\chi$EFT interactions up to next-to-next-to-next-to-leading
order (N3LO). Both calculations include all two-, three-, and four-neutron
interactions up to this order. The results from Ref. Drischler  _et al._
(2019), which we label MBPT${}^{\text{(2019)}}_{\text{N${}^{3}$LO}}$, include
contributions up to higher order in MBPT as well as EFT truncation
uncertainties (for two cutoffs: 450 and 500 MeV), while the results from Ref.
Tews  _et al._ (2013), labeled MBPT${}^{\text{(2013)}}_{\text{N${}^{3}$LO}}$,
are lower order in MBPT but include other uncertainties in two- and three-
nucleon interactions as well. Therefore, we find it useful to explore both EOS
bands here. We note that the combined 450 and 500 MeV N3LO bands from Ref.
Drischler  _et al._ (2019) overlap very closely with the recent GP uncertainty
bands (GP-B) from Ref. Drischler  _et al._ (2020a), labeled
$\mathrm{MBPT}^{(2020\mathrm{\,GP})}_{\mathrm{N^{3}LO}}$in the following (see
also Ref. Huth  _et al._ (2021)). Finally, we also consider the MBPT
calculations with two-nucleon interactions at N3LO and three-nucleon
interactions at N2LO, labeled MBPT${}^{\text{(2010)}}_{\text{mixed}}$, based
on a broader range of three-nucleon couplings Hebeler and Schwenk (2010);
Hebeler  _et al._ (2013). Exploring these four bands allows us to account for
different nuclear interactions and many-body approaches, increasing the
robustness of our results.

##  III Extraction of Nuclear Parameters from Nonparametric EOS Realizations

The nuclear EOS can be described by the nucleonic energy per particle, $E_{\rm
nuc}/A(n,x)$, which depends on the density $n$ and the proton fraction
$x=n_{p}/n$ with $n_{p}$ being the proton density. The symmetry energy $S(n)$
is encoded in the $x$ dependence of $E_{\rm nuc}/A(n,x)$. In our approach, we
approximate the $x$ dependence of the nucleonic energy per particle with the
standard quadratic expansion,

$\frac{E_{\rm nuc}}{A}(n,x)=\frac{E_{\rm SNM}}{A}+S(n)(1-2x)^{2}\,,$ (8)

where higher-order terms beyond $\mathcal{O}(x^{2})$ are expected to be small
around $n_{0}$, and can be safely neglected given current EOS uncertainties
Drischler  _et al._ (2014); Somasundaram  _et al._ (2021). For example, Ref.
Goriely  _et al._ (2013) suggested systematic shifts of
$\mathcal{O}(3\,\mathrm{MeV})$ in $L$ when higher-order terms are included in
Eq. (8) (compare $L$ and $\tilde{L}$ in Table V), but these are much smaller
than the statistical uncertainty in all our priors (Table 1). $S(n)$ can be
computed as

$\displaystyle S(n)$ $\displaystyle=\frac{E_{\rm
nuc}}{A}\left(n,0\right)-\frac{E_{\rm nuc}}{A}\left(n,1/2\right)$
$\displaystyle=\frac{E_{\rm PNM}}{A}-\frac{E_{\rm SNM}}{A}\,.$ (9)

In our nonparametric EOS inference, each EOS realization is represented in
terms of the baryon density $n$, the energy density $\varepsilon_{\beta}$, and
the pressure $p_{\beta}$ in $\beta$-equilibrium. These quantities are related
to the energy per particle $E/A$ through

$\displaystyle\varepsilon$ $\displaystyle=n\cdot\left(\frac{E}{A}+m_{\rm
N}\right)\,,$ (10) $\displaystyle p$ $\displaystyle=n^{2}\frac{\partial
E/A}{\partial n}\,,$ (11)

where $m_{N}$ is the average nucleon mass and we use units with $\hbar=c=1$.
We need to correct the total energy density in $\beta$-equilibrium for the
contribution of electrons:

$\frac{E_{\rm nuc}}{A}(n,x)=\frac{\varepsilon_{\beta}(n)-\varepsilon_{\rm
e}(n,x)}{n}-m_{\rm N}\,.$ (12)

In this work, we describe the electron contribution using the relations for a
relativistic Fermi gas Chamel and Haensel (2008):

$\varepsilon_{e}(n_{\rm e})=\frac{m_{\rm
e}^{4}}{8\pi^{2}}\left(x_{r}(2x_{r}^{2}+1)\sqrt{x_{r}^{2}+1}\right.\\\
\left.-\ln(x_{r}+\sqrt{x_{r}^{2}+1})\right)\,.$ (13)

where $n_{e}$ is the electron density, and $x_{r}=k_{F}/m_{\rm
e}=(3\pi^{2}n_{\rm e})^{1/3}/m_{e}$ with the electron mass $m_{\rm e}=0.511$
MeV. We neglect the contribution from muons because their effect on the EOS
around nuclear saturation density is small. Then, due to charge neutrality,
the electron density in $\beta$-equilibrium equals the proton density, $n_{\rm
e}=x(n)\cdot n$.

The proton fraction $x(n)$ is unknown for each EOS draw but it can be
constrained by enforcing the $\beta$-equilibrium condition,

$\mu_{\rm n}(n,x)=\mu_{\rm p}(n,x)+\mu_{\rm e}(n,x)\,,$ (14)

where $\mu_{\rm i}(n,x)$ is the chemical potential for particle species $i$.
The electron chemical potential is given by

$\mu_{\rm e}(n_{\rm e})=\sqrt{(3\pi^{2}n_{\rm e})^{2/3}+m_{\rm e}^{2}}\,,$
(15)

and the neutron and proton chemical potentials $\mu_{\rm n}$ and $\mu_{\rm p}$
in asymmetric nuclear matter are given by

$\displaystyle\mu_{\rm p}(n,x)$
$\displaystyle=\frac{d\varepsilon_{\mathrm{nuc}}}{dn_{p}}=n\frac{\partial\left(E_{\rm
nuc}/A\right)}{\partial n}+\frac{\partial\left(E_{\rm nuc}/A\right)}{\partial
x}(1-x)+\frac{E_{\rm nuc}}{A}+m_{\rm p}\,,$ (16) $\displaystyle\mu_{\rm
n}(n,x)$
$\displaystyle=\frac{d\varepsilon_{\mathrm{nuc}}}{dn_{n}}=n\frac{\partial\left(E_{\rm
nuc}/A\right)}{\partial n}-\frac{\partial\left(E_{\rm nuc}/A\right)}{\partial
x}x+\frac{E_{\rm nuc}}{A}+m_{\rm n}\,,$ (17)

with the neutron and proton masses $m_{\rm n}$ and $m_{\rm p}$, respectively.
Hence, the $\beta$-equilibrium condition is given by

$m_{\rm n}-m_{\rm p}-\frac{\partial\left(E_{\rm nuc}/A\right)}{\partial
x}-\mu_{\rm e}(n,x)=0\,.$ (18)

From Eqs. (8) and (9), the derivative of the nucleonic energy per particle
with respect to $x$ is given by

$\frac{\partial\left(E_{\rm nuc}/A\right)}{\partial x}=-4\left(\frac{E_{\rm
PNM}}{A}-\frac{E_{\rm SNM}}{A}\right)(1-2x)\,.$ (19)

For the energy per particle of SNM, we can employ the standard Taylor
expansion about $n_{0}$,

$\frac{E_{\rm
SNM}}{A}(n)=E_{0}+\frac{1}{2}K_{0}\left(\frac{n-n_{0}}{3n_{0}}\right)^{2}+\cdots\,,$
(20)

where $n_{0}$, the saturation energy $E_{0}$, and the incompressibility
$K_{0}$ are constrained empirically. Higher-order terms beyond $K_{0}$ can be
neglected because we determine the symmetry energy only around $n_{0}$. See
the Supplemental Material in Ref. Essick  _et al._ (2021) for a quantification
of the effect of higher-order terms in $n$ and the presence of muons near
saturation density. For the parameters $n_{0}$, $E_{0}$, and $K_{0}$, we use
the ranges from Ref. Huth  _et al._ (2021) (means $\pm$ standard deviations of
Gaussian distributions):

$\displaystyle n_{0}$ $\displaystyle=0.164\pm 0.007\,\text{fm}^{-3}\,,$
$\displaystyle E_{0}$ $\displaystyle=-15.86\pm 0.57\,\text{MeV}\,,$ (21)
$\displaystyle K_{0}$ $\displaystyle=215\pm 40\,\text{MeV}\,.$

Putting all of this together, $\beta$-equilibrium must satisfy

$\frac{1-2x_{\beta}}{4}\Bigl{(}m_{p}-m_{n}+\mu_{\mathrm{e}}(n,x_{\beta})\Bigr{)}=\\\
\left(\frac{\varepsilon_{\beta}-\varepsilon_{\mathrm{e}}(n,x_{\beta})}{n}-m_{\mathrm{N}}-\frac{E_{\rm
SNM}}{A}(n)\right)\,.$ (22)

We self-consistently reconstruct the proton fraction for each EOS realization
by solving Eq. (22) for $x_{\beta}$ as a function of $n$ around $n_{0}$. For
this, we draw the parameters $E_{0}$, $K_{0}$, and $n_{0}$ from their
empirical distributions in Eq. (21) separately for each EOS, thereby
marginalizing over their uncertainty within our Monte-Carlo sums over EOS
realizations. We then calculate the PNM energy per particle $E_{\rm
PNM}/A(n)$, the symmetry energy $S_{0}$, its derivative $L$, and its curvature
$K_{\rm sym}$ as a function of baryon density $n$ in the vicinity of $n_{0}$
and report their values at the reference density,
$n_{0}^{\mathrm{ref}}=0.16\,\mathrm{fm}^{-3}$. In the following we use $n_{0}$
to denote this reference density, but note again that the uncertainty in the
empirical saturation point, Eq. (21), is included when extracting $S_{0}$,
$L$, and $K_{\mathrm{sym}}$ from EOS samples.

Figure 1:  Uncertainty relation between
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ and $L$ modeled on the 31 models
from Ref. Mondal  _et al._ (2016) (_red circles_) compared with 47 models from
Ref. Viñas  _et al._ (2014) (_blue squares_). (_left_) We model the
theoretical uncertainty with a conditional probability
$P(R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}|L)$ using a normal distribution
with mean given by Eq. (24). Shaded bands correspond to 1, 2, and 3-$\sigma$
uncertainties for $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ at each $L$.
(_bottom right_) Predicted cumulative distribution of residuals and empirical
distribution based on the fit to Ref. Mondal  _et al._ (2016), showing good
quantitative agreement between our model and the scatter between the
theoretical calculations. We note that the models from Ref. Viñas  _et al._
(2014) are systematically shifted compared to Ref. Mondal  _et al._ (2016),
but they are well represented by our uncertainty model.  Figure 2:  Analogous
to Fig. 1, but showing the conditional uncertainty for
$P(\alpha_{D}^{{}^{208}\mathrm{Pb}}S_{0}|L)$, modeled as a Gaussian with mean
given by Eq. (25), based on Ref. Roca-Maza  _et al._ (2013). Shaded bands
represent 1, 2, and 3-$\sigma$ uncertainty within our model. We again obtain
good quantitative agreement between our uncertainty model and the observed
scatter of the theoretical models.

With the mapping between the EOS and the parameters $E_{\mathrm{PNM}}/A$,
$S_{0}$, $L$, and $K_{\mathrm{sym}}$ established, we calculate a posterior
distribution

$P(E_{\mathrm{PNM}}/A,S_{0},L,K_{\rm sym}|\\{d\\})=\\\
\int\mathcal{D}\varepsilon_{\beta}\,P(\varepsilon_{\beta}|\\{d\\})P(E_{\mathrm{PNM}}/A,S_{0},L,K_{\rm
sym}|\varepsilon_{\beta})$ (23)

over the nuclear physics properties by conditioning on the astrophysical
observations and marginalizing over many EOS realizations.

To extract the neutron-skin thickness of 208Pb, we use an empirical fit
between $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ and $L$ based on the data in
Ref. Mondal  _et al._ (2016):

$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}\,[\mathrm{fm}]=0.0724+0.0019\times(L\,[\mathrm{MeV}]).$
(24)

This fit is calculated from a range of nonrelativistic Skyrme and relativistic
energy-density functionals. To model the uncertainty in this empirical
relation, we fit the distribution of
($R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$, $L$) from Ref. Mondal  _et al._
(2016) to a Gaussian with a mean given by Eq. (24), obtaining a standard
deviation of $0.0143\,\mathrm{fm}$. This uncertainty model and the residuals
of the fit are shown in Fig. 1. We also compare this fit with the density
functionals used in Ref. Viñas  _et al._ (2014). Our fit provides a good
representation of the spread between all these models.

Similarly, to connect our results to the electric dipole polarizability of
208Pb, $\alpha_{D}^{{}^{208}\mathrm{Pb}}$, we use an empirical fit between
$\alpha_{D}^{{}^{208}\mathrm{Pb}}\cdot S_{0}$ and $L$ based on Ref. Roca-Maza
_et al._ (2013), finding

$\alpha_{D}^{{}^{208}\mathrm{Pb}}\cdot
S_{0}\,[\mathrm{fm}^{3}\mathrm{MeV}]=493.5+3.08\times(L\,[\mathrm{MeV}])\,.$
(25)

We again model the conditional distribution
$P(\alpha_{D}^{{}^{208}\mathrm{Pb}}S_{0}|L)$ as a Gaussian with mean given by
Eq. (25) and a standard deviation of $27.6\,\mathrm{fm}^{3}\mathrm{MeV}$. This
uncertainty model is shown in Fig. 2.

##  IV Results

PREX-II Figure 3:  Correlations between $S_{0}$, $L$, $K_{\rm sym}$, and $R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$ within our nonparametric prior (_unshaded yellow_) and Astro-only posterior (_shaded green_) as well as the $\chi\mathrm{EFT}$-marginalized (_shaded blue_), $\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$, $\mathrm{MBPT}^{(2010)}_{\mathrm{mixed}}$, $\mathrm{MBPT}^{(2013)}_{\mathrm{N^{3}LO}}$, and $\mathrm{MBPT}^{(2019)}_{\mathrm{N^{3}LO}}$ Astro-only posteriors (_unshaded greys_ , ordered from lighter to darker with increasing $L$, see Table 1). Joint distributions show 90% credible regions, and the horizontal bands (_pink_) represent PREX-II 90% credible regions, with dashed lines the corresponding 68% (1-$\sigma$) regions. The expanded ($S_{0}$, $L$) panel (_upper right_) compares our nonparametric prior, Astro-only posterior, and $\chi\mathrm{EFT}$+Astro posterior to other constraints: (_white region_) Lattimer and Prakash Lattimer and Prakash (2016) (overlap region of various nuclear experimental constraints), the unitary-gas (UG) bound from Ref. Tews  _et al._ (2017), and the values reported by Reed et al. Reed  _et al._ (2021) based on the PREX-II results. In addition, we show the correlation obtained from the experimental $\alpha_{D}^{{}^{208}\mathrm{Pb}}$ Tamii  _et al._ (2011) using Eq. (25).  |  | $\frac{E_{\mathrm{PNM}}}{A}(n_{0})\ [\mathrm{MeV}]$ | $S_{0}\ [\mathrm{MeV}]$ | $L\ [\mathrm{MeV}]$ | $K_{\mathrm{sym}}\ [\mathrm{MeV}]$ | $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}\ [\mathrm{fm}]$ | $\alpha_{D}^{{}^{208}\mathrm{Pb}}\ [\mathrm{fm}^{3}]$
---|---|---|---|---|---|---|---
Nonparametric | Prior | $17.5^{+14.6}_{-7.7}$ | $33.3^{+14.7}_{-8.2}$ | $38^{+109}_{-41}$ | $-255^{+853}_{-566}$ | $0.14^{+0.19}_{-0.09}$ | $18.9^{+4.1}_{-4.7}$
Astro Posterior | $19.3^{+11.7}_{-8.5}$ | $35.1^{+11.6}_{-8.9}$ | $58^{+61}_{-56}$ | $-240^{+559}_{-503}$ | $0.19^{+0.12}_{-0.11}$ | $19.0^{+3.8}_{-3.9}$
Astro+PREX-II Post. | $21.5^{+10.8}_{-8.3}$ | $37.3^{+11.8}_{-7.5}$ | $80^{+51}_{-46}$ | $-223^{+608}_{-565}$ | $0.23^{+0.10}_{-0.10}$ | $19.6^{+3.9}_{-4.4}$
Astro+$\alpha_{D}^{{}^{208}\mathrm{Pb}}$ Post. | $18.4^{+7.4}_{-7.8}$ | $34.2^{+7.4}_{-7.9}$ | $61^{+49}_{-57}$ | $-172^{+483}_{-388}$ | $0.19^{+0.10}_{-0.12}$ | $19.8^{+2.0}_{-2.0}$
$\chi\mathrm{EFT}$-marginalized | Prior | $16.7^{+1.5}_{-1.3}$ | $32.5^{+1.9}_{-1.8}$ | $47^{+15}_{-15}$ | $-119^{+129}_{-133}$ | $0.16^{+0.04}_{-0.04}$ | $19.6^{+1.7}_{-2.0}$
Astro Posterior | $16.9^{+1.5}_{-1.4}$ | $32.7^{+1.9}_{-1.8}$ | $49^{+14}_{-15}$ | $-107^{+124}_{-128}$ | $0.17^{+0.04}_{-0.04}$ | $19.6^{+1.9}_{-1.7}$
Astro+PREX-II Post. | $17.1^{+1.5}_{-1.5}$ | $33.0^{+2.0}_{-1.8}$ | $53^{+14}_{-15}$ | $-91^{+118}_{-130}$ | $0.17^{+0.04}_{-0.04}$ | $19.8^{+1.7}_{-1.9}$
Astro+$\alpha_{D}^{{}^{208}\mathrm{Pb}}$ Post. | $16.9^{+1.5}_{-1.4}$ | $32.7^{+1.9}_{-1.8}$ | $51^{+13}_{-14}$ | $-98^{+117}_{-124}$ | $0.17^{+0.04}_{-0.03}$ | $19.8^{+1.5}_{-1.9}$
$\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$ Lynn  _et al._ (2016) | Original Work | [14.2, 18.8] | [28.6, 36.2] | [23.8, 58.2] | - | - | -
Prior | $16.4^{+1.0}_{-0.9}$ | $32.2^{+1.5}_{-1.5}$ | $39^{+11}_{-100}$ | $-179^{+111}_{-112}$ | $0.15^{+0.03}_{-0.03}$ | $19.1^{+1.7}_{-1.7}$
Astro Posterior | $16.5^{+1.1}_{-0.9}$ | $32.4^{+1.5}_{-1.5}$ | $41^{+11}_{-11}$ | $-165^{+114}_{-112}$ | $0.15^{+0.03}_{-0.03}$ | $19.2^{+1.6}_{-1.9}$
Astro+PREX-II Post. | $16.7^{+1.1}_{-1.0}$ | $32.5^{+1.7}_{-1.4}$ | $44^{+12}_{-12}$ | $-151^{+124}_{-108}$ | $0.16^{+0.03}_{-0.03}$ | $19.3^{+1.6}_{-1.9}$
Astro+$\alpha_{D}^{{}^{208}\mathrm{Pb}}$ Post. | $16.5^{+1.1}_{-0.9}$ | $32.2^{+1.5}_{-1.5}$ | $43^{+11}_{-10}$ | $-153^{+111}_{-107}$ | $0.16^{+0.03}_{-0.03}$ | $19.4^{+1.5}_{-1.8}$
$\mathrm{MBPT}^{(2010)}_{\mathrm{mixed}}$ Hebeler and Schwenk (2010); Hebeler  _et al._ (2013) | Original Work | [14.3, 18.4] | [29.7, 33.2] | [32.5, 57.0] | - | [0.14, 0.20] | -
Prior | $16.6^{+1.2}_{-1.2}$ | $32.4^{+1.7}_{-1.6}$ | $43^{+11}_{-11}$ | $-149^{+104}_{-100}$ | $0.16^{+0.03}_{-0.03}$ | $19.3^{+1.7}_{-1.7}$
Astro Posterior | $16.7^{+1.3}_{-1.2}$ | $32.6^{+1.7}_{-1.7}$ | $44^{+12}_{-11}$ | $-145^{+101}_{-103}$ | $0.16^{+0.03}_{-0.03}$ | $19.3^{+1.7}_{-1.7}$
Astro+PREX-II Post. | $16.9^{+1.3}_{-1.3}$ | $32.8^{+1.8}_{-1.7}$ | $47^{+12}_{-12}$ | $-138^{+100}_{-102}$ | $0.16^{+0.03}_{-0.04}$ | $19.4^{+1.6}_{-1.8}$
Astro+$\alpha_{D}^{{}^{208}\mathrm{Pb}}$ Post. | $16.7^{+1.3}_{-1.3}$ | $32.5^{+1.7}_{-1.7}$ | $46^{+12}_{-11}$ | $-138^{+97}_{-101}$ | $0.16^{+0.03}_{-0.03}$ | $19.5^{+1.5}_{-1.8}$
$\mathrm{MBPT}^{(2013)}_{\mathrm{N^{3}LO}}$ Tews  _et al._ (2013) | Original Work | [13.4, 20.1] | [28.9, 34.9] | [43.0, 66.6] | - | - | -
Prior | $16.9^{+1.9}_{-1.9}$ | $32.8^{+2.2}_{-2.2}$ | $52^{+13}_{-13}$ | $-86^{+94}_{-103}$ | $0.17^{+0.04}_{-0.03}$ | $19.9^{+1.6}_{-1.8}$
Astro Posterior | $17.1^{+1.8}_{-1.9}$ | $32.9^{+2.2}_{-2.1}$ | $53^{+13}_{-12}$ | $-86^{+96}_{-101}$ | $0.18^{+0.03}_{-0.04}$ | $19.9^{+1.6}_{-1.8}$
Astro+PREX-II Post. | $17.4^{+1.9}_{-1.9}$ | $33.2^{+2.2}_{-2.2}$ | $55^{+13}_{-12}$ | $-80^{+99}_{-93}$ | $0.18^{+0.03}_{-0.03}$ | $19.9^{+1.6}_{-1.8}$
Astro+$\alpha_{D}^{{}^{208}\mathrm{Pb}}$ Post. | $17.1^{+1.8}_{-1.9}$ | $32.9^{+2.1}_{-2.0}$ | $54^{+13}_{-12}$ | $-84^{+102}_{-92}$ | $0.18^{+0.03}_{-0.04}$ | $19.9^{+1.5}_{-1.8}$
$\mathrm{MBPT}^{(2019)}_{\mathrm{N^{3}LO}}$ Drischler  _et al._ (2019) | Original Work | [15.3, 18.7] | - | - | - | - | -
Prior | $17.0^{+1.4}_{-1.4}$ | $32.8^{+1.8}_{-1.8}$ | $53^{+12}_{-12}$ | $-63^{+117}_{-113}$ | $0.18^{+0.03}_{-0.03}$ | $20.0^{+1.6}_{-1.9}$
Astro Posterior | $17.1^{+1.3}_{-1.2}$ | $32.9^{+1.8}_{-1.7}$ | $54^{+11}_{-11}$ | $-63^{+114}_{-117}$ | $0.18^{+0.03}_{-0.03}$ | $20.0^{+1.6}_{-1.9}$
Astro+PREX-II Post. | $17.2^{+1.3}_{-1.3}$ | $33.1^{+1.7}_{-1.8}$ | $56^{+11}_{-12}$ | $-53^{+115}_{-116}$ | $0.18^{+0.03}_{-0.03}$ | $20.1^{+1.5}_{-2.0}$
Astro+$\alpha_{D}^{{}^{208}\mathrm{Pb}}$ Post. | $17.1^{+1.3}_{-1.3}$ | $32.9^{+1.7}_{-1.6}$ | $54^{+11}_{-11}$ | $-61^{+111}_{-114}$ | $0.18^{+0.03}_{-0.03}$ | $20.0^{+1.5}_{-1.8}$
| Prior | $16.9^{+1.2}_{-1.2}$ | $32.8^{+1.7}_{-1.7}$ | $53^{+10}_{-10}$ | $-87^{+99}_{-101}$ | $0.17^{+0.03}_{-0.03}$ | $20.0^{+1.5}_{-1.9}$
$\mathrm{MBPT}^{(2020\mathrm{\,GP})}_{\mathrm{N^{3}LO}}$ | Astro Posterior | $17.0^{+1.3}_{-1.1}$ | $32.8^{+1.7}_{-1.5}$ | $53^{+9}_{-10}$ | $-86^{+95}_{-104}$ | $0.18^{+0.03}_{-0.03}$ | $20.0^{+1.5}_{-1.9}$
Drischler  _et al._ (2020b, c) | Astro+PREX-II Post. | $17.1^{+1.2}_{-1.1}$ | $32.9^{+1.7}_{-1.6}$ | $54^{+10}_{-9}$ | $-81^{+98}_{-97}$ | $0.18^{+0.03}_{-0.03}$ | $20.0^{+1.5}_{-1.9}$
| Astro+$\alpha_{D}^{{}^{208}\mathrm{Pb}}$ Post. | $17.0^{+1.3}_{-1.1}$ | $32.8^{+1.7}_{-1.4}$ | $53^{+10}_{-9}$ | $-85^{+93}_{-103}$ | $0.18^{+0.03}_{-0.03}$ | $20.0^{+1.4}_{-1.9}$
Table 1:  Medians and 90% highest-probability-density credible regions for
selected nuclear properties. All $\chi\mathrm{EFT}$ results trust the
theoretical prediction up to $p_{\mathrm{max}}/c^{2}=4.3\times
10^{12}\,\mathrm{g}/\mathrm{cm}^{3}$, corresponding to
$n(p_{\mathrm{max}})\sim n_{0}$. $\chi\mathrm{EFT}$-marginalized results
combine results from $\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$ Lynn  _et al._
(2016), $\mathrm{MBPT}^{(2010)}_{\mathrm{mixed}}$ Hebeler and Schwenk (2010);
Hebeler  _et al._ (2013), $\mathrm{MBPT}^{(2013)}_{\mathrm{N^{3}LO}}$ Tews
_et al._ (2013), and $\mathrm{MBPT}^{(2019)}_{\mathrm{N^{3}LO}}$ Drischler
_et al._ (2019) with equal weight a priori. We also tabulate results from each
of these 4 $\chi\mathrm{EFT}$ predictions separately. In addition, we provide
results from $\mathrm{MBPT}^{(2020\mathrm{\,GP})}_{\mathrm{N^{3}LO}}$
Drischler  _et al._ (2020b, c) for comparison with
$\mathrm{MBPT}^{(2019)}_{\mathrm{N^{3}LO}}$, both of which use the same
microscopic calculations. Where possible, we also provide bounds quoted for
the original studies, given by envelopes containing all models considered
within the original studies. As such, they do not have an immediate
statistical interpretation and are wider than our 90% credible regions.

We first summarize our conclusions about
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ in Sec. IV.1 before comparing
constraints on broader sets of nuclear properties near $n_{0}$ in Sec. IV.2.
Section IV.3 summarizes what we can learn about NS properties from current
experimental constraints and possible future improvements.

###  IV.1 Symmetry-Energy Parameters and Neutron-Skin Thickness in Lead

We begin by discussing our findings for $S_{0}$, $L$, $K_{\rm sym}$, and
$R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$, shown in Fig. 3. We plot the
nonparametric prior, the posterior constrained only by astrophysical data, and
the posterior additionally constrained by $\chi$EFT calculations up to
$n\approx n_{0}$. Our GPs are conditioned on $\chi$EFT up to a maximum
pressure, $p_{\mathrm{max}}$. To translate this into a density, we report the
median density at $p_{\mathrm{max}}$ a priori; the exact density at
$p_{\mathrm{max}}$ varies due to uncertainty in the EOS from $\chi$EFT. In
addition to the constraints obtained by marginalizing over the four separate
$\chi\mathrm{EFT}$ calculations, we also show the posteriors for each
individual $\chi$EFT calculation. Finally, we also compare our results with
the recent constraints on $R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$ and $L$ from
the PREX-II experiment Adhikari  _et al._ (2021), where we have translated
from $R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$ to $L$ using our model of the
theoretical uncertainty in the correlation between these two quantities. Prior
and posterior credible regions are also provided in Table 1.

The priors and Astro-only posteriors for the nonparametric inference are very
broad, and we find large ranges for $S_{0}$, $L$, $K_{\rm sym}$, and
$R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$ (see Table 1). The astrophysical data
slightly informs our uncertainty in $S_{0}$ and $L$, shifting the median
values of their distributions, but the 90% confidence intervals are less
impacted. The astrophysical data does not strongly constrain $K_{\rm sym}$,
but suggests that it is negative. Taken together, this highlights the fact
that astrophysical information alone is not sufficient to pin down properties
of the EOS around nuclear saturation density.

When we additionally constrain the nonparametric EOSs using the four $\chi$EFT
calculations, we obtain much narrower posteriors. It is noteworthy that the
$\chi\mathrm{EFT}$ posteriors fall near the maximum of the Astro-only
nonparametric posterior. We stress that this need not have been the case,
because the nonparametric Astro-only posterior does not know anything about
$\chi\mathrm{EFT}$. While the four individual calculations result in slightly
different values for $L$ and, hence, $R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$,
overall all four $\chi$EFT calculations are very consistent.

When we compare our findings with the recent PREX-II results, we find that the
nonparametric Astro-only posterior prefers lower values for $L$ and
$R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$, in good agreement with the result that
includes $\chi\mathrm{EFT}$. Both posteriors peak at similar values of $L$, on
the order of $50$–$60$ MeV, and of $R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$, on
the order of $0.15$–$0.20$ fm. However, uncertainties are large and
nonparametric Astro-only results remain compatible with both the $\chi$EFT
prediction and the PREX-II results. Nonetheless, when we additionally
condition on $\chi$EFT calculations, we find that the PREX-II result for
$R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$ and the associated range for $L$
($69$–$143$ MeV at $1\sigma$ Reed  _et al._ (2021)), are only in mild tension
with the $\chi\mathrm{EFT}$ predictions.

Finally, we compare our findings for $S_{0}$ and $L$ with other constraints in
the upper-right panel of Fig. 3. Our $\chi\mathrm{EFT}$+Astro posterior is
very consistent with the overlap region from various experimental constraints
from Lattimer and Prakash Lattimer and Prakash (2016) and lies fully witin the
bounds of the unitary gas conjecture Tews  _et al._ (2017). While the
extraction of $S_{0}$ and $L$ from PREX-II by Reed et al. Reed  _et al._
(2021) leads to significantly larger central values, it also has large 90%
credible regions, which overlap with our $\chi\mathrm{EFT}$+Astro posterior.
In addition, we show here the correlation obtained from the experimental value
of the dipole polarizability $\alpha_{D}^{{}^{208}\mathrm{Pb}}$ Tamii  _et
al._ (2011) with our uncertainty model Eq. (25) assuming uninformative priors
for $S_{0}$ and $L$. This overlaps nicely with all extractions.

###  IV.2 Compatibility of Astrophysical, Experimental, and Theoretical
Results for Nuclear Properties

Figure 4:  Priors (_grey, unshaded_), Astro-only posteriors (_left side of
violins, green unshaded_), Astro+PREX-II posteriors (_right side of violins,
red shaded_), Astro+$\alpha_{D}^{{}^{208}\mathrm{Pb}}$ posteriors (_right side
of violins, blue shaded+hatched_), and Astro+PREX-
II+$\alpha_{D}^{{}^{208}\mathrm{Pb}}$ posteriors (_left side of violins, grey
shaded+dots_) for $\chi\mathrm{EFT}$-marginalized results as a function of the
maximum pressure up to which we trust $\chi\mathrm{EFT}$. The left-most curves
(median $n\sim 0.3n_{0}$) are equivalent to the nonparametric results in Fig.
3. Horizontal bands (dashed lines) correspond to 90% (1-$\sigma$) credible
regions from PREX-II Adhikari  _et al._ (2021)
($R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$; _pink_) and the electric dipole
polarizability Tamii  _et al._ (2011) ($\alpha_{D}^{{}^{208}\mathrm{Pb}}$;
_orange_). When translating experimental data to their correlated properties
in this figure (e.g., horizontal $\alpha_{D}^{{}^{208}\mathrm{Pb}}$ bands for
$L$ and $R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$), we employ our uncertainty
relations in the theoretical correlations (Eqs. (24) and (25), assuming
$S_{0}=32.5\,\mathrm{MeV}$ for the latter).

In Fig. 4, we show the evolution of our constraints on $L$,
$R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$, and $\alpha_{D}^{{}^{208}\mathrm{Pb}}$
as a function of the maximum density up to which we condition our prior on
$\chi$EFT. In addition to the posterior conditioned only on astrophysical
data, we show results for three cases that are additionally conditioned on
either the PREX-II $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ data Adhikari
_et al._ (2021), the $\alpha_{D}^{{}^{208}\mathrm{Pb}}$ data from Ref. Tamii
_et al._ (2011), or both.

Figure 5:  Probability of PREX-II disagreeing with posteriors conditioned on
$\chi\mathrm{EFT}$ up to $p_{\mathrm{max}}$ by at least the measured
difference given experimental uncertainties ($p$-values, _solid lines_). We
also show the $p$-values for a hypothetical experiment producing the same mean
as PREX-II with half the uncertainty (_dashed lines_). Results are given for
nonparametric Astro-only posteriors (_black horizontal lines_),
$\chi\mathrm{EFT}$-marginalized (_blue_),
$\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$ (_yellow_),
$\mathrm{MBPT}^{(2010)}_{\mathrm{mixed}}$ (_orange_),
$\mathrm{MBPT}^{(2013)}_{\mathrm{N^{3}LO}}$ (_purple_), and
$\mathrm{MBPT}^{(2019)}_{\mathrm{N^{3}LO}}$ (_red_).

Figure 6:  Bayes factors between priors conditioned on $\chi\mathrm{EFT}$
calculations up to different $p_{\mathrm{max}}$ vs. the priors not conditioned
on $\chi\mathrm{EFT}$ at all for (_top_) astrophysical data, (_middle_)
Astro+$\alpha_{D}^{{}^{208}\mathrm{Pb}}$, and (_bottom_) Astro+PREX-II data.
We show results for the $\chi\mathrm{EFT}$-marginalized calculations (_blue_)
as well as the $\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$ (_yellow_),
$\mathrm{MBPT}^{(2010)}_{\mathrm{mixed}}$ (_orange_),
$\mathrm{MBPT}^{(2013)}_{\mathrm{N^{3}LO}}$ (_purple_), and
$\mathrm{MBPT}^{(2019)}_{\mathrm{N^{3}LO}}$ (_red_) calculations separately.

Theory Favored
Theory Disfavored

Figure 7:  Bayes factors between priors conditioned on $\chi\mathrm{EFT}$ vs.
priors not conditioned on $\chi\mathrm{EFT}$ at all for different nuclear data
when we first condition on the astrophysical observations (include them as
part of the prior). We show the result for (_top_)
$\alpha_{D}^{{}^{208}\mathrm{Pb}}$ and (_bottom_) PREX-II data.

If we do not condition the prior on $\chi\mathrm{EFT}$ (left-most violins,
where we match directly to the crust at $0.3n_{0}$), the Astro-only posterior
retains large uncertainties for all three quantities. As stated before,
astrophysical data inform our knowledge of $L$ and
$R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$ to some degree, but they do not add
further information about $\alpha_{D}^{{}^{208}\mathrm{Pb}}$ because $S_{0}$
is not strongly constrained. When we additionally condition on the recent
PREX-II result, uncertainties remain large, but the posteriors for $L$ and
$R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$ are pushed to higher values.
Alternatively, conditioning instead on the $\alpha_{D}^{{}^{208}\mathrm{Pb}}$
measurement, the posteriors for $L$ and $R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$
agree very well with the Astro-only result, highlighting the consistency of
this experiment and neutron-star observations; see also Table 1. In this case,
as expected, the posterior for $\alpha_{D}^{{}^{208}\mathrm{Pb}}$ is much
narrower. Conditioning on astrophysical observations and both PREX-II and
$\alpha_{D}^{{}^{208}\mathrm{Pb}}$ produces posteriors for $L$ and
$R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$ similar to those obtained by only
conditioning on astrophysical observations and PREX-II because there is enough
additional freedom in $S_{0}$ to accomodate the
$\alpha_{D}^{{}^{208}\mathrm{Pb}}$ measurements for almost any $L$ (see also
Fig. 9).

When conditioning the priors on $\chi$EFT constraints to higher densities, all
posteriors start to overlap more. They agree with each other very closely if
we condition up to $n_{0}$, where the $\chi\mathrm{EFT}$ constraints dominate.
In this case, the tension of our process with the PREX-II results is maximized
but nonetheless remains mild due to the large PREX-II uncertainties. On the
other hand, the agreement with the $\alpha_{D}^{{}^{208}\mathrm{Pb}}$ result
improves the more we trust the $\chi$EFT constraints.

Prior

Astro Posterior

Astro+PREX-II Posterior

Figure 8:  Median and 90% symmetric credible regions for the prior (_left_),
Astro-only posterior (_middle_), and Astro+PREX-II posterior (_right_) for all
EOS and all values of $L$ (_green_), EOS with
$30\,\,\text{MeV}<L<70\,\,\text{MeV}$ (_hatched blue_), and EOS with
$100\,\,\text{MeV}<L$ (_purple_). The main effect of the PREX-II data is to
rule out some of the very soft EOS at low densities ($L\lesssim
30\,\,\text{MeV}$).

Figure 5 shows how the probability ($p$-value) that the true
$R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$ differs from the PREX-II mean at least
as much as the Astro+$\chi\mathrm{EFT}$ posterior suggests, given the
uncertainty in PREX-II’s measurement. The $p$-values decrease as we trust
$\chi\mathrm{EFT}$ up to higher densities, and we estimate a p-value of
$12.3\%$ when trusting $\chi\mathrm{EFT}$ up to $n\sim n_{0}$ (c.f., $25.3\%$
for the nonparametric Astro-only posterior). However, if a hypothetical
experiment confirmed the PREX-II mean value with half the uncertainty, this
$p$-value would be reduced to $0.6\%$. In fact, a hypothetical
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ measurement with half the
uncertainty has a smaller $p$-value under the nonparametric Astro-only
posterior than the $\chi\mathrm{EFT}$-marginalized posterior has with the
current $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ measurement uncertainties.

To investigate this further, we compute Bayes factors between the processes
conditioned on $\chi\mathrm{EFT}$ up to various pressures vs. processes not
conditioned on $\chi$EFT at all (Figs. 6 and 7) for different sets of data:
Astro-only, Astro+$\alpha_{D}^{{}^{208}\mathrm{Pb}}$, and Astro+PREX-II (Fig.
6) and when astrophysical data is already included in the prior (Fig. 7). In
addition to the posteriors marginalized over all four $\chi$EFT results, we
also show the Bayes factors for the individual
$\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$,
$\mathrm{MBPT}^{(2010)}_{\mathrm{mixed}}$,
$\mathrm{MBPT}^{(2013)}_{\mathrm{N^{3}LO}}$, and
$\mathrm{MBPT}^{(2019)}_{\mathrm{N^{3}LO}}$ results. These Bayes factors
quantify the relative likelihood of obtaining the observed data under
different models, specifically whether $\chi\mathrm{EFT}$=informed priors are
more ($\mathcal{B}^{\mathrm{theory}}_{\mathrm{agnostic}}>1$) or less
($\mathcal{B}^{\mathrm{theory}}_{\mathrm{agnostic}}<1$) likely to have
produced the observed data compared to our completely nonparametric prior.

Considering only astrophysical data, we find that $\chi$EFT is preferred over
the theory-agnostic result up to at least nuclear saturation density. This is
also true for the individual calculations, although we find that the Bayes
factor in favor of $\mathrm{MBPT}^{(2013)}_{\mathrm{N^{3}LO}}$ and
$\mathrm{MBPT}^{(2019)}_{\mathrm{N^{3}LO}}$ are a factor of two larger than
for $\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$. This agrees with previous
results Essick  _et al._ (2020b) and could be associated with the higher-order
$\chi\mathrm{EFT}$ interactions included in
$\mathrm{MBPT}^{(2013)}_{\mathrm{N^{3}LO}}$ and
$\mathrm{MBPT}^{(2019)}_{\mathrm{N^{3}LO}}$ that tend to increase the pressure
and are not included in $\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$. It could
also be associated with the different regularization schemes employed in these
calculations. However, this preference may be due to different widths of the
theoretical uncertainty bands within different $\chi\mathrm{EFT}$
calculations. These Bayes factors are likely driven partly by Occam factors
where a wider prior is penalized even though all models may achieve similar
maximum likelihoods. For example, $\chi\mathrm{EFT}$ yields a narrower prior
which penalizes the freedom in the nonparametric model without
$\chi\mathrm{EFT}$. Similarly, the $\mathrm{MBPT}^{(2013)}_{\mathrm{N^{3}LO}}$
and $\mathrm{MBPT}^{(2019)}_{\mathrm{N^{3}LO}}$ priors predict higher median
pressures with smaller uncertainties than
$\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$, and both effects will tend to
increase the relative Bayes factor. We also find that the astrophysical
observations can only distinguish between individual $\chi\mathrm{EFT}$
calculations if we trust them up to $\gtrsim 0.75n_{0}$.

When additionally including $\alpha_{D}^{{}^{208}\mathrm{Pb}}$, the Bayes
factors in favor of $\chi\mathrm{EFT}$ increase by a factor of two. In
contrast, including the PREX-II information decreases the Bayes factors by a
factor of $\lesssim 2$. Figure 7 shows this behavior explicitly by first
conditioning on the astrophysical observations, thereby isolating the new
information obtained from the inclusion of each nuclear experiment.
Nonetheless, in all cases, models conditioned on $\chi\mathrm{EFT}$
information are favored when we consider all nuclear experiments and
astrophysical observations simultaneously (i.e., Bayes factors remain larger
than $1$ in Fig. 6). We find that the Bayes factors are largest for
$\mathrm{MBPT}^{(2013)}_{\mathrm{N^{3}LO}}$ and
$\mathrm{MBPT}^{(2019)}_{\mathrm{N^{3}LO}}$ and smallest for
$\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$. Again, this is likely due to a
combination of high-order interactions only present in some calculations,
choices of the regulator scheme, and the widths of prior uncertainty bands.

Nonparametric
Astro+PREX-II
Posterior
$\chi\mathrm{EFT}$
Astro+PREX-II
Posterior

Figure 9:  Correlations between $S_{0}$ and $L$ when we model the PREX-II
estimate with different uncertainties: (_top_) the actual measurement
uncertainty, (_middle_) a hypothetical measurement with half the PREX-II
uncertainty, and (_bottom_) a hypothetical measurement with vanishingly small
uncertainty for $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$. We show the
nonparametric prior (_unshaded yellow_), Astro-only posterior (_unshaded
green_), and Astro+PREX-II posterior (_shaded red_) as well as the
$\chi\mathrm{EFT}$-marginalized Astro-only posterior (_unshaded blue_) and
Astro+PREX-II posterior (_shaded light blue_). As in Fig. 3, (_pink_) shaded
vertical bands represent (real and hypothetical) PREX-II 90% credible regions
and dashed lines show the 1-$\sigma$ credible regions uncertainty. Improved
measurements of $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ are still consistent
with a wide range of $S_{0}$ within the nonparametric inference.

Given the mild tension between the PREX-II value for
$R_{\mathrm{skin}}^{{}^{208}{\rm Pb}}$ and that inferred from the
astrophysical inference with $\chi$EFT information, we investigate what kind
of EOS behavior is required to satisfy both the PREX-II and astrophysical
constraints. In Fig. 8, we show the pressure and the speed of sound $c_{s}$ as
a function of density for the nonparametric process conditioned only on
astrophysical data for all values of $L$, for $30\,\text{MeV}<L\leq
70\,\text{MeV}$, and for $L>100\,\text{MeV}$. Note that this is a stricter
requirement than the nominal PREX-II observations suggest at 1-$\sigma$. We
find that the speed of sound generally increases with density. However, if we
assume $L>100$ MeV, we find a local maximum in the median $c_{s}(n)$ just
below $n_{0}$, although the uncertainties in $c_{s}$ are large. The reason for
this feature is that EOSs that are stiff at low densities (large $L$) need to
soften beyond $n_{0}$ to remain consistent with astrophysical data (small
tidal deformabilities from GWs). Should the PREX-II constraints be confirmed
with smaller uncertainty in the future, this might favor the existence of a
phase transition between $1$–$2n_{0}$. However, given current uncertainties,
there is no strong preference for such exotic EOS phenomenology based on the
data.

Figure 10:  Correlations of $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ and the
radii of NSs with M=1.0, 1.4, and 1.8 $M_{\odot}$ with $L$. Colors and shading
match those in Fig. 9.

Finally, we can ask what would happen to our uncertainty in $S_{0}$ and $L$ if
a series of hypothetical future experiments confirmed the mean of
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ from PREX-II but with smaller
uncertainties. In Fig. 5, we already showed the $p$-values for such a case,
which highlight the increased tension with $\chi\mathrm{EFT}$ calculations. In
Fig. 9, we show the joint posteriors on $S_{0}$ and $L$ with the current PREX-
II uncertainty, half the current uncertainty, and with a perfect
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ measurement with vanishing
uncertainty, where the remaining uncertainty in $L$ is due purely to the
uncertainty in the theoretical correlation in Eq. (24). An increased
hypothetical precision for $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ could
change our knowledge of $L$ dramatically, possibly rendering it incompatible
with the $\chi\mathrm{EFT}$ predictions when using Eq. (24). However, although
the nonparametric Astro+PREX-II posteriors shift compared to the Astro-only
posteriors, we never find any significant disagreement. Indeed, the width of
our posterior for $S_{0}$ is nearly unchanged, even if we assume vanishingly
small measurement uncertainty for $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$.
This is another demonstration that current astrophysical data from NSs in the
observed mass range cannot strongly constrain nuclear interactions around
$n_{0}$ without further assumptions about the EOS. The agnostic priors do not
closely follow any particular theory (which would generically predict stronger
correlations between $S_{0}$ and $L$).

###  IV.3 Comparisons between PREX-II, $\chi\mathrm{EFT}$, and Astrophysical
Data for NS Observables

Having shown that current strophysical observations of NSs carry only limited
information about densities below nuclear saturation, we demonstrate that the
inverse is true as well. Improved measurements of
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$, or even hypothetical direct
measurements of $L$, will not significantly improve our knowledge of the
macroscopic properties of NSs with masses of $\gtrsim 1.2\,M_{\odot}$, without
additional theory input for the EOS. Fundamentally, this is because the
central densities of astrophysical NSs are above 2$n_{0}$ (see Ref. Legred
_et al._ (2021) for a recent inference of the relation between NS masses and
central densities), while the neutron-skin thickness and the symmetry energy
parameters describe matter around $n_{0}$. Constraints at nuclear saturation
density, then, must be extrapolated to higher densities to inform the
properties of NSs. In the nonparametric priors used here, there is enough
freedom that such extrapolations only introduce weak correlations between $L$
and, e.g., the radius of NSs. Strong correlations, like those in Ref. Reed
_et al._ (2021), thus also depend on the model used to describe the EOS above
nuclear densities.

Figure 11:  Pearson correlation coefficients between $L$ and $R(M)$ for
marginalized $\chi\mathrm{EFT}$ results (_blue circles, solid lines_) and
$\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$ (_yellow squares, dashed lies_). In
order from lightest to darkest lines (_top to bottom_), we plot the
correlation between $L$ and $R(1.0\,M_{\odot})$, $R(1.4\,M_{\odot})$, and
$R(1.8\,M_{\odot})$.

Figure 12:  Correlations between (_left_) $R_{1.4}$, (_right_) $\Lambda_{1.4}$
and $L$ when we model the PREX-II measurement with different uncertainties:
the actual measurement uncertainty (_top_), a hypothetical measurement with
half the PREX-II uncertainty (_middle_), and a hypothetical measurement with
vanishingly small uncertainty for $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$
(_bottom_). Colors and shading match those in Fig. 9. We see that even a
perfect measurement of $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ does not
significantly alter our knowledge of the macroscopic properties of typical
astrophysical NSs within the nonparametric inference.

We summarize the impact of current $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$
constraints from PREX-II on NSs observables in Fig. 10. As in Figs. 4 and 9,
we see that the PREX-II observations do increase the inferred value of $L$
when we do not condition on $\chi\mathrm{EFT}$. However, this translates only
into a modest shift in the radius of $1.0\,M_{\odot}$ stars ($R_{1.0}$) and
virtually no change for the radii of $1.4$ or $1.8\,M_{\odot}$ stars
($R_{1.4}$ and $R_{1.8}$, respectively) when we condition on existing
astrophysical data. While we observe correlations between $L$ and $R(M)$ a
priori, these are intrinsically broad (broader than is often assumed Lattimer
and Prakash (2001) and _not_ one-to-one) and weaken for NSs with higher
masses. These broad correlations are loose enough that astrophysical
observations are able to constrain the NS properties while remaining
consistent with a wide range of $L$ values. We find this behavior also in our
processes which are conditioned on $\chi\mathrm{EFT}$ calculations up to
$n_{0}$.

We also consider whether the inclusion of nuclear theory predictions up to
higher densities induces stronger correlations between $L$ and $R(M)$ in Fig.
11. Specifically, we show the Pearson correlation coefficient between $L$ and
$R(M)$ under the Astro-only posteriors as a function of the maximum density up
to which we trust $\chi\mathrm{EFT}$. Generally, we see an increase in the
correlation as we trust $\chi\mathrm{EFT}$ up to higher densities, as
expected, although the rate of increase slows at higher densities ($\gtrsim
1.3n_{0}$) for the $\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$ calculation.
This is likely due to the increase of the theoretical uncertainty band from
$\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$ with density, and therefore
conditioning on this theoretical prediction imposes a looser constraint. Taken
to an extreme (high $p_{\mathrm{max}}$ and small theoretical uncertainties),
one sees how trusting a particular theoretical extrapolation to high densities
will introduce a strong correlation between $L$ and $R(M)$. However, we note
that the theoretical uncertainties in current $\chi\mathrm{EFT}$ calculations
naturally limit the strength of such correlations, reaching a maximum
correlation coefficient of only $\simeq 0.5$ between $L$ and $R_{1.4}$, even
when we trust $\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$ up to $>1.8n_{0}$.
This may be refined with improved nuclear theory calculations at higher
densities. As expected, the correlation with $L$ is weaker for heavier NSs,
see for $R_{1.8}$ in Figs. 10 and 11, which is why the recent NICER+XMM
observations of J0740+662 ($M=2.08\pm 0.07\,M_{\odot}$) Miller  _et al._
(2021); Riley  _et al._ (2021) will not constrain the EOS substanially at
$n_{0}$.

Fig. 12 further demonstrates that improved constraints on $L$ will only
significantly change our knowledge of $R_{1.4}$ with improved nuclear theory
calculations to higher densities (Fig. 12 trusts $\chi\mathrm{EFT}$ up to
$n_{0}$). Figure 12 demonstrates this explicitly, where $\chi\mathrm{EFT}$
input is used only up to $n_{0}$. Similar to Fig. 9, we present current
constraints on $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ along with
hypothetical measurements with half the uncertainty and with vanishingly small
uncertainty for $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$. Again, while our
knowledge of $L$ improves with better measurements of
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$, the inferred posteriors for
$R_{1.4}$ and $\Lambda_{1.4}$ are nearly unaffected.111Reference Biswas (2021)
finds that improved measurements of $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$
can reduce the uncertainty in $R_{1.4}$. We attribute the apparent improvement
to correlations introduced by modeling choices made in Ref. Biswas (2021)
(e.g., the extent of the “low-density” nuclear parameterization and the
polytropic extension to higher densities) that are not introduced within our
nonparametric analysis. As elsewhere (e.g., Reed  _et al._ (2021)), reduced
uncertainty in NS observables from improved measurements at densities at or
below nuclear saturation are contingent upon specific model assumptions that
may not be correct. In fact, $L$ seems to be particularly uncorrelated with
$\Lambda_{1.4}$ within nonparametric extensions, implying that even a perfect
measurement of $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ additionally requires
reliable nuclear theory calculations to higher densities to impact our
expectations for future GW observations.

##  V Further Discussion

Finally, we discuss possible future areas of improvement and their expected
impact, from the assumptions made about the crust EOS, the different neutron
matter calculations, translations from pure neutron matter to matter in
$\beta$-equilibrium, and the likelihood modeling. We also briefly discuss
additional experimental probes of $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$.

Although we follow the uncertainty of individual $\chi\mathrm{EFT}$
calculations down to very low densities $n\leq 0.3n_{0}$, we match all EOS
draws to a single BPS crust model Baym  _et al._ (1971) below that. Previous
work suggested that the uncertainty in the crust at densities below $\simeq
10^{14}\mathrm{g}/\mathrm{cm}^{3}=0.36n_{0}$ can lead to a $\leq
0.3\,\mathrm{km}$ change in the radii of typical NSs Gamba  _et al._ (2019).
This effect is smaller than our current uncertainty in, e.g., $R_{1.4}$ at the
90% level ($12.39^{+1.02}_{-1.46}$ km), but it may not be negligible. However,
our results are qualitatively and quantitatively similar to the results of
Ref. Essick  _et al._ (2020b), which used $\chi\mathrm{EFT}$ uncertainties
down to similarly low densities but connected to a different crust model (SLy
Douchin and Haensel (2001)), as well as Refs. Essick  _et al._ (2020a); Landry
_et al._ (2020), which directly marginalized over 3 different crust EOSs (from
SLy, ENG Engvik  _et al._ (1996) and HQC18 Baym  _et al._ (2019)). Therefore,
any uncertainty within the crust model appears to have a minimal impact on our
results.

In our work, we explore 4 different $\chi\mathrm{EFT}$ calculations. These
explore $\chi\mathrm{EFT}$ interactions at different orders, employ different
local and nonlocal regularization schemes, and use different many-body methods
for the calculation of neutron matter. The PNM results are then extended to
matter in $\beta$-equilibrium, containing a small fraction of protons and
electrons around saturation densities. We emphasize here that for our
inference of nuclear matter properties, we focus on densities around $n_{0}$.
This enables the use of expansions around the empirical saturation point.
These expansions need to be truncated, but this approximation has a negligible
effect for the density expansion, again due to the focus on properties at or
around $n_{0}$. In the asymmetry expansion, the truncated higher-order terms
beyond $\mathcal{O}(x^{2})$ are estimated to be sub-MeV corrections around
$n_{0}$, and can be safely neglected given current EOS uncertainties Drischler
_et al._ (2014); Somasundaram  _et al._ (2021). Nonetheless, this could be
improved by future calculations of asymmetric matter around saturation
density.

We also note that several approaches to neutron matter calculations and their
associated uncertainties exist (see, e.g., discussion in Ref. Drischler  _et
al._ (2021a)). Our goal in this work was to span a range of different
$\chi\mathrm{EFT}$ calculations instead of attempting to quantify the errors
or term-by-term convergence within each individual calculation; thus our
choice to marginalize over separate $\chi\mathrm{EFT}$ estimates. As such, we
took the “best” constraint from each calculation instead of, e.g., considering
multiple orders within the same calculation (as, e.g., in Ref. Drischler  _et
al._ (2021b)). While our marginalization renders our conclusions robust and
tends to emphasize general trends, future work searching for astrophysical
evidence for, e.g. the breakdown scale within $\chi\mathrm{EFT}$ calculations,
will benefit from explicitly checking term-by-term convergence within
individual calculations against astrophysical data and further exploring the
effects of regulator artefacts.

While our results suggest that higher-order chiral interactions might be
important (compare N2LO $\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$
calculations with all other calculations that employ some N3LO contributions)
and that locally regularized interactions are less favored (again, compare
$\mathrm{QMC}^{(2016)}_{\mathrm{N^{2}LO}}$ to other calculations) we stress
that all $\chi\mathrm{EFT}$ calculations are consistent with each other and
that our conclusions about consistency with nuclear experiment and
astrophysical observations apply equally to all four $\chi\mathrm{EFT}$
calculations. This highlights the robustness of our findings.

Additionally, one may be concerned with the single-event likelihood models
constructed within our hierarchical inference. We use optimized Gaussian KDEs
(see Sec. II.1), which have previously been shown to robustly model the
associated likelihoods (see, e.g., discussion within Ref. Essick  _et al._
(2020a)). Indeed, while KDEs are known to be biased approximations to
probability densities, these effects are small given the current sample sizes
available within public posterior samples for each astrophysical observation
we consider. As Ref. Essick  _et al._ (2020a) discussed, we primarily expect
these to impact our estimate of the evidence that a particular object was a BH
rather than a NS (due to the sharp boundary at $\Lambda=0$ within GW
likelihoods). We do not consider such an inference here, and therefore expect
our KDE models to suffice for the task at hand. Similar to Refs. Landry and
Essick (2019); Essick  _et al._ (2020a), we also confirm that we retain large
effective numbers of samples throughout all stages of our Monte Carlo
inference scheme (typically, $\gtrsim\mathcal{O}(10^{4})$ effective EOS
samples for our nonparametric and $\chi\mathrm{EFT}$-marginalized results).
Nevertheless, it is worth noting that other approaches to modeling single-
event likelihoods exist in the literature (e.g., Ref. Hernandez Vivanco  _et
al._ (2019)) which may be of increasing importance with larger numbers of
astrophysical observations.

Similarly, marginal likelihoods from astrophysical observations implicitly
depend on the mass distributions assumed. Although the impact of our current
assumptions is expected to be small for the existing set of events, larger
sample sizes may require simultaneous inference of the NS mass distribution
and the EOS, e.g. Wysocki  _et al._ (2020); Golomb and Talbot (2021).

Finally, in addition to the approach using weak probes employed by PREX, and
the strong correlation with the dipole polarizability from $(p,p^{\prime})$
scattering, there are other experiments sensitive to
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ that rely on strong probes, see,
e.g., the reviews Thiel  _et al._ (2019) and Tsang  _et al._ (2012). While
here we have focused on the recent PREX result, and also explored
$\alpha_{D}^{{}^{208}\mathrm{Pb}}$ due to its well studied strong correlation
with $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$, we note that many of the
measurements of $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ that employ strong
probes tend to agree more closely with our $\chi\mathrm{EFT}$ priors, similar
to the $\alpha_{D}^{{}^{208}\mathrm{Pb}}$ results we consider. For example,
Ref. Tarbert  _et al._ (2014) estimates
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}=0.15\pm
0.03\,(\mathrm{stat.})^{+0.01}_{-0.03}\,(\mathrm{sys.})\,\mathrm{fm}$ based on
coherent pion production, and Ref. Jastrzȩbski  _et al._ (2004) estimates
$0.15\pm 0.02\,(\mathrm{stat.})\,\mathrm{fm}$ based on analyses of
antiprotonic atoms. While we do not explicitly consider these in our analysis
because of the difficulty in estimating the associated model systematics,
future analyses may include them if the model dependence implicit within the
experimental results is better understood.

##  VI Summary

In summary, we used nonparametric EOS inference to constrain the symmetry
energy, its density dependence, and $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$
directly from astrophysical data, leading to
$S_{0}=35.1^{+11.6}_{-8.9}\,\,\text{MeV}$, $L=58^{+61}_{-56}\,\,\text{MeV}$,
and
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}=0.19^{+0.12}_{-0.11}\,\mathrm{fm}$.
Folding in $\chi\mathrm{EFT}$ constraints reduces these ranges to
$S_{0}=32.7^{+1.9}_{-1.8}\,\,\text{MeV}$, $L=49^{+14}_{-15}\,\,\text{MeV}$,
and
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}=0.17^{+0.04}_{-0.04}\,\mathrm{fm}$.
While these results prefer values below the ones that PREX-II recently
reported Adhikari  _et al._ (2021); Reed  _et al._ (2021), the PREX-II
uncertainties are still broad and any tension is very mild. Furthermore, our
findings are in good agreement with other nuclear physics information. Our
analysis suggests that a future measurement of
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ with an uncertainty of $\pm
0.04\,\mathrm{fm}$ (a factor of $\simeq 2$ smaller than the current
uncertainty) could challenge current $\chi\mathrm{EFT}$ calculations, although
the tension with astrophysical data would still be relatively mild ($p$-value
of $11.5\%$). However, we also note that the formation of light clusters at
the surface of heavy nuclei could affect the extracted L value Tanaka  _et
al._ (2021).

Finally, our results demonstrate that the correlation between $R_{1.4}$ and
$L$ (or $R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$) is looser than suggested by
analyses based on a specific class of EOS models. In fact, even a
hypothetically perfect measurements of
$R_{\mathrm{skin}}^{{}^{208}\mathrm{Pb}}$ will not strongly impact our
knowledge of the radius and tidal deformability of $1.4\,M_{\odot}$ NSs when
using nonparametric EOS representations. The inverse is also true for such
EOSs: observations of NSs at astrophysically relevant masses will carry only
limited information about nuclear interactions at or below nuclear saturation
density. Extrapolating neutron-skin thickness measurements to NS scales thus
requires a careful treatment of systematic EOS model uncertainties to
distinguish implicit modeling assumptions from the data’s impact. In
particular, we find that the PREX-II data does not require NSs to have large
radii. However, if the high $L$ values of PREX-II persist, this may suggest a
peak in the sound speed around saturation density in order to accommodate both
the moderate radii inferred from astrophysical data and the large $L$ observed
in terrestrial experiments. Although tantalizing, it remains to be seen
whether astrophysical observations of low-mass NSs or future nuclear
experiments will bear this out.

Finally, we note that a confirmation of high values for $S_{0}$ and $L$
implied by the central PREX-II results would challenge all available
microscopic models for nuclear interactions (see, e.g., Refs. Lattimer and Lim
(2013); Tews  _et al._ (2017); Huth  _et al._ (2021); Drischler  _et al._
(2021a)). This affects both phenomenological two- and three-nucleon potentials
as well as interactions derived from $\chi\mathrm{EFT}$, and would require a
significant increase of the repulsion between neutrons at densities of the
order of $n_{0}$. This would have direct implications for studies of the
structure of medium-mass to heavy nuclei.

###### Acknowledgements.

R.E. was supported by the Perimeter Institute for Theoretical Physics and the
Kavli Institute for Cosmological Physics. R.E. also 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. The
Kavli Institute for Cosmological Physics at the University of Chicago is
supported by an endowment from the Kavli Foundation and its founder Fred
Kavli. P.L. is supported by National Science Foundation award PHY-1836734 and
by a gift from the Dan Black Family Foundation to the Nicholas & Lee Begovich
Center for Gravitational-Wave Physics & Astronomy. The work of A.S. was
supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) – Project-ID 279384907 – SFB 1245. 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 numbers 20190617PRD1 and 20190021DR, 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. This
work benefited from discussions within IReNA, which is supported in part by
the National Science Foundation under Grant No. OISE-1927130. The authors also
gratefully acknowledge the computational resources provided by the LIGO
Laboratory and supported by NSF grants PHY-0757058 and PHY-0823459.
Computational resources have also been provided by the Los Alamos National
Laboratory Institutional Computing Program, which is supported by the U.S.
Department of Energy National Nuclear Security Administration under Contract
No. 89233218CNA000001, and by the National Energy Research Scientific
Computing Center (NERSC), which is supported by the U.S. Department of Energy,
Office of Science, under contract No. DE-AC02-05CH11231.

## References

  * Lattimer and Lim (2013) J. M. Lattimer and Y. Lim, Astrophys. J. 771, 51 (2013), arXiv:1203.4286 [nucl-th] .
  * Tsang  _et al._ (2012) M. B. Tsang, J. R. Stone, F. Camera, P. Danielewicz, S. Gandolfi, K. Hebeler, C. J. Horowitz, J. Lee, W. G. Lynch, Z. Kohley, _et al._ , Phys. Rev. C 86, 015803 (2012), arXiv:1204.0466 [nucl-ex] .
  * Huth  _et al._ (2021) S. Huth, C. Wellenhofer, and A. Schwenk, Phys. Rev. C 103, 025803 (2021), arXiv:2009.08885 [nucl-th] .
  * Typel and Brown (2001) S. Typel and B. A. Brown, Phys. Rev. C 64, 027302 (2001).
  * Viñas  _et al._ (2014) X. Viñas, M. Centelles, X. Roca-Maza, and M. Warda, Eur. Phys. J. A 50, 27 (2014), arXiv:1308.1008 [nucl-th] .
  * Reinhard and Nazarewicz (2016) P. G. Reinhard and W. Nazarewicz, Phys. Rev. C 93, 051303 (2016), arXiv:1601.06324 [nucl-th] .
  * Mondal  _et al._ (2016) C. Mondal, B. K. Agrawal, M. Centelles, G. Colò, X. Roca-Maza, N. Paar, X. Viñas, S. K. Singh, and S. K. Patra, Phys. Rev. C 93, 064303 (2016), arXiv:1605.05048 [nucl-th] .
  * Tamii  _et al._ (2011) A. Tamii, I. Poltoratska, P. von Neumann-Cosel, Y. Fujita, T. Adachi, C. A. Bertulani, J. Carter, M. Dozono, H. Fujita, K. Fujita, _et al._ , Phys. Rev. Lett. 107, 062502 (2011), arXiv:1104.5431 [nucl-ex] .
  * Piekarewicz  _et al._ (2012) J. Piekarewicz, B. K. Agrawal, G. Colo, W. Nazarewicz, N. Paar, P. G. Reinhard, X. Roca-Maza, and D. Vretenar, Phys. Rev. C 85, 041302 (2012), arXiv:1201.3807 [nucl-th] .
  * Roca-Maza  _et al._ (2013) X. Roca-Maza, M. Brenna, G. Colò, M. Centelles, X. Viñas, B. K. Agrawal, N. Paar, D. Vretenar, and J. Piekarewicz, Phys. Rev. C 88, 024316 (2013), arXiv:1307.4806 [nucl-th] .
  * Roca-Maza  _et al._ (2015) X. Roca-Maza, X. Viñas, M. Centelles, B. K. Agrawal, G. Colò, N. Paar, J. Piekarewicz, and D. Vretenar, Phys. Rev. C 92, 064304 (2015), arXiv:1510.01874 [nucl-th] .
  * Lattimer and Prakash (2001) J. M. Lattimer and M. Prakash, Astrophys. J. 550, 426 (2001), arXiv:astro-ph/0002232 .
  * Steiner and Gandolfi (2012) A. W. Steiner and S. Gandolfi, Phys. Rev. Lett. 108, 081102 (2012), arXiv:1110.4142 [nucl-th] .
  * Neill  _et al._ (2021) D. Neill, W. G. Newton, and D. Tsang, Mon. Not. Roy. Astron. Soc. 504, 1129 (2021), arXiv:2012.10322 [astro-ph.HE] .
  * Russotto  _et al._ (2016) P. Russotto, S. Gannon, S. Kupny, P. Lasko, L. Acosta, M. Adamczyk, A. Al-Ajlan, M. Al-Garawi, S. Al-Homaidhi, F. Amorini, _et al._ , Phys. Rev. C 94, 034608 (2016), arXiv:1608.04332 [nucl-ex] .
  * Tews  _et al._ (2013) I. Tews, T. Krüger, K. Hebeler, and A. Schwenk, Phys. Rev. Lett. 110, 032504 (2013), arXiv:1206.0025 [nucl-th] .
  * Hebeler  _et al._ (2013) K. Hebeler, J. Lattimer, C. Pethick, and A. Schwenk, Astrophys. J. 773, 11 (2013), arXiv:1303.4662 [astro-ph.SR] .
  * Drischler  _et al._ (2016) C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev. C 93, 054314 (2016), arXiv:1510.06728 [nucl-th] .
  * Lonardoni  _et al._ (2020) D. Lonardoni, I. Tews, S. Gandolfi, and J. Carlson, Phys. Rev. Res. 2, 022033 (2020), arXiv:1912.09411 [nucl-th] .
  * Drischler  _et al._ (2020a) C. Drischler, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, Phys. Rev. Lett. 125, 202702 (2020a), arXiv:2004.07232 [nucl-th] .
  * Adhikari  _et al._ (2021) D. Adhikari, H. Albataineh, D. Androic, K. Aniol, D. S. Armstrong, T. Averett, C. Ayerbe Gayoso, S. Barcus, V. Bellini, R. S. Beminiwattha, _et al._ (PREX), Phys. Rev. Lett. 126, 172502 (2021), arXiv:2102.10767 [nucl-ex] .
  * 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] .
  * Reinhard  _et al._ (2021) P.-G. Reinhard, X. Roca-Maza, and W. Nazarewicz, (2021), arXiv:2105.15050 [nucl-th] .
  * Essick  _et al._ (2021) R. Essick, I. Tews, P. Landry, and A. Schwenk, (2021), arXiv:2102.10074 [nucl-th] .
  * 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] .
  * Alford  _et al._ (2013) M. G. Alford, S. Han, and M. Prakash, Phys. Rev. D 88, 083013 (2013), arXiv:1302.4732 [astro-ph.SR] .
  * Tews  _et al._ (2018) I. Tews, J. Margueron, and S. Reddy, Phys. Rev. C 98, 045804 (2018), arXiv:1804.02783 [nucl-th] .
  * Greif  _et al._ (2019) S. Greif, G. Raaijmakers, K. Hebeler, A. Schwenk, and A. Watts, Mon. Not. Roy. Astron. Soc. 485, 5363 (2019), arXiv:1812.08188 [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] .
  * Lindblom (2010) L. Lindblom, Phys. Rev. D 82, 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] .
  * Alam  _et al._ (2016) N. Alam, B. K. Agrawal, M. Fortin, H. Pais, C. Providência, A. R. Raduta, and A. Sulaksono, Phys. Rev. C 94, 052801 (2016), arXiv:1610.06344 [nucl-th] .
  * Carson  _et al._ (2019) Z. Carson, A. W. Steiner, and K. Yagi, Phys. Rev. D 99, 043010 (2019), arXiv:1812.08910 [gr-qc] .
  * Biswas  _et al._ (2020) B. Biswas, P. Char, R. Nandi, and S. Bose, (2020), arXiv:2008.01582 [astro-ph.HE] .
  * Yue  _et al._ (2021) T.-G. Yue, L.-W. Chen, Z. Zhang, and Y. Zhou, (2021), arXiv:2102.05267 [nucl-th] .
  * Baym  _et al._ (1971) G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. 170, 299 (1971).
  * Landry  _et al._ (2020) P. Landry, R. Essick, and K. Chatziioannou, Phys. Rev. D 101, 123007 (2020), arXiv:2003.04880 [astro-ph.HE] .
  * 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, _et al._ , Science 340, 448 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Cromartie  _et al._ (2020) H. T. Cromartie, E. Fonseca, S. M. Ransom, P. B. Demorest, Z. Arzoumanian, H. Blumer, P. R. Brook, M. E. DeCesar, T. Dolch, J. A. Ellis, _et al._ , Nature Astronomy 4, 72 (2020), arXiv:1904.06759 [astro-ph.HE] .
  * Abbott  _et al._ (2017a) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017a).
  * Abbott  _et al._ (2019) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, _et al._ , Phys. Rev. X 9, 011001 (2019), arXiv:1805.11579 [gr-qc] .
  * Aasi  _et al._ (2015) J. Aasi, B. P. Abbott, R. Abbott, T. Abbott, M. R. Abernathy, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, _et al._ (LIGO Scientific Collaboration), Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ (2015) F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, A. Allocca, J. Amarni, P. Astone, G. Balestri, G. Ballardin, _et al._ (Virgo Collaboration), Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Miller  _et al._ (2019) 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, _et al._ , Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * 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, _et al._ , Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [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, _et al._ , (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Miller  _et al._ (2019) M. C. Miller, C. Chirenti, and F. K. Lamb, Astrophysical J. 888, 12 (2019), 1904.08907 [astro-ph.HE] .
  * Abbott  _et al._ (2017b) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, _et al._ , Astrophysical J. 850, L39 (2017b), arXiv:1710.05836 [astro-ph.HE] .
  * Coughlin  _et al._ (2018) M. W. Coughlin, T. Dietrich, Z. Doctor, D. Kasen, S. Coughlin, A. Jerkstrand, G. Leloudas, O. McBrien, B. D. Metzger, R. O’Shaughnessy, _et al._ , Mon. Not. Roy. Astron. Soc. 480, 3871 (2018), arXiv:1805.09371 [astro-ph.HE] .
  * Kasen and Barnes (2019) D. Kasen and J. Barnes, Astrophys. J. 876, 128 (2019), arXiv:1807.03319 [astro-ph.HE] .
  * Siegel (2019) D. M. Siegel, Eur. Phys. J. A 55, 203 (2019), arXiv:1901.09044 [astro-ph.HE] .
  * Metzger (2020) B. D. Metzger, Living Rev. Rel. 23, 1 (2020), arXiv:1910.01617 [astro-ph.HE] .
  * Dietrich  _et al._ (2020) T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Science 370, 1450 (2020), arXiv:2002.11355 [astro-ph.HE] .
  * Stachie  _et al._ (2021) C. Stachie, M. W. Coughlin, T. Dietrich, S. Antier, M. Bulla, N. Christensen, R. Essick, P. Landry, B. Mours, F. Schianchi, _et al._ , Mon. Not. Roy. Astron. Soc. 505, 4235 (2021), arXiv:2103.01733 [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, Astrophys. J. Lett. 918, L29 (2021), arXiv:2105.06981 [astro-ph.HE] .
  * Miller  _et al._ (2021) M. C. Miller _et al._ , (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Riley  _et al._ (2021) T. E. Riley _et al._ , (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Pang  _et al._ (2021) P. T. H. Pang, I. Tews, M. W. Coughlin, M. Bulla, C. Van Den Broeck, and T. Dietrich, (2021), arXiv:2105.08688 [astro-ph.HE] .
  * Legred  _et al._ (2021) I. Legred, K. Chatziioannou, R. Essick, S. Han, and P. Landry, (2021), arXiv:2106.05313 [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] .
  * 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] .
  * Drischler  _et al._ (2019) C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev. Lett. 122, 042501 (2019), arXiv:1710.08220 [nucl-th] .
  * Hebeler and Schwenk (2010) K. Hebeler and A. Schwenk, Phys. Rev. C 82, 014314 (2010), arXiv:0911.0483 [nucl-th] .
  * Drischler  _et al._ (2014) C. Drischler, V. Soma, and A. Schwenk, Phys. Rev. C 89, 025806 (2014), arXiv:1310.5627 [nucl-th] .
  * Somasundaram  _et al._ (2021) R. Somasundaram, C. Drischler, I. Tews, and J. Margueron, Phys. Rev. C 103, 045803 (2021), arXiv:2009.04737 [nucl-th] .
  * Goriely  _et al._ (2013) S. Goriely, N. Chamel, and J. M. Pearson, Phys. Rev. C 88, 024308 (2013).
  * Chamel and Haensel (2008) N. Chamel and P. Haensel, Living Rev. Rel. 11, 10 (2008), arXiv:0812.3955 [astro-ph] .
  * Lattimer and Prakash (2016) J. M. Lattimer and M. Prakash, Phys. Rept. 621, 127 (2016), arXiv:1512.07820 [astro-ph.SR] .
  * Tews  _et al._ (2017) I. Tews, J. M. Lattimer, A. Ohnishi, and E. E. Kolomeitsev, Astrophys. J. 848, 105 (2017), arXiv:1611.07133 [nucl-th] .
  * Drischler  _et al._ (2020b) C. Drischler, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, Phys. Rev. Lett. 125, 202702 (2020b).
  * Drischler  _et al._ (2020c) C. Drischler, J. A. Melendez, R. J. Furnstahl, and D. R. Phillips, Phys. Rev. C 102, 054315 (2020c), arXiv:2004.07805 [nucl-th] .
  * Biswas (2021) B. Biswas, (2021), arXiv:2105.02886 [astro-ph.HE] .
  * Gamba  _et al._ (2019) R. Gamba, J. S. Read, and L. E. Wade, Class. Quant. Grav. 37, 025008 (2019), arXiv:1902.04616 [gr-qc] .
  * Douchin and Haensel (2001) F. Douchin and P. Haensel, Astronomy & Astrophysics 380, 151 (2001), arXiv:astro-ph/0111092 [astro-ph] .
  * Engvik  _et al._ (1996) L. Engvik, E. Osnes, M. Hjorth-Jensen, G. Bao, and E. Ostgaard, Astrophys. J. 469, 794 (1996), arXiv:nucl-th/9509016 [nucl-th] .
  * Baym  _et al._ (2019) G. Baym, S. Furusawa, T. Hatsuda, T. Kojo, and H. Togashi, Astrophys. J. 885, 42 (2019), arXiv:1903.08963 [astro-ph.HE] .
  * Drischler  _et al._ (2021a) C. Drischler, J. W. Holt, and C. Wellenhofer, Annu. Rev. Nucl. Part. Sci. 71, 403 (2021a), arXiv:2101.01709 .
  * 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] .
  * Hernandez Vivanco  _et al._ (2019) F. Hernandez Vivanco, R. Smith, E. Thrane, P. D. Lasky, C. Talbot, and V. Raymond, Phys. Rev. D 100, 103009 (2019), arXiv:1909.02698 [gr-qc] .
  * Wysocki  _et al._ (2020) D. Wysocki, R. O’Shaughnessy, L. Wade, and J. Lange, (2020), arXiv:2001.01747 [gr-qc] .
  * Golomb and Talbot (2021) J. Golomb and C. Talbot, (2021), arXiv:2106.15745 [astro-ph.HE] .
  * Thiel  _et al._ (2019) M. Thiel, C. Sfienti, J. Piekarewicz, C. J. Horowitz, and M. Vanderhaeghen, J. Phys. G 46, 093003 (2019).
  * Tarbert  _et al._ (2014) C. M. Tarbert, D. P. Watts, D. I. Glazier, P. Aguar, J. Ahrens, J. R. M. Annand, H. J. Arends, R. Beck, V. Bekrenev, B. Boillat, A. Braghieri, D. Branford, W. J. Briscoe, J. Brudvik, S. Cherepnya, R. Codling, E. J. Downie, K. Foehl, P. Grabmayr, R. Gregor, E. Heid, D. Hornidge, O. Jahn, V. L. Kashevarov, A. Knezevic, R. Kondratiev, M. Korolija, M. Kotulla, D. Krambrich, B. Krusche, M. Lang, V. Lisin, K. Livingston, S. Lugert, I. J. D. MacGregor, D. M. Manley, M. Martinez, J. C. McGeorge, D. Mekterovic, V. Metag, B. M. K. Nefkens, A. Nikolaev, R. Novotny, R. O. Owens, P. Pedroni, A. Polonski, S. N. Prakhov, J. W. Price, G. Rosner, M. Rost, T. Rostomyan, S. Schadmand, S. Schumann, D. Sober, A. Starostin, I. Supek, A. Thomas, M. Unverzagt, T. Walcher, L. Zana, and F. Zehr (Crystal Ball at MAMI and A2 Collaboration), Phys. Rev. Lett. 112, 242502 (2014).
  * Jastrzȩbski  _et al._ (2004) J. Jastrzȩbski, A. Trzcińska, P. Lubiński, B. Kłos, F. J. Hartmann, T. von Egidy, and S. Wycech, Int. J. Mod. Phys. E 13, 343 (2004).
  * Tanaka  _et al._ (2021) J. Tanaka, Z. Yang, S. Typel, S. Adachi, S. Bai, P. van Beek, D. Beaumel, Y. Fujikawa, J. Han, S. Heil, _et al._ , Science 371, 260 (2021).