# Needle in a Bayes Stack: a Hierarchical Bayesian Method for Constraining the
Neutron Star Equation of State with an Ensemble of Binary Neutron Star Post-
merger Remnants

Alexander W. Criswell<EMAIL_ADDRESS>Minnesota Institute for Astrophysics,
University of Minnesota, Minneapolis, Minnesota 55455, USA  Jesse Miller
School of Statistics, University of Minnesota, Minneapolis, Minnesota 55455,
USA  Noah Woldemariam  Department of Computer Science and Engineering,
University of Minnesota, Minneapolis, Minnesota 55455, USA  Theodoros
Soultanis  GSI Helmholtzzentrum für Schwerionenforschung, Planckstraße 1,
64291 Darmstadt, Germany  Heidelberg Institute for Theoretical Studies,
Schloss-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany  Max-Planck-Institut für
Astronomie, Königstuhl 17, 69117 Heidelberg, Germany  Andreas Bauswein  GSI
Helmholtzzentrum für Schwerionenforschung, Planckstraße 1, 64291 Darmstadt,
Germany  Helmholtz Research Academy Hesse for FAIR (HFHF), GSI Helmholtz
Center for Heavy Ion Research, Campus Darmstadt, 64291 Darmstadt, Germany
Katerina Chatziioannou  Department of Physics, California Institute of
Technology, Pasadena, California 91125, USA  LIGO Laboratory, California
Institute of Technology, Pasadena, California 91125, USA  Michael W. Coughlin
Minnesota Institute for Astrophysics, University of Minnesota, Minneapolis,
Minnesota 55455, USA  Galin Jones  School of Statistics, University of
Minnesota, Minneapolis, Minnesota 55455, USA  Vuk Mandic  Minnesota Institute
for Astrophysics, University of Minnesota, Minneapolis, Minnesota 55455, USA

###### Abstract

Binary neutron star (BNS) post-merger gravitational-wave emission can occur in
the aftermath of a BNS merger — provided the system avoids prompt collapse to
a black hole — as a quasistable hypermassive remnant experiences quadrupolar
oscillations and non-axisymmetric deformations. The post-merger gravitational-
wave spectrum possesses a characteristic peak frequency that has been shown to
be dependent on the binary chirp mass and the neutron star equation of state
(EoS), rendering post-merger gravitational waves a powerful tool for
constraining neutron star composition. Unfortunately, the BNS post-merger
signal is emitted at high ($\gtrsim 1.5$ kHz) frequencies, where ground-based
gravitational wave detectors suffer from reduced sensitivity. It is therefore
unlikely that post-merger signals will be detected with sufficient signal-to-
noise ratio (SNR) until the advent of next-generation detectors. However, by
employing empirical relations derived from numerical relativity simulations,
we can combine information across an ensemble of BNS mergers, allowing us to
obtain EoS constraints with many low-SNR signals. We present a hierarchical
Bayesian method for deriving constraints on $R_{1.6}$, the radius of a
1.6$\mathrm{M_{\odot}}$ neutron star, through an ensemble analysis of binary
neutron star mergers. We apply this method to simulations of the next two
LIGO-Virgo-KAGRA observing runs, O4 and O5, as well as an extended 4-year run
at A+ sensitivity, demonstrating the potential of our approach to yield EoS
information from the post-merger signal with current-generation detectors. The
A+ 4-year scenario is predicted to improve the constraint on $R_{1.6}$ from
the currently available multimessenger-based 95% credible interval (C.I.)
uncertainty of $R_{1.6}=12.07^{+0.98}_{-0.77}$ km to
$R_{1.6}=11.91^{+0.80}_{-0.56}$ km, a 22% reduction of the 95% C.I. width.

##  I Introduction

###  I.1 On Neutron Star Interiors

The nature of matter at the temperatures and densities present in neutron star
(NS) interiors is an open question. This is, in some ways, unsurprising; NSs
are among the most extreme material objects in the universe. Highly
relativistic and immensely dense, with masses/radii of order 1-2
$\mathrm{M_{\odot}}$/10-15 km, NSs reside in the final stable bastion of
matter before collapse into a black hole and are the sole known occupants of
the cold, high-density regime of dense nuclear matter [1]. As such, NSs
provide a powerful tool for exploring dense nuclear theory.

In practice, one wishes to constrain the NS equation of state (EoS), which
describes pressure as a function of mass-energy density (or, equivalently,
mass as a function of radius). Such constraints are realized through a variety
of NS observables, including direct mass-radius measurements of individual
pulsars by NICER (e.g., [2, 3, 4, 5]), observations of individual high-mass
NSs [6, 7, 8, 9, 10], measurements of the NS tidal deformability $\Lambda$
from binary neutron star (BNS) inspiral gravitational wave (GW) signals by the
LIGO-Virgo-KAGRA (LVK) detectors [11, 12, 13], and studies of the kilonova
(AT2017gfo) [14, 15, 16, 17, 18, 19, 20] and gamma ray burst (GRB170817) [14,
21, 22] counterparts to the sole BNS merger detected in both GWs and
electromagnetic radiation, GW170817. By combining the information contained in
each of these measurements in a cohesive multimessenger framework, one can
arrive at the most stringent available constraints. At the time of writing, a
number of complementary studies have presented such collective constraints,
reaching similar conclusions: a (95% credible interval (C.I.)) constraint on
the radius of a $1.6\mathrm{M_{\odot}}$ NS of $R_{1.6}=12.07^{+0.98}_{-0.77}$
km [23, 24, 3, 25, 26, 27, 28, 29, 30, 31, 32, 5, 33, 34, 35].

###  I.2 Binary Neutron Star Mergers and Their Remnants

BNS mergers consist of two phases, both of which produce GWs: inspiral and
post-merger. The inspiral signal encodes a wealth of information about the
merging binary, including the binary masses, spins, distance, redshift, tidal
deformation, and so on, although some of these parameters are degenerate with
each other. Importantly for this work, the inspiral frequency evolution is
dependent (to leading order) on a specific combination of the binary masses,
known as the chirp mass:

$\mathcal{M}=\frac{(m_{1}m_{2})^{3/5}}{(m_{1}+m_{2})^{1/5}}\,,$ (1)

where $m_{1}$ and $m_{2}$ are the individual NS masses. The majority of the
frequency content and frequency evolution of a GW binary inspiral signal
occurs in the most sensitive band of the LIGO [36], Virgo [37], and KAGRA [38,
39] detectors, allowing inspiral analyses to measure the chirp mass
$\mathcal{M}$ to extremely high precision. For example, the chirp mass of
GW170817 was measured to be $1.188^{+0.004}_{-0.002}\mathrm{M_{\odot}}$ (90%
C.I.) [11].

After the inspiral phase, there are four possible outcomes for the post-merger
remnant object: direct collapse to a black hole, the formation of a short-
lived hypermassive neutron star (HMNS), the formation of a long-lived
supramassive neutron star (SMNS), and the formation of a stable NS. Given our
understanding of the BNS mass distribution, it is likely that a large fraction
of BNS events will result in the second scenario: formation of a HMNS that
collapses to a black hole within a few seconds [40, 41]. Indeed, it has been
argued that this was the nature of GW170817’s remnant (e.g., [42, 43, 44, 45,
46, 47, 48, 49, 50].) This scenario is of particular interest to us, as it is
accompanied by a post-merger GW signal with frequencies of 1.5-4 kHz and
properties that are strongly tied to the NS EoS.

HMNSs exceed the maximum allowed mass of a cold, uniformly-rotating NS;111Thus
surpassing the limitations of the SMNS scenario, which exceeds the maximum
allowed mass of a cold, non-spinning NS through some degree of uniform
rotation. such a configuration should, in principle, be unstable and undergo
gravitational collapse. However, as shown in Baumgarte  _et al._ [40], a
differentially-rotating NS featuring a higher rotation rate in its core than
its envelope — of the kind formed in many BNS mergers [51] — can temporarily
maintain support for much higher222An increase of 30-70%; see Bauswein  _et
al._ [52]. masses against gravitational collapse into a black hole. This
support can be compromised by neutrino losses, gravitational-wave emission,
and angular momentum transport within the NS reducing the degree of
differential rotation present [53, 41]. Precisely modelling these quenching
processes is challenging, and as such the exact timescale for loss of support
from differential rotation is difficult to specify. However, it is generally
agreed that this process happens on relatively short timescales [53, 41].

Numerical relativity simulations of hypermassive BNS remnants have shown that
the dominant frequency of GW emission in such remnants is driven by $f$-mode
pulsations333The spherical harmonic $l=m=2$ mode of oscillation. See
Stergioulas  _et al._ [54] for details concerning the generation of GWs in
post-merger remnants through $f$-mode oscillations and non-axisymmetric
deformations. A review of spherical harmonic modes as applied to
asteroseismology of compact objects can be found in Stergioulas [55]. and the
rotating bar-shaped structure of the remnant that, due to their quadrupolar
nature, are efficient GW emitters; additionally, the post-merger spectrum
possesses a number of subdominant peaks arising from other oscillation modes
or non-linear dynamical features [56, 57, 54, 58, 59, 60]. We follow the
convention of referring to this characteristic dominant frequency as the post-
merger peak frequency $f_{\mathrm{peak}}$.444Sometimes referred to elsewhere
as $f_{2}$. In Fig. 1, we plot the spectrum of a simulated BNS gravitational
waveform, drawing attention to the post-merger peak frequency.

Figure 1: GW power spectral density (PSD) of a simulated merger + post-merger
waveform for the SFHX EoS (see §II.6) with $m_{1}=1.03\mathrm{M_{\odot}}$ and
$m_{2}=1.09\mathrm{M_{\odot}}$, injected at a distance of 50 Mpc. The current,
projected, and proposed detector noise PSDs discussed in this work are
provided for reference: that of the LIGO Livingston detector during the 3rd
LVK observing run (O3 L1); the projected LIGO sensitivities for the 4th and
5th LVK observing runs (O4 and O5), assuming LIGO A+ design sensitivity for O5
[61]; and the proposed NEMO and Einstein Telescope detectors (see § IV.1) [62,
63]. The post-merger peak frequency of 2.67 kHz is marked.

Interestingly, $f_{\mathrm{peak}}$ has been shown to be strongly dependent on
the NS EoS and, for fixed mass, shows significant correlations with a number
of EoS-dependent quantities. This EoS dependence has been for instance
described in terms of the tidal coupling constant $\kappa_{2}^{T}$ [64, 58]
and the radii of cold, non-rotating NS of a fixed mass (e.g., that of a
$1.6\mathrm{M_{\odot}}$ NS, $R_{1.6}$) [65, 66]. These correlations imply that
the post-merger GW signal (or, at least, its peak frequency) is primarily
determined by the merger mass and the NS EoS, and can serve as a powerful tool
for constraining the nuclear EoS. Additionally, if EoS constraints arising
from post-merger studies differ significantly from those arising from inspiral
measurements or other studies of cold, stable NSs, this could be indicative of
phase transitions occurring during the merger process and provide evidence for
more exotic matter in NS cores such as deconfined quark matter [67].

Unfortunately, the post-merger GW signal is predicted to be difficult to
detect [68, 69]; while its total emitted power is in fact larger than that of
the inspiral signal, the post-merger signal is emitted at high (1.5-4 kHz)
frequencies, where the LVK detectors are less sensitive. Unambiguous
individual detection of post-merger GW by current-generation detectors is
unlikely unless a suitable BNS merger occurs within $\lesssim 20$ Mpc,
corresponding to a post-merger signal-to-noise ratio (SNR)$~{}\gtrsim 4-8$
[68, 70].

Indeed, while there have been two detected BNS mergers in the first three
LIGO/Virgo observing runs, no post-merger GW emission has been detected. For
GW170817, searches were performed for post-merger GWs of short ($\lesssim$ 1 s
- the variety considered in this work), intermediate ($\lesssim$ 500 s), and
long ($\lesssim$ 8.5 days) durations.555We limit this discussion to searches
for short-duration signals from HMNS remnants; see Abbott  _et al._ [71] and
Abbott  _et al._ [72] for details on the searches for intermediate- and long-
duration signals, respectively. The search for short-duration signals employed
two primary approaches. The first was a maximum-likelihood evaluation of
coherent excess power using the Coherent WaveBurst (cWB) algorithm [73]. The
second was a BayesWave search of the kind outlined in Chatziioannou  _et al._
[74] (see §II.3). Neither search yielded a detection of the post-merger
signal; the cWB search placed upper limits on the root-sum-square GW strain
from 1-4 kHz of $h_{rss}^{50\%}=8.4\times 10^{-22}\text{ Hz}^{-1/2}$ [71] and
the BayesWave analysis placed 90% credible upper limits on the network SNR
from 1-4 kHz of $\rho_{\text{net}}\leq 6.7$, as well as frequency-dependent
90% upper limits on the strain amplitude spectral density and the spectral
energy density [75]. These upper limits are larger than any expected strain
produced by the post-merger remnant of a GW170817-like event; the LIGO/Virgo
high-frequency sensitivity was simply insufficient to detect the post-merger
signal [71, 75].

For GW190425, only a BayesWave search was performed. No post-merger signal was
detected and 90% credible upper limits were placed on the strain amplitude
spectral density of $1.1\times 10^{-22}\text{ Hz}^{-1/2}$ for a frequency of
2.5 kHz, higher than those placed for GW170817 due to increased source
distance [13]. It is worth noting in this case that not only is this strain
upper limit higher than would be expected for a post-merger signal, but the
inferred binary masses for GW190425 are so large that the likely post-merger
scenario for this binary is direct collapse to a black hole (and therefore no
post-merger GW signal emitted in the first place) [52, 76, 13].

Taken together, this indicates that unambiguous individual post-merger
detections will need to wait until next-generation detectors [69]. With recent
improvements to analytic post-merger models [77, 78, 79, 60], Einstein
Telescope and Cosmic Explorer are predicted to observe $\lesssim 2-4$ and
$\lesssim 10$ post-merger signals with sufficient SNR per year, respectively
[80, 81]; these rates can be elevated to $\lesssim 40-120$ if one considers a
current + next-generation network, multiple next-generation detectors, or
tuning Cosmic Explorer to optimally detect the post-merger signal [82, 83].
Given that the prospects for unambiguous individual detection of post-merger
GWs in the near future seem bleak — and thus their constraining power for the
NS EoS appears, for now, out of reach — it is worth considering how we might
make use of the ambiguous and marginal detections we may be able to achieve
with current-generation detectors.

One promising approach is to combine many such signals in a coherent fashion,
allowing us to learn what the ensemble of BNS mergers can tell us in concert,
rather than individually, about the NS EoS. Previous studies have proposed
coherent mode stacking using empirical relations and inspiral mass
measurements [84] or phase information from the inspiral signal [85], although
this latter method is predicted to only produce meaningful results for next-
generation detectors. There remains, however, another potential tactic for
aggregating information about the NS EoS across many marginal post-merger
signals: that of hierarchical Bayesian inference.

