# Constraints on the dense matter equation of state and neutron star
properties from
NICER’s mass-radius estimate of PSR J0740$+$6620 and multimessenger
observations

G. Raaijmakers GRAPPA, Anton Pannekoek Institute for Astronomy and Institute
of High-Energy Physics, University of Amsterdam, Science Park 904, 1098 XH
Amsterdam, Netherlands  S. K. Greif Technische Universität Darmstadt,
Department of Physics, 64289 Darmstadt, Germany  ExtreMe Matter Institute
EMMI, GSI Helmholtzzentrum für Schwerionenforschung GmbH, 64291 Darmstadt,
Germany  K. Hebeler Technische Universität Darmstadt, Department of Physics,
64289 Darmstadt, Germany  ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum
für Schwerionenforschung GmbH, 64291 Darmstadt, Germany  T. Hinderer Institute
for Theoretical Physics, Utrecht University, Princetonplein 5, 3584CC Utrecht,
The Netherlands  S. Nissanke GRAPPA, Anton Pannekoek Institute for Astronomy
and Institute of High-Energy Physics, University of Amsterdam, Science Park
904, 1098 XH Amsterdam, Netherlands  Nikhef, Science Park 105, 1098 XG
Amsterdam, The Netherlands  A. Schwenk 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  T. E. Riley Anton Pannekoek Institute for Astronomy,
University of Amsterdam, Science Park 904, 1090GE Amsterdam, the Netherlands
A. L. Watts Anton Pannekoek Institute for Astronomy, University of Amsterdam,
Science Park 904, 1090GE Amsterdam, the Netherlands  J. M. Lattimer
Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY
11794-3800, USA  W. C. G. Ho Department of Physics and Astronomy, Haverford
College, 370 Lancaster Avenue, Haverford, PA 19041, USA

###### Abstract

In recent years our understanding of the dense matter equation of state (EOS)
of neutron stars has significantly improved by analyzing multimessenger data
from radio/X-ray pulsars, gravitational wave events, and from nuclear physics
constraints. Here we study the additional impact on the EOS from the jointly
estimated mass and radius of PSR J0740$+$6620, presented in Riley et al.
(2021) by analyzing a combined dataset from X-ray telescopes NICER and XMM-
Newton. We employ two different high-density EOS parameterizations: a
piecewise-polytropic (PP) model and a model based on the speed of sound in a
neutron star (CS). At nuclear densities these are connected to microscopic
calculations of neutron matter based on chiral effective field theory
interactions. In addition to the new NICER data for this heavy neutron star,
we separately study constraints from the radio timing mass measurement of PSR
J0740$+$6620, the gravitational wave events of binary neutron stars GW190425
and GW170817, and for the latter the associated kilonova AT2017gfo. By
combining all these, and the NICER mass-radius estimate of PSR J0030$+$0451,
we find the radius of a $1.4\,$M⊙ neutron star to be constrained to the 95%
credible ranges $12.33^{+0.76}_{-0.81}\,$km (PP model) and
$12.18^{+0.56}_{-0.79}\,$km (CS model). In addition, we explore different
chiral effective field theory calculations and show that the new NICER results
provide tight constraints for the pressure of neutron star matter at around
twice saturation density, which shows the power of these observations to
constrain dense matter interactions at intermediate densities.

dense matter — equation of state — stars: neutron — X-rays: stars —
gravitational waves

††software: Python/C language (Oliphant, 2007), GNU Scientific Library (GSL;
Gough, 2009), NumPy (van der Walt et al., 2011), Cython (Behnel et al., 2011),
SciPy (Virtanen et al., 2020), MPI (Forum, 1994), MPI for Python (Dalcín et
al., 2008), Matplotlib (Hunter, 2007; Droettboom et al., 2018), IPython (Perez
& Granger, 2007), Jupyter (Kluyver et al., 2016), MultiNest (Feroz et al.,
2009), PyMultiNest (Buchner et al., 2014), kalepy (Kelley, 2021).

##  1 Introduction

Our understanding of the dense matter equation of state (EOS) of neutron stars
has made significant progress over the last few years due to the arrival of
new avenues to measure observables like mass, radius and tidal deformability,
that connect to the behavior of matter at supranuclear densities. Recently
NASA’s X-ray timing telescope, the Neutron Star Interior Composition Explorer
(NICER), has delivered the first joint measurement of mass and radius through
pulse profile modeling of the millisecond pulsar PSR J0030$+$0451 (Riley et
al., 2019; Miller et al., 2019). The impact of this measurement on the dense
matter EOS has been extensively studied in various EOS frameworks (see, e.g.,
Raaijmakers et al., 2019; Miller et al., 2019; Raaijmakers et al., 2020;
Essick et al., 2020b; Landry et al., 2020; Dietrich et al., 2020; Jiang et
al., 2020; Al-Mamun et al., 2021), including EOS with phase transitions to
quark matter (see, e.g., Xie & Li, 2021; Li et al., 2020; Tang et al., 2021;
Blaschke et al., 2020; Alvarez-Castillo et al., 2020) and models that explore
the possibility of there being two stable neutron star branches (Christian &
Schaffner-Bielich, 2020).

Concurrently, the second and third observing runs of LIGO/Virgo have so far
resulted in the confirmed gravitational wave detections of two (most-likely)
binary neutron star mergers: GW170817 (Abbott et al., 2017c, 2019a) and
GW190425 (Abbott et al., 2020a). By accurately measuring the gravitational
wave phase, limits can be put on the EOS-dependent tidal deformability of the
neutron stars (Flanagan & Hinderer, 2008; Hinderer et al., 2010). While for
GW170817 the tidal deformability could be measured within a $90\%$ highest
posterior density interval when adopting low spin priors (see, e.g. Abbott et
al., 2018, 2019b), the low signal-to-noise ratio (SNR) of GW190425 resulted in
only weak upper limits on the tidal deformability even when assuming low spins
(Abbott et al., 2020a). We consider the $\sim 2.6$ M⊙ secondary object in
GW190814 (Abbott et al., 2020d) to be a black hole (Nathanail et al., 2021),
and will therefore not use this third event in our analysis.

At nuclear densities, the EOS is well constrained by nuclear theory and
experiments (see, e.g., Tsang et al., 2012; Lattimer & Lim, 2013; Huth et al.,
2021). In particular, many-body calculations based on chiral effective field
theory (EFT) interactions have enabled systematic predictions for the neutron
matter EOS up to nuclear saturation density including theoretical
uncertainties (see, e.g., Hebeler et al., 2013; Tews et al., 2013; Lynn et
al., 2016; Drischler et al., 2019, 2020). Up to saturation density, the
resulting symmetry energy and pressure of neutron matter are also consistent
with extractions from nuclear experiments (Lattimer & Lim, 2013), including
from measurements of the dipole polarizability of neutron-rich nuclei (Roca-
Maza et al., 2015; Birkhan et al., 2017; Kaufmann et al., 2020). Taking these
results at nuclear densities, combined with standard crust EOS, different
extrapolations to high densities have been found to lead to NS radii
consistent with all multimessenger observations (see, e.g., Raaijmakers et
al., 2020; Essick et al., 2020b; Annala et al., 2020; Dietrich et al., 2020;
Biswas et al., 2021). Recently, the results of PREX-II have pointed to higher
pressures (Adhikari et al., 2021; Reed et al., 2021), but with very large
uncertainties, so that in a combined analysis with astrophysical and chiral
EFT constraints, the overall consistency still persists (Essick et al., 2021).

NICER data has now enabled a joint estimate of the mass and radius of the
high-mass rotation-powered millisecond pulsar PSR J0740$+$6620\. Since PSR
J0740$+$6620 (unlike PSR J0030$+$0451) is in a binary with an inclination that
allows measurement of the Shapiro delay, its mass can be measured
independently via radio timing. Cromartie et al. (2020) reported a mass of
$2.14^{+0.10}_{-0.09}$ M⊙, and a joint campaign by the North American
Nanohertz Observatory for Gravitational Waves (NANOGrav) and the Canadian
Hydrogen Intensity Mapping Experiment (CHIME)/Pulsar collaborations has now
resulted in an updated mass of $2.08\pm 0.07$ M⊙ (Fonseca et al., 2021).

Riley et al. (2021) have used this mass measurement as an informative prior
for pulse-profile modeling analysis that is joint over the phase-resolved
spectroscopic data from NICER and phase-averaged data from the XMM-Newton
European Photon Imaging Camera (EPIC). The inclusion of the smaller XMM-Newton
(hereafter XMM) data set allows for better constraints on the proportion of
the X-ray emission that is attributable to background rather than PSR
J0740$+$6620, ultimately acting to cut out solutions with high compactness.
This results in an inferred radius of $12.39_{-0.98}^{+1.30}$ km, and a mass
of $2.072_{-0.066}^{+0.067}$ M⊙ that is little changed from the radio prior.
For a full description of the methodology employed in the mass-radius
inference we refer the reader to Riley et al. (2021).

In this Letter, we use the mass and radius from Riley et al. (2021) for PSR
J0740$+$6620 as input for inferring the dense matter EOS, combining it with
other constraints from nuclear theory and multi-messenger observations. It
should be considered as a follow-up to our previous work that built on NICER’s
results for PSR J0030$+$0451 (Raaijmakers et al., 2019, 2020), where in this
work we explore also a broader range of multi-messenger constraints. As the
high-density constraints from astrophysical observations get more precise,
with the new NICER results and future LIGO/Virgo measurements, it will be
intriguing to see them play out with the present nuclear constraints. In this
Letter, we also explore this for the new NICER results and how they constrain
the EOS above nuclear densities starting from different chiral EFT
calculations111The posterior samples and scripts to make the plots in this
Letter are available in a Zenodo repository at Raaijmakers et al. (2021a)..

##  2 Inference framework

In this work we will closely follow the analysis framework developed
previously in Greif et al. (2019), Raaijmakers et al. (2019) and Raaijmakers
et al. (2020). Below, we summarize this method and highlight several updates
to the framework.

We consider two EOS parameterizations: i) a piecewise polytropic (PP) model
with three segments between varying transitions densities (Hebeler et al.,
2013), and ii) a speed-of-sound (CS) model first introduced in Greif et al.
(2019). To capture the uncertainty in the EOS around nuclear saturation
density ($n_{0}=0.16\,\text{fm}^{-3}$), both parameterizations are matched to
a power law fit of a range of EOS calculated from chiral effective field
theory interactions (Hebeler & Schwenk, 2010; Hebeler et al., 2013) below
$1.1n_{0}$. At densities below $0.5n_{0}$ this power law fit is connected to
the BPS crust EOS (Baym et al., 1971).

Figure 1: Constraints on the mass-radius relation of neutron stars, given the
posterior distribution on EOS parameters $\bm{\theta}$ using the PP model
(left) and CS model (right panel). The constraints from the updated radio
timing mass of PSR J0740$+$6620 from Fonseca et al. (2021) (present work,
green) are compared to the mass from Cromartie et al. (2020) used in our
previous works (Raaijmakers et al., 2019, 2020) (orange, dashed-dotted),
showing both the $68\%$ and $95\%$ credible regions. The black dashed lines
indicate the $95\%$ credible region of the prior distribution. Note that the
slightly lower mass measurement does not have a significant impact on the EOS
posterior.

To constrain these EOS parameterizations, governed by the EOS parameters
$\bm{\theta}$, we employ Bayes’ theorem and write the posterior distributions
of the EOS parameters and central energy densities $\bm{\varepsilon}$ as

$p(\bm{\theta},\bm{\varepsilon}\,|\,\bm{d},\mathbb{M})\propto
p(\bm{\theta}\,|\,\mathbb{M})~{}p(\bm{\varepsilon}\,|\,\bm{\theta},\mathbb{M})~{}p(\bm{d}\,|\,\bm{\theta},\mathbb{M})\,,$
(1)

where $\mathbb{M}$ denotes the model including all assumed physics and
$\bm{d}$ the dataset used to constrain the EOS, consisting of, e.g., radio-,
X-ray and gravitational wave data. When assuming each of these datasets to be
independent of each other, we can separate the likelihoods and write

$\displaystyle p(\bm{\theta},\bm{\varepsilon}\,|\,$
$\displaystyle\bm{d},\mathbb{M})\propto
p(\bm{\theta}\,|\,\mathbb{M})~{}p(\bm{\varepsilon}\,|\,\bm{\theta},\mathbb{M})$
$\displaystyle\times\prod_{i}p(\Lambda_{1,i},\Lambda_{2,i},M_{1,i},M_{2,i}\,|\,\bm{d}_{\rm
GW,i}(,\bm{d}_{\rm EM,i}))$
$\displaystyle\times\prod_{j}p(M_{j},R_{j}\,|\,\bm{d}_{\rm NICER,j})$
$\displaystyle\times\prod_{k}p(M_{k}\,|\,\bm{d}_{\rm radio,k})\,.$ (2)

Here the products run over the number of different observed stars, or mergers,
in the case of the gravitational wave data. Furthermore, in Equation (2) we
have equated the nuisance-marginalized likelihoods to the nuisance-
marginalized posterior distributions derived in Riley et al. (2019); Fonseca
et al. (2021); Riley et al. (2021); Abbott et al. (2019a, 2020a). This
approximation is justifiable when the priors used in estimating these
nuisance-marginalized posterior distributions are uninformative, which for
simplicity we will assume to be a uniform prior in this case. The posterior
distributions derived by Riley et al. (2019) and Riley et al. (2021) already
use a jointly uniform prior in mass and radius. The posterior distributions
derived by Abbott et al. (2019a) and Abbott et al. (2020a) use a jointly
uniform prior in the tidal deformabilities of the two components $\Lambda_{i}$
within the range $\Lambda_{i}\subset[0,5000]$ (for GW190425 the upper bound of
$\Lambda_{2}$ was set to $10^{4}$.). The prior on the detector frame masses,
which are redshifted with respect to the source frame masses ($M_{\rm
det}=M_{i}(1+z)$), is uniform within the range $M_{\rm det}\subset[0.5,7.7]$
and $M_{\rm det}\subset[1,5.31]$ for GW170817 and GW190425 respectively.
However, the posterior distribution on component masses from gravitational
waves is highly degenerate because of the accurately measured chirp mass
$\mathcal{M}_{c}=(M_{1}M_{2})^{3/5}/(M_{1}+M_{2})^{1/5}$. To speed up the
convergence of our parameter estimation, we therefore transform the
gravitational wave posterior distributions to include the two tidal
deformabilities, chirp mass and mass ratio $q$, while reweighing such that the
prior distribution on these parameters is uniform. Further, we also fix the
chirp mass to its median value, since the small uncertainty in this parameter
does not affect the EOS parameter estimation (see Raaijmakers et al., 2020),
and thus have:

$\displaystyle p(\bm{\theta},\bm{\varepsilon}\,|\,$
$\displaystyle\bm{d},\mathbb{M})\propto
p(\bm{\theta}\,|\,\mathbb{M})~{}p(\bm{\varepsilon}\,|\,\bm{\theta},\mathbb{M})$
$\displaystyle\times\prod_{i}p(\Lambda_{1,i},\Lambda_{2,i},q_{i}\,|\,\mathcal{M}_{c},\bm{d}_{\rm
GW,i}(,\bm{d}_{\rm EM,i}))$
$\displaystyle\times\prod_{j}p(M_{j},R_{j}\,|\,\bm{d}_{\rm NICER,j})$
$\displaystyle\times\prod_{k}p(M_{k}\,|\,\bm{d}_{\rm radio,k})\,.$ (3)

Fixing the chirp mass means that the vector $\bm{\varepsilon}$ only contains
one central density per merger, where the tidal deformability of the second
component is now set by $\Lambda_{2}=\Lambda_{2}(\bm{\theta};q)$. If a
gravitational wave event has an associated electromagnetic (EM) counterpart,
the likelihood for that event becomes a product of the nuisance-marginalized
posterior distribution from the gravitational wave data and the nuisance-
marginalized posterior distribution from the EM analysis, such that:

$\displaystyle p(\Lambda_{1},\Lambda_{2},q\,|\,\mathcal{M}_{c},\bm{d}_{\rm
GW},$ $\displaystyle\bm{d}_{\rm EM})\propto
p(\Lambda_{1},\Lambda_{2},q\,|\,\mathcal{M}_{c},\bm{d}_{\rm GW})$
$\displaystyle\times
p(\Lambda_{1},\Lambda_{2},q\,|\,\mathcal{M}_{c},\bm{d}_{\rm EM})\,.$ (4)

Obtaining the posterior distribution
$p(\Lambda_{1},\Lambda_{2},q\,|\,\mathcal{M}_{c},\bm{d}_{\rm EM})$ is
discussed in Section 3.2.1 for the specific case of AT2017gfo, the kilonova
associated with GW170817 (see, e.g., Abbott et al., 2017a, b; Arcavi et al.,
2017; Coulter et al., 2017; Chornock et al., 2017; Cowperthwaite et al., 2017;
Kasliwal et al., 2017; Nicholl et al., 2017; Tanvir et al., 2017).

We then sample from the posterior distribution
$p(\bm{\theta},\bm{\varepsilon}\,|\,\bm{d},\mathbb{M})$, compute the
corresponding $M$, $R$, and $\Lambda$, and then evaluate the likelihood by
applying a kernel density estimation to the posterior distributions from Riley
et al. (2019, 2021); Abbott et al. (2019a, 2020a) using the nested sampling
software MultiNest. The same prior distribution
$p(\bm{\theta}\,|\,\mathbb{M})$ is used as in previous work; we refer the
reader to Section 2.3 of Raaijmakers et al. (2020) and references therein for
a more detailed description.

##  3 EOS constraints

In this Section we investigate the impact of the Riley et al. (2021) mass-
radius measurement for PSR J0740$+$6620 on the dense matter EOS, both
separately and when combined with previous constraints.

###  3.1 Radio mass measurement of PSR J0740$+$6620

Firstly, we constrain the EOS using the updated mass measurement of $2.08\pm
0.07$ M⊙ for PSR J0740$+$6620 derived using radio timing (Fonseca et al.,
2021), and compare this to the constraints from the previously published mass
of $2.14^{+0.1}_{-0.09}$ M⊙ (Cromartie et al., 2020). In Figure 1 we show the
posterior distribution on EOS parameters $\bm{\theta}$ when transformed to the
mass-radius parameter space. We note that, as expected, the slightly lower
updated mass measurement shifts the posterior distributions to lower maximum
neutron star masses and lower radii, although the effect is almost negligible.
Since the radio timing mass measurement is already incorporated in the joint
mass-radius estimate from NICER we will not use this measurement in the
remainder of this work.

###  3.2 GW170817 and GW190425