###  I.3 Hierarchical Bayesian Inference

Bayesian inference is a common framework for parameter estimation in GW
astronomy. In its most basic application, we wish to calculate
$p(\vec{\theta}|d)$, the posterior probability density of a set of model
parameters $\vec{\theta}=\langle\theta_{1}...\theta_{M}\rangle$, given some GW
data $d$. We use Bayes’ theorem to do so:

$p(\vec{\theta}|d)=\frac{p(d|\vec{\theta})}{p(d)}p(\vec{\theta})\propto
p(d|\vec{\theta})p(\vec{\theta}),$ (2)

where $p(d|\vec{\theta})$ is the likelihood function of the data given a
signal model parameterized by $\vec{\theta}$, $p(d)$ is the evidence, and
$p(\vec{\theta})$ is the prior probability of the parameters $\vec{\theta}$.
If we have multiple independent events $i=1...N$, the joint likelihood is

$p(d_{1}...d_{N}|\vec{\theta})=\prod_{i=1}^{N}p(d_{i}|\vec{\theta}),$ (3)

and the posterior is

$\begin{split}p(\vec{\theta}|d_{1}...d_{N})&=\frac{1}{p(d)}\prod_{i=1}^{N}\left[p(d_{i}|\vec{\theta})\right]p(\vec{\theta})\\\
&\propto\prod_{i=1}^{N}\left[p(d_{i}|\vec{\theta})\right]p(\vec{\theta}),\end{split}$
(4)

where $p(d)\equiv p(d_{1}...d_{N})$, provided that we expect our model
parameters $\vec{\theta}$ to be identical for all N events. If this is not the
case — if each event is expected to be described by its own individual
$\vec{\theta_{i}}$, where each $\vec{\theta_{i}}$ is informed by a shared
process or characteristic across all events — then we must turn to a more
general case of hierarchical Bayesian inference.

In hierarchical Bayesian inference, each event’s individual model parameters
$\vec{\theta}_{i}$ are drawn from some common distribution described by a
shared parameter (or parameters) $\phi$. This creates a hierarchy where the
population-level parameter $\phi$ informs the individual event parameters
$\vec{\theta}_{i}$, which in turn describe the data $d$. Using Bayes’ theorem
and the definition of conditional probability
$p(\vec{\theta},\phi)=p(\vec{\theta}|\phi)p(\phi)$, we can write for a single
event

$\begin{split}p(\vec{\theta},\phi|d)&=\frac{p(d|\vec{\theta},\phi)p(\vec{\theta},\phi)}{p(d)}\\\
&=\frac{p(d|\vec{\theta})p(\vec{\theta}|\phi)p(\phi)}{p(d)}\\\ &\propto
p(d|\vec{\theta})p(\vec{\theta}|\phi)p(\phi).\end{split}$ (5)

This hierarchy is often described in shorthand using the following notation:

$\begin{split}d|\vec{\theta}&\sim p(d|\vec{\theta}),\\\ \vec{\theta}|\phi&\sim
p(\vec{\theta}|\phi),\\\ \phi&\sim p(\phi).\end{split}$ (6)

To generalize to $N$ events, we can follow the same procedure as we did for
non-hierarchical inference, taking a product of the individual event
likelihoods:

$\begin{split}p(\vec{\theta}_{1...N},\phi|d_{1...N})&=\frac{1}{p(d)}\prod_{i=1}^{N}\mathrlap{\left[p(d_{i}|\vec{\theta}_{i})p(\vec{\theta}_{i}|\phi)\right]p(\phi)}\\\
&\propto\prod_{i=1}^{N}\left[p(d_{i}|\vec{\theta}_{i})p(\vec{\theta}_{i}|\phi)\right]p(\phi).\end{split}$
(7)

where again $p(d)\equiv p(d_{1}...d_{N})$. This returns a population-informed
posterior on both the individual event parameters $\vec{\theta}_{i}$ and the
population-level parameter, $\phi$. If we are ultimately only interested in
$\phi$, we can marginalize over $\vec{\theta}_{i}$, treating them as
“nuisance”, or latent parameters and recovering a marginalized posterior
distribution on $\phi$ only:

$p(\phi|d_{1...N})\propto\prod_{i=1}^{N}\left[\int
p(d_{i}|\vec{\theta}_{i})p(\vec{\theta}_{i}|\phi)d\vec{\theta}_{i}\right]p(\phi).$
(8)

Hierarchical Bayesian methods are often applied in astronomy to situations
where a set of astrophysical objects or events have differing individual
properties which all stem from the same underlying physical processes.
Examples include the spatial distribution of trans-Neptunian objects [86],
exoplanet populations with a shared mass-radius relationship [87], galactic
satellites whose orbital dynamics are dominated by their host galaxy [88],
ensembles of BNS inspiral signals with tidal signatures determined by
$\Lambda(m)$ [89], observations of compact binary mergers arising from an
underlying population distribution and cosmology [90, 91, 92] or, more
pertinently, BNS post-merger remnants with various individual values of
$\mathcal{M}$ and $f_{\mathrm{peak}}$ but a shared EoS.

In this work, we present a method for combining post-merger information from
an ensemble of BNS mergers through a hierarchical Bayesian schema, in order to
constrain $R_{1.6}$. We utilize quasi-universal, EoS-agnostic empirical
relations based on those of Vretinaris  _et al._ [93] to connect posterior
distributions for $\mathcal{M}$ and $f_{\mathrm{peak}}$ to the EoS proxy
$R_{1.6}$, allowing for collective inference of the latter quantity across the
ensemble of BNS events. This approach is presented in full in §II and its
efficacy is explored in §III, including application to simulations of the next
two LVK observing runs (O4 and O5) and an extended run at LIGO A+ sensitivity.

##  II Methods

###  II.1 Empirical Relations

It has been shown via large numbers of numerical relativity simulations that
there are quasi-universal, EoS-agnostic empirical relations connecting
$f_{\mathrm{peak}}$, $\mathcal{M}$, and $R_{1.6}$ [66, 65, 93]. While
Vretinaris  _et al._ [93] have calculated a wide variety of empirical
relations for $R_{1.\text{x}}\;(\text{x}\in[2,4,6,8])$ and several different
post-merger signal characteristics, we use a relation of the form
$f_{\mathrm{peak}}=f(R_{1.6},\mathcal{M})$, based on Equation 4 of Vretinaris
_et al._ [93],666Recently, it has been speculated that these relations may be
less tight, especially for EoSs with significant positive M-R slope [94].
However, a re-analysis of previous results does not support this conclusion
for microphysical models with similar variations of the slope in the M-R
curves [65]. Also, using the large sample of microphysical EoSs of Vretinaris
_et al._ [93], we do not find evidence that the introduction of an
$R_{1.4}/R_{1.8}$ term as proposed in Raithel and Most [94] would improve the
fits (for an optimally chosen radius) nor do we find the same magnitude of
scatter in the relations for $f_{\mathrm{peak}}$ (see also Lioutas  _et al._
[95] for a discussion of frequency deviations). which possessed the lowest
residuals of the relations considered in Vretinaris  _et al._ [93]. However,
any of the relations for $R_{1.\text{x}}$ could be used with this procedure.
The original relation of Vretinaris  _et al._ [93] specifies a model of the
chirp-mass-scaled peak frequency $f_{\mathrm{peak}}/\mathcal{M}$ as a function
of $\mathcal{M}$ and $R_{1.6}$; that is,
$f_{\mathrm{peak}}/\mathcal{M}=f(\mathcal{M},R_{1.6})+\epsilon$, where
$\epsilon\sim\mathcal{N}(0,\sigma^{2})$ is the error term. This exact form
would suffice if we only needed a point estimate of $f_{\mathrm{peak}}$ as a
function of $\mathcal{M}$ and $R_{1.6}$, as we could simply multiply both
sides by $\mathcal{M}$. However, we must account for the uncertainty in the
fitted values we obtain from this relation, and multiplying both sides by
$\mathcal{M}$ produces a relation of the form
$f_{\mathrm{peak}}=g(\mathcal{M},R_{1.6})+\mathcal{M}\times\epsilon$, which
has an error term dependent on $\mathcal{M}$. To avoid this, we consider a
modified relation, based on the original Vretinaris  _et al._ [93] relations,
of the form

$\begin{split}f_{\mathrm{peak}}=\beta_{0}\
&+\beta_{1}\mathcal{M}+\beta_{2}\mathcal{M}^{2}+\beta_{3}R_{1.6}\mathcal{M}\\\
&+\beta_{4}R_{1.6}\mathcal{M}^{2}+\beta_{5}R_{1.6}^{2}\mathcal{M},\end{split}$
(9)

with an additional error term $\epsilon\sim\mathcal{N}(0,\sigma^{2})$, where
$\beta_{0}$ through $\beta_{5}$ are our fit coefficients, and
$f_{\mathrm{peak}}$, $\mathcal{M}$, and $R_{1.6}$ are in kHz,
$\mathrm{M_{\odot}}$, and km, respectively. Fitting this model to the data in
Tables I and II in Vretinaris  _et al._ [93] using least squares, we obtain
the fit coefficients in Table 1 with a maximum residual of $0.185$ kHz and a
residual standard error $0.061$ kHz. We obtain a set of relations that is
consistent with those of Vretinaris  _et al._ [93] to within the intrinsic
error of their relations. We also obtain the covariance matrix $\hat{\Sigma}$
for the vector of estimated coefficients
$\hat{\beta}=(\hat{\beta}_{0},\ldots,\hat{\beta}_{5})$. This enables us to
characterize the intrinsic empirical relation uncertainty by sampling
$N_{\mathrm{R}}=1000$ relations from the multivariate normal distribution
$\mathcal{N}(\hat{\beta},\hat{\Sigma})$. These fits and their corresponding
sampled relations are shown in Fig. 2.

Coefficient | Value | Unit
---|---|---
$\beta_{0}$ | $1.5220$ | kHz
$\beta_{1}$ | $8.4021$ | kHz/$\mathrm{M_{\odot}}$
$\beta_{2}$ | $2.3876$ | kHz/$\mathrm{M}_{\odot}^{2}$
$\beta_{3}$ | $-1.1133$ | kHz/(km $\mathrm{M_{\odot}}$)
$\beta_{4}$ | $-0.1291$ | kHz/(km $\mathrm{M}_{\odot}^{2}$)
$\beta_{5}$ | $0.0366$ | kHz/($\mathrm{km}^{2}$ $\mathrm{M_{\odot}}$)
Table 1: Fit coefficients for the relation of Eq. (9). Figure 2: Distribution
of $N_{\mathrm{R}}=1000$ sampled relations, evaluated at a range of values of
$R_{1.6}$, alongside our fitted empirical relations, likewise evaluated over
the same range. The sampled distribution is consistent with the fit relation,
and characterizes the relation’s intrinsic uncertainty.

It is worth noting that these relations are derived for $\mathcal{M}$ and
$f_{\mathrm{peak}}$ in the source rest frame. These quantities as measured by
GW detectors include some amount of cosmological redshift, which must be
accounted for in order to accurately employ the empirical relations. Inspiral
source frame $\mathcal{M}$ posteriors are a standard data product, but we need
to manually correct for redshift in our $f_{\mathrm{peak}}$ posteriors. This
correction is of order 5-10% in $f_{\mathrm{peak}}$ for the most distant (and
therefore least informative) BNS events we consider in this study. For the
purposes of this work, we make the simplifying assumption that the redshift of
each event is known. However, we note that the extension of our formalism to
include inspiral redshift posteriors is straightforward and should be
incorporated in the future.777Not only will next-generation detectors be able
to observe more events at appreciable redshifts, but Haster  _et al._ [96]
also showed that high-redshift post-merger signals contribute more information
due to the peak frequency being redshifted towards the detectors’ more
sensitive band.

###  II.2 Derivation of Hierarchical Framework

We wish to glean whatever information on $R_{1.6}$ is contained in the GW
observations of an ensemble of BNS mergers. We can do this by considering
measurements of the chirp mass (via the inspiral signal; see §II.4) and the
post-merger peak frequency (see §II.3), connected through the empirical
relations of Vretinaris  _et al._ [93] (see §II.1). We construct a statistical
hierarchy such that

$\begin{split}D_{i}|\mathcal{M}_{i},f_{\mathrm{peak},i}&\sim
p(D_{\mathrm{IN},i}|\mathcal{M}_{i})p(D_{\mathrm{PM},i}|f_{\mathrm{peak},i}),\\\
\mathcal{M}_{i},f_{\mathrm{peak},i}|R_{1.6}&\sim
p(\mathcal{M}_{i},f_{\mathrm{peak},i}|R_{1.6}),\\\
R_{1.6}&\sim\pi(R_{1.6}),\\\ \end{split}$ (10)

where $D_{i}$ is the collected GW data of the $i^{th}$ BNS merger, consisting
of its inspiral data $D_{\mathrm{IN},i}$ and post-merger data
$D_{\mathrm{PM},i}$. The first line incorporates the inspiral and post-merger
analyses, and makes an important assumption: that the chirp mass as measured
from the inspiral signal is conditionally independent of the post-merger peak
frequency from the post-merger signal. Despite the fact that both measured
quantities arise from the same BNS merger, this is a reasonable assumption.
While $\mathcal{M}$ and $f_{\mathrm{peak}}$ are related physically, the data
analysis approaches leading to their posterior distributions are independent
of each other (see §II.3 and §II.4). Parameter estimation performed on the
inspiral signal does not account for post-merger GWs, and is not biased by the
presence/absence or morphology of a post-merger signal [97]. Similarly, the
morphology-independent analysis used to reconstruct the post-merger signal and
compute the peak frequency posterior is entirely agnostic of all system
parameters, including the chirp mass. It would be inaccurate to say that
$\mathcal{M}$ and $f_{\mathrm{peak}}$ are unrelated, but their measurements
can be approximated as statistically independent in the low-SNR limit.888For
louder signals, this issue can be handled by simultaneous analysis of the
inspiral and post-merger signal [98].

We additionally neglect selection effects arising from the increased
sensitivity of our detectors to BNS mergers with higher values of
$\mathcal{M}$. In practice, this effect will need to be accounted for to avoid
bias; however, the means of doing so is known [90] and will be incorporated in
the future.

Ultimately, we want to calculate the posterior for $R_{1.6}$ given $N$ events’
post-merger and inspiral data ($D_{\mathrm{PM}}$ and $D_{\mathrm{IN}}$,
respectively):

$p(R_{1.6}|D)=p(R_{1.6}|(D_{\mathrm{PM}},D_{\mathrm{IN}})_{i...N}).$ (11)

This is expanded into a likelihood and a prior:

$p(R_{1.6}|D)=\frac{p((D_{\mathrm{PM}},D_{\mathrm{IN}})_{i...N}|R_{1.6})}{p(D)}\pi(R_{1.6}).$
(12)