The gravitational wave events GW170817 and GW190425 have so far been the only
confirmed neutron star binary mergers during the recent observing runs of the
LIGO/Virgo collaboration (Abbott et al., 2020c). Although both events have a
non-negligible chance of being neutron star-black hole mergers (see, e.g.,
Yang et al. (2018); Ascenzi et al. (2019); Coughlin & Dietrich (2019);
Hinderer et al. (2019) for GW170817 and e.g., Kyutoku et al. (2020); Han et
al. (2020) for GW190425), in the following we will assume both objects to be
neutron stars. We use the low-spin222The low-spin assumption is chosen to be
consistent with measurements of spins in Galactic neutron star binaries that
merge within a Hubble time. posterior distributions on tidal deformability and
mass ratio with the IMRPhenomPv2_NRTidal333See Table 1 of Abbott et al.
(2019a) for a description of the waveform model. waveform model (Hannam et
al., 2014; Khan et al., 2016; Dietrich et al., 2019) for GW170817 and
GW190425. Furthermore we use the median chirp mass values of
$\mathcal{M}_{c}=1.186$ M⊙ for
GW170817444https://dcc.ligo.org/LIGO-P1800370/public/ and
$\mathcal{M}_{c}=1.44$ M⊙ for
GW190425555https://dcc.ligo.org/LIGO-P2000223/public/.

Figure 2: In blue we show the bolometric luminosity of GW170817 from the data
compiled in Kasliwal et al. (2017). The red band contains $95\%$ of the light
curves of the posterior distribution when fitted with the model described in
Section 3.2.1.

Figure 3: Upper panels: Constraints on the mass-radius relation of neutron
stars, given the posterior distribution on EOS parameters $\bm{\theta}$ using
the PP model (left) and CS model (right) when analyzing the gravitational wave
events GW170817 (Abbott et al., 2017c) and GW190425 (Abbott et al., 2020a),
both separately and combined. The estimated tidal deformability from GW170817
offers more posterior support for softer EOS, and thus lower radii. For
GW190425 only weak upper limits could be set on the tidal deformability, but
the relatively high estimated mass of the primary object disfavors softer EOS,
as we are not considering any high-mass information from radio pulsars here.
Lower panels: The change in the posterior distribution on the EOS when
including information from the kilonova associated with GW170817, AT2017gfo
(Kasliwal et al., 2017). The estimated mass that was ejected during the merger
favors higher tidal deformabilities, and thus constrains the mass-radius space
at low radii.

The upper panels of Figure 3 show the posterior distributions on the EOS for
both events in the mass-radius space. We note that the constraints on tidal
deformability from GW170817 give more support to softer EOS, although the
$95\%$ credible region spans a relatively large range of radii. GW190425 only
led to weak upper limits on the tidal deformability due to its low SNR and
single-detector detection. The EOS is however constrained as a result of the
high mass of the primary component (with $95\%$ credible range $1.60-1.87$
M⊙), excluding EOS that do not support these masses.

Figure 4: Upper panels: The $68\%$ and $95\%$ credible regions of the EOS
given the mass-radius estimate of PSR J0740$+$6620 by Riley et al. (2021),
using the PP model (left) and CS model (right). The black dashed lines and
orange dashed-dotted lines indicate the $95\%$ credible region of the prior
and the constraints given the radio mass measurement of PSR J0740$+$6620 by
Fonseca et al. (2021), respectively. The red contour shows the posterior
distribution on central energy density and pressure for this source, and in
the inset we plot the KL-divergence as a function of energy density. Lower
panels: Same as upper panels but for the mass-radius space. Also shown in
blue, dotted lines is the $95\%$ credible region of the EOS posterior
distribution, when analyzing the result from Riley et al. (2021) without the
inclusion of the XMM-dataset (so NICER only). In addition, we show the mass-
radius posterior for PSR J0740$+$6620 by Riley et al. (2021) as dark-green
contours (68% and 95%). Note that when considering both NICER and XMM data,
the posterior distribution (green shaded) is very close to the constraints
obtained from the radio mass measurement of PSR J0740$+$6620 (orange), due to
this mass-radius posterior (dark green) showing support over an extended range
of radii.

Figure 5: Upper panels: Constraints on the mass-radius space of neutron stars,
given the posterior distribution of EOS parameters $\bm{\theta}$ using the PP
model (left) and CS model (right). Shown are the $68\%$ and $95\%$ credible
regions when analyzing PSR J0030$+$0451, PSR J0740$+$6620 and the combination
of the two pulsars. Note that the distribution of PSR J0030$+$0451 is
different than in Raaijmakers et al. (2019), because here we have not included
any high-mass pulsar information. Lower panels: Similar to upper panels, but
when analyzing jointly mass-radius estimates from PSR J0740$+$6620 (Riley et
al., 2021), PSR J0030$+$0451 (Riley et al., 2019), mass-tidal deformability
estimates from GW170817 (Abbott et al., 2019a) and GW190425 (Abbott et al.,
2020a) and the kilonova data of Kasliwal et al. (2017) as described in Section
3.2.1. Combined, we find the radius of a $1.4\,$M⊙ neutron star to be
constrained to the 95% credible ranges $12.33^{+0.76}_{-0.81}\,$km (PP model)
and $12.18^{+0.56}_{-0.79}\,$km (CS model). To show the impact of the radius
measurement of PSR J0740$+$6620 we also plot the posterior distribution when
analyzing combined constraints with only the $2.08$M⊙ mass measurement of PSR
J0740$+$6620 (orange dashed-dotted lines).

####  3.2.1 AT2017gfo

Following the detection of GW170817 an EM counterpart was observed across the
frequency spectrum (see, e.g., Abbott et al. (2017a, b) and references
therein; Coulter et al. (2017); Chornock et al. (2017); Drout et al. (2017);
Hallinan et al. (2017); Kasliwal et al. (2017, 2019); Margutti et al. (2017);
Pian et al. (2017); Smartt et al. (2017); Troja et al. (2017)). Of particular
interest here is the thermal infrared-optical-ultraviolet transient powered by
radioactive decay of r-process nucleosynthesis in the neutron-rich material
ejected during merger; the so-called kilonova or macronova (e.g. Li &
Paczyński, 1998; Kulkarni, 2005; Metzger et al., 2010). The kilonova
properties depend on the mass, velocity, and composition of the ejected
material, which in turn depend on the binary progenitor parameters such as the
tidal deformability of the neutron stars. Using this connection it is possible
to constrain the EOS from the kilonova light curve (see, e.g., Coughlin et
al., 2018; Radice & Dai, 2019; Hinderer et al., 2019; Capano et al., 2020;
Dietrich et al., 2020).

Here we analyze the bolometric luminosity of GW170817 (as compiled in Kasliwal
et al., 2017) via the new Bayesian framework outlined in Raaijmakers et al.
(2021b). We consider a two-component kilonova model, where the first
component, the dynamical ejecta, is associated with material ejected through
tidal forces and the shock interface between the two neutron stars (see, e.g.,
Radice et al., 2018, and references therein). The second component is
associated with neutrino-driven winds or material ejected through viscous
forces. We connect the outflow properties of these components to the binary
progenitor properties by using the formulae presented in Krüger & Foucart
(2020) for dynamical ejecta and disk mass, which are fitted to numerical
simulations of compact mergers. The velocity of the dynamical ejecta is
calculated using the formula in Coughlin et al. (2019), while the velocity of
the disk wind ejecta is left as a free parameter. The dynamical ejecta
includes both material ejected through tidal forces and material ejected
through shocks on the contact interface between the stars (see, e.g.,
Sekiguchi et al., 2016; Dietrich & Ujevic, 2017; Tanaka et al., 2020; Nedora
et al., 2021). To distinguish these we consider two different opacities in the
dynamical ejecta, corresponding to tidal tail and shock ejecta , where the
latter is less neutron-rich compared to the tidal tail and thus has a lower
opacity (see Table 1). For simplicity we take a single opacity for the disk
wind ejecta. The outflow properties are then connected to a bolometric
luminosity through the semi-analytic light curve model by Hotokezaka & Nakar
(2020). The priors on all parameters are shown in Table 1.

The fit to the bolometric luminosity of AT2017gfo using the data compiled in
Kasliwal et al. (2017) is shown in Figure 2, showing all datapoints to be
contained within the $95\%$ credible region of the posterior distribution. In
the lower panels of Figure 3 we show the updated prior distribution for the
EOS with GW170817 and with the inclusion of AT2017gfo. The EM data gives more
posterior support to stiffer over softer EOS, due to the estimated ejected
mass requiring a neutron star with larger tidal deformability. The estimated
radius of a $1.4$ M⊙ neutron star for the PP and CS model is
$12.12^{+1.10}_{-1.44}\,$km and $11.53^{+1.16}_{-1.15}\,$km, respectively,
which is broadly consistent with multimessenger constraints obtained by other
works (see, e.g., Coughlin et al., 2019; Dietrich et al., 2020; Capano et al.,
2020; Breschi et al., 2021; Nicholl et al., 2021). Important to note is that
the EM modeling of the kilonova here is simplified and relies on a few
assumptions that are known to affect results, such as spherical ejecta
geometry (see, e.g., Heinzel et al., 2021; Korobkin et al., 2021), fixed
nuclear heating rate (see, e.g., Barnes et al., 2020), and an incomplete
mapping between properties of the binary system and the ejecta outflows. It is
also dependent on the choice of light curve modeling, where the distinction
can be made between semi-analytic modeling (such as in this work and, e.g.,
Breschi et al., 2021; Nicholl et al., 2021) and interpolating between
radiative transfer simulations (e.g., Coughlin et al., 2018; Dietrich et al.,
2020). We use a semi-analytical model from Hotokezaka & Nakar (2020), which
for the current statistical uncertainty in gravitational wave parameter
estimation and uncertainty in light curve observations produces consistent
results to full radiative transport models (see, e.g., Coughlin et al., 2020a,
b), although this will change in the future with improved gravitational wave
detectors and optical telescopes.