The likelihood for $N$ events can be written as

$\begin{split}p(D|R_{1.6})&=\prod_{i=1}^{N}p((D_{\mathrm{PM}},D_{\mathrm{IN}})_{i}|R_{1.6})\\\
&=\prod_{i=1}^{N}\bigg{[}\int\int
p(D_{\mathrm{IN},i},D_{\mathrm{PM},i}|\mathcal{M}_{i},f_{\mathrm{peak},i})p(\mathcal{M}_{i},f_{\mathrm{peak},i}|R_{1.6})df_{\mathrm{peak},i}d\mathcal{M}_{i}\bigg{]}\\\
&=\prod_{i=1}^{N}\bigg{[}\int\int
p(D_{\mathrm{IN},i}|\mathcal{M}_{i})p(D_{\mathrm{PM},i}|f_{\mathrm{peak},i})p(\mathcal{M}_{i},f_{\mathrm{peak},i}|R_{1.6})df_{\mathrm{peak},i}d\mathcal{M}_{i}\bigg{]},\end{split}$
(13)

where we have introduced a marginalization over $f_{\mathrm{peak}}$ and
$\mathcal{M}$ (as we are only interested in $R_{1.6}$) and incorporated the
assumption that the inspiral and post-merger measurements are conditionally
independent. As the distributions available on $f_{\mathrm{peak}}$ and
$\mathcal{M}$ are posteriors arising from established GW analyses (see §II.3,
§II.4), we use Bayes’ theorem to expand each of our parameter likelihoods
$p(D_{\mathrm{PM},i}|f_{\mathrm{peak},i})$ and
$p(D_{\mathrm{IN},i}|\mathcal{M}_{i})$ as follows:

$p(D|R_{1.6})=\prod_{i=1}^{N}\bigg{[}\int\int\frac{p(\mathcal{M}_{i}|D_{\mathrm{IN},i})p(D_{\mathrm{IN},i})}{p_{0}(\mathcal{M}_{i})}\frac{p(f_{\mathrm{peak},i}|D_{\mathrm{PM},i})p(D_{\mathrm{PM},i})}{p(f_{\mathrm{peak},i})}p(f_{\mathrm{peak},i},\mathcal{M}_{i}|R_{1.6})df_{\mathrm{peak},i}d\mathcal{M}_{i}\bigg{]}.$
(14)

To expand the likelihoods for $f_{\mathrm{peak}}$ and $\mathcal{M}$ in this
way, the $p_{0}(\mathcal{M})$ and $p(f_{\mathrm{peak}})$ priors we use must be
identical to those used in the initial parameter estimation process. Finally,
we expand
$p(f_{\mathrm{peak},i},\mathcal{M}_{i}|R_{1.6})=p(f_{\mathrm{peak},i}|\mathcal{M}_{i},R_{1.6})p(\mathcal{M}_{i})$,
assuming that the distribution of $\mathcal{M}$ is independent of $R_{1.6}$;
that is, $p(\mathcal{M}|R_{1.6})=p(\mathcal{M})$.999The astrophysical
distribution of NS masses is currently uncertain [99, 92], including whether
the maximum NS mass is determined by the EoS through $M_{\mathrm{TOV}}$ (the
Tolman–Oppenheimer–Volkoff limit for the maximum mass of a non-rotating NS) or
astrophysical processes. At current statistical uncertainties these
assumptions have a small impact on the inferred $M_{\mathrm{TOV}}$; however,
they are negligible for the radius [33]. Though the impact on
$M_{\mathrm{TOV}}$ will become more prominent with additional data, it has
been shown that — for the pertinent case of generic, flexible analyses —
$R_{1.6}$ is largely independent of both $M_{\mathrm{TOV}}$ and such
astrophysical considerations [100], allowing for the assertion that
$p(\mathcal{M}|R_{1.6})=p(\mathcal{M})$. In practice, one should parameterize
$p(\mathcal{M})$ and simultaneously infer $R_{1.6}$ alongside the
astrophysical $\mathcal{M}$ distribution; however, we assume for simplicity
that $p(\mathcal{M})$ is known. At this juncture we must note the distinction
between $p_{0}(\mathcal{M})$, the prior for $\mathcal{M}$ used in the original
inspiral analyses from which we derive our $\mathcal{M}$ posteriors, and
$p(\mathcal{M})$, our own choice of prior for $\mathcal{M}$. While
$p_{0}(\mathcal{M})$ and $p(\mathcal{M})$ can be identical, in a realistic
case $p_{0}(\mathcal{M})$ will have been chosen by others, and may or may not
align with our own choice of prior.101010In the realistic case of a future
post-observing-run analysis, we would use the published $\mathcal{M}$
posterior distributions for each BNS event. This would then necessitate that
we use the $\mathcal{M}$ prior $p_{0}(\mathcal{M})$ chosen for these analyses
— likely the standard uniform priors used in previous LVK catalogues (see,
e.g., GWTC-3 [101]) — to ensure that we correctly recover the likelihood from
the posterior samples available to us. However, our choice of
$p(\mathcal{M})$is not constrained in this fashion and should instead describe
our expectations of the overall population distribution of $\mathcal{M}$. The
specific distributions for $p_{0}(\mathcal{M})$ and $p(\mathcal{M})$ used in
this work are discussed in §II.5.

We can now write the posterior for $R_{1.6}$, factoring out the individual
event evidences in favor of the usual proportionality:

$\begin{split}p(R_{1.6}|D)&=\pi(R_{1.6})\frac{\prod_{k=1}^{N}\big{[}p(D_{\mathrm{IN},k})p(D_{\mathrm{PM},k})\big{]}}{p(D)}\\\
&\quad\times\prod_{i=1}^{N}\bigg{[}\int\int\frac{p(\mathcal{M}_{i}|D_{\mathrm{IN},i})}{p_{0}(\mathcal{M}_{i})}\frac{p(f_{\mathrm{peak},i}|D_{\mathrm{PM},i})}{p(f_{\mathrm{peak},i})}p(f_{\mathrm{peak},i}|\mathcal{M}_{i},R_{1.6})p(\mathcal{M}_{i})df_{\mathrm{peak},i}d\mathcal{M}_{i}\bigg{]}\\\
&\propto\pi(R_{1.6})\prod_{i=1}^{N}\bigg{[}\int\int\frac{p(\mathcal{M}_{i}|D_{\mathrm{IN},i})}{p_{0}(\mathcal{M}_{i})}\frac{p(f_{\mathrm{peak},i}|D_{\mathrm{PM},i})}{p(f_{\mathrm{peak},i})}p(f_{\mathrm{peak},i}|\mathcal{M}_{i},R_{1.6})p(\mathcal{M}_{i})df_{\mathrm{peak},i}d\mathcal{M}_{i}\bigg{]}.\end{split}$
(15)

The term $p(f_{\mathrm{peak},i}|\mathcal{M}_{i},R_{1.6})$ serves to tie our
$f_{\mathrm{peak}}$ and $\mathcal{M}$ posteriors to $R_{1.6}$ while accounting
for the intrinsic scatter of the empirical relations. As discussed in §II.1,
we model the uncertainty of the empirical relation by assuming that the
coefficients in Eq. (9) are well modeled by a multivariate normal distribution
$\mathcal{N}(\hat{\beta},\hat{\Sigma})$ and sampling from this distribution.
Formally, let

$p(f_{\mathrm{peak},i}|\mathcal{M}_{i},R_{1.6})=\int
p(f_{\mathrm{peak},i}|\mathcal{M}_{i},R_{1.6},F_{E})p(F_{E})dF_{E},$ (16)

where the integral is a marginalization over the intrinsic uncertainty of the
empirical relations, represented by a distribution of potential empirical
relations $F_{E}$. We approximate this using our empirical relation samples
$F_{E,j}$ such that

$p(f_{\mathrm{peak},i}|\mathcal{M}_{i},R_{1.6})=\frac{1}{N_{\mathrm{R}}}\sum_{j}p(f_{\mathrm{peak},i}|\mathcal{M}_{i},R_{1.6},F_{E,j}),$
(17)

where $N_{\mathrm{R}}$ is the number of empirical relation samples. We take an
empirical Bayes approach and assume that the prediction of each sampled
relation is exact; that is,

$p(f_{\mathrm{peak},i}|\mathcal{M}_{i},R_{1.6},F_{E,j})=\delta(f_{\mathrm{peak},i}-\hat{f}_{\mathrm{peak},ij}),$
(18)

where the $\delta$ function is a Dirac delta and
$\hat{f}_{\mathrm{peak},ij}=F_{E,j}(\mathcal{M}_{i},R_{1.6})\times(1+z)^{-1}$,
the $j^{th}$ sampled relation’s redshift-corrected111111We assume the redshift
is known; see §II.1. prediction for $f_{\mathrm{peak},i}$ as a function of
$\mathcal{M}$ for a given value of $R_{1.6}$. The posterior then becomes

$\begin{split}p(R_{1.6}|D)\propto&\
\pi(R_{1.6})\prod_{i=1}^{N}\bigg{[}\int\int\frac{p(\mathcal{M}_{i}|D_{\mathrm{IN},i})}{p_{0}(\mathcal{M}_{i})}\frac{p(f_{\mathrm{peak},i}|D_{\mathrm{PM},i})}{p(f_{\mathrm{peak},i})}\frac{1}{N_{\mathrm{R}}}\sum_{j}\big{[}\delta(f_{\mathrm{peak},i}-\hat{f}_{\mathrm{peak},ij})\big{]}p(\mathcal{M}_{i})df_{\mathrm{peak},i}d\mathcal{M}_{i}\bigg{]}\\\
\propto&\
\pi(R_{1.6})\prod_{i=1}^{N}\bigg{[}\int\frac{1}{N_{\mathrm{R}}}\sum_{j}\bigg{(}\frac{p(\mathcal{M}_{i}|D_{\mathrm{IN},i})}{p_{0}(\mathcal{M}_{i})}\frac{p(\hat{f}_{\mathrm{peak},ij}|D_{\mathrm{PM},i})}{p(\hat{f}_{\mathrm{peak},ij})}\bigg{)}p(\mathcal{M}_{i})d\mathcal{M}_{i}\bigg{]}.\\\
\end{split}$ (19)

Finally, we must account for an additional subtlety: we only consider data on
a [1.5, 4] kHz frequency band. This choice is partially driven by the fact
that, at frequencies above 4 kHz, the rapid degradation of the detector
sensitivity becomes a severely limiting factor. More importantly, all
numerical relativity simulations used to fit the Vretinaris  _et al._ [93]
relations have peak frequencies within this band; relying on the empirical
relations’ predictions outside of this regime would necessarily comprise
extrapolation. Thankfully, given the broad coverage of the Vretinaris  _et
al._ [93] data, we can be assured that few (if any) systems will fall outside
of this frequency band.

Despite this choice, some combinations of $R_{1.6}$ and $\mathcal{M}$ —
especially for $R_{1.6}\lesssim 10$ km — yield an empirical-relation-predicted
$\hat{f}_{\mathrm{peak}}$ that lies outside this range. As we do not consider
data from these out-of-band frequencies, any $(\mathcal{M}_{i},R_{1.6})$ pair
that results in such a prediction should not encode any additional
probability, instead returning only the prior probability on
$R_{1.6}$.121212This can be considered an “analysis selection effect” (as
opposed to the more typical “observational selection effect”) wherein some
data is unanalyzed/unobserved depending on what the _observed_ (rather than
the true) parameters are. Accordingly, we split the marginalization over
$\mathcal{M}$ into two integrals: one for each of the in-band ($\mathcal{I}$)
and out-of-band ($\mathcal{O}$) cases. We can then write

$\begin{split}p(R_{1.6}|D)\propto\pi(R_{1.6})\prod_{i=1}^{N}&\bigg{[}\int_{\mathcal{I}}\frac{1}{N_{\mathrm{R}}}\sum_{j}\bigg{(}\frac{p(\mathcal{M}_{i}|D_{\mathrm{IN},i})}{p_{0}(\mathcal{M}_{i})}\frac{p(\hat{f}_{\mathrm{peak},ij}|D_{\mathrm{PM},i})}{p(\hat{f}_{\mathrm{peak},ij})}\bigg{)}p(\mathcal{M}_{i})d\mathcal{M}_{i}\\\
&\qquad\qquad\qquad\qquad\qquad\ +\
\int_{\mathcal{O}}\frac{1}{N_{\mathrm{R}}}\sum_{j}\bigg{(}\frac{p(\mathcal{M}_{i}|D_{\mathrm{IN},i})}{p_{0}(\mathcal{M}_{i})}\frac{p(\hat{f}_{\mathrm{peak},ij}|D_{\mathrm{PM},i})}{p(\hat{f}_{\mathrm{peak},ij})}\bigg{)}p(\mathcal{M}_{i})d\mathcal{M}_{i}\bigg{]},\\\
\end{split}$ (20)

where the integrals over $\mathcal{I}$ and $\mathcal{O}$ integrate over only
$\mathcal{M}$ values associated with in-band and out-of-band values of
$\hat{f}_{\mathrm{peak},ij}(\mathcal{M}_{i},R_{1.6})$, respectively. Finally,
we note that if $\hat{f}_{\mathrm{peak},ij}$ is out-of-band,
$p(\hat{f}_{\mathrm{peak},ij}|D_{\mathrm{PM},i})$ is independent of the post-
merger signal. That is,
$p(\hat{f}_{\mathrm{peak},ij}|D_{\mathrm{PM},i})=p(\hat{f}_{\mathrm{peak},ij})$,
so

$\sum_{j}\frac{p(\hat{f}_{\mathrm{peak},ij}|D_{\mathrm{PM},i})}{p(\hat{f}_{\mathrm{peak},ij})}=\sum_{j}\frac{p(\hat{f}_{\mathrm{peak},ij})}{p(\hat{f}_{\mathrm{peak},ij})}=N_{\mathrm{R}}\times
1.$ (21)

This then (at long last!) allows us to write the final form of our posterior:

$p(R_{1.6}|D)\propto\pi(R_{1.6})\prod_{i=1}^{N}\bigg{[}\int_{\mathcal{I}}\frac{1}{N_{\mathrm{R}}}\sum_{j}\bigg{(}\frac{p(\mathcal{M}_{i}|D_{\mathrm{IN},i})}{p_{0}(\mathcal{M}_{i})}\frac{p(\hat{f}_{\mathrm{peak},ij}|D_{\mathrm{PM},i})}{p(\hat{f}_{\mathrm{peak},ij})}\bigg{)}p(\mathcal{M}_{i})d\mathcal{M}_{i}+\int_{\mathcal{O}}\bigg{(}\frac{p(\mathcal{M}_{i}|D_{\mathrm{IN},i})}{p_{0}(\mathcal{M}_{i})}\bigg{)}p(\mathcal{M}_{i})d\mathcal{M}_{i}\bigg{]}.$
(22)