Table 1: The parameters used in the model described in Section 3.2.1 and their prior support in the analysis of AT2017gfo. The notation $U(a,b)$ here means uniformly drawn between boundaries $a$ and $b$. Parameters | Prior density and support
---|---
Binary properties
$\mathcal{M}_{c}$ [M⊙] |  $\sim$ U(1.18, 1.2)
q |  $\sim$ U(0.2, 1)
$\Lambda_{1}$ |  $\sim$ U(0, 2500)
$\Lambda_{2}$ |  $\sim$ U(0, 2500)
Ejecta and light curve properties
Mdyn [M⊙] | Eq. (2) Raaijmakers et al. (2021b)
vdyn [c] | Eq. (D5) Coughlin et al. (2019)
$v_{\rm{min,dyn}}$ [c] |  $\sim$ U(0.1, 1.0) vdyn
$v_{\rm{max,dyn}}$ [c] |  $\sim$ U(1.5, 2.5) vdyn
vκ [c] |  $\sim$ U(vmin,dyn, vmax,dyn)
$\kappa_{\rm{low}}$ [cm2 g-1] |  $\sim$ U(0.1, 5)
$\kappa_{\rm{high}}$ [cm2 g-1] |  $\sim$ U(5, 30)
Mwind [M⊙] | Eq. (10) Raaijmakers et al. (2021b)
vwind [c] |  $\sim$ U(0.03, 0.15)
$v_{\rm{min,wind}}$ [c] |  $\sim$ U(0.1, 1.0) vwind
$v_{\rm{max,wind}}$ [c] |  $\sim$ U(1.5, 2.0) vwind
$\kappa_{\rm{wind}}$ [cm2 g-1] |  $\sim$ U(0.1, 5)

###  3.3 NICER mass-radius and multimessenger constraints

Next we study the constraints on the EOS from the new mass-radius estimate of
PSR J0740$+$6620 using data from NICER and XMM, presented in Riley et al.
(2021). They find a radius of $12.39_{-0.98}^{+1.30}$ km, and a mass of
$2.072_{-0.066}^{+0.067}$ M⊙, where the upper and lower limit bound the $68\%$
credible regions. The EOS results are shown in Figure 4, both in energy
density-pressure and mass-radius space. From the Kullback-Leibler
(KL)-divergence (Kullback & Leibler, 1951) plotted as a function of energy
density in the upper insets, we find that especially at higher energy
densities there is a significant information gain from prior-to-posterior.
Note that similar but, especially for the CS model, broader constraints are
found for the posterior distribution when only using the radio mass
measurement of PSR J0740$+$6620, as indicated by the orange dashed-dotted
lines. This is a result of the mass-radius estimate of PSR J0740$+$6620 being
very consistent with our prior ranges informed by low-density chiral EFT
calculations. The chiral EFT calculations do exclude however stiffer EOS with
radii $>14$ km, where the mass-radius posterior of PSR J0740$+$6620 has non-
negligible posterior support. For the CS model this effect is stronger as
additional constraints on the speed of sound at $1.5n_{0}$ in the CS model
lead to overall smaller radii than in the PP model (see Section 2.3 of
Raaijmakers et al. 2020).

In the mass-radius space we also plot the EOS constraints given the joint
NICER mass-radius estimate excluding the XMM data. For this analysis Riley et
al. (2021) report a value of $11.29_{-0.81}^{+1.20}$ km for the radius, and
$2.078_{-0.063}^{+0.066}$ M⊙ for the mass. As this joint mass-radius estimate
has slightly more posterior support for lower radii, the corresponding EOS
constraints suggest a softening of the EOS at high densities. These results
should be interpreted with caution however, because the NICER-only analysis
leads to an under-prediction of the background (the contribution from
instrumental or astrophysical background to the unpulsed component of the
pulse profile). This results in more of the unpulsed component being
attributed to the hot regions via high compactness solutions. The XMM data
show that a larger component of the unpulsed emission must come from true
background, eliminating these high compactness solutions and increasing the
inferred radius in the joint NICER-XMM analysis (see also Section 4.2 in Riley
et al. (2021)).

Finally, in Figure 5 we show the constraints on the EOS from PSR J0740$+$6620,
PSR J0030$+$0451 (first derived in Raaijmakers et al. 2019, but here no
information on high-mass pulsars is included) and the combination of the two
pulsars. Note that for the combined constraints, most of the information comes
from PSR J0740$+$6620, since the $68\%$ credible region of the mass-radius
posterior of PSR J0030$+$0451 covers a broad range in radii that are
consistent with the EOS constraints from PSR J0740$+$6620.

In the lower panels of Figure 5 we show the combined constraints on the EOS
including mass-radius estimates from PSR J0740$+$6620 (Riley et al., 2021),
PSR J0030$+$0451 (Riley et al., 2019) and mass-tidal deformability estimates
from GW170817 (Abbott et al., 2019a) and GW190425 (Abbott et al., 2020a) and
the kilonova AT2017gfo (Kasliwal et al., 2017). We find that especially the
pulsar mass-radius estimates by NICER favor stiffer EOS, as well as GW170817
when the associated kilonova AT2017gfo (Kasliwal et al., 2017) is included.
The weak constraints from GW190425 on the tidal deformability are also broadly
consistent with the constraints coming from the other sources. As a comparison
we show the posterior distribution when combining all analyses excluding the
mass-radius estimate of PSR J0740$+$6620, but with the radio mass measurement
of Fonseca et al. (2021). We note that the additional radius information on
PSR J0740$+$6620 constrains the softer EOS, especially for the CS model.

Figure 6: Different chiral EFT bands for the pressure of neutron star matter
at nuclear densities, $n/n_{0}$ in units of saturation density $n_{0}=0.16$
fm-3, and their matching to the BPS crust EOS at $0.5n_{0}$. The different
bands are based on microscopic calculations of neutron matter from Hebeler et
al. (2013), Tews et al. (2013), Lynn et al. (2016) and Drischler et al. (2019)
and include beta equilibrium (with protons and electrons) following the
construction in Hebeler et al. (2013). The four chiral EFT calculations are
considered between $0.5n_{0}$ and $1.1n_{0}$ in the analyses presented in
Section 4. Also shown are examples of the fit we use to approximate the EOS
within these uncertainty bands, see Eq. (5), and connect to the BPS crust EOS.
For a comparison of the chiral EFT bands in pure neutron matter, see Figure 1
in Huth et al. (2021).

Figure 7: Upper panels: $95\%$ credible region for the mass-radius space given
the mass-radius estimate of PSR J0740$+$6620 by Riley et al. (2021), using the
PP model (left) and CS model (right). The different results correspond to
using the four different chiral EFT calculations between $0.5$ and $1.1n_{0}$
as shown in Fig. 6. Moroever, the red, dashed lines correspond to the $95\%$
credible region, if the PP or CS parameterization is used down to $0.5n_{0}$,
i.e., immediately following the BPS crust, so that no information from chiral
EFT is used. Lower panels: Marginalized posterior distributions for the
pressure $P$ above saturation density, at density $n=1.5n_{0}$ (left) and
$n=2n_{0}$ (right) above the chiral EFT bands.

##  4 Sensitivity of posteriors to nuclear constraints at low densities

To investigate the impact of the EOS constraints from nuclear physics we
compare our analysis of PSR J0740$+$6620 using four different chiral EFT
uncertainty bands. All bands are based on microscopic calculations for pure
neutron matter, which are then extended to neutron star matter in beta-
equilibrium using the formalism discussed in Hebeler et al. (2013). In order
to improve the description of all employed EOSs, we generalized the density
dependence of the energy-density functional [see Eq. (2) in Hebeler et al.
(2013)] by enlarging the range of the exponent $\gamma$ to
$\gamma\in\left[1.2,2.5\right]$.

The results from Hebeler et al. (2013) formed the basis of our previous
studies (Raaijmakers et al., 2019, 2020). The calculations for pure neutron
matter were initially performed in Hebeler & Schwenk (2010) using many-body
perturbation theory, while the uncertainty band results mainly from variations
of the couplings involved in three-nucleon interactions. Second, in Tews et
al. (2013) the calculations for neutron matter were improved by including for
the first time all two-, three-, and four-neutron interactions to next-to-
next-to-next-to-leading order (N3LO), which are predicted in a parameter-free
way for neutron matter (see, e.g., Hebeler et al. 2015; Hebeler 2021 for
reviews). Third, in Drischler et al. (2019) the calculations were further
optimized by improving the treatment of three-nucleon interactions and
extending the many-body expansion to higher orders. In addition, the EOS
uncertainty bands also include effects from variations of regulator scales in
state-of-the-art nucleon-nucleon and three-nucleon interactions. In this work,
we use the combined 450 MeV and 500 MeV N3LO uncertainty bands from Drischler
et al. (2019). Finally, we include results of Lynn et al. (2016). These were
obtained by nonperturbative quantum Monte-Carlo simulations of neutron matter
at next-to-next-to-leading order (N2LO). This represents a completely
different many-body method than those used for the other three bands, and the
results of Lynn et al. (2016) are also based on a different set of local two-
and three-nucleon interactions derived from chiral EFT.