Each component of Eq. (22) is explored in detail below. The process of
calculating the peak frequency prior and posterior via BayesWave [102] is
described is §II.3; the approximation we employ to simulate the chirp mass
posteriors from inspiral analyses is discussed in §II.4; the choice of priors
is considered in §II.5. Finally, the implementation of this formalism in the
bayestack package is detailed in §II.7.

###  II.3 BayesWave: Recovering the Post-merger Peak Frequency

To obtain posterior distributions for the post-merger peak frequency, we use
the BayesWave package [102, 103, 104] and follow the methods of Chatziioannou
_et al._ [74] and Wijngaarden  _et al._ [98]. In the BayesWave algorithm, each
potential signal is constructed from an arbitrary number of Morlet-Gabor
wavelets with associated wavelet parameters, accounting for sky location,
polarization angle, and ellipticity. BayesWave uses a reversible jump Markov
chain Monte Carlo algorithm [105] to sample the joint and variable-dimension
posterior of these parameters. The resulting distribution of signal
reconstructions maps to a posterior distribution of peak frequencies, where
each $f_{\mathrm{peak}}$ sample is the spectrum maximum on [1.5, 4] kHz of its
corresponding signal reconstruction sample131313Power contribution from the
inspiral and merger is suppressed via windowing in time-domain; see
Chatziioannou  _et al._ [74] for further detail. [74]. If a given signal
reconstruction sample has no such maximum, a peak frequency sample is instead
drawn from the $f_{\mathrm{peak}}$ prior. Sample BayesWave waveform
reconstructions and their corresponding peak frequency posteriors are shown in
Figs. 3 and 4 for strong and weak post-merger signals, respectively, using the
sensitivity of the proposed Neutron Star Extreme Matter Observatory (NEMO)
detector.

We model each BayesWave peak frequency posterior by fitting the posterior
samples with a PESummary [106] reflection-bounded kernel density estimator
(KDE). The use of KDEs to model the underlying probability distribution
associated with a set of marginal posterior samples introduces an inherent
approximation, but allows for rapid post-processing of data products from
otherwise time-intensive parameter estimation processes. Accordingly, this
approximation is often employed in hierarchical gravitational-wave data
analyses (e.g., Lackey and Wade [107], Golomb and Talbot [108]), although
other approximate representations have been proposed, including interpolation
of the individual marginal likelihoods with random forest regressors [89] or
Gaussian processes [109]. Employing a reflection-bounded KDE allows us to
mitigate the edge effects otherwise characteristic of standard KDEs at the
edges of bounded distributions such as those of our $f_{\mathrm{peak}}$
posteriors.

Figure 3: BayesWave signal reconstruction and corresponding peak frequency
posterior for a strong BNS post-merger signal. This BayesWave run was
performed on a SFHX post-merger waveform with $m_{1}=1.03\mathrm{M_{\odot}}$
and $m_{2}=1.09\mathrm{M_{\odot}}$ injected at 30 Mpc at NEMO sensitivity. In
the top panel, the injected waveform is in teal and the median BayesWave
reconstruction with 50% (dark shading) and 90% (light shading) C.I. is shown
in purple. The post-merger signal is reconstructed nicely, and the injected
waveform’s peak frequency $f_{\mathrm{peak}}=2.67$ kHz is recovered to within
$\sim 0.01$ kHz. Figure 4: BayesWave signal reconstruction and corresponding
peak frequency posterior for a marginal BNS post-merger signal. This BayesWave
run was performed on a post-merger waveform with
$m_{1}=1.03\mathrm{M_{\odot}}$ and $m_{2}=1.09\mathrm{M_{\odot}}$ injected at
70 Mpc at NEMO sensitivity. In the top panel, the injected waveform is in teal
and the 50% (dark shading) and 90% (light shading) C.I. of the BayesWave
reconstruction are shown in purple. No median reconstruction is shown because
the reconstructed signal is consistent with zero. Despite this, the spectrum
and corresponding $f_{\mathrm{peak}}$ posterior show a slight preference for
the injected waveform’s peak frequency of 2.67 kHz. Taken by itself, this
recovery of the peak frequency is unconvincing — and would be considered a
nondetection — but nonetheless contains information about the post-merger
signal.

###  II.4 Inspiral Chirp Mass Posteriors

We employ a convenient approximation for the chirp mass posteriors of each
event obtained from the inspiral signal. As shown in Farr  _et al._ [110], the
relative width of the chirp mass posterior distribution scales inversely with
SNR. That is, $\sigma_{\mathcal{M}}/\mathcal{M}\propto(\textrm{SNR})^{-1}$. We
therefore approximate our chirp mass posteriors as draws from
$\mathcal{N}(\mu_{\mathcal{M}},\sigma_{\mathcal{M}})$, where
$\mu_{\mathcal{M}}$ is itself drawn from
$\mathcal{N}(\mathcal{M}_{i,\mathrm{true}},\sigma_{\mathcal{M}})$ to represent
the scatter in the $\mathcal{M}$ posterior mean due to detector noise and
$\sigma_{\mathcal{M}}$ is defined using the chirp mass posterior width
$\sigma_{\textrm{170817}}=0.004\mathrm{M_{\odot}}$,141414This value of
$\sigma_{\textrm{170817}}=0.004\mathrm{M_{\odot}}$ is the upper 90% C.I.
bound; we are being conservative in our estimate of the chirp mass precision.
chirp mass $\mathcal{M}_{\textrm{170817}}=1.188\mathrm{M_{\odot}}$, and
network SNR $\textrm{SNR}_{170817}=32.4$ of the BNS merger GW170817 as our
point of reference [11] such that

$\sigma_{\mathcal{M}_{i}}=\frac{\mathcal{M}_{i}}{\mathcal{M}_{\textrm{170817}}}\frac{\textrm{SNR}_{170817}}{\textrm{SNR}_{i}}\sigma_{\textrm{170817}}.$
(23)

We simulate $N_{\mathrm{R}}=1000$ samples from the resulting distribution for
each event and fit a KDE to the samples as discussed in §II.3. Simulating our
chirp mass posteriors in this way allows us to maintain the expedience
afforded by such an approximation while retaining the ability to generalize to
real data products in the future.

###  II.5 Choice of Priors

For $R_{1.6}$, we consider both uniform and astrophysically motivated priors.
The former will allow us to consider our results in isolation, whereas the
latter will demonstrate what we may be able to contribute within a larger
multimessenger framework. The uniform prior is on $[9,15]$ km. The
astrophysical prior uses the posterior distribution from the multimessenger
constraints on $R_{1.6}$ found in Huth  _et al._ [35]. This distribution is
consistent with current constraints on the NS EoS and was derived through a
combination of maximum NS mass constraints from observed pulsars and the
remnant classification of GW170817, NICER NS mass measurements, GW analyses of
the BNS mergers GW170817 and GW190425 along with their electromagnetic (EM)
counterpart (AT2017gfo in the case of GW170817) or lack thereof (in the case
of GW190425), and data from terrestrial heavy-ion collision experiments [35].
While the results presented in Huth  _et al._ [35] are for $\text{R}_{1.4}$
and $\text{R}_{2.0}$ only, the authors provided all considered EoSs and their
associated posterior probabilities, from which any $\text{R}_{1.\textrm{x}}$
can be calculated, and from which we obtain a distribution on $R_{1.6}$
specifically. Both priors are shown in Fig. 5.

Figure 5: $R_{1.6}$ priors used in this work. One is a uniform prior on
$[9,15]$ km, the other is an astrophysically motivated $R_{1.6}$ prior,
derived from the multimessenger analysis of Huth  _et al._ [35].

Eq. (22) further includes two priors on $\mathcal{M}$ — the original analysis
prior $p_{0}(\mathcal{M})$ and our own prior $p(\mathcal{M})$ — and one on
$f_{\mathrm{peak}}$. We are free to choose $p(\mathcal{M})$, and let
$p(\mathcal{M})=\mathcal{N}(\mu=1.33\mathrm{M_{\odot}},\sigma=0.09\mathrm{M_{\odot}})$
to match the distributions used in Abbott  _et al._ [111] and Petrov  _et al._
[112] as we use the latter to simulate our event catalogues in §III.2. We note
that the choice of $p(\mathcal{M})$ has minimal impact on the final results of
our analysis due to the precision of inspiral $\mathcal{M}$ posteriors.
However, we emphasize that $p_{0}(\mathcal{M})$ and $p(f_{\mathrm{peak}})$ are
intrinsically tied to the choices made in the original analyses for
$\mathcal{M}$ and $f_{\mathrm{peak}}$, respectively, solely serving to cancel
out the original priors used to produce our posterior samples. Determining the
original $p_{0}(\mathcal{M})$ is trivial as we are simulating our
$\mathcal{M}$ posteriors directly with a uniform prior on
$[0.4,4.4]\mathrm{M_{\odot}}$. In principle, $p(f_{\mathrm{peak}})$ can be
obtained from BayesWave by sampling from the prior (i.e., setting the
likelihood to a constant) and then following the same procedure as above to
extract $f_{\mathrm{peak}}$ from the signal reconstruction samples. However,
there is a subtlety here: due to correlations between $f_{\mathrm{peak}}$ and
other parameters $\alpha$ in BayesWave, the marginal prior will not
generically be the same as the conditional prior; that is,
$p(f_{\mathrm{peak}})\neq p(f_{\mathrm{peak}}|\alpha)$. In such cases, the
formally correct approach would be to extend the framework of Eq. (22) to
account for the other parameters (as for example should be done for the joint
inference between EoS and mass distribution discussed above). We defer such an
approach to future work and here instead continue with Eq. (22) where
$p(f_{\mathrm{peak}})$ is the conditional prior, which we obtain through
successive BayesWave analyses of pure detector noise. The resultant prior
distribution is shown in Fig. 6.

Figure 6: Peak frequency prior samples, with a reflection-bounded KDE
approximation of its probability distribution. These samples were computed
from an ensemble of N=50 BayesWave analyses over simulated instrumental noise.

###  II.6 Simulated Post-merger Waveforms

Figure 7: Neutron star mass-radius relationship for each of the equations of
state used in this work, along with each relationship’s intersection with
$\text{M}=1.6\mathrm{M_{\odot}}$. The corresponding values of $R_{1.6}$ are
11.54 km (SLy4), 11.98 km (SFHX), and 13.26 km (DD2). Figure 8: Simulation
parameters (post-merger peak frequency $f_{\mathrm{peak}}$ and chirp mass
$\mathcal{M}$) for each EoS, with the corresponding Vretinaris  _et al._ [93]
empirical relation with $3\sigma$ bounds for that EoS’s value of $R_{1.6}$.

In this study we use a library of simulated post-merger signals created by
relativistic hydrodynamical simulations based on smooth particle hydrodynamics
(SPH). For code details we refer to [113, 114, 115], and note that in this
work we employ a Wendland kernel function (see [116, 117]) as compared to a
cubic spline kernel [118] used in previous studies with this code. The new SPH
kernel features less noise and numerical damping but does not strongly affect
the GW frequency. It leads to slightly tighter peaks in the post-merger
spectrum.

This library is comprised of simulations for three different EoSs: DD2 [119],
SLy4 [120], and SFHX [121], representative of stiffer, softer, and medium
EoSs, respectively. Softer EoSs result in smaller NS radii and in higher post-
merger GW frequencies at the same total binary mass. The DD2, SLy4, and SFHX
EoSs were chosen accordingly to explore the effectiveness of our analysis
throughout the EoS parameter space. The DD2 [119] and SFHX [121] EoSs are
derived using the relativistic mean field approximation, whereas SLy4 [120] is
based on a Skyrme interaction. Mass-radius relationships for each of these
EoSs can be found in Fig. 7. While the DD2 and SFHX models provide the
temperature dependence of the EoS, the SLy4 is only available at zero
temperature assuming neutrinoless beta-equilibrium. Thus, the latter EoS model
is supplemented by an approximate thermal treatment choosing
$\Gamma_{\mathrm{th}}=1.75$ (see Bauswein  _et al._ [122] for details and a
justification of this choice).

For each EoS, we perform simulations for a grid of total binary masses
$M_{\mathrm{tot}}$ and mass ratios $q$. We vary the total binary mass
$M_{\mathrm{tot}}$ in a broad range, but we do not consider systems which
promptly form a black hole after merging. Hence, we only study binaries with a
total mass well below the threshold mass $M_{\mathrm{thres}}$ for direct black
hole formation. In all of these systems the remnants survive for at least 15
ms before collapsing to a black hole. Generally, $M_{\mathrm{thres}}$ depends
on the EoS and exhibits an additional (mild) dependence on the mass ratio $q$
(see e.g. [123]). Table 2 provides the values of $M_{\mathrm{thres}}$ for the
respective EoSs included in this work for mass ratios $q\in\\{0.7,0.85,1.0\\}$
[123]. We consider six binary mass ratios
$q\in\\{0.7,0.8,0.85,0.9,0.95,1.0\\}$. For the SFHX models, we use a finer
$M_{\mathrm{tot}}$ grid, allowing for a deeper exploration of the analysis
using a large number of events for one EoS. In total, the library contains 41
DD2 waveforms, 132 SFHX waveforms, and 41 SLy4 waveforms. The $\mathcal{M}$
and $f_{\mathrm{peak}}$ values for all simulations can be found on Zenodo at
10.5281/zenodo.7007630, and are shown in Fig. 8 along with the empirical
relation $f_{\mathrm{peak}}({\mathcal{M}})$ of Vretinaris  _et al._ [93] for a
constant $R_{1.6}$ fixed for the respective EOS.

EoS |  | $M_{\textrm{thres}}\ (\mathrm{M_{\odot}})$ |
---|---|---|---
| $q=0.7$ | $q=0.85$ | $q=1.0$
SFHX | 2.925 | 2.975 | 2.975
DD2 | 3.275 | 3.325 | 3.325
SLy4 | 2.775 | 2.825 | 2.825
Table 2: Threshold masses for direct collapse to a black hole at multiple
values of the mass ratio $q$ for each EoS in this work.

###  II.7 Implementation: the BAYESTACK Package

The final $R_{1.6}$ posterior is one-dimensional due to our marginalization
over $f_{\mathrm{peak}}$ and $\mathcal{M}$. As we utilize KDEs to model the
$f_{\mathrm{peak}}$ and $\mathcal{M}$ posterior samples for each event, the
low dimensionality of this problem allows us to directly compute the final
$N$-event posterior on a grid, as opposed to more time-consuming stochastic
sampling techniques. Accordingly, we present the hierarchical BAYEsian
STACKing (bayestack) Python package as a direct numerical implementation of
the formalism described in §II.2. This code is publicly available and can be
found on GitHub at https://github.com/criswellalexander/bayestack. The central
algorithm of the bayestack package as employed in this work is described in
brief below.

  1. 1.

A dictionary is created to store the data for $N$ events. Each event is then
assigned its own dictionary to hold its associated data and distributions.

  2. 2.

For each event, the peak frequency posterior samples are loaded. We then use
the PESummary [106] reflection-bounded KDE adaptation of the SciPy [124]
Gaussian KDE scipy.stats.gaussian_kde with a KDE bandwidth of 0.15 to estimate
the $f_{\mathrm{peak}}$ posterior probability density on [1.5,4] kHz.

  3. 3.

A similar procedure is followed for the chirp mass, simulating samples from a
normal distribution as described in §II.4 — or loading pre-computed
$\mathcal{M}$ samples if desired — and using a KDE to estimate the posterior
probability density.

  4. 4.

The pre-computed $f_{\mathrm{peak}}$ prior samples are loaded (see §II.3,
§II.5, and Fig. 6). The prior is fit with a PESummary reflection-bounded KDE.

  5. 5.

The individual event likelihoods are then computed. For every point in
$R_{1.6}$, the integrands of Eq. (22) are calculated on a grid in
$\mathcal{M}$ with $\Delta\mathcal{M}=0.005\mathrm{M_{\odot}}$. For each grid
point:

     * •

The original prior $p_{0}(\mathcal{M})$, prior $p(\mathcal{M})$, and posterior
$p(\mathcal{M}|D_{\textrm{IN}})$ probabilities for $\mathcal{M}$ are computed
and stored.

     * •

The corresponding value of $\hat{f}_{\mathrm{peak},ij}$ is computed for each
of $N_{\mathrm{R}}=1000$ sampled $F_{E,j}$ empirical relations.

     * •

If $\hat{f}_{\mathrm{peak},ij}$ is in-band (see §II.2), the estimated
$f_{\mathrm{peak}}$ prior and posterior probabilities at that value are
stored.

     * •

If $\hat{f}_{\mathrm{peak},ij}$ is out-of-band, the $f_{\mathrm{peak}}$
likelihood is taken to be 1 as discussed in §II.2.

  6. 6.

The grid is then summed over $F_{E,j}$ and $\mathcal{M}$ to yield the
individual marginal likelihood in $R_{1.6}$, $p(D_{i}|R_{1.6})$.

  7. 7.

We then calculate the $N$-event likelihood $\prod_{i=1}^{N}p(D_{i}|R_{1.6})$
given in Eq. (22). In addition to calculating the product across all event
likelihoods, each $p(D_{i}|R_{1.6})$ is also stored individually.

  8. 8.

This final likelihood is then multiplied by an $R_{1.6}$ prior of our choice
to yield the final posterior $p(R_{1.6}|D_{i...N})$. We use either a uniform
prior on [9,15] km or the multimessenger prior of Huth  _et al._ [35] as
described in §II.5.

##  III Results

###  III.1 Model Validation

We perform a set of tests under ideal conditions to ensure the analysis works
as intended. For each of these tests, we follow the above procedure, save that
we directly simulate the $f_{\mathrm{peak}}$ posteriors, so as to have fine
control over the embedded signal strength. Additionally, this choice allows us
the flexibility to consider a broad range of merger masses and EoSs beyond the
bounds of our library of simulated post-merger waveforms. For each event,
$\mathcal{M}_{i,\mathrm{true}}$ is determined by drawing $m_{1}$ and $m_{2}$
from the distributions used in Abbott  _et al._ [111] and Petrov  _et al._
[112] (as we use the events of the latter work in our simulations of future
observing runs; see §III.2).

To simulate $f_{\mathrm{peak}}$ posterior samples, we draw a mixture of
$N_{s}$ signal and $N_{n}$ noise samples. The noise samples are drawn from our
$f_{\mathrm{peak}}$ prior samples. To generate signal samples, we first use
the associated $\mathcal{M}_{i}$ of each event to determine the corresponding
value of $f_{\mathrm{peak},i,\mathrm{true}}$, the intrinsic true peak
frequency of the event. We do so by using the relation of Vretinaris  _et al._
[93], with dispersion reflective of the relation’s intrinsic scatter,
$\sigma_{\mathrm{ER}}=0.053$ kHz. That is, we draw each event’s
$f_{\mathrm{peak},i,\mathrm{true}}$ from
$\mathcal{N}(\mu=\hat{f}_{\mathrm{peak}},\sigma=\sigma_{\mathrm{ER}})$, with
$\hat{f}_{\mathrm{peak}}=F_{E}(R_{\mathrm{true}},\mathcal{M}_{i})$. Unless
otherwise stated, we have (arbitrarily) chosen $R_{1.6,\mathrm{true}}=12.5$ km
for all tests in this section. With the intrinsic
$f_{\mathrm{peak},i,\mathrm{true}}$ in hand, we then draw signal samples for
each event from $\mathcal{N}(\mu=\mu_{f,i},\sigma=\sigma_{f})$, where
$\sigma_{f}=0.1$ kHz and $\mu_{f,i}$ is itself drawn
$\mathcal{N}(\mu=f_{\mathrm{peak},i,\mathrm{true}},\sigma=\sigma_{f})$ to
represent scatter due to detector noise. The posterior width $\sigma_{f}$ is
reflective of low-SNR BayesWave $f_{\mathrm{peak}}$ recoveries (see, e.g., the
lower panel of Fig. 4.)

The result of this process is a set of simulated $f_{\mathrm{peak}}$
posteriors for which the ratio of signal samples to noise samples
$N_{s}:N_{n}$ acts as a loose proxy for signal SNR. A posterior with $N_{s}=0$
contains no hint of a signal whatsoever; conversely, one with $N_{n}=0$ would
be emblematic of a tight, normally-distributed posterior about the injected
$\mu_{f}$. One can then move along the spectrum created by these end-points to
simulate a variety of ambiguous $f_{\mathrm{peak}}$ posteriors of differing
strengths, representative of the $f_{\mathrm{peak}}$ posteriors produced by
BayesWave for different post-merger signal SNRs. However, we stress that
$N_{s}:N_{n}$ is only a proxy for SNR insofar as increasing $N_{s}:N_{n}$
correlates with increasing SNR; the two are not interchangeable, nor can one
directly convert between them.

We begin by considering the case of uninformative $f_{\mathrm{peak}}$ data,
with $N_{s}=0$. We apply the analysis to a set of 100 synthetic ensembles,
each containing 100 events, whose $f_{\mathrm{peak}}$ posteriors consist
solely of noise samples ($N_{s}=0,N_{n}=5000$). As no EoS information is
encoded in the $f_{\mathrm{peak}}$ data, we expect to recover results
consistent with our prior on $R_{1.6}$, with some fluctuations due to random
alignment across events of noise samples with the empirical relations. As is
shown in Fig. 9, this is what occurred. Moreover, this set of $R_{1.6}$
posteriors from uninformative data provides us with a reference dataset with
which to compare our recoveries of informative data.

Figure 9: 100 realizations of $R_{1.6}$ posteriors arising from 100 events
consisting solely of noise draws. These distributions are representative of
what can be produced from random fluctuations in the $f_{\mathrm{peak}}$
posteriors, absent any impact from an embedded signal. As expected, while we
do observe some level of fluctuation in these $R_{1.6}$ posteriors, the set as
a whole is consistent with the prior. The flattening at low $R_{1.6}$ arises
from our band-limited search; as discussed in §II.2, most predicted values of
$\hat{f}_{\mathrm{peak}}$ for $R_{1.6}\lesssim 10$ km fall out-of-band and do
not contribute information to our final posterior. As a result, most analyses
without a significant signal present in the data will tend to cleave to the
prior for these values of $R_{1.6}$.

We then turn to the most pertinent test case: an ensemble of weak signals.
With our same 100-event ensemble, we slowly increase the ratio of
$N_{s}:N_{n}$ and investigate how the analysis responds to increasing signal
strength in the weak-signal regime. In Fig. 10 we show sample simulated
$f_{\mathrm{peak}}$ posteriors and the corresponding final $R_{1.6}$ posterior
for the 100-event ensemble for each ratio of $N_{s}:N_{n}$.

Figure 10: Sample synthetic $f_{\mathrm{peak}}$ posteriors and corresponding
100-event ensemble $R_{1.6}$ posteriors for increasing ratio of signal samples
to noise samples (S:N). Note that in all cases, $N_{s}$ is sufficiently small
that the corresponding sample $f_{\mathrm{peak}}$ posterior is not strongly
peaked (and in fact by eye appears to be entirely uninformative). For
reference, we have also plotted the set of noise-only $R_{1.6}$ posteriors
described in §III.1. From top to bottom, these posteriors correspond to
$R_{1.6}$ means and 95% C.I. of $12.3^{+2.5}_{-3.2}$ km, $12.1^{+2.1}_{-2.9}$
km, $12.5^{+0.9}_{-2.6}$ km, and $12.59^{+0.33}_{-0.31}$ km, respectively.

At low $N_{s}:N_{n}$, we obtain an $R_{1.6}$ posterior consistent with random
fluctuations due to noise. As $N_{s}:N_{n}$ increases, we rapidly accrue
posterior density around the injected value of $R_{1.6}=12.5$ km. By
$N_{s}=25$, the posterior peak begins to emerge from the noise-only
distributions, by $N_{s}=50$, the peak is unmistakable, and by $N_{s}=100$, we
are able to exclude wide swathes of the parameter space, corresponding to a
mean and 95% C.I. of $R_{1.6}=12.59^{+0.33}_{-0.31}$ km.

We then explore the efficacy of our analysis across the $R_{1.6}$ parameter
space. We repeat the 100-event, ($N_{s}=50,N_{n}=5000$) test, changing the
injected value of $R_{1.6}$ for each iteration. As shown in Fig. 11, we are
able to recover the injected value in each case, but the strength and
precision of the recovery has some dependence on $R_{1.6}$. This behavior is
largely driven by two factors. At low $R_{1.6}$, $f_{\mathrm{peak}}$ tends to
be higher, and more frequently falls outside of our analysis frequency band.
As fewer signals are contained within our data, less overall information is
available and the resulting posterior is correspondingly weaker. At high
$R_{1.6}$, however, the increased posterior width is simply driven by a larger
intrinsic scatter in the empirical relation. As can be seen in Fig. 2, as we
approach higher values of $R_{1.6}$, the empirical relation curves in the
$f_{\mathrm{peak}}-\mathcal{M}$ plane become shallower and closer together for
equal intervals in $R_{1.6}$. Accordingly, the ability of a
($f_{\mathrm{peak}},\mathcal{M}$) pair to differentiate between nearby values
of $R_{1.6}$ is reduced.

Figure 11: 100-event ensemble $R_{1.6}$ posteriors with $N_{s}=50,N_{n}=5000$
for different injected values of $R_{1.6}$. For reference, we have also
plotted the set of noise-only $R_{1.6}$ posteriors shown in Fig. 9. We accrue
posterior density around the injected value in each case, although the
recovery strength is reduced by mitigating factors at low and high $R_{1.6}$.

Finally, we investigate the efficacy scaling of our analysis with the number
of events included. As in Fig. 10, we consider $N_{s}:N_{n}$ ratios of
$10:5000$, $25:5000$, $50:5000$, and $100:5000$ with $R_{1.6}=12.5$ km and
analyze a 2000-event ensemble generated as described above. Fig. 12 shows the
evolution with increasing $N_{\mathrm{events}}$ of the ensemble posterior 68%
C.I. width evolution as a fraction of the original prior C.I.. In each case,
we can see three stages of posterior evolution (although different stages
dominate for different signal strengths). For low $N_{\mathrm{events}}$, few
gains are made. However, as the ensemble becomes more robust and covers more
of the parameter space, the power of our analysis comes to the forefront, and
we enter a second phase where the posterior drastically and rapidly improves.
Ultimately, the diversity of events saturates, and most additional events do
not explore a new area of the parameter space; accordingly, the evolution
moves into its final stage, wherein its slope matches that of the “standard”
$1/\sqrt{N}$ evolution.

Figure 12: Scaling with $N_{\mathrm{events}}$ of the ensemble $R_{1.6}$
posterior 68% C.I. width as a fraction of the prior C.I. width for a
2000-event ensemble with $R_{1.6}=12.5$ km and four different injected signal
strengths. Note that the same behavior occurs for other C.I. percentages, but
viewing the entire trend at low signal strengths for e.g. 95% C.I. requires
many more events.

We emphasize that, for each case considered, the individual
$f_{\mathrm{peak}}$ posteriors remain ambiguous, but our analysis is able to
combine the trace amounts of information embedded therein to produce
informative constraints on $R_{1.6}$. This showcases the power of our
hierarchical approach under ideal, controlled conditions. We now turn to
investigate how our analysis performs on its intended input: the peak
frequency posteriors of BayesWave.

###  III.2 Simulations of Future Observing Runs

We perform simulations of the next two LVK observing runs (O4 and O5), using
the catalogues of Petrov  _et al._ [112]. Originally created to accurately
explore the prospects for EM followup of GW events in O4/O5, these catalogues
comprise realistic simulations of compact binary merger detections in O4 and
O5, using astrophysical distributions of binary parameters and SNR detection
thresholds reflecting those utilized in O3. For BNS mergers specifically, NS
component masses $m_{1}$ and $m_{2}$ were normally distributed as
$\mathcal{N}(1.33\mathrm{M_{\odot}},0.09\mathrm{M_{\odot}})$ and component
spins were uniformly distributed on [0,0.05], following Abbott  _et al._
[111]. Events were isotropically distributed in space, given redshifts
uniformly distributed in comoving rate density, and generated with the rate
measured in the latest O3 results for the BNS astrophysical rate density of
$320^{+490}_{-240}\;\textrm{Gpc}^{-3}\textrm{yr}^{-1}$ (see Petrov  _et al._
[112] for further details).