Table 2: Key quantities from the posterior distributions obtained in Sections 3 and 4: The radius of a $1.4$ M⊙, $1.6$ M⊙, and $1.8$ M⊙ neutron star, as well as $\Delta R=R_{2}-R_{1.4}$, and the maximum mass of a non-rotating neutron star MTOV. For the analyses of Section 4, we also show the inferred central energy density $\varepsilon_{c}$, the corresponding central pressure Pc, and the Bayes’ factor $K$ comparing with the model using the chiral EFT band from Hebeler et al. (2013). The first four column results are for the different chiral EFT bands from Hebeler et al. (2013) (Heb 13), Tews et al. (2013) (Tews 13), Lynn et al. (2016) (Lynn 16), and Drischler et al. (2019) (Dri 19), while all other results are for the baseline inference using Heb 13. The column “Combined with” refers to the NICER $\times$ XMM analysis of PSR J0740$+$6620, the NICER analysis of PSR J0030$+$0451 and multimessenger constraints combined, while in the column “Combined without” the NICER $\times$ XMM analysis of PSR J0740$+$6620 is replaced with just the radio mass measurement by Fonseca et al. (2021). The radii are given in km, MTOV in M⊙, and $\varepsilon_{c}$ and Pc in g/cm3 and dyn/cm2, respectively. The upper and lower values correspond to the $95\%$ credible interval. | PSR J0740$+$6620, NICER x XMM | PSR J0030 | GW170817 | GW170817 | Combined | Combined
---|---|---|---|---|---|---
| Heb 13 | Tews 13 | Lynn 16 | Dri 19 | +0451 | \+ GW190425 | \+ AT2017gfo | without | with
PP model |
$R_{1.4}$ | $12.56^{+0.80}_{-0.91}$ | $12.85^{+0.77}_{-0.95}$ | $12.35^{+0.96}_{-0.98}$ | $12.87^{+0.85}_{-0.98}$ | $12.35^{+0.99}_{-1.99}$ | $11.51^{+1.51}_{-1.47}$ | $12.12^{+1.10}_{-1.44}$ | $12.30^{+0.72}_{-0.76}$ | $12.33^{+0.76}_{-0.81}$
R1.6 | $12.60^{+0.87}_{-1.00}$ | $12.87^{+0.87}_{-1.05}$ | $12.40^{+1.03}_{-1.04}$ | $12.90^{+0.94}_{-1.08}$ | $12.50^{+0.96}_{-2.08}$ | $11.43^{+1.68}_{-1.53}$ | $12.10^{+1.23}_{-1.69}$ | $12.32^{+0.92}_{-0.99}$ | $12.35^{+0.83}_{-0.90}$
R1.8 | $12.62^{+0.98}_{-1.19}$ | $12.86^{+1.00}_{-1.26}$ | $12.42^{+1.13}_{-1.19}$ | $12.89^{+1.08}_{-1.29}$ | $12.68^{+0.94}_{-1.99}$ | $11.65^{+1.64}_{-1.80}$ | $12.22^{+1.26}_{-1.91}$ | $12.29^{+1.06}_{-1.19}$ | $12.33^{+0.97}_{-1.06}$
$\Delta R$ | $-0.24^{+0.65}_{-1.04}$ | $-0.17^{+1.26}_{-1.21}$ | $-0.22^{+0.60}_{-1.05}$ | $-0.45^{+1.14}_{-1.07}$ | $-0.13^{+0.76}_{-1.02}$ | $-0.35^{+0.80}_{-1.09}$ | $-0.26^{+0.77}_{-1.14}$ | $-0.30^{+0.64}_{-1.06}$ | $-0.29^{+0.61}_{-0.98}$
MTOV | $2.26^{+0.15}_{-0.23}$ | $2.33^{+0.14}_{-0.30}$ | $2.22^{+0.19}_{-0.21}$ | $2.33^{+0.18}_{-0.31}$ | $1.74^{+0.66}_{-0.57}$ | $1.84^{+0.51}_{-0.17}$ | $1.96^{+0.42}_{-0.44}$ | $2.23^{+0.15}_{-0.23}$ | $2.23^{+0.14}_{-0.23}$
$\log_{10}(\varepsilon_{c})$ | $14.99^{+0.27}_{-0.16}$ | $14.99^{+0.28}_{-0.18}$ | $15.00^{+0.26}_{-0.16}$ | $14.99^{+0.28}_{-0.19}$ | $14.86^{+0.28}_{-0.13}$ | - | - | - |
$\log_{10}$(P${}_{c})$ | $35.39^{+0.39}_{-0.24}$ | $35.37^{+0.41}_{-0.26}$ | $35.41^{+0.37}_{-0.25}$ | $35.37^{+0.43}_{-0.28}$ | $34.92^{+0.30}_{-0.21}$ | - | - | - |
K | 1.00 | 0.89 | 1.00 | 0.85 | - | - | - | - |
CS model |
R1.4 | $12.27^{+0.54}_{-0.90}$ | $12.49^{+0.49}_{-0.87}$ | $12.16^{+0.63}_{-0.97}$ | $12.56^{+0.51}_{-0.92}$ | $11.51^{+1.12}_{-1.90}$ | $11.18^{+1.33}_{-1.51}$ | $11.53^{+1.16}_{-1.15}$ | $11.98^{+0.63}_{-0.71}$ | $12.18^{+0.56}_{-0.79}$
R1.6 | $12.25^{+0.59}_{-0.94}$ | $12.43^{+0.55}_{-0.92}$ | $12.16^{+0.66}_{-0.99}$ | $12.50^{+0.55}_{-0.96}$ | $11.48^{+1.20}_{-1.92}$ | $10.92^{+1.58}_{-1.56}$ | $11.33^{+1.38}_{-1.41}$ | $11.91^{+0.78}_{-0.94}$ | $12.14^{+0.61}_{-0.84}$
R1.8 | $12.14^{+0.69}_{-1.05}$ | $12.27^{+0.66}_{-1.04}$ | $12.08^{+0.74}_{-1.07}$ | $12.33^{+0.65}_{-1.06}$ | $11.52^{+1.23}_{-1.77}$ | $10.85^{+1.73}_{-1.47}$ | $11.34^{+1.44}_{-1.65}$ | $11.72^{+0.96}_{-1.06}$ | $12.00^{+0.74}_{-0.96}$
$\Delta R$ | $-0.69^{+1.10}_{-1.02}$ | $-0.72^{+1.12}_{-1.08}$ | $-0.58^{+1.03}_{-1.08}$ | $-1.06^{+1.46}_{-0.83}$ | $-0.93^{+1.31}_{-0.86}$ | $-0.93^{+1.36}_{-0.83}$ | $-0.81^{+1.22}_{-0.92}$ | $-0.91^{+1.15}_{-0.85}$ | $-0.74^{+1.09}_{-0.95}$
MTOV | $2.13^{+0.33}_{-0.16}$ | $2.13^{+0.29}_{-0.18}$ | $2.14^{+0.34}_{-0.17}$ | $2.12^{+0.31}_{-0.16}$ | $1.46^{+0.82}_{-0.42}$ | $1.81^{+0.45}_{-0.15}$ | $1.85^{+0.56}_{-0.30}$ | $2.09^{+0.26}_{-0.15}$ | $2.11^{+0.29}_{-0.16}$
$\log_{10}(\varepsilon_{c})$ | $15.19^{+0.21}_{-0.20}$ | $15.19^{+0.20}_{-0.20}$ | $15.18^{+0.21}_{-0.21}$ | $15.20^{+0.19}_{-0.20}$ | $15.03^{+0.33}_{-0.21}$ | - | - | - |
$\log_{10}$(P${}_{c})$ | $35.61^{+0.30}_{-0.27}$ | $35.62^{+0.30}_{-0.28}$ | $35.60^{+0.30}_{-0.28}$ | $35.63^{+0.30}_{-0.29}$ | $35.05^{+0.30}_{-0.24}$ | - | - | - |
K | 1.00 | 1.04 | 0.92 | 1.05 | - | - | - | - |

Similar to Raaijmakers et al. (2020) we approximate the EOS within these bands
with a single polytrope $P=Nn^{\Gamma}$. However, to obtain a better fit to
the additional bands considered here, we vary the polytropic index $\Gamma$ as
a function of the normalization $N$,

$\Gamma(N)=\frac{(N-N_{\rm min})}{(N_{\rm max}-N_{\rm min})}(\Gamma_{\rm
max}-\Gamma_{\rm min})+\Gamma_{\rm min}\,,$ (5)

where $N_{\rm min/max}$ and $\Gamma_{\rm min/max}$ are determined by fitting a
polytrope to the lower and upper bound of the band. In Figure 6 we show the
four different bands for the pressure of neutron star matter with an example
of the fit through each band. This shows the consistency of these different
chiral EFT calculations, with different methods, interactions, and
approximations. The first point of the band where $n/n_{0}>0.5$ is matched to
the BPS crust EOS at $0.5n_{0}$ via a linear interpolation.

We study the dependence of the EOS constraints on the different chiral EFT
bands by inferring the EOS from the mass-radius estimate of PSR J0740$+$6620
using each band and both high-density parameterizations. The results are shown
in Figure 7. We also show the $95\%$ credible region of the updated prior
distribution when directly joining the PP or CS high-density parameterization
to the crust EOS at $0.5\rho_{ns}$. As expected the chiral EFT calculations
mostly exclude stiffer EOS. While the different chiral EFT bands yield very
good agreement on the upper bound of the radius estimates, the lower bound on
the radius does slightly depend on the chiral EFT band used, especially at
lower neutron star masses, depending on how soft the chiral EFT band is (see
Figure 6).

In the lower panels of Figure 7 we also show the posterior distributions on
the pressure at densities $n=1.5n_{0}$ and $n=2n_{0}$ above the chiral EFT
bands. These results demonstrate that the PSR J0740$+$6620 mass-radius
measurement systematically prefers higher pressures at these densities
compared to the corresponding prior distributions of each chiral EFT band.
Furthermore, the posteriors at $n=2n_{0}$ agree very well for all chiral EFT
bands and are peaked around $P\sim 10^{34.5}\text{dyn}/\text{cm}^{2}\sim
20\,\text{MeV}/\text{fm}^{3}$.

##  5 Discussion