We use a one year duration for both O4 and O5. Petrov  _et al._ [112] predict
a median annual number of detected BNS mergers of 34 for O4 and 190 for O5
with an SNR cut of inspiral SNR$>$8, reflective of the alerts sent out during
O3; see Petrov  _et al._ [112] for justification. We therefore take 34 and 190
draws from the Petrov  _et al._ [112] distributions to simulate O4 and O5,
respectively, ensuring that both realizations are representative of the median
number of events in the nearby universe (which will be the most important for
post-merger searches). The Petrov  _et al._ [112] simulations contain many BNS
mergers fated to direct collapse, but we are only interested in the HMNS post-
merger scenario. Accordingly, we discard all events with
$M_{\mathrm{tot}}>3\mathrm{M_{\odot}}$. Such a hard cut is somewhat inelegant;
for example, $M_{\mathrm{thres}}$ for the DD2 EoS lies above this cut by $\sim
0.3\mathrm{M_{\odot}}$ whereas that of SLy4 lies below it by $\sim
0.15\mathrm{M_{\odot}}$, leading us to discard events we should keep in the
former case and vice versa for the latter. Ultimately, this approach is a
necessary evil; it allows for identical treatment and therefore direct
comparison of each EoS with respect towards the others, and as such suits our
present purposes. A more refined treatment of $M_{\mathrm{thres}}$ and the
direct collapse scenario is left to future work (see § IV.3). After this cut,
we are left with 26 (83) events in O4 (O5). For each selected event in the
catalogue, we analyze the post-merger waveform from our library
closest151515Specifically, we select the waveform with the smallest 2-norm
distance in $m_{1}-m_{2}$ space from the original catalogue event’s
parameters. to the catalogue parameters, injected and recovered using the
event’s distance as given in Petrov  _et al._ [112]. The final set of events
with associated $f_{\mathrm{peak}}$ posteriors can be found on Zenodo at
10.5281/zenodo.7007630. We use the detector strain noise spectral density
released in LIGO-T2000012 [61]. For the purposes of this study, we assume each
event was optimally-inclined and detected by LIGO Hanford, LIGO Livingston,
and VIRGO interferometers. For each post-merger signal, we simulate a
corresponding chirp mass posterior as discussed in §II.4. We perform our
simulations thrice, once for each EoS (SLy4, SFHX, and DD2). In Fig. 13, we
present the full O4+O5 simulation posterior for each EoS with a uniform prior
in $R_{1.6}$. For SLy4, SFHX, and DD2, the full O4+O5 posteriors correspond to
means and 95% C.I. of $R_{1.6}=11.8^{+3.0}_{-2.6}$ km,
$R_{1.6}=12.0^{+2.8}_{-2.8}$ km, and $R_{1.6}=11.8^{+3.0}_{-2.7}$ km,
respectively, compared to the original uniform prior mean and 95% C.I. of
$R_{1.6}=12.0^{+2.9}_{-2.9}$ km. Additionally, we present the O4, O5, and
combined O4+O5 posteriors for the SFHX EoS in Fig. 14 with the multimessenger
prior of Huth  _et al._ [35]. The full O4+O5 ensemble posterior equates to a
mean and 95% C.I. of $R_{1.6}=12.02^{+0.93}_{-0.68}$ km, a reduction of 0.13
km (7.5%) compared to the original prior 95% C.I. of
$R_{1.6}=12.07^{+0.98}_{-0.77}$ km. As would be expected, the final posterior
is primarily driven by O5 events, reflecting both the increased detector
sensitivity and event count. However, the combined O4+O5 posterior does
improve slightly when compared with O5 alone, indicating that the O4 data is
in fact providing some information.

Figure 13: $R_{1.6}$ posteriors for simulations of LVK’s upcoming O4 and O5
observing runs for the SLy4, SFHX, and DD2 equations of state. These
posteriors correspond to means and 95% C.I. of $R_{1.6}=11.8^{+3.0}_{-2.6}$
km, $R_{1.6}=12.0^{+2.8}_{-2.8}$ km, and $R_{1.6}=11.8^{+3.0}_{-2.7}$ km,
respectively, as compared to the original uniform prior mean and 95% C.I. of
$R_{1.6}=12.0^{+2.9}_{-2.9}$ km. While there appear to be hints of a signal in
the SLy4 and SFHX cases, the final posterior is not conclusive and could arise
from noise fluctuations, as seen in the case of DD2. Figure 14: $R_{1.6}$
posteriors for simulations of LVK’s upcoming O4 and O5 observing runs with the
SFHX EoS and the multimessenger prior of Huth  _et al._ [35]. The individual
O4 and O5 posteriors are presented along with the combined O4+O5 posterior.
For O4, O5, and O4+O5, these posteriors correspond to means and 95% C.I. of
$R_{1.6}=12.07^{+0.95}_{-0.75}$ km, $R_{1.6}=12.02^{+0.95}_{-0.70}$ km, and
$R_{1.6}=12.02^{+0.93}_{-0.68}$ km, respectively. For comparison, the prior
mean and 95% C.I. are $R_{1.6}=12.07^{+0.98}_{-0.77}$ km. The combined O4+O5
posterior provides a reduction of the 95% C.I. width by $0.13$ km (7.5%) over
the prior. Figure 15: $R_{1.6}$ posteriors arising from a simulation of 4
years of operation at A+ sensitivity for the SFHX EoS, presented with both a
uniform prior and the multimessenger prior of Huth  _et al._ [35]. For the
multimessenger (uniform) prior, we attain posterior means and 95% C.I. of
$R_{1.6}=11.91^{+0.80}_{-0.56}$ km ($R_{1.6}=11.9^{+2.7}_{-2.6}$ km). For
reference, the multimessenger (uniform) prior means and 95% C.I. are
$R_{1.6}=12.07^{+0.98}_{-0.77}$ km ($R_{1.6}=12.0^{+2.9}_{-2.9}$ km). In the
case of the multimessenger prior, we achieve a reduction in the 95% C.I. width
of $0.39$ km (22%).

While these posteriors do not comprise stringent constraints on $R_{1.6}$,
there is an accumulation of posterior probability about the simulated values
of $R_{1.6}$ corresponding to the SFHX and SLy4 EoSs. However, the result for
the DD2 EoS is less informative and does not show accumulation about the true
value as seen with the other two equations of state. That we see such behavior
in the case of DD2 is not entirely unsurprising; as demonstrated in Fig. 11,
our analysis is most sensitive for middling values of $R_{1.6}$, and decreases
in statistical power towards the high-$R_{1.6}$ regime occupied by the DD2
EoS. The nature of the O4+O5 posteriors is consistent with the behavior of our
analysis as seen in §III.1 for the case of low signal strength and
(relatively) small $N_{\mathrm{events}}$; it is possible that the embedded
information contained in the O4+O5 simulation data is sufficiently weak that
this high-$R_{1.6}$ reduction in sensitivity proved insurmountable. Repetition
of the DD2 O4+O5 study with all events injected at $20$ Mpc — corresponding to
a mix of low- and mid-SNR signals — yielded posterior accumulation about the
injected value of $R_{1.6}=13.26$ for DD2, lending additional credence to this
explanation. We further investigate methods for characterizing the probability
of the observed posterior accumulation arising from noise alone in §III.3.

As seen in §III.1, the primary source of our resolving power in the low-SNR
regime is the inclusion of additional events. If the posterior accumulation
observed in e.g. the SFHX O4+O5 simulation is reflective of real information
contained within the post-merger data, we would expect to see continued gains
as we increase $N_{\mathrm{events}}$. Accordingly, we consider what may be
possible with 4 years of data from the LIGO/Virgo network operating at O5 (A+)
sensitivity. We follow the procedure described for the simulations of O4 and
O5, save that we draw 760 events (4 years at a median annual detection count
of 190 events), of which 357 are retained after our cut in $M_{\mathrm{tot}}$.
We consider only the SFHX EoS, as the finer $q-M_{\mathrm{tot}}$ grid used to
compute the SFHX waveforms lends itself to deep studies with larger
$N_{\mathrm{events}}$. The resulting combined posterior can be found in Fig.
15, and is presented both with a uniform prior and with the astrophysical
prior of Huth  _et al._ [35] as a demonstration of what this method may be
able to contribute to our knowledge of the NS EoS in concert with current
constraints. For the multimessenger prior, the final $R_{1.6}$ posterior
equates to a mean and 95% C.I. of $R_{1.6}=11.91^{+0.80}_{-0.56}$ km, a
reduction in the 95% C.I. width of $0.39$ km (22%) over the prior 95% C.I. of
$R_{1.6}=12.07^{+0.98}_{-0.77}$ km. We emphasize that this posterior arises
solely from post-merger signals whose individual $f_{\mathrm{peak}}$
posteriors are consistent with nondetection; this result demonstrates the
potential of our approach to produce informative constraints on the NS EoS
from post-merger GWs with current-generation detector sensitivities.

###  III.3 Validation

As the post-merger signals considered in this work are individually consistent
with nondetection, it is important that we be able to assess if the final
posterior for $R_{1.6}$ has arisen from running our analysis over detector
noise, as opposed to containing useful information about the NS EoS. To this
end, we repeat the simulations of O4, O5, and the A+ 4 year dataset 1000 times
each, following the procedure of §III.2, save that the $f_{\mathrm{peak}}$
posteriors for each event do not contain any information about the post-merger
signal. Instead, each $f_{\mathrm{peak}}$ posterior consists only of an
equivalent number of prior draws (which in turn are draws from noise-only
BayesWave posteriors; see §II.5). The result of this process is a set of
“null” posterior densities for $R_{1.6}$ for each of our observing run
simulations.

This null set of posterior densities is then used to investigate whether we
can differentiate our simulation results from the noise-only case. Standard
model comparison techniques, such as Bayes factors, are not well-suited to our
setting; we therefore present the following method, based on the usual
statistical hypothesis testing framework, to characterize the difference
between the set of null densities and the final $R_{1.6}$ posterior
distributions we obtain from the original simulations. For each null $R_{1.6}$
posterior density, we compute two metrics: the Wasserstein distance161616Also
called the earth-mover’s distance. from the prior to the posterior and the
standard 95% C.I. posterior width. The Wasserstein distance between two
probability density functions $f$ and $g$ is defined as

$W(f,g)=\int_{\mathbb{R}}|F(x)-G(x)|\,dx,$

where $F$ and $G$ are the cumulative distribution functions associated with
$f$ and $g$, respectively. The distribution of values for a given metric
associated with the collection of null densities forms a simulated null
distribution; we can then use this distribution to estimate the p-value of the
same metric as computed for the final $R_{1.6}$ posterior for each of our
simulated observing runs. For the Wasserstein distance, we estimate the
p-value as the proportion of the null distribution greater than the observed
density distance, whereas the equivalent estimate for the 95% C.I. width is
the proportion of null C.I. widths less than the C.I. width of the observed
density. The metric values and their associated p-values for each EoS and
simulation are shown in Tables 3 and 4 for analyses with the uniform and
multimessenger priors, respectively.

Run | EoS | WD | p (WD) | CI95 | p (CI95)
---|---|---|---|---|---
O4+O5 | SLy4 | 0.21 | 0.48 | 5.64 km | 0.17
O4+O5 | SFHX | 0.10 | 0.88 | 5.70 km | 0.42
O4+O5 | DD2 | 0.24 | 0.38 | 5.64 km | 0.17
A+ 4yr | SFHX | 0.62 | 0.23 | 5.37 km | 0.14
Table 3: Estimated p-values for the Wasserstein distance (WD) and 95% credible interval width (CI95) metrics for simulation results with a uniform $R_{1.6}$ prior. Run | EoS | WD | p (WD) | CI95 | p (CI95)
---|---|---|---|---|---
O4+O5 | SLy4 | 0.02 | 0.97 | 1.78 km | 0.67
O4+O5 | SFHX | 0.05 | 0.68 | 1.60 km | 0.20
O4+O5 | DD2 | 0.07 | 0.49 | 1.63 km | 0.24
A+ 4yr | SFHX | 0.16 | 0.40 | 1.36 km | 0.12
Table 4: Estimated p-values for the Wasserstein distance (WD) and 95% credible
interval width (CI95) metrics for simulation results with the Huth  _et al._
[35] multimessenger $R_{1.6}$ prior.

Despite some posteriors hinting at the presence of a signal by eye, the
corresponding p-values for our simulated results are insufficient to claim a
decisive constraint. This is reflected in the fact that none of our posteriors
confidently exclude large swathes of the $R_{1.6}$ parameter space, instead
displaying in most cases posterior accumulation near the injected value of
$R_{1.6}$. We note that the A+ 4yr dataset has the highest significance for
all configurations. The case of the A+ 4yr dataset with the multimessenger
prior has the highest overall significance at $p=0.12$ for the 95% C.I.
metric. Both the corresponding analysis with a uniform prior and that of the
multimessenger prior and Wasserstein distance metric yield lower significance
estimates. The increased significance of the multimessenger result is not
entirely surprising; as discussed in §III.1, our analysis is less effective at
constraining the far tails of the $R_{1.6}$ distribution, and gains power
accordingly when combined with other data that rules out those tails. The
difference in significance for the A+ 4yr result between the 95% C.I. and
Wasserstein distance metrics may indicate a need for further exploration; in
general, this kind of approach to model comparison is an interesting
statistical problem and warrants further work beyond the scope of this paper.

##  IV Discussion and Conclusions

We present a hierarchical Bayesian approach for constraining the NS EoS
through BNS post-merger GW emission. The analysis is able to combine an
ensemble of barely-informative, marginal detections of the BNS post-merger
signal in a coherent fashion to achieve informative constraints on the NS EoS.
Additionally, we provide results of the analysis as applied to simulations of
the upcoming LVK observing runs O4 and O5 for three different NS equations of
state, alongside an investigation of what could be achieved with an extended
run of 4 years at A+ sensitivity. The results of these simulations indicate
that this method has the potential to produce informative posteriors on
$R_{1.6}$ from BNS post-merger GW with current generation detectors. If
reflected in the GW observations of coming years, the resulting constraints
would comprise a novel contribution to our understanding of the NS EoS from a
previously inaccessible source.

###  IV.1 Prospects for Future Detection

The results we achieve in simulations at current-generation sensitivity
indicate a promising future as detector sensitivity improves. In the near
future, the proposed LIGO Voyager upgrades would greatly improve detection and
characterization of BNS inspiral signals [125]. In addition to Voyager’s own
increased high-frequency sensitivity, a promising project by OzGrav, NEMO
[126], is explicitly designed for sensitivity in the high-frequency regime. A
2.5-generation detector network combining BNS inspiral detections by Voyager
and post-merger data from NEMO could enable this approach to yield greatly
improved constraints on the NS EoS.

The advent of 3rd-generation detectors such as Einstein Telescope [127] and
Cosmic Explorer [128] will bring BNS post-merger GW well within the detection
horizon. As we enter the fully-operational era of 3rd-generation detectors, it
is expected that we will observe a handful of unambiguous post-merger signals
each year. However, it is likely that in the early days of these detectors, a
large fraction of any post-merger detections will be marginal. Moreover, we
expect to see many more marginal signals of all types than so-called “gold-
plated” detections. This will be no less true of the BNS post-merger signal
[96], providing utility for our approach well into future generations of GW
detectors.

###  IV.2 Potential to Probe Exotic Equations of State

It is worth noting that the large number of BNS detections required to achieve
competitive constraints on $R_{1.6}$ using this method will also yield greatly
improved constraints on the NS EoS through GW inspiral measurements of
$\Lambda(m)$ and EM observations of kilonovae [26, 129, 130, 131]. However,
even if the constraints placed by this method are outpaced by the
aforementioned inspiral- and EM-counterpart-derived constraints, there is
still significant utility in considering NS EoS constraints specifically
arising from post-merger GW. The reasons for this are twofold. First, any
additional information is useful. Even mild constraints produced through our
method will improve the potentially more accurate inspiral measurements. More
importantly, though, there is a possibility that some post-merger remnants may
undergo a phase transition occurring during the merger. If this is the case,
our post-merger analysis will yield EoS constraints that differ from inspiral
constraints arising from the same BNS population [67, 98]. If the former and
latter distributions were to exclude each other to a statistically significant
extent, our approach could provide strong evidence for exotic matter in BNS
post-merger remnants. This possibility renders crucial the concurrent,
comparative consideration of both inspiral and post-merger constraints.