In this Letter, we have investigated the constraints on the EOS posed by the
new joint mass-radius estimate from NICER $\times$ XMM data (Riley et al.,
2021), and compared and combined with multimessenger EOS constraints from
radio timing, gravitational wave mergers and their counterparts, and the
previous PSR J0030$+$0451 mass-radius estimate by NICER. In Table 2 we
summarize the results obtained in Sections 3 and 4 for the constraints on the
radius of a $1.4,1.6$ and $1.8$ M⊙ neutron star, as well as $\Delta
R=R_{2}-R_{1.4}$, and the maximum mass of a non-rotating neutron star MTOV, as
well as the constraints on the central energy density and pressure for PSR
J0740$+$6620.

###  5.1 Implications for nuclear physics

We have studied the sensitivity of the EOS constraints from PSR J0740$+$6620
using four different low-density EOS calculations from chiral EFT (see Section
4). From the results presented in Figure 7 and Table 2 we conclude that the
constraints on the EOS are only weakly dependent on the choice of low-density
calculations, although small differences exist at lower radii. Assuming all
four low-density calculations to be equally probable, we can compute the
Bayes’ factor $K$ by taking the ratio of the evidence of each MultiNest run,
and assess whether one model is preferred over another by the data of PSR
J0740$+$6620\. We list the Bayes’ factors in Table 2, where each model is
compared to using the chiral EFT band from Hebeler et al. (2013). All values
are close to one, indicating that there is no substantial support for one
model over the other, based on the mass-radius estimate of PSR J0740$+$6620\.
These results are consistent with the observation that predictions for pure
neutron matter are well constrained by modern nuclear forces derived within
chiral EFT (Huth et al., 2021; Hebeler, 2021).

Also shown in Table 2 are the values of $\Delta R=R_{2}-R_{1.4}$, the
difference in radius of a $2$ M⊙ and $1.4$ M⊙ neutron star. As pointed out by
Drischler et al. (2021), the value of $\Delta R$, if positive, can give an
indication that possibly unusual stiffening happens at high densities. We find
however all values to be consistent with the mean $\Delta R$ being negative,
but due to the broad uncertainty no conclusive statements can be made.

###  5.2 Implications for maximum mass

An important quantity relating to the EOS is the maximum stable mass of a non-
rotating neutron star, MTOV. Accurate knowledge of MTOV can aid in classifying
compact mergers and merger remnants. In Figure 8 we show posterior
distributions on MTOV when analyzing the updated radio mass measurement of PSR
J0740$+$6620, the joint mass-radius estimate of PSR J0740$+$6620 and combining
GW170817, GW190425, AT2017gfo, PSR J0740$+$6620 and PSR J0030$+$0451\. The
latter results in 95% credible ranges for M${}_{\rm TOV}=2.23^{+0.15}_{-0.25}$
M⊙ and M${}_{\rm TOV}=2.11^{+0.28}_{-0.16}$ M⊙ for the PP and CS model,
respectively. This is in agreement with values previous found (see, e.g.,
Nathanail et al., 2021, and references therein) when assuming the secondary
component in GW190814 was a black hole (Abbott et al., 2020d). Note that the
higher end of the distribution in Figure 8 is very dependent on our choice of
parameterization, as no information is included from sources with masses above
$2.08$ M⊙. One could use information on the merger remnant of GW170817 to put
an upper bound on MTOV (see, e.g., Margalit & Metzger, 2017; Shibata et al.,
2017; Ruiz et al., 2018), but that is beyond the scope of this Letter. The
lower end of the distribution on the other hand is strongly correlated with
the radio mass measurement of PSR J0740$+$6620\. The recently lowered mass
distribution presented in Fonseca et al. (2021) results in slightly lower
values for MTOV compared to the distributions found in Raaijmakers et al.
(2020).

Figure 8: Posterior distribution of the maximum mass of a non-rotating neutron
star MTOV for the PP model (left) and CS model (right) when considering only
the radio mass measurement of PSR J0740$+$6620, the joint mass-radius estimate
of PSR J0740$+$6620 (NICER $\times$ XMM), and when combining NICER’s results
on PSR J0740$+$6620 and PSR J0030$+$0451 with GW170817 and GW190425, and
AT2017gfo. For the latter (“Combined”) we find a 95% credible range for
M${}_{\rm TOV}=2.23^{+0.14}_{-0.23}\,M_{\odot}$ and M${}_{\rm
TOV}=2.11^{+0.29}_{-0.16}\,M_{\odot}$ for the PP and CS model, respectively.
Also shown in pink is the radio mass measurement of PSR J0740$+$6620 from
Fonseca et al. (2021), as the heaviest pulsar measured to date.

###  5.3 Systematic uncertainties and framework comparisons

The analysis presented in this Letter is conditional on both the modeling
choices of the dense matter EOS and on modeling choices within each analysis
of the multimessenger sources considered here. The sensitivity to the EOS
modeling is explored here by employing two different high-density
parameterizations and four different low-density chiral EFT calculations (see
Section 4). From Table 2 we conclude that the CS model systematically predicts
lower radii, as a result of the additional constraints on the speed of sound
that are not considered in the PP model. The discrepancy between the two
models increases with increasing neutron star mass, as high-mass stars depend
more sensitively on the choice of high-density parameterization. The two
models considered here are however not exhaustive as many more high-density
parameterizations exist (see, e.g., Lindblom, 2018; Capano et al., 2020;
O’Boyle et al., 2020).

Furthermore, we do not consider the impact of any systematic effects present
in estimating the posterior distributions on $M$, $R$ and $\Lambda$. For
example, the uncertainty in modeling the hot regions in pulse profile modeling
and the effect on the EOS has been studied in Raaijmakers et al. (2020) by
using two different models to fit PSR J0030$+$0451, which led to slightly
different constraints. For PSR J0740$+$6620, different assumptions and priors
lead to a higher estimated radius in the independent analysis of Miller et al.
(2021) (see the extensive discussion of this issue in Section 4.4 of Riley et
al. 2021), and we refer the reader to that paper for an EOS analysis using
those results. 666Note however that one of the main reasons for the higher
inferred radius reported by Miller et al. (2021) is that they do not truncate
the prior on radius during the pulse profile modelling step, which Riley et
al. (2021) do (truncating above 16 km, reflecting the lack of EOS models
predicting higher radii, and thereby lowering the computational cost by
reducing the parameter space). In the analysis by Miller et al. (2021) the
lack of prior support for high radii is effectively incorporated at a later
stage, in the EOS analysis..

Measurements in $\Lambda$ from gravitational wave data are also sensitive to
choice of priors and gravitational waveform models (see, e.g., Kastaun & Ohme,
2019; Gamba et al., 2020). Lastly, many different kilonova models exist (see,
e.g., Dietrich et al., 2020; Nicholl et al., 2021; Breschi et al., 2021, for
recent analyses) that derive slightly different constraints on the EOS due to
differences in modeling assumptions on, e.g., geometry, composition and the
connection between binary properties and outflow properties.

The inference framework employed in this Letter was first discussed in Riley
et al. (2018) and subsequently developed in Greif et al. (2019); Raaijmakers
et al. (2019, 2020), which also introduced the chiral EFT constraints.
Although an exhaustive comparison with other frameworks is out of the scope of
this work, we will briefly mention similarities and differences with some
commonly used frameworks in the field. Firstly, we make use of two particular
high-density EOS parameterizations. Besides many different existing choices in
these parameterizations, a completely different approach is to use non-
parametric inference involving Gaussian Processes (see, e.g., Landry & Essick,
2019; Essick et al., 2020a; Han et al., 2020), or discretely sampling a set of
pre-computed EOS (see, e.g., Capano et al., 2020; Dietrich et al., 2020).
Secondly, we compute likelihoods by performing kernel density estimation on
posterior samples of neutron star properties such as mass, radius and tidal
deformability (see also, e.g, Miller et al., 2019; Al-Mamun et al., 2021). It
is also possible to directly infer EOS properties from the observational data,
for example X-ray or gravitational wave data. For the first, Riley et al.
(2018) argue that this approach would be computationally too expensive, while
for the latter this has been done by, e.g., Capano et al. (2020); Dietrich et
al. (2020). A slightly different approach is used by Hernandez Vivanco et al.
(2020), where the likelihood is computed by interpolating marginalized
likelihoods using machine learning.

###  5.4 Summary and future prospects

In summary, the new joint mass-radius estimate of PSR J0740$+$6620
significantly constrains the EOS. For the PP model the information gain is
mostly a result of the high mass of the pulsar, as the $68\%$ credible range
of the radius estimate exactly encompasses our prior distribution, informed by
chiral EFT calculations, in that mass range. For the CS model the relatively
high radius estimate does constrain the model at lower radii on top of
constraints coming from the mass estimate. Combined with other current
observational data from gravitational waves and kilonova light curves, as well
as the NICER mass-radius estimate of PSR J0030$+$0451, we find the 95%
credible ranges $12.38^{+0.70}_{-0.97}\,$km (PP model) and
$12.23^{+0.48}_{-0.97}\,$km (CS model) for the radius of a $1.4$ M⊙ neutron
star.

In the near future, the detailed analysis of gravitational wave events
observed during the second part of the third observing run of LIGO/Virgo are
expected to be published, among them a few candidate events which, in an
initial rapid classification, were identified as containing at least one
neutron star. Any measured tidal deformability from these gravitational waves
events will help constrain the EOS further. There were unfortunately no EM
counterparts for the potential binary neutron star or black hole-neutron star
events during this observing run. The fourth observing run is planned to start
next year, with the LIGO and Virgo detectors close to their design sensitivity
and KAGRA fully joining the network (Abbott et al., 2020b). At design
sensitivity, GW170817-like signals will have signal-to-noise ratios of 100 and
enable measurements of tidal deformability with more than three times better
accuracy (Capano et al., 2020). Subsequent further detector improvements are
already planned for the mid to late 2020s (Abbott et al., 2020b), and an
ongoing worldwide effort is paving the way for next decade’s third generation
detectors. These will improve current measurements of tidal deformability by a
factor of $\sim 10$ and observe the population of tens to hundreds of
thousands of neutron star binaries, with EM counterparts detectable for a
fraction of them (Maggiore et al., 2020; Sathyaprakash et al., 2019a, b).