###  IV.3 Future Work

Looking forward, this work has several natural extensions. First, Vretinaris
_et al._ [93] provide empirical relations for several values of
$\text{R}_{\text{1.x}}$ beyond $R_{1.6}$. Extending this analysis to provide
constraints throughout the mass-radius relation space would greatly improve
its ability to differentiate between proposed equations of state.

Additionally, given this work’s results for current-generation detectors, it
is certainly worth investigating what can be accomplished with our analysis
when applied to 2.5- and 3rd-generation detectors. Realistic simulations of a
2.5-generation network consisting of Voyager and NEMO may provide a new
science case for the implementation of such a network, and a similar treatment
of Cosmic Explorer/Einstein Telescope would provide important insight into
what can be accomplished in the future — and how constraints derived from
marginal detections can compare to those from unambiguous ones. With the
extended reach of these detectors, it will be necessary to directly
incorporate inspiral redshift posteriors into our statistical hierarchy.
Additionally, a more refined treatment of remnant classification would allow
us to reduce the potential impact of contamination from direct collapse
scenarios as well as avoid discarding true HMNS remnants with masses near the
(unknown) mass threshold for direct collapse, $M_{\text{thres}}$. This could
be accomplished through determination of the probability of direct collapse —
marginalized over inspiral $M_{\text{tot}}$ posteriors and current EoS
uncertainties, as performed in Abbott  _et al._ [13] — in conjunction with the
addition of a remnant-classification mixture model to our formalism and/or
simultaneous inference of $R_{1.6}$ and $M_{\text{thres}}$ using (e.g.) the
relation between $M_{\text{thres}}$ and $f_{\mathrm{peak}}^{max}$ [132].

Finally, as for instance discussed in Sarin and Lasky [133], Bauswein  _et
al._ [67], Chatziioannou and Han [134] and mentioned in §IV.2, there is a
possibility that during the merger or post-merger phase, the BNS remnant may
experience a hadron-quark phase transition. As shown in Bauswein  _et al._
[67] and Blacker  _et al._ [135], this phase transition can manifest as a
mass-dependent offset in $f_{\mathrm{peak}}$ from what is predicted by the
empirical relations discussed in this work, dependent on the strength of the
phase transition. This phenomenon could enable us to extend our hierarchical
formalism to additionally infer the presence/absence of a strong phase
transition, and — if present — the transition strength and nature of the high-
mass, post-transition regime of the EoS. Such an extension would lend
significant utility to our analysis well into the era of next-generation
detectors, even if unambiguous detections of post-merger signals become
common.

## Acknowledgements

The authors would like to thank James Clark and Tim Dietrich for many helpful
conversations about BayesWave and multimessenger considerations, respectively.
This research makes use of the Numpy [136], SciPy [124], Pandas [137],
Matplotlib [138] scientific computing packages, the BayesWave data analysis
pipeline [102, 139, 104], and utilized computing resources provided by the
Minnesota Supercomputing Institute (MSI) at the University of Minnesota. The
authors are grateful for computational resources provided by the LIGO
Laboratory and supported by National Science Foundation Grants PHY-0757058 and
PHY-0823459. This material is based upon work supported by NSF’s LIGO
Laboratory which is a major facility fully funded by the National Science
Foundation. A.W.C. and J.M. acknowledge support by the National Science
Foundation Research Traineeship Program ‘Data Science in Multimessenger
Astrophysics’ under grant No. 1922512. M.W.C. acknowledges support from the
National Science Foundation with grant numbers PHY-2010970 and OAC-2117997.
A.B. acknowledges support by the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme under grant
agreement No. 759253, and support by Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) - Project-ID 279384907 - SFB 1245 and DFG -
Project-ID 138713538 - SFB 881 (“The Milky Way System”, subproject A10) and
support by the State of Hesse within the Cluster Project ELEMENTS. The work of
T.S. is supported by the State of Hesse within the Cluster Project ELEMENTS,
and the Klaus Tschira Foundation. T.S. is Fellow of the International Max
Planck Research School for Astronomy and Cosmic Physics at the University of
Heidelberg (IMPRS-HD) and acknowledges financial support from IMPRS-HD. T.S.
acknowledges support by the High Performance and Cloud Computing Group at the
Zentrum für Datenverarbeitung of the University of Tübingen, the state of
Baden-Württemberg through bwHPC and the German Research Foundation (DFG)
through grant no INST 37/935-1 FUGG. A.W.C. and J.M. contributed equally to
this work.