In the coming months, NICER is expected to deliver mass-radius measurements
for three additional pulsars: two for which independent mass constraints exist
(the $\sim 1.4$ M⊙ pulsar PSR J0437-4715 and the $\sim 1.9$ M⊙ pulsar PSR
J1614-2230); and the pulsar PSR J1231-1411, which has no independently known
constraint on the mass. There will be an update to the inferred mass and
radius of PSR J0030+0451, using a larger data set, taking into account
improvements to our understanding of the NICER instrument response, and
including XMM data in a joint analysis (as done for PSR J0740+6620). There are
also good prospects for narrowing the mass-radius measurements for PSR
J0740+6620, using models of the NICER background. All of these promise further
improvements to our understanding of the dense matter EOS.

## Acknowledgments

This work was supported in part by NASA through the NICER mission and the
Astrophysics Explorers Program. G.R., T.H., S.N. and T.E.R. are grateful for
support from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)
through the VIDI and Projectruimte grants (PI Nissanke). T.H. also
acknowledges support from the NWO sectorplan. T.E.R. and A.L.W. acknowledge
support from ERC Consolidator Grant No. 865768 AEONS (PI: Watts). S.K.G.,
K.H., and A.S. acknowledge support from the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) – Project-ID 279384907 – SFB 1245. W.C.G.H.
acknowledges support through grant 80NSSC20K0278 from NASA. This work was
sponsored by NWO Exact and Natural Sciences for the use of supercomputer
facilities, and was carried out on the Dutch national e-infrastructure with
the support of SURF Cooperative. We thank Sharon Morsink and Renee Ludlam for
comments on a draft manuscript. We thank Mansi Kasliwal for providing us with
the data for AT2017gfo. We thank the anonynomous referee for providing useful
feedback on the manuscript.

## References

  * Abbott et al. (2017a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, ApJ, 848, L12
  * Abbott et al. (2017b) —. 2017b, ApJ, 848, L13
  * Abbott et al. (2017c) —. 2017c, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2018) —. 2018, Phys. Rev. Lett., 121, 161101
  * Abbott et al. (2019a) —. 2019a, PhRvX, 9, 011001
  * Abbott et al. (2019b) —. 2019b, PhRvX, 9, 031040
  * Abbott et al. (2020a) —. 2020a, ApJ, 892, L3
  * Abbott et al. (2020b) —. 2020b, Living Reviews in Relativity, 23, 3
  * Abbott et al. (2020c) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020c, arXiv e-prints, arXiv:2010.14527
  * Abbott et al. (2020d) —. 2020d, ApJ, 896, L44
  * Adhikari et al. (2021) Adhikari, D., Albataineh, H., Androic, D., et al. 2021, Phys. Rev. Lett., 126, 172502
  * Al-Mamun et al. (2021) Al-Mamun, M., Steiner, A. W., Nättilä, J., et al. 2021, Phys. Rev. Lett., 126, 061101
  * Alvarez-Castillo et al. (2020) Alvarez-Castillo, D., Ayriyan, A., Barnaföldi, G. G., & Pósfay, P. 2020, Physics of Particles and Nuclei, 51, 725
  * Annala et al. (2020) Annala, E., Gorda, T., Kurkela, A., Nättilä, J., & Vuorinen, A. 2020, Nature Phys., 16, 907
  * Arcavi et al. (2017) Arcavi, I., Hosseinzadeh, G., Howell, D. A., et al. 2017, Nature, 551, 64
  * Ascenzi et al. (2019) Ascenzi, S., De Lillo, N., Haster, C.-J., Ohme, F., & Pannarale, F. 2019, ApJ, 877, 94
  * Barnes et al. (2020) Barnes, J., Zhu, Y. L., Lund, K. A., et al. 2020, arXiv e-prints, arXiv:2010.11182
  * Baym et al. (1971) Baym, G., Pethick, C., & Sutherland, P. 1971, ApJ, 170, 299
  * Behnel et al. (2011) Behnel, S., Bradshaw, R., Citro, C., et al. 2011, Computing in Science Engineering, 13, 31
  * Birkhan et al. (2017) Birkhan, J., Miorelli, M., Bacca, S., et al. 2017, Phys. Rev. Lett., 118, 252501
  * Biswas et al. (2021) Biswas, B., Char, P., Nandi, R., & Bose, S. 2021, Phys. Rev. D, 103, 103015
  * Blaschke et al. (2020) Blaschke, D., Ayriyan, A., Alvarez-Castillo, D. E., & Grigorian, H. 2020, Universe, 6, 81
  * Breschi et al. (2021) Breschi, M., Perego, A., Bernuzzi, S., et al. 2021, MNRAS, arXiv:2101.01201
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125
  * Capano et al. (2020) Capano, C. D., Tews, I., Brown, S. M., et al. 2020, Nature Astronomy, 4, 625
  * Chornock et al. (2017) Chornock, R., Berger, E., Kasen, D., et al. 2017, ApJ, 848, L19
  * Christian & Schaffner-Bielich (2020) Christian, J.-E., & Schaffner-Bielich, J. 2020, ApJ, 894, L8
  * Coughlin & Dietrich (2019) Coughlin, M. W., & Dietrich, T. 2019, Phys. Rev. D, 100, 043011
  * Coughlin et al. (2019) Coughlin, M. W., Dietrich, T., Margalit, B., & Metzger, B. D. 2019, MNRAS, 489, L91
  * Coughlin et al. (2018) Coughlin, M. W., Dietrich, T., Doctor, Z., et al. 2018, MNRAS, 480, 3871
  * Coughlin et al. (2020a) Coughlin, M. W., Dietrich, T., Antier, S., et al. 2020a, MNRAS, 492, 863
  * Coughlin et al. (2020b) —. 2020b, MNRAS, 497, 1181
  * Coulter et al. (2017) Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017, Science, 358, 1556
  * Cowperthwaite et al. (2017) Cowperthwaite, P. S., Berger, E., Villar, V. A., et al. 2017, ApJ, 848, L17
  * Cromartie et al. (2020) Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2020, Nature Astronomy, 4, 72
  * Dalcín et al. (2008) Dalcín, L., Paz, R., Storti, M., & D’Elía, J. 2008, Journal of Parallel and Distributed Computing, 68, 655
  * Dietrich et al. (2020) Dietrich, T., Coughlin, M. W., Pang, P. T. H., et al. 2020, Science, 370, 1450
  * Dietrich & Ujevic (2017) Dietrich, T., & Ujevic, M. 2017, Classical and Quantum Gravity, 34, 105014
  * Dietrich et al. (2019) Dietrich, T., Khan, S., Dudi, R., et al. 2019, Phys. Rev. D, 99, 024029
  * Drischler et al. (2020) Drischler, C., Furnstahl, R. J., Melendez, J. A., & Phillips, D. R. 2020, Phys. Rev. Lett., 125, 202702
  * Drischler et al. (2021) Drischler, C., Han, S., Lattimer, J. M., et al. 2021, Phys. Rev. C, 103, 045808
  * Drischler et al. (2019) Drischler, C., Hebeler, K., & Schwenk, A. 2019, Phys. Rev. Lett., 122, 042501
  * Droettboom et al. (2018) Droettboom, M., Caswell, T. A., Hunter, J., et al. 2018, matplotlib/matplotlib v2.2.2, , , doi:10.5281/zenodo.1202077
  * Drout et al. (2017) Drout, M. R., Piro, A. L., Shappee, B. J., et al. 2017, Science, 358, 1570
  * Essick et al. (2020a) Essick, R., Landry, P., & Holz, D. E. 2020a, Phys. Rev. D, 101, 063007
  * Essick et al. (2020b) Essick, R., Tews, I., Landry, P., Reddy, S., & Holz, D. E. 2020b, Phys. Rev. C, 102, 055803
  * Essick et al. (2021) Essick, R., Tews, I., Landry, P., & Schwenk, A. 2021, arXiv e-prints, arXiv:2102.10074
  * Feroz et al. (2009) Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601
  * Flanagan & Hinderer (2008) Flanagan, E. E., & Hinderer, T. 2008, Phys. Rev. D, 77, 021502
  * Fonseca et al. (2021) Fonseca, E., Cromartie, H. T., Pennucci, T. T., et al. 2021, arXiv e-prints, arXiv:2104.00880
  * Forum (1994) Forum, M. P. 1994, MPI: A Message-Passing Interface Standard, Tech. rep., Knoxville, TN, USA
  * Gamba et al. (2020) Gamba, R., Breschi, M., Bernuzzi, S., Agathos, M., & Nagar, A. 2020, arXiv e-prints, arXiv:2009.08467
  * Gough (2009) Gough, B. 2009, GNU Scientific Library Reference Manual - Third Edition, 3rd edn. (Network Theory Ltd.)
  * Greif et al. (2019) Greif, S. K., Raaijmakers, G., Hebeler, K., Schwenk, A., & Watts, A. L. 2019, MNRAS, 485, 5363
  * Hallinan et al. (2017) Hallinan, G., Corsi, A., Mooley, K. P., et al. 2017, Science, 358, 1579
  * Han et al. (2020) Han, M.-Z., Tang, S.-P., Hu, Y.-M., et al. 2020, ApJ, 891, L5
  * Hannam et al. (2014) Hannam, M., Schmidt, P., Bohé, A., et al. 2014, Phys. Rev. Lett., 113, 151101
  * Hebeler (2021) Hebeler, K. 2021, Phys. Rep., 890, 1
  * Hebeler et al. (2015) Hebeler, K., Holt, J. D., Menéndez, J., & Schwenk, A. 2015, Annu. Rev. Nucl. Part. Sci., 65, 457
  * Hebeler et al. (2013) Hebeler, K., Lattimer, J. M., Pethick, C. J., & Schwenk, A. 2013, ApJ, 773, 11
  * Hebeler & Schwenk (2010) Hebeler, K., & Schwenk, A. 2010, Phys. Rev. C, 82, 014314
  * Heinzel et al. (2021) Heinzel, J., Coughlin, M. W., Dietrich, T., et al. 2021, MNRAS, 502, 3057
  * Hernandez Vivanco et al. (2020) Hernandez Vivanco, F., Smith, R., Thrane, E., & Lasky, P. D. 2020, MNRAS, 499, 5972
  * Hinderer et al. (2010) Hinderer, T., Lackey, B. D., Lang, R. N., & Read, J. S. 2010, Phys. Rev. D, 81, 123016
  * Hinderer et al. (2019) Hinderer, T., Nissanke, S., Foucart, F., et al. 2019, Phys. Rev. D, 100, 063021
  * Hotokezaka & Nakar (2020) Hotokezaka, K., & Nakar, E. 2020, ApJ, 891, 152
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90
  * Huth et al. (2021) Huth, S., Wellenhofer, C., & Schwenk, A. 2021, Phys. Rev. C, 103, 025803
  * Jiang et al. (2020) Jiang, J.-L., Tang, S.-P., Wang, Y.-Z., Fan, Y.-Z., & Wei, D.-M. 2020, ApJ, 892, 55
  * Kasliwal et al. (2017) Kasliwal, M. M., Nakar, E., Singer, L. P., et al. 2017, Science, 358, 1559
  * Kasliwal et al. (2019) Kasliwal, M. M., Kasen, D., Lau, R. M., et al. 2019, MNRAS, arXiv:1812.08708
  * Kastaun & Ohme (2019) Kastaun, W., & Ohme, F. 2019, Phys. Rev. D, 100, 103023
  * Kaufmann et al. (2020) Kaufmann, S., et al. 2020, Phys. Rev. Lett., 124, 132502
  * Kelley (2021) Kelley, L. Z. 2021, Journal of Open Source Software, 6, 2784
  * Khan et al. (2016) Khan, S., Husa, S., Hannam, M., et al. 2016, Phys. Rev. D, 93, 044007
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt, IOS Press, 87 – 90
  * Korobkin et al. (2021) Korobkin, O., Wollaeger, R. T., Fryer, C. L., et al. 2021, ApJ, 910, 116
  * Krüger & Foucart (2020) Krüger, C. J., & Foucart, F. 2020, Phys. Rev. D, 101, 103002
  * Kulkarni (2005) Kulkarni, S. R. 2005, arXiv:astro-ph/0510256
  * Kullback & Leibler (1951) Kullback, S., & Leibler, R. A. 1951, Ann. Math. Statist., 22, 79
  * Kyutoku et al. (2020) Kyutoku, K., Fujibayashi, S., Hayashi, K., et al. 2020, ApJ, 890, L4
  * Landry & Essick (2019) Landry, P., & Essick, R. 2019, Phys. Rev. D, 99, 084049
  * Landry et al. (2020) Landry, P., Essick, R., & Chatziioannou, K. 2020, Phys. Rev. D, 101, 123007
  * Lattimer & Lim (2013) Lattimer, J. M., & Lim, Y. 2013, ApJ, 771, 51
  * Li et al. (2020) Li, A., Jiang, J.-L., Tang, S.-P., et al. 2020, arXiv e-prints, arXiv:2009.12571
  * Li & Paczyński (1998) Li, L.-X., & Paczyński, B. 1998, ApJ, 507, L59
  * Lindblom (2018) Lindblom, L. 2018, Phys. Rev. D, 97, 123019
  * Lynn et al. (2016) Lynn, J. E., Tews, I., Carlson, J., et al. 2016, Phys. Rev. Lett., 116, 062501
  * Maggiore et al. (2020) Maggiore, M., et al. 2020, JCAP, 03, 050
  * Margalit & Metzger (2017) Margalit, B., & Metzger, B. D. 2017, ApJ, 850, L19
  * Margutti et al. (2017) Margutti, R., Berger, E., Fong, W., et al. 2017, ApJ, 848, L20
  * Metzger et al. (2010) Metzger, B. D., Martínez-Pinedo, G., Darbha, S., et al. 2010, MNRAS, 406, 2650
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24
  * Miller et al. (2021) —. 2021, arXiv e-prints, arXiv:2105.06979
  * Nathanail et al. (2021) Nathanail, A., Most, E. R., & Rezzolla, L. 2021, ApJ, 908, L28
  * Nedora et al. (2021) Nedora, V., Bernuzzi, S., Radice, D., et al. 2021, ApJ, 906, 98
  * Nicholl et al. (2021) Nicholl, M., Margalit, B., Schmidt, P., et al. 2021, arXiv e-prints, arXiv:2102.02229
  * Nicholl et al. (2017) Nicholl, M., Berger, E., Kasen, D., et al. 2017, ApJ, 848, L18
  * O’Boyle et al. (2020) O’Boyle, M. F., Markakis, C., Stergioulas, N., & Read, J. S. 2020, Phys. Rev. D, 102, 083027
  * Oliphant (2007) Oliphant, T. E. 2007, Computing in Science Engineering, 9, 10
  * Perez & Granger (2007) Perez, F., & Granger, B. E. 2007, Computing in Science Engineering, 9, 21
  * Pian et al. (2017) Pian, E., D’Avanzo, P., Benetti, S., et al. 2017, Nature, 551, 67
  * Raaijmakers et al. (2019) Raaijmakers, G., Riley, T. E., Watts, A. L., et al. 2019, ApJ, 887, L22
  * Raaijmakers et al. (2020) Raaijmakers, G., Greif, S. K., Riley, T. E., et al. 2020, ApJ, 893, L21
  * Raaijmakers et al. (2021a) Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021a, Constraints on the dense matter equation of state and neutron star properties from NICER’s mass-radius estimate of PSR J0740+6620 and multimessenger observations: posterior samples and scripts for generating plots, Zenodo, doi:10.5281/zenodo.4696232
  * Raaijmakers et al. (2021b) Raaijmakers, G., Nissanke, S., Foucart, F., et al. 2021b, arXiv e-prints, arXiv:2102.11569
  * Radice & Dai (2019) Radice, D., & Dai, L. 2019, Eur. Phys. J., A55, 50
  * Radice et al. (2018) Radice, D., Perego, A., Hotokezaka, K., et al. 2018, ApJ, 869, L35
  * Reed et al. (2021) Reed, B. T., Fattoyev, F. J., Horowitz, C. J., & Piekarewicz, J. 2021, Phys. Rev. Lett., 126, 172503
  * Riley et al. (2018) Riley, T. E., Raaijmakers, G., & Watts, A. L. 2018, MNRAS, 478, 1093
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJ, 887, L21
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, arXiv e-prints, arXiv:2105.06980
  * Roca-Maza et al. (2015) Roca-Maza, X., Viñas, X., Centelles, M., et al. 2015, PhRvC, 92, 064304
  * Ruiz et al. (2018) Ruiz, M., Shapiro, S. L., & Tsokaros, A. 2018, Phys. Rev. D, 97, 021501
  * Sathyaprakash et al. (2019a) Sathyaprakash, B., Buonanno, A., Lehner, L., et al. 2019a, BAAS, 51, 251
  * Sathyaprakash et al. (2019b) Sathyaprakash, B., Bailes, M., Kasliwal, M. M., et al. 2019b, BAAS, 51, 276
  * Sekiguchi et al. (2016) Sekiguchi, Y., Kiuchi, K., Kyutoku, K., Shibata, M., & Taniguchi, K. 2016, Phys. Rev. D, 93, 124046
  * Shibata et al. (2017) Shibata, M., Fujibayashi, S., Hotokezaka, K., et al. 2017, Phys. Rev. D, 96, 123012
  * Smartt et al. (2017) Smartt, S. J., Chen, T. W., Jerkstrand, A., et al. 2017, Nature, 551, 75
  * Tanaka et al. (2020) Tanaka, M., Kato, D., Gaigalas, G., & Kawaguchi, K. 2020, MNRAS, 496, 1369
  * Tang et al. (2021) Tang, S.-P., Jiang, J.-L., Gao, W.-H., Fan, Y.-Z., & Wei, D.-M. 2021, Phys. Rev. D, 103, 063026
  * Tanvir et al. (2017) Tanvir, N. R., Levan, A. J., González-Fernández, C., et al. 2017, ApJ, 848, L27
  * Tews et al. (2013) Tews, I., Krüger, T., Hebeler, K., & Schwenk, A. 2013, Phys. Rev. Lett., 110, 032504
  * Troja et al. (2017) Troja, E., Piro, L., van Eerten, H., et al. 2017, Nature, 551, 71
  * Tsang et al. (2012) Tsang, M. B., Stone, J. R., Camera, F., et al. 2012, Phys. Rev. C, 86, 015803
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science Engineering, 13, 22
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Xie & Li (2021) Xie, W.-J., & Li, B.-A. 2021, Phys. Rev. C, 103, 035802
  * Yang et al. (2018) Yang, H., East, W. E., & Lehner, L. 2018, ApJ, 856, 110