## References

  * Baym  _et al._ [2018] G. Baym, T. Hatsuda, T. Kojo, _et al._ , Reports on Progress in Physics 81, 056902 (2018).
  * Riley  _et al._ [2019] T. E. Riley, A. L. Watts, S. Bogdanov, _et al._ , The Astrophysical Journal 887, L21 (2019), arXiv:1912.05702 .
  * Miller  _et al._ [2019] M. C. Miller, F. K. Lamb, A. J. Dittmann, _et al._ , The Astrophysical Journal 887, L24 (2019), arXiv:1912.05705 .
  * Riley  _et al._ [2021] T. E. Riley, A. L. Watts, P. S. Ray, _et al._ , arXiv:2105.06980 [astro-ph, physics:nucl-th] (2021), arXiv:2105.06980 [astro-ph, physics:nucl-th] .
  * Miller  _et al._ [2021] M. C. Miller, F. K. Lamb, A. J. Dittmann, _et al._ , The Astrophysical Journal Letters 918, L28 (2021), arXiv:2105.06979 .
  * Demorest  _et al._ [2010] P. B. Demorest, T. Pennucci, S. M. Ransom, _et al._ , Nature 467, 1081 (2010).
  * Antoniadis  _et al._ [2013] J. Antoniadis, P. C. C. Freire, N. Wex, _et al._ , Science 340, 448 (2013).
  * Cromartie  _et al._ [2020] H. T. Cromartie, E. Fonseca, S. M. Ransom, _et al._ , Nature Astronomy 4, 72 (2020), arXiv:1904.06759 .
  * Fonseca  _et al._ [2021] E. Fonseca, H. T. Cromartie, T. T. Pennucci, _et al._ , The Astrophysical Journal Letters 915, L12 (2021), arXiv:2104.00880 .
  * Romani  _et al._ [2022] R. W. Romani, D. Kandel, A. V. Filippenko, _et al._ , The Astrophysical Journal Letters 934, L17 (2022).
  * Abbott  _et al._ [2017a] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ , Physical Review Letters 119, 161101 (2017a).
  * Abbott  _et al._ [2018] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ , Physical Review Letters 121, 161101 (2018).
  * Abbott  _et al._ [2020a] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ , The Astrophysical Journal 892, L3 (2020a).
  * Abbott  _et al._ [2017b] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ , The Astrophysical Journal 848, L12 (2017b).
  * Soares-Santos  _et al._ [2017] M. Soares-Santos, D. E. Holz, J. Annis, _et al._ , The Astrophysical Journal 848, L16 (2017).
  * Cowperthwaite  _et al._ [2017] P. S. Cowperthwaite, E. Berger, V. A. Villar, _et al._ , The Astrophysical Journal 848, L17 (2017).
  * Nicholl  _et al._ [2017] M. Nicholl, E. Berger, D. Kasen, _et al._ , The Astrophysical Journal 848, L18 (2017).
  * Chornock  _et al._ [2017] R. Chornock, E. Berger, D. Kasen, _et al._ , The Astrophysical Journal 848, L19 (2017).
  * Margutti  _et al._ [2017] R. Margutti, E. Berger, W. Fong, _et al._ , The Astrophysical Journal 848, L20 (2017).
  * Alexander  _et al._ [2017] K. D. Alexander, E. Berger, W. Fong, _et al._ , The Astrophysical Journal 848, L21 (2017).
  * Goldstein  _et al._ [2017] A. Goldstein, P. Veres, E. Burns, _et al._ , The Astrophysical Journal 848, L14 (2017).
  * Savchenko  _et al._ [2017] V. Savchenko, C. Ferrigno, E. Kuulkers, _et al._ , The Astrophysical Journal 848, L15 (2017).
  * Coughlin  _et al._ [2019] M. W. Coughlin, T. Dietrich, B. Margalit, and B. D. Metzger, Monthly Notices of the Royal Astronomical Society: Letters 489, L91 (2019), arXiv:1812.04803 .
  * Raaijmakers  _et al._ [2020] G. Raaijmakers, S. K. Greif, T. E. Riley, _et al._ , The Astrophysical Journal 893, L21 (2020), arXiv:1912.11031 .
  * Dietrich  _et al._ [2020] T. Dietrich, M. W. Coughlin, P. T. H. Pang, _et al._ , Science 370, 1450 (2020).
  * Landry  _et al._ [2020] P. Landry, R. Essick, and K. Chatziioannou, Physical Review D 101, 123007 (2020), arXiv:2003.04880 .
  * Capano  _et al._ [2020] C. D. Capano, I. Tews, S. M. Brown, _et al._ , Nature Astronomy 4, 625 (2020).
  * Raaijmakers  _et al._ [2021] G. Raaijmakers, S. K. Greif, K. Hebeler, _et al._ , arXiv:2105.06981 [astro-ph, physics:nucl-ex, physics:nucl-th] (2021), arXiv:2105.06981 [astro-ph, physics:nucl-ex, physics:nucl-th] .
  * Al-Mamun  _et al._ [2021] M. Al-Mamun, A. W. Steiner, J. Nättilä, _et al._ , Physical Review Letters 126, 061101 (2021), arXiv:2008.12817 .
  * Biswas [2021] B. Biswas, The Astrophysical Journal 921, 63 (2021), arXiv:2105.02886 .
  * Pang  _et al._ [2021] P. T. H. Pang, I. Tews, M. W. Coughlin, _et al._ , The Astrophysical Journal 922, 14 (2021), arXiv:2105.08688 .
  * Breschi  _et al._ [2021] M. Breschi, A. Perego, S. Bernuzzi, _et al._ , Monthly Notices of the Royal Astronomical Society 505, 1661 (2021).
  * Legred  _et al._ [2021] I. Legred, K. Chatziioannou, R. Essick, _et al._ , Physical Review D 104, 063003 (2021).
  * Nicholl  _et al._ [2021] M. Nicholl, B. Margalit, P. Schmidt, _et al._ , Monthly Notices of the Royal Astronomical Society 505, 3016 (2021).
  * Huth  _et al._ [2022] S. Huth, P. T. H. Pang, I. Tews, _et al._ , Nature 606, 276 (2022), arXiv:2107.06229 [astro-ph, physics:gr-qc, physics:nucl-ex, physics:nucl-th] .
  * Aasi  _et al._ [2015] a. J. Aasi, B. P. Abbott, R. Abbott, _et al._ , Classical and Quantum Gravity 32, 074001 (2015).
  * Acernese  _et al._ [2015] F. Acernese, M. Agathos, K. Agatsuma, _et al._ , Classical and Quantum Gravity 32, 024001 (2015), arXiv:1408.3978 .
  * Somiya [2012] K. Somiya, Classical and Quantum Gravity 29, 124007 (2012).
  * Aso  _et al._ [2013] Y. Aso, Y. Michimura, K. Somiya, _et al._ , Physical Review D 88, 043007 (2013).
  * Baumgarte  _et al._ [1999] T. W. Baumgarte, S. L. Shapiro, and M. Shibata, The Astrophysical Journal 528, L29 (1999).
  * Hotokezaka  _et al._ [2013] K. Hotokezaka, K. Kiuchi, K. Kyutoku, _et al._ , Physical Review D 88, 044026 (2013), arXiv:1307.5888 .
  * Margalit and Metzger [2017] B. Margalit and B. Metzger, The Astrophysical Journal 850, L19 (2017), arXiv:1710.05938 [astro-ph] .
  * Granot  _et al._ [2017] J. Granot, D. Guetta, and R. Gill, The Astrophysical Journal 850, L24 (2017), arXiv:1710.06407 [astro-ph] .
  * Perego  _et al._ [2017] A. Perego, D. Radice, and S. Bernuzzi, The Astrophysical Journal 850, L37 (2017).
  * Ma  _et al._ [2018] P.-X. Ma, J.-L. Jiang, H. Wang, _et al._ , The Astrophysical Journal 858, 74 (2018), arXiv:1711.05565 [astro-ph, physics:gr-qc] .
  * Ruiz  _et al._ [2018] M. Ruiz, S. L. Shapiro, and A. Tsokaros, Physical Review D 97, 021501 (2018).
  * Pooley  _et al._ [2018] D. Pooley, P. Kumar, J. C. Wheeler, and B. Grossan, The Astrophysical Journal 859, L23 (2018), arXiv:1712.03240 [astro-ph] .
  * Rezzolla  _et al._ [2018] L. Rezzolla, E. R. Most, and L. R. Weih, The Astrophysical Journal 852, L25 (2018), arXiv:1711.00314 [astro-ph, physics:gr-qc] .
  * Metzger  _et al._ [2018] B. D. Metzger, T. A. Thompson, and E. Quataert, The Astrophysical Journal 856, 101 (2018), arXiv:1801.04286 [astro-ph] .
  * Gill  _et al._ [2019] R. Gill, A. Nathanail, and L. Rezzolla, The Astrophysical Journal 876, 139 (2019).
  * Shibata and Uryu [2000] M. Shibata and K. Uryu, Physical Review D 61, 064001 (2000), arXiv:gr-qc/9911058 .
  * Bauswein  _et al._ [2013] A. Bauswein, T. W. Baumgarte, and H.-T. Janka, Physical Review Letters 111, 131101 (2013), arXiv:1307.5191 .
  * Shapiro [2000] S. L. Shapiro, The Astrophysical Journal 544, 397 (2000).
  * Stergioulas  _et al._ [2011] N. Stergioulas, A. Bauswein, K. Zagkouris, and H.-T. Janka, Monthly Notices of the Royal Astronomical Society 418, 427 (2011).
  * Stergioulas [2003] N. Stergioulas, Living Reviews in Relativity 6, 3 (2003).
  * Zhuge  _et al._ [1994] X. Zhuge, J. M. Centrella, and S. L. W. McMillan, Physical Review D 50, 6247 (1994).
  * Shibata  _et al._ [2005] M. Shibata, K. Taniguchi, and K. Uryū, Physical Review D 71, 084021 (2005).
  * Takami  _et al._ [2015] K. Takami, L. Rezzolla, and L. Baiotti, Physical Review D 91, 064001 (2015), arXiv:1412.3240 .
  * Bauswein and Stergioulas [2015] A. Bauswein and N. Stergioulas, Physical Review D 91, 124056 (2015).
  * Soultanis  _et al._ [2022] T. Soultanis, A. Bauswein, and N. Stergioulas, Physical Review D 105, 043020 (2022), arXiv:2111.08353 [astro-ph, physics:gr-qc] .
  * OReilly  _et al._ [2020] B. OReilly, M. Branchesi, S. Haino, and G. Gamme, LIGO-T2000012-v1: Noise curves used for Simulations in the update of the Observing Scenarios Paper, https://dcc.ligo.org/LIGO-T2000012/public (2020).
  * Evans  _et al._ [2020] M. Evans, R. Sturani, S. Vitale, and E. Hall, LIGO-T1500293-v13: Unofficial sensitivity curves (ASD) for aLIGO, Kagra, Virgo, Voyager, Cosmic Explorer, and Einstein Telescope, https://dcc.ligo.org/LIGO-T1500293/public (2020).
  * Adya  _et al._ [2020] V. Adya, D. Brown, J. Eichholz, _et al._ , LIGO-T2000062-v1: NEMO noise budget, https://dcc.ligo.org/LIGO-T2000062 (2020).
  * Bernuzzi  _et al._ [2015] S. Bernuzzi, T. Dietrich, and A. Nagar, Physical Review Letters 115, 091101 (2015), arXiv:1504.01764 .
  * Bauswein  _et al._ [2012] A. Bauswein, H.-T. Janka, K. Hebeler, and A. Schwenk, Physical Review D 86, 063001 (2012).
  * Bauswein and Janka [2012] A. Bauswein and H.-T. Janka, Physical Review Letters 108, 011101 (2012), arXiv:1106.1616 .
  * Bauswein  _et al._ [2019] A. Bauswein, N.-U. F. Bastian, D. B. Blaschke, _et al._ , Physical Review Letters 122, 061102 (2019), arXiv:1809.01116 .
  * Clark  _et al._ [2014] J. Clark, A. Bauswein, L. Cadonati, _et al._ , Physical Review D 90, 062004 (2014), arXiv:1406.5444 .
  * Torres-Rivas  _et al._ [2019] A. Torres-Rivas, K. Chatziioannou, A. Bauswein, and J. A. Clark, Physical Review D 99, 044014 (2019), arXiv:1811.08931 .
  * Tsang  _et al._ [2019] K. W. Tsang, T. Dietrich, and C. V. D. Broeck, Physical Review D 100, 044047 (2019), arXiv:1907.02424 .
  * Abbott  _et al._ [2017c] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ , The Astrophysical Journal 851, L16 (2017c).
  * Abbott  _et al._ [2019a] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ , The Astrophysical Journal 875, 160 (2019a).
  * Klimenko  _et al._ [2016] S. Klimenko, G. Vedovato, M. Drago, _et al._ , Physical Review D 93, 042004 (2016).
  * Chatziioannou  _et al._ [2017] K. Chatziioannou, J. A. Clark, A. Bauswein, _et al._ , Physical Review D 96, 124035 (2017), arXiv:1711.00040 .
  * Abbott  _et al._ [2019b] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ , Physical Review X 9, 011001 (2019b).
  * Agathos  _et al._ [2020] M. Agathos, F. Zappa, S. Bernuzzi, _et al._ , Physical Review D 101, 044006 (2020).
  * Breschi  _et al._ [2019] M. Breschi, S. Bernuzzi, F. Zappa, _et al._ , Physical Review D 100, 104029 (2019), arXiv:1908.11418 .
  * Breschi  _et al._ [2022a] M. Breschi, R. Gamba, S. Borhanian, _et al._ , Kilohertz Gravitational Waves from Binary Neutron Star Mergers: Inference of Postmerger Signals with the Einstein Telescope (2022a), arXiv:2205.09979 [astro-ph, physics:gr-qc] .
  * Easter  _et al._ [2020] P. J. Easter, S. Ghonge, P. D. Lasky, _et al._ , arXiv:2006.04396 [astro-ph, physics:gr-qc] (2020), arXiv:2006.04396 [astro-ph, physics:gr-qc] .
  * Martynov  _et al._ [2019] D. Martynov, H. Miao, H. Yang, _et al._ , Physical Review D 99, 102004 (2019), arXiv:1901.03885 .
  * Breschi  _et al._ [2022b] M. Breschi, S. Bernuzzi, K. Chakravarti, _et al._ , Kilohertz Gravitational Waves From Binary Neutron Star Mergers: Numerical-relativity Informed Postmerger Model (2022b), arXiv:2205.09112 [astro-ph, physics:gr-qc] .
  * Evans  _et al._ [2021] M. Evans, R. X. Adhikari, C. Afle, _et al._ , A Horizon Study for Cosmic Explorer: Science, Observatories, and Community (2021), arXiv:2109.09882 [astro-ph, physics:gr-qc] .
  * Srivastava  _et al._ [2022] V. Srivastava, D. Davis, K. Kuns, _et al._ , The Astrophysical Journal 931, 22 (2022), arXiv:2201.10668 [astro-ph, physics:gr-qc, physics:physics] .
  * Bose  _et al._ [2018] S. Bose, K. Chakravarti, L. Rezzolla, _et al._ , Physical Review Letters 120, 031102 (2018), arXiv:1705.10850 .
  * Yang  _et al._ [2018] H. Yang, V. Paschalidis, K. Yagi, _et al._ , Physical Review D 97, 024049 (2018), arXiv:1707.00207 .
  * Loredo [2004] T. J. Loredo, AIP Conference Proceedings 735, 195 (2004).
  * Wolfgang  _et al._ [2016] A. Wolfgang, L. A. Rogers, and E. B. Ford, The Astrophysical Journal 825, 19 (2016), arXiv:1504.07557 .
  * Martinez [2015] G. D. Martinez, Monthly Notices of the Royal Astronomical Society 451, 2524 (2015), arXiv:1309.2641 .
  * Hernandez Vivanco  _et al._ [2019] F. Hernandez Vivanco, R. Smith, E. Thrane, _et al._ , Physical Review D 100, 103009 (2019).
  * Mandel  _et al._ [2019] I. Mandel, W. M. Farr, and J. R. Gair, Monthly Notices of the Royal Astronomical Society 486, 1086 (2019), arXiv:1809.02063 [astro-ph, physics:physics] .
  * Abbott  _et al._ [2021a] R. Abbott, H. Abe, F. Acernese, _et al._ , arXiv:2111.03604 [astro-ph, physics:gr-qc] (2021a), arXiv:2111.03604 [astro-ph, physics:gr-qc] .
  * Abbott  _et al._ [2022] R. Abbott, T. D. Abbott, F. Acernese, _et al._ , arXiv:2111.03634 [astro-ph, physics:gr-qc] (2022), arXiv:2111.03634 [astro-ph, physics:gr-qc] .
  * Vretinaris  _et al._ [2020] S. Vretinaris, N. Stergioulas, and A. Bauswein, Physical Review D 101, 084039 (2020), arXiv:1910.10856 .
  * Raithel and Most [2022] C. A. Raithel and E. R. Most, Characterizing the breakdown of quasi-universality in the post-merger gravitational waves from binary neutron star mergers (2022), arXiv:2201.03594 [astro-ph, physics:gr-qc] .
  * Lioutas  _et al._ [2021] G. Lioutas, A. Bauswein, and N. Stergioulas, Physical Review D 104, 043011 (2021), arXiv:2102.12455 [astro-ph, physics:gr-qc] .
  * Haster  _et al._ [2020] C.-J. Haster, K. Chatziioannou, A. Bauswein, and J. A. Clark, Physical Review Letters 125, 261101 (2020), arXiv:2004.11334 .
  * Dudi  _et al._ [2018] R. Dudi, F. Pannarale, T. Dietrich, _et al._ , Physical Review D 98, 084061 (2018), arXiv:1808.09749 [gr-qc] .
  * Wijngaarden  _et al._ [2022] M. Wijngaarden, K. Chatziioannou, A. Bauswein, _et al._ , Physical Review D 105, 104019 (2022), arXiv:2202.09382 [astro-ph, physics:gr-qc] .
  * Landry and Read [2021] P. Landry and J. S. Read, The Astrophysical Journal 921, L25 (2021).
  * Legred  _et al._ [2022] I. Legred, K. Chatziioannou, R. Essick, and P. Landry, arXiv:2201.06791 [astro-ph, physics:gr-qc, physics:nucl-th] (2022), arXiv:2201.06791 [astro-ph, physics:gr-qc, physics:nucl-th] .
  * Abbott  _et al._ [2021b] R. Abbott, T. D. Abbott, F. Acernese, _et al._ , GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run (2021b), arXiv:2111.03606 [astro-ph, physics:gr-qc] .
  * Cornish and Littenberg [2015] N. J. Cornish and T. B. Littenberg, Classical and Quantum Gravity 32, 135012 (2015).
  * Littenberg and Cornish [2015a] T. B. Littenberg and N. J. Cornish, Physical Review D 91, 084034 (2015a), arXiv:1410.3852 .
  * Cornish  _et al._ [2021] N. J. Cornish, T. B. Littenberg, B. Bécsy, _et al._ , Physical Review D 103, 044006 (2021).
  * Green [1995] P. Green, J., Biometrika 82, 711 (1995).
  * Hoy  _et al._ [2021] C. Hoy, V. Raymond, A. Vijaykumar, _et al._ , Pesummary/pesummary: 0.13.0 release, Zenodo (2021).
  * Lackey and Wade [2015] B. D. Lackey and L. Wade, Physical Review D 91, 043002 (2015), arXiv:1410.8866 .
  * Golomb and Talbot [2022] J. Golomb and C. Talbot, The Astrophysical Journal 926, 79 (2022).
  * Lange  _et al._ [2018] J. Lange, R. O’Shaughnessy, and M. Rizzo, Rapid and accurate parameter inference for coalescing, precessing compact binaries (2018), arXiv:1805.10457 [gr-qc] .
  * Farr  _et al._ [2016] B. Farr, C. P. L. Berry, W. M. Farr, _et al._ , The Astrophysical Journal 825, 116 (2016), arXiv:1508.05336 .
  * Abbott  _et al._ [2020b] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ , Living Reviews in Relativity 23, 3 (2020b).
  * Petrov  _et al._ [2022] P. Petrov, L. P. Singer, M. W. Coughlin, _et al._ , The Astrophysical Journal 924, 54 (2022).
  * Oechslin  _et al._ [2002] R. Oechslin, S. Rosswog, and F.-K. Thielemann, Physical Review D 65, 103005 (2002).
  * Oechslin  _et al._ [2007] R. Oechslin, H.-T. Janka, and A. Marek, Astronomy & Astrophysics 467, 395 (2007), arXiv:astro-ph/0611047 .
  * Bauswein  _et al._ [2010a] A. Bauswein, R. Oechslin, and H.-T. Janka, Physical Review D 81, 024012 (2010a), arXiv:0910.5169 [astro-ph] .
  * Schaback and Wendland [2006] R. Schaback and H. Wendland, Acta Numerica 15, 543 (2006).
  * Rosswog [2015] S. Rosswog, SPH Methods in the Modelling of Compact Objects (2015), arXiv:1406.4224 [astro-ph] .
  * Monaghan [1992] J. J. Monaghan, Annual Review of Astronomy and Astrophysics 30, 543 (1992).
  * Hempel and Schaffner-Bielich [2010] M. Hempel and J. Schaffner-Bielich, Nuclear Physics A 837, 210 (2010), arXiv:0911.4073 .
  * Douchin and Haensel [2001] F. Douchin and P. Haensel, Astronomy & Astrophysics 380, 151 (2001), arXiv:astro-ph/0111092 .
  * Steiner  _et al._ [2013] A. W. Steiner, M. Hempel, and T. Fischer, The Astrophysical Journal 774, 17 (2013), arXiv:1207.2184 .
  * Bauswein  _et al._ [2010b] A. Bauswein, H.-T. Janka, and R. Oechslin, Physical Review D 82, 084043 (2010b).
  * Bauswein  _et al._ [2021] A. Bauswein, S. Blacker, G. Lioutas, _et al._ , Physical Review D 103, 123004 (2021), arXiv:2010.04461 [astro-ph] .
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, _et al._ , Nature Methods 17, 261 (2020), arXiv:1907.10121 .
  * Adhikari  _et al._ [2020] R. X. Adhikari, O. Aguiar, K. Arai, _et al._ , Classical and Quantum Gravity 37, 165003 (2020), arXiv:2001.11173 .
  * Ackley  _et al._ [2020] K. Ackley, V. B. Adya, P. Agrawal, _et al._ , Publications of the Astronomical Society of Australia 37, e047 (2020), arXiv:2007.03128 .
  * Maggiore  _et al._ [2020] M. Maggiore, C. van den Broeck, N. Bartolo, _et al._ , Journal of Cosmology and Astroparticle Physics 2020 (03), 050, arXiv:1912.02622 .
  * Reitze  _et al._ [2019] D. Reitze, R. X. Adhikari, S. Ballmer, _et al._ , arXiv:1907.04833 [astro-ph, physics:gr-qc] (2019), arXiv:1907.04833 [astro-ph, physics:gr-qc] .
  * Gamba  _et al._ [2021] R. Gamba, M. Breschi, S. Bernuzzi, _et al._ , Physical Review D 103, 124015 (2021).
  * Pratten  _et al._ [2021] G. Pratten, P. Schmidt, and N. Williams, Impact of Dynamical Tides on the Reconstruction of the Neutron Star Equation of State (2021), arXiv:2109.07566 [astro-ph, physics:gr-qc] .
  * Chatziioannou [2022] K. Chatziioannou, Physical Review D 105, 084021 (2022).
  * Bauswein  _et al._ [2014] A. Bauswein, N. Stergioulas, and H.-T. Janka, Physical Review D 90, 023002 (2014).
  * Sarin and Lasky [2021] N. Sarin and P. D. Lasky, General Relativity and Gravitation 53, 59 (2021), arXiv:2012.08172 .
  * Chatziioannou and Han [2020] K. Chatziioannou and S. Han, Physical Review D 101, 044019 (2020), arXiv:1911.07091 .
  * Blacker  _et al._ [2020] S. Blacker, N.-U. F. Bastian, A. Bauswein, _et al._ , Physical Review D 102, 123023 (2020), arXiv:2006.03789 [astro-ph, physics:hep-ph, physics:nucl-th] .
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, _et al._ , Nature 585, 357 (2020), arXiv:2006.10256 .
  * McKinney [2010] W. McKinney, Proceedings of the 9th Python in Science Conference , 56 (2010).
  * Hunter [2007] J. D. Hunter, Computing in Science and Engineering 9, 90 (2007).
  * Littenberg and Cornish [2015b] T. B. Littenberg and N. J. Cornish, Physical Review D 91, 084034 (2015b).