# Probing neutron stars with the full premerger and postmerger gravitational
wave signal from binary coalescences

Marcella Wijngaarden  Center for Computational Astrophysics, Flatiron
Institute, 162 5th Ave, New York, NY 10010, USA  Mathematical Sciences and
STAG Research Centre, University of Southampton, SO17 1BJ, Southampton, UK
Katerina Chatziioannou  Department of Physics, California Institute of
Technology, Pasadena, California 91125, USA  LIGO Laboratory, California
Institute of Technology, Pasadena, CA 91125, USA  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, Planckstraße 1, 64291
Darmstadt, Germany  James A. Clark  LIGO Laboratory, California Institute of
Technology, Pasadena, CA 91125, USA  Center for Relativistic Astrophysics and
School of Physics, Georgia Institute of Technology, Atlanta, GA 30332, USA
Neil J. Cornish  eXtreme Gravity Institute, Department of Physics, Montana
State University, Bozeman, Montana 59717, USA

###### Abstract

The gravitational wave signal emitted during the coalescence of two neutron
stars carries information about the stars’ internal structure. During the long
inspiral phase the main matter observable is the tidal interaction between the
binary components, an effect that can be parametrically modeled with compact-
binary solutions to general relativity. After the binary merger the main
observable is frequency modes of the remnant, most commonly giving rise to a
short-duration signal accessible only through numerical simulations. The
complicated morphology and the decreasing detector sensitivity in the relevant
frequencies currently hinder detection of the postmerger signal and motivate
separate analyses for the premerger and postmerger data. However, planned and
ongoing detector improvements could soon put the postmerger signal within
reach. In this study we target the whole premerger and postmerger signal
without an artificial separation at the binary merger. We construct a hybrid
analysis that models the inspiral with templates based on analytical
calculations and calibrated to numerical relativity and the postmerger signal
with a flexible morphology-independent analysis. Applying this analysis to
GW170817 we find, as expected, that the postmerger signal remains undetected.
We further study simulated signals and find that we can reconstruct the full
signal and simultaneously estimate both the premerger tidal deformation and
the postmerger signal frequency content. Our analysis allows us to study
neutron star physics using all the data available and directly test the
premerger and postmerger signal for consistency thus probing effects such as
the onset of the hadron-quark phase transition.

##  I Introduction

Understanding of properties of supranuclear matter in the core of neutron
stars (NSs), commonly encoded through their Equation of State (EoS) is an
ongoing challenge [1, 2, 3, 4]. Electromagnetic and gravitational wave (GW)
observations of the coalescence of two NSs provide an additional tool to study
the ultra-dense interiors of NSs, as they affect the emitted signals. Two such
GW events, GW170817 [5] and GW190425 [6], have already been detected by LIGO
[7] and Virgo [8] and planned, improved GW detectors promise more and louder
detections in the future [9]. These data will add to the interdisciplinary
effort to study the EoS, including observations of macroscopic NS properties
such as heavy pulsar observations [10, 11, 12], X-ray pulse-profiling with
NICER [13, 14, 15, 16], and nuclear theory and experiment, e.g. [17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28].

The GW signal from a binary NS (BNS) coalescence consists of two parts: a
premerger and a postmerger. In the premerger phase the binary components
inspiral toward each other as they lose orbital energy to GWs and eventually
collide [29]. In the case of GW170817, the detectable premerger signal
increased in frequency from $\sim 23$Hz to many hundreds of Hz as the orbital
separation between the two NSs rapidly decayed over $\sim 2$ minutes [30]
before finally merging at a frequency of $\sim 1500$Hz [31]. During the late
stages of the inspiral, tidal interactions accelerate the binary evolution and
leave an imprint on the signal that depends on the NS matter properties and
can thus be used to infer the EoS [32]. The premerger GW signal is typically
modeled with compact binary coalescence (CBC) templates that are based on
approximate solutions to the general relativity field equations and numerical
relativity (NR) simulations, see [33] for a recent review.

After the inevitable merger, the remnant star evolves in a way that depends on
the system mass and the EoS, see [34, 35, 36, 37] for reviews. For most EoSs
and NS masses, a hypermassive NS supported by differential rotation and
thermal effects is expected to be formed and sustained for $\sim 10-100$ms
[38]. During this time, the NS remnant emits a short-duration postmerger GW
signal with a characteristic peak at $1500-4000$Hz [39, 40, 41, 42, 43, 44,
45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57] depending on the mass and
EoS, a promising frequency range for upcoming GW detectors111Depending on the
NS mass and EoS, other possibilities for the postmerger remnant include direct
collapse to a black hole (BH) whose ringdown signal is at $\sim 6$kHz, and a
long-lived NS that emits a signal that can last for minutes or more.. In what
follows this is the case we consider. While the exact physics governing the
postmerger phase is not fully understood and NR simulations are complicated by
factors such as thermal effects, turbulence, magnetohydrodynamical
instabilites, neutrinos, and possible phase transitions, simulations make a
robust empirical prediction for the frequency of the peak of the postmerger
spectrum given the NS mass and EoS.

Information extracted from a postmerger signal about NS properties, is
complementary to the premerger. Analysis of a postmerger signal would probe
different density [58] and temperature [59, 60] regimes of the EoS, as the
more massive, hot remnant will have a higher maximum density than the
premerger NSs [40]. This complementarity can lead to insights about high-
density phenomena such as phase-transitions between the densities probed by
the premerger and postmerger signals [61, 62, 63, 64, 65, 66, 67, 68, 69].
Finally, detection of a postmerger signal would allow us to determine the
nature of the merger remnant, namely NS or BH, which has implications for a
potential electromagnetic counterpart and its interpretation [70] as well as
the determination of the threshold mass for prompt collapse [71].

In the case of GW1708017 and GW190425 only piecewise analyses of the GW signal
have been performed focusing either on the premerger or on the postmerger part
of the signal [30, 72, 73, 6]. This reduces the computational cost
significantly, as the postmerger signal requires a large sampling rate
(typically $8194$Hz) while the premerger analysis involves long duration data
segments (typically $128$s). As no postmerger signal is detectable yet [30,
72], excluding it from the premerger analysis did not bias the results for the
binary parameters [74]. However, looking forward and in the case of a
postmerger signal detection, separate analyses will not be able to exploit
phase coherence through merger or parameter relations, such as those between
the tidal deformability from the premerger signal and the frequency peak of
the postmerger spectrum [52, 58, 75].

Another reason for separate premerger and postmerger analyses is that the
morphology and details of the postmerger signal are not well understood,
making it difficult to construct a first-principles physical model such as the
premerger CBC templates. While approximate analytic models for the postmerger
signal have been proposed with the aid of NR simulations [76, 77, 78, 79, 80,
81, 82], these are still limited in accuracy as well as by uncertainties of
the simulations they are based on [80]. An alternative are morphology-
independent analyses that are not limited to a specific signal type [83, 84,
85]. Specifically, a model-agnostic approach with the BayesWave algorithm [86,
87, 88] has been shown to accurately determine the main features of the
postmerger signal [85, 31], and in some cases do so more accurately than
tailored models [81].

BayesWave models a GW signal with a sum of sine-Gaussian wavelets. Both the
number of wavelets and their parameters are marginalized over, resulting in a
flexible analysis that has been applied to a variety of signals [89, 90, 91,
92, 93, 94, 95, 96, 97, 98]. BayesWave’s sine-Gaussian wavelets are
particularly adept for postmerger signals that are dominated by distinct
frequency components, as each component can be approximately modeled by a
wavelet. In this context, BayesWave has been used to search for a short-
duration, high-frequency signal after both GW170817 [30] and GW190425 [6],
returning null results and upper limits on the energy content [30]. A detailed
study of BayesWave’s sensitivity to postmerger signals is presented in [85].

In this paper we construct a hybrid analysis of the full BNS GW signal that
targets both the premerger and the postmerger data. We simultaneously analyze
the full signal using 1) a CBC template to describe the well-modeled premerger
phase and 2) sine-Gaussian wavelets to capture the less well understood
postmerger. The parameters of the CBC template and the wavelets are
simultaneously sampled over to obtain the combined multidimensional posterior
for all components of the model. Those parameters include the premerger tidal
deformability that quantifies the inspiral tidal deformation and the peak
frequency of the postmerger spectrum. Our analysis extracts both
simultaneously and allows for direct consistency comparisons.

After validating the hybrid modeled/unmodeled analysis on a toy model based on
GW150914 data, we apply it to GW170817 and find that the postmerger signal
remains undetected, as expected from the upper limit estimates of [30]. We
further analyze simulated signals of high signal-to-noise ratio (SNR) for
which the postmerger signal is detectable. We show that our analysis can
reconstruct the premerger signal as well as the main components of the
postmerger signal such as the dominant frequency mode. We simulate signals for
which the premerger and postmerger parts are consistent (corresponding to a
hadronic EoS) and inconsistent (corresponding to an EoS with a strong phase
transition in the relevant density scales) and demonstrate that our analysis
can extract either behavior.

The rest of the paper is organized as follows. In Sec. II we describe in
detail our analysis and algorithm. In Sec. III we present a proof-of-principle
analysis on GW150914. In Sec. IV we reanalyze the GW170817 data. In Sec. V we
apply our analysis to simulated data. In Sec. VI we conclude.

##  II Methodology and Models

In order to analyze a GW signal that contains both modeled and unmodeled
features we employ a hybrid analysis that makes use of both CBC templates and
flexible models for the signal. We base our analysis on the morphology-
agnostic data analysis algorithm BayesWave [86, 87] by extending it to account
for CBC templates similar to [99]. In its core functionality BayesWave
simultaneously models a GW signal (referred to as the “signal model” in
BayesWave literature), instrumental glitches (the “glitch model”), and the
gaussian detector noise (the “PSD model”). BayesWave uses minimal assumptions,
i.e., no physical model, to describe these components. Instead, the signal and
the glitch are modeled as the sum of a variable number of sine-Gaussian
wavelets, whereas the noise model describes the power spectral density (PSD)
of the Gaussian noise with a variable number of spline points and Lorentzians.
For this work and [99] a fourth component has been added, which models GW
signals through CBC templates obtained by solving the two-body problem in
general relativity; the CBC model.

###  II.1 Overview of models

BayesWave uses sampling methods to characterize the multi- and variable-
dimension posterior distribution for independent models that target different
features of the GW data. The goal is to draw samples from the posterior
distribution function $p(M|d)$, the probability that a model $M$ describes the
data $d$ and a prior $p(M)$, defined as

$\displaystyle p(M|d)=\frac{p(M)p(d|M)}{p(d)},$ (1)

where $p(d)$ is the evidence, and $p(d|M)$ expresses the likelihood of the
data for a given $M$. We model the data in each interferometric detector
$d_{I}$ as a linear combination of multiple components

$\displaystyle d_{I}=g_{I}+h^{\mathrm{cbc}}_{I}+h^{\mathrm{w}}_{I}+n_{I}\equiv
M_{I}+n_{I},$ (2)

where $g_{I}$ denotes any detector glitches222For the remainder of the paper
we will work with data without instrumental glitches and thus ignore the
glitch model $g_{I}$. As we will not be making use of the glitch model, we
will also refer to the wavelet signal model simply as the wavelet model to
avoid confusion with the CBC signal model., the full GW signal
($h^{\mathrm{cbc}}_{I}+h^{\mathrm{w}}_{I}$) consists of the CBC waveform model
$h^{\mathrm{cbc}}_{I}$ and any signal (such as the postmerger) that is not
included in the CBC template and is captured by wavelets $h^{\mathrm{w}}_{I}$,
and $n_{I}$ is the Gaussian component of the noise. The combination of a
signal, glitch, and noise defines the full model $M$. Given such a model, the
likelihood function $p(d_{I}|M_{I})$, in short $L_{I}$, can be computed by
noting that the residual $r_{I}=d_{I}-M_{I}$ is consistent with Gaussian noise
(see also Eq. (3) of [87]):

$\displaystyle\log{L_{I}}=-\frac{1}{2}(r_{I}|r_{I})+C_{I},$ (3)

where $C_{I}$ is a constant that depends on the PSD of the detector noise
$S_{n,I}(f)$, and the noise weighted inner product is defined as

$\displaystyle(a|b)\equiv
2\int\frac{a^{*}(f)b(f)+b^{*}(f)a(f)}{S_{n,I}(f)}df,$ (4)

where an asterisk denotes a complex conjugate. The remaining ingredients of
the analysis are a specific model for each data component, the corresponding
prior, and the sampling procedure.

We express the GW signal at geocenter and then project it onto each detector
in the network to compute the detector response in the frequency domain

$\displaystyle h^{{\mathrm{cbc/w}}}_{I}(f)$ $\displaystyle=e^{2\pi i\Delta
t(\alpha,\delta)}$ (5)
$\displaystyle\left[F^{+}(\alpha,\delta,\psi)h^{{\mathrm{cbc/w}}}_{+}+F^{\times}(\alpha,\delta,\psi)h^{{\mathrm{cbc/w}}}_{\times}\right].$

The intrinsic part of the signal for each GW polarization mode is given by
$h^{{\mathrm{cbc/w}}}_{+}$ and $h^{{\mathrm{cbc/w}}}_{\times}$ and it is
different under the CBC or the wavelet model. The projection then involves a
time delay in the arrival of the signal in the different detectors $\Delta
t(\alpha,\delta)$, the detector antenna patterns $F^{+}(\alpha,\delta,\psi)$,
$F^{\times}(\alpha,\delta,\psi)$ for each GW polarization mode, parametrized
through the sky location $\alpha$, $\delta$ corresponding to the right
ascension and declination respectively, and the polarization angle $\psi$. The
projection is independent of whether the signal is expressed through CBC
templates or wavelets so the functional form of Eq. (5) and the parameters
$(\alpha,\delta,\psi)$ are common. The ratio
$h^{{\mathrm{cbc/w}}}_{+}/h^{{\mathrm{cbc/w}}}_{\times}$ is related to the
binary inclination $\iota$ and is also shared between the CBC and wavelet
models.

By construction BayesWave breaks down the full parameter space defined by all
models in different blocks and uses a blocked Gibbs sampler to explore the
posterior. Each block of parameters is sampled over with different (reversible
jump [100]) Markov Chain Monte Carlo (RJMCMC) samplers. Parallel tempering is
employed both for efficient posterior exploration and evidence calculation
[101], so all samplers share the same number and temperature of chains and
exchanges between the chains proceed for all models simultaneously. In order
for this procedure to lead to efficient exploration of the posterior
distribution, each block of parameters needs to be selected to be as
independent from and uncorrelated to other blocks as possible.

The sampler ensemble consists of (1) a CBC MCMC, Sec. II.2, (2) an RJMCMC
sampler for the wavelet intrinsic parameters, Sec. II.3, (3) an MCMC sampler
for the common extrinsic parameters for the CBC and wavelet models, Sec. II.4,
and (4) an RJMCMC noise model sampler, Sec. II.5. The structure is given in
Fig. 1 and it is similar to the code used in [99] but instead of
nonoverlapping parameter blocks, we now allow for a set of common parameters
to be updated in multiple blocks as detailed in the subsequent subsections.
Sampling then proceeds iteratively both within each block (${\cal{O}}(10^{2})$
iterations at a time) and between different blocks (${\cal{O}}(10^{4})$
iterations). In the rest of this section, we discuss each (RJ)MCMC, which
parameters it updates, and how it relates to the other samplers. The
discussion is structured to be as self-contained and detailed as possible;
readers primarily interested on the results can find them starting in Sec.
III.

Figure 1: General code workflow. After some preconditioning, the data are used
to obtain a quick fit to the noise PSD $S_{n}(f)$ and CBC
$\theta^{\mathrm{cbc}}$ and extrinsic $\theta_{\mathrm{ext}}$ parameters.
These are then used as a starting point for the blocked Gibbs sampler (red
box) which consists of $4$ independent samplers. Each iteration of the Gibbs
sampler consists of a wavelet RJMCMC (blue box) that updates the wavelet
parameters $\theta^{\mathrm{w}}$, an extrinsic MCMC (green box) that updates
the extrinsic parameters $\theta_{\mathrm{ext}}$, a noise RJMCMC (gray box)
that updates $S_{n}(f)$, and a CBC MCMC (pink box) that updates
$\theta^{\mathrm{cbc}}$. Orange rectangles enclose the parameter that is being
updated from the preceding box. Each block consists of ${\cal{O}}(10^{2})$
iterations, while the entire red Gibbs box consists of ${\cal{O}}(10^{4})$
iterations.

###  II.2 CBC-specific parameters: CBC sampler

A nonprecessing quasicircular inspiral of two compact objects is characterized
by a set of $6$ intrinsic source parameters and a set of $7$ extrinsic
parameters that describe the binary’s location and orientation relative to the
detectors. The CBC sampler updates the parameters that are specific to
$h^{\mathrm{cbc}}_{+}$ and $h^{\mathrm{cbc}}_{\times}$, namely masses, spins,
tidal parameters (for BNSs), distance, time, and phase of coalescence. In Fig.
1 the CBC sampler is depicted by the pink box “CBC MCMC” and the relevant
parameters are collectively denoted by $\theta^{\mathrm{cbc}}$.

####  II.2.1 Waveform model and Parameters

The CBC sampler is integrated with the LALSimulation [102] waveforms suite and
can make use of any non-precessing waveform model; we here focus on PhenomD
[103, 104] and PhenomDNRT [105, 106, 107]. The former includes the full
inspiral-merger-ringdown (IMR) signal expected from a BBH coalescence. The
latter targets the inspiral of BNS systems as it includes the effect of tidal
interactions between the NSs, manifesting as a change in the waveform phase.
PhenomDNRT models the long inspiral phase and instead of terminating abruptly
at merger, it decays rapidly through a window function after the merger
frequency [106, 107]. PhenomDNRT is well suited for our purposes regarding BNS
signals as the waveform remains smooth and continuous, naturally allowing the
resulting CBC (inspiral) + wavelet (postmerger) model to be smooth without any
forced stitching.

A total of $9$ CBC parameters are updated by the CBC sampler
$\theta^{\mathrm{cbc}}\equiv(m_{1},m_{2},\chi_{1},\chi_{2},\Lambda_{1},\Lambda_{2},D_{L},t_{c},\phi_{c})$.

  * •

The masses of the binary components are $m_{1}$ and $m_{2}$ with the
convention $m_{1}\geq m_{2}$. We also parametrize the masses through the total
mass $M=m_{1}+m_{2}$, chirp mass
${\cal{M}}=\left(m_{1}m_{2}\right)^{3/5}/\left(m_{1}+m_{2}\right)^{1/5}$, mass
ratio $q=m_{2}/m_{1}<1$, and symmetric mass ratio
$\eta=({\cal{M}}/M)^{5/3}<1/4$.

  * •

The dimensionless spin components aligned with the orbital angular momentum
are $\chi_{1}$ and $\chi_{2}$. Current detector sensitivities make individual
spin components hard to measure [108, 109] so a commonly used spin combination
is the mass-weighted effective spin parameter,

$\displaystyle\chi_{\text{eff}}=\frac{m_{1}\chi_{1}+m_{2}\chi_{2}}{m_{1}+m_{2}},$
(6)

which is conserved under spin-precession to at least the second post-Newtonian
order [110].

  * •

Tidal interactions in BNS signals are encoded by the NS dimensionless tidal
deformabilities

$\displaystyle\Lambda_{1,2}=\frac{2}{3}k_{1,2}\left(\frac{R_{1,2}}{m_{1,2}}\right)^{5},$
(7)

where $k_{1,2}$ is the tidal Love number, $R_{1,2}$ is the radius, and
$m_{1,2}$ is the mass of each binary component. To leading order, the tidal
effects are imprinted in the waveform through the binary tidal deformability
[111, 112]

$\displaystyle\tilde{\Lambda}=\frac{16}{13}\frac{(12q+1)\Lambda_{1}+(12+q)q^{4}\Lambda_{2}}{(1+q)^{5}},$
(8)

with subleading terms introducing another tidal parameter that is generally
not measurable, $\delta\tilde{\Lambda}$ [112]. For BBH systems,
$\Lambda_{1}=\Lambda_{2}=0$ [113, 114].

  * •

The CBC sampler further handles $3$ extrinsic parameters either because they
do not explicitly appear in the wavelet model (time of coalescence $t_{c}$,
luminosity distance $D_{L}$) or because they are correlated with intrinsic
parameters and need to be simultaneously sampled for efficiency (phase of
coalescence $\phi_{c}$). Despite the well-known correlation between the
distance and the inclination [115], we do _not_ update the latter here as the
inclination affects the wavelet model as well.

####  II.2.2 Priors

The priors on the binary intrinsic parameters depend on whether we consider
BNSs or BBHs. We sample the detector frame chirp mass
${\cal{M}}\in\left[0.23,2.6\right]M_{\odot}$ and total mass
$M\in\left[1.0,6.0\right]M_{\odot}$ for BNSs and
${\cal{M}}\in\left[0.23,174\right]M_{\odot}$,
$M\in\left[1.0,400.0\right]M_{\odot}$ for BBHs and apply a Jacobian correction
that makes the prior flat for the individual component masses $m_{1}$ and
$m_{2}$ [116]. We use flat priors for the dimensionless spin components
$\chi_{1},\chi_{2}\in\left[-0.05,0.05\right]$ for BNSs and
$\chi_{1},\chi_{2}\in\left[-1,1\right]$ for BBHs. In the BNS case we sample
the effective dimensionless tidal parameters uniform
$\tilde{\Lambda}\in\left[0,1000\right]$ and
$\delta\tilde{\Lambda}\in\left[-500,500\right]$ with the constraint that the
corresponding individual component tidal deformabilities cannot be negative.
We can optionally directly sample the individual component tidal parameters
uniformly with the default prior range $\Lambda_{1,2}$ $\in$
$\left[0,1000\right]$ though all results presented here use the former tidal
prior. In all cases, the individual component tidal deformabilities are
assumed to be independent of each other and not linked under any assumption
about the NS EoS [117, 118, 119].

We sample the merger time relative to the GPS trigger time
($t_{\textrm{trig}}$) with a default window of $t_{c}$ $\in$
$\left[t_{\textrm{trig}}-0.5\text{s},t_{\textrm{trig}}+1.5\text{s}\right]$.
Alternatively, it is possible to use a symmetric window of arbitrary size
around the trigger time. We also sample in the logarithmic luminosity distance
with an adjustment to the Jacobian such that the prior is uniform in volume
(i.e., uniform in $~{}\text{D}_{L}^{3}$) with $\text{D}_{L}$ $\in$
$\left[1,10000\right]\text{Mpc}$. Finally, the coalescence phase prior is
uniform in $\phi_{c}\in\left[0,2\pi\right]$.

####  II.2.3 Data

Since the CBC and the wavelet models are active simultaneously, the overall
likelihood has to take into account both CBCs and wavelets. The blocked nature
of the code means that the wavelet model is not updated in the CBC sampler, so
the wavelet-subtracted data $d^{\mathrm{cbc}}\equiv d-h^{\mathrm{w}}$ are
constant as the CBC parameters are being updated. We compute
$d^{\mathrm{cbc}}=d-h^{\mathrm{w}}$ once at the beginning of the CBC sampler
block given the current wavelet parameters. The wavelet-subtracted data
$d^{\mathrm{cbc}}$ then play the role of the data in the sense that
$d^{\mathrm{cbc}}-h^{\mathrm{cbc}}$ is the residual during the likelihood
computation.

####  II.2.4 Proposals

We employ the CBC sampler described in detail in [120] within the BayesWave
blocked Gibbs sampler, see Fig. 1. The sampler performs by default 300
iterations that update the CBC parameters described above using a mixture of
proposals:
Fisher jumps: Gaussian jumps along the eigenvectors of the Fisher information
matrix. The Fisher information matrix is recomputed at the start of each CBC
sampler block for each parallel chain and using each chain’s current location.
Differential evolution [121]: Two samples from each chain’s history propose a
new location. The initial history array is drawn from the prior and then
updated every 20 iterations. As each CBC sampler block does 300 iterations by
default, the history array contains information about previous blocks.
Jiggle: Small Gaussian jumps in individual parameter directions. The size of
the proposed jump is scaled with the chain temperature such that high
temperature chains have larger jumps and can explore a larger area of the
parameter space. The jiggle proposal allows for exploring areas outside those
covered by Fisher matrix jumps and are especially useful early on.
Mixed: Uniform proposal for the tidal parameters and differential evolution
(see above) for the other parameters. This proposal is employed $5\%$ of the
time and aids convergence for cases where the tidal parameters cannot be
tightly constrained. The proposal facilitates taking larger steps for the
tidal parameters whereas the other better constrained parameters evolve with
smaller steps, thus leading to a higher acceptance rate than a fully uniform
proposal.

####  II.2.5 Heterodyned likelihood

We integrate the heterodyned likelihood technique [122, 120] in the BayesWave
CBC (and extrinsic, to be described later) sampler in order to speed up the
likelihood evaluation. The technique is based on the fact that all CBC
waveforms that contribute to the posterior (i.e. have a likelihood high enough
that the MCMC proposal can be accepted) are similar. The likelihood of a model
$h^{\mathrm{cbc}}$ can then be computed based on a fixed and precomputed
reference model $h^{\mathrm{cbc}}_{\mathrm{ref}}$ that is similar to
$h^{\mathrm{cbc}}$, and given the noise PSD and the data.

We can rewrite the residuals $r$ using the reference heterodyne, omitting the
detector subscript $I$ for clarity, as:

$\displaystyle
r=d^{\mathrm{cbc}}-h^{\mathrm{cbc}}=d^{\mathrm{cbc}}-(h^{\mathrm{cbc}}_{\mathrm{ref}}-\Delta
h^{\mathrm{cbc}})=r_{\mathrm{ref}}+\Delta h^{\mathrm{cbc}},$ (9)

where $\Delta h^{\mathrm{cbc}}$ is the difference between the reference
waveform and the exact waveform at a given set of CBC parameters,
$d^{\mathrm{cbc}}$ is the wavelet-subtracted data as given to the CBC sampler,
$d^{\mathrm{cbc}}=d-h^{\mathrm{w}}$, and $r_{\mathrm{ref}}\equiv
d^{\mathrm{cbc}}-h^{\mathrm{cbc}}_{\mathrm{ref}}$. The likelihood from Eq. (3)
then becomes

$\displaystyle\log{L}\sim-\frac{1}{2}(r_{\mathrm{ref}}+\Delta
h^{\mathrm{cbc}}|r_{\mathrm{ref}}+\Delta h^{\mathrm{cbc}}),$ (10)

where we omit the second term of Eq. (3) as the PSD is held fixed during the
CBC updates and the term does not affect the relative likelihoods. This form
of the likelihood enables several computational advantages as we split the
inner product in its components

$\displaystyle-2\log{L}\sim$
$\displaystyle(r_{\mathrm{ref}}|r_{\mathrm{ref}})+(\Delta
h^{\mathrm{cbc}}|\Delta h^{\mathrm{cbc}})$ (11)
$\displaystyle+\left[(r_{\mathrm{ref}}|\Delta h^{\mathrm{cbc}})+(\Delta
h^{\mathrm{cbc}}|r_{\mathrm{ref}})\right].$

A detailed description of how each components of the heterodyned likelihood
above is calculated is available in [122, 120]; we here briefly discuss how
this is adapted and implemented within the BayesWave blocked Gibbs sampler.
The first term of the inner product $(r_{\mathrm{ref}}|r_{\mathrm{ref}})$ is
the most straightforward as it is computed once at the beginning of each CBC
sampler block using the current wavelet-subtracted data and CBC parameters for
each chain.

The key performance improvement of the heterodyne likelihood is achieved by
recognizing that the latter two terms in Eq. (11) can be further split up into
the product of two terms: a slowly varying term that is computed at each
iteration, and a rapidly varying term that is computed once. Both terms are
evaluated through a sum over frequencies with bin width $\Delta f$ using a
Legendre polynomial expansion [120]. The rapidly-varying term requires a small
$\Delta f$ but this is computed once. The slowly-varying term that needs to be
calculated on each likelihood evaluation can instead make use of a coarse
frequency grid of $\Delta f\leq$ 4 Hz thus significantly reducing
computational cost.

The term $(\Delta h^{\mathrm{cbc}}|\Delta h^{\mathrm{cbc}})$ in Eq. (11)
introduces a slowly varying phase difference as a function of frequency
between the reference waveform and the proposed waveform. However, in the
noise weighted inner product, Eq. (4) the denominator might not be slowly
varying if the PSD includes sharp spectral lines. Following the approach
discussed above, we split the PSD in a smooth broadband component that is
slowly varying as a function of frequency, and line features which are only
evaluated at discrete intervals. If the PSD is held fixed throughout the
analysis, the smooth part is computed once at the beginning. When using a
variable PSD model instead, we update the smooth PSD component at the
beginning of the CBC sampler block.

The final term of Eq. (11) combining $r_{\mathrm{ref}}$ and $\Delta
h^{\mathrm{cbc}}$ is split up into the product of a slowly and rapidly varying
part by introducing the smooth component of the PSD in the numerator and
denominator of, respectively, the rapidly and slowly varying part. Due to
this, the heterodyned difference in the waveforms $\Delta h^{\mathrm{cbc}}$ is
effectively whitened using only the smooth PSD component on a coarse frequency
grid. The rapidly varying part consisting of the heterodyned whitened
residuals $r_{\mathrm{ref}}$ now includes an extra product with the smooth
part of the PSD, but since this is only computed at the start of the CBC
sampler block, it adds little computational cost.

The heterodyned likelihood speeds up the likelihood calculation by a factor of
$\sim T_{\mathrm{obs}}\Delta f$, where $T_{\mathrm{obs}}$ is the duration of
the analyzed data, thus making it especially useful for long duration signals
such as BNSs. The heterodyne implementation in BayesWave is more complex than
the one in [120] as it has to account for a potentially changing PSD and the
wavelet model. We thus expect more modest overall computational improvements,
though the speed-up of each individual likelihood calculation is consistent
with [120].

###  II.3 Wavelet parameters: Wavelet sampler

BayesWave has a common RJMCMC sampler for its signal and glitch models where
the non-Gaussian features (either a GW signal or a detector glitch) are
modeled as a sum of a variable number of sine-Gaussian wavelets which can be
expressed in the time-domain as

$\displaystyle\psi(t;A,f_{0},\tau,t_{0},\phi_{0})=Ae^{-(t-t_{0})^{2}/\tau^{2}}\cos\left(2\pi
f_{0}(t-t_{0})+\phi_{0}\right),$ (12)

where $\tau\equiv Q/(2\pi f_{0})$ is the decay time, $Q$ is the quality
factor, $A$ is the wavelet amplitude, $f_{0}$ is the wavelet central
frequency, $t_{0}$ is the wavelet central time and $\phi_{0}$ is the phase
offset. The glitch model uses independent sums of wavelets for the different
detectors, whereas the signal models use a set of wavelets at geocenter and
projects them onto detectors through Eq. (5) with

$\displaystyle h^{\mathrm{w}}_{+}(f)$
$\displaystyle=\sum_{i=0}^{N_{s}}\psi(f;A_{i},f_{0,i},Q_{i},t_{0,i},\phi_{0,i})$
(13) $\displaystyle h^{\mathrm{w}}_{\times}(f)$ $\displaystyle=\epsilon
h^{\mathrm{w}}_{+}(f)e^{i\pi/2},$ (14)

where $\epsilon$ is the ellipticity parameter and $N_{s}$ is the number of
wavelets in the signal model. Details about the wavelet sampler and
priors/proposals for each wavelet parameter and the number of wavelets are
discussed in [86, 87] and remain mostly unaltered in our analysis. The prior
on the number of wavelets has an upper limit of $100$. We extend the standard
settings by adding an option to limit the central frequency range of the
wavelets to
$f_{0}\in\left[f^{\mathrm{w}}_{\mathrm{min}},f^{\mathrm{w}}_{\mathrm{max}}\right]$
which defaults to the analysis bandwidth. Since the ellipticity $\epsilon$ is
related to the remaining extrinsic parameters and to the binary inclination,
it is not sampled by the wavelet sampler but the extrinsic sampler discussed
next.

In Fig. 1 the wavelet sampler is depicted by the blue box termed “wavelet
RJMCMC” where $\theta^{\mathrm{w}}$ corresponds to the relevant parameters:
the number of wavelets and the amplitude, quality factor frequency, time, and
phase of each wavelet. Similarly to the CBC sampler, the wavelet sampler makes
use of the CBC-subtracted data $d^{\mathrm{w}}\equiv d-h^{\mathrm{cbc}}$,
which are computed once at the beginning of each wavelet sampler block using
the current CBC sample.

###  II.4 Common parameters: Extrinsic sampler

Since both the CBC and wavelet models target the same GW source, they share a
number of extrinsic source parameters that govern the signal projection onto
the detector network, Eq. (5). These common parameters $\theta_{\mathrm{ext}}$
are updated by the extrinsic sampler block, denoted as “extrinsic MCMC” in
Fig. 1.

####  II.4.1 Parameters and priors

The extrinsic parameters shared between the wavelet and CBC models are
$\theta_{\mathrm{ext}}\equiv(\alpha,\delta,\psi,\epsilon,\varphi_{0},A^{\mathrm{cbc}})$.

  * •

The source sky location is given by the right ascension
$\alpha\in\left[0,2\pi\right]$ and the declination
$\sin{\delta}\in\left[0,1\right]$, and the corresponding prior is flat on the
sphere.

  * •

The polarization angle $\psi$ is sampled with a uniform prior
$\psi\in\left[0,2\pi\right]$.

  * •

The degree of elliptical polarization in the signal is encoded in
$\epsilon\in\left[-1,1\right]$, Eq. (14). For elliptically polarized signals,
the ellipticity parameter can be related to the orbital inclination with
respect to the line of sight to the binary $\iota$ as

$\displaystyle\cos{\iota}=\frac{1-\sqrt{1-\epsilon^{2}}}{\epsilon}.$ (15)

We use a uniform prior in $\cos{\iota}\in\left[-1,1\right]$, which is a
deviation from previous works where BayesWave samples uniformly in $\epsilon$
instead.

  * •

An overall phase shift $\varphi_{0}$ is sampled uniformly as
$\varphi_{0}\in\left[0,2\pi\right]$ and it is equivalent to a phase shift in
all wavelets simultaneously (see Eq. 12). Despite the degeneracy, introducing
this parameter leads to more efficient convergence in the wavelet model.

  * •

We also introduce an overall amplitude scaling parameter $A^{\mathrm{cbc}}$,
which is sampled uniformly and is degenerate with the CBC luminosity distance
$D_{L}$. We restrict $A^{\mathrm{cbc}}$ to correspond to values for $D_{L}$
within the $D_{L}$ prior range and add a Jacobian factor such that the prior
remains uniform in volume, i.e., $(D_{L}/A^{\mathrm{cbc}})^{3}$.

The extrinsic parameters can be efficiently sampled independently from the CBC
intrinsic parameters by noting that a change in the extrinsic parameters
induces an amplitude, time, and phase shift in the signal seen in each
detector, while leaving its phase evolution unaltered. We therefore compute
the geocenter CBC waveform given the fixed intrinsic parameters once and then
apply the shifts accordingly to evaluate the waveform at each new proposed set
of extrinsic parameters. This has the advantage that we can sample the
extrinsic parameters efficiently while not having to recompute the full,
expensive CBC waveform phase at each iteration.

In practice we write the CBC response at a given detector $I$ using the
geocenter waveform $h^{\mathrm{cbc}}_{+}$ and applying the shifts as

$\displaystyle
h^{\mathrm{cbc}}_{I}(f)=h^{\mathrm{cbc}}_{+}(f)F^{\mathrm{cbc}}_{I}e^{i\lambda_{I}}e^{2\pi
if\Delta t_{I}},$ (16)

where $\Delta t_{I}$ is the arrival time at the detector relative to the
geocenter, the phase shift is

$\displaystyle\lambda_{I}\equiv\arctan[\epsilon
F^{\times}(\alpha,\delta,\psi)/F^{+}(\alpha,\delta,\psi)]+\phi_{c},$ (17)

and $F^{\mathrm{cbc}}_{I}$ is the magnitude of the detector response

$\displaystyle F^{\mathrm{cbc}}_{I}\equiv
A^{\mathrm{cbc}}\sqrt{\left[\frac{\cos{\iota}}{\epsilon}F^{+}(\alpha,\delta,\psi)\right]^{2}+\left[\cos{\iota}\,F^{\times}(\alpha,\delta,\psi)\right]^{2}},$
(18)

Even though $A^{\mathrm{cbc}}$ is completely degenerate with $D_{L}$, we find
that it increases sampling efficiency in the extrinsic sampler due to the
distance-inclination degeneracy [115].

At the end of the block of extrinsic updates, the scaling and shift parameters
are reset and propagated to the corresponding CBC parameters before proceeding
to the CBC sampler

$\displaystyle D_{L}$ $\displaystyle\rightarrow\frac{D_{L}}{A^{\mathrm{cbc}}}$
(19) $\displaystyle t_{c}$ $\displaystyle\rightarrow t_{c}+\Delta t_{a}$ (20)
$\displaystyle\phi_{c}$ $\displaystyle\rightarrow\phi_{c}-\Delta\varphi_{0}$
(21)

In these equations $\Delta t_{a}$ and $\Delta\varphi_{0}$ correspond to the
overall time and phase shift induced on the waveform by the change in
extrinsic parameters. The time of arrival of the merger at geocenter $t_{c}$
changes due to the new sky location and this shift is encoded in $\Delta
t_{a}$.

####  II.4.2 Proposals

To update the parameters $\theta_{\mathrm{ext}}$ the extrinsic sampler uses a
mixture of proposals, which are described in detail in [86, 87]. In the
following we focus on some CBC-specific updates made to the existing
proposals.

The most commonly used proposal is the Fisher proposal, where jumps are
proposed along the eigenvectors of the Fisher information matrix. The matrix
elements, eigenvalues, and eigenvectors are computed once at the start of the
extrinsic sampler block by finite differencing using the current chain
position. We extend the existing BayesWave Fisher calculation by adding the
CBC model such that the final Fisher matrix captures information from the full
CBC+wavelet signal. When the sky location is known and fixed (for example for
sources with an electromagnetic counterpart), the Fisher matrix is reduced to
the parameters that are varied.

The BayesWave extrinsic sampler also employs a two-part proposal that exploits
a degeneracy of the extrinsic parameters when multiple detectors are used
(also proposed by [123]). The first part consists of a sky-ring proposal that
finds a new sky location ($\alpha$, $\sin{\delta}$) such that the time delays
between the detectors are preserved [124, 86]. The second part of the proposal
updates the remaining parameters ($\psi$, $\cos{\iota}$, $A^{\mathrm{CBC}}$
and $\varphi_{0}$) either uniformly in their prior range or deterministically
such that the waveform is identical at the new sky location compared to the
waveform at the previous sky location, see [120] for details. Since the
waveform is conserved, so is the likelihood, though this does not guarantee
that the proposed sky location will be accepted due the necessary Jacobian
factor [120].

Finally, the BayesWave extrinsic sampler includes a uniform proposal that
draws from the prior of each parameter. These jumps are particularly useful
for higher temperature chains to explore the entire space available as well as
early on in the sampling if a signal has not been found.

####  II.4.3 Heterodyned likelihood

A change in the extrinsic parameters $\theta_{\mathrm{ext}}$ will affect both
the CBC and the wavelet models, therefore both need to be taken into account
when computing the likelihood. However, we find that the CBC part of the model
dominates the computational cost so we again use the heterodyned likelihood to
speed up the calculation. In this case, there are some differences compared to
the implementation for the CBC sampler as discussed in Sec. II.2.5 because we
can no longer assume that the wavelet-subtracted data are fixed. Indeed, a new
set of extrinsic parameters changes the full signal projection and the full
residual $r=d-h^{\mathrm{w}}-h^{\mathrm{cbc}}$ has to be recomputed on each
iteration.

Despite this, we can achieve a computational improvement by splitting up the
heterodyne computation even further to avoid recomputing all components of the
heterodyne when the residuals change, i.e., at every iteration. All
subcomponents of the heterodyne that do not involve the residuals are computed
and stored once at the start of each extrinsic MCMC sampler block, whereas the
remaining components are computed at each iteration. If the signal wavelet
model is not employed (i.e. an analysis with only the CBC model), the
heterodyne likelihood implementation falls back to the same, more efficient,
version as for the CBC sampler.

Since the likelihood calculation is more involved and parts of the heterodyne
procedure need to be recomputed in each iteration, we expect more modest
computational improvements compared to those of Sec. II.2.5.

###  II.5 Noise parameters: PSD sampler with BayesLine

The PSD of the Gaussian noise is modeled as a combination of splines and
Lorenzians targeting the broadband behavior and spectral lines respectively
and marginalized over through an RJMCMC in BayesLine [125, 126, 87]. The
analysis presented here makes no updates on this noise marginalization
procedure and can optionally include BayesLine as one of the blocks,
corresponding to the “Noise RJMCMC” gray box in Fig. 1. Alternatively, the
analysis can skip the noise marginalization and use a predetermined fixed PSD.
If the PSD is modeled by BayesLine, the heterodyne procedure needs to be
repeated at the end of the BayesLine sampling block and before moving to CBC
sampling with the latest PSD sample. However, this coincides with a
recomputing of the heterodyne elements due to a change of the wavelet model
that is subtracted from the data so no further computations are necessary. In
the current implementation [125], the wavelet model is subtracted from the
data before entering the BayesLine sampling block and we extend this
subtraction to include the CBC model.

##  III Toy model: GW150914

To illustrate the wavelet+CBC model, we begin with a toy model by analyzing
the first detected BBH merger, GW150914 [127]. In reality, the full BBH signal
can be modeled with existing CBC templates, something confirmed in the case of
GW150914 through comparisons between different models and physical effects
[128]. In our case, and in order to test our hybrid analysis on what is
probably the best studied signal, we pretend we lacked efficient modeling of a
BBH merger and waveform templates terminated shortly before merger.

We analyze 4s of data from LIGO Hanford and LIGO Livingston available through
the Gravitational Wave Open Science Center (GWOSC) [129, 130] in the frequency
band from 16Hz to 2048Hz. The noise PSD is marginalized over as discussed in
Sec. II.5. For the CBC model, we use either the PhenomD waveform or the
PhenomDNRT waveform. PhenomD is appropriate for BBH systems as it models the
full IMR signal. On the other hand, we do not expect PhenomDNRT to accurately
model the BBH data as it lacks an accurate merger and ringdown portion due to
tapering [106]. Regardless, we use it to demonstrate how any leftover signal
that is not captured by the CBC waveform can be captured by the wavelets.

Label | GPS time (s) | Model | Waveform |  | Wavelet
---
window (s)
| Wavelet
---
bandwidth (Hz)
CBCIMR | 1126259462.391 | CBC/noise | PhenomD | N/A | N/A
CBCInsp | 1126259462.391 | CBC/noise | PhenomDNRT | N/A | N/A
CBCInsp+wavelet | 1126259462.416 | CBC/wavelet/noise | PhenomDNRT | (-0.025, 0.025) | (150, 2048)
Table 1:  Settings for the runs of Sec. III on GW150914. All runs use a
segment length of 4s, a low frequency cut off of 16Hz, and a sampling rate of
4096Hz. From left to right, columns correspond to the run label in Fig. 2, the
GPS time ($t_{\textrm{trig}}$), the models active, the CBC waveform, the
wavelet window around the trigger time, and the wavelet bandwidth.

Figure 2: Results on the GW150914 data with different waveforms and analyses.
The left panel presents 90% credible intervals for the whitened signal
reconstructions in each detector while the right panel shows one- and two-
dimensional marginalized posterior distributions for selected CBC parameters:
mass ratio $q$, detector-frame chirp mass ${\cal{M}}$, and the effective spin
$\chi_{\text{eff}}$. In both panels red corresponds to an analysis with the
CBC model and PhenomD (CBCIMR), yellow corresponds again to a CBC model
analysis but with PhenomDNRT (CBCInsp), and black corresponds to an analysis
with the CBC+wavelet model with PhenomDNRT (CBCInsp+wavelet). In the left
panel the gray dashed line indicates the data whitened with a fair draw PSD
from the noise model posterior.

We perform three analyses with different models and CBC templates:

  1. 1.

A full “CBCIMR” analysis where the data are described only by the CBC model
with PhenomD.

  2. 2.

A “CBCInsp” analysis where the data are again described only by the CBC model
but with PhenomDNRT.

  3. 3.

A “CBCInsp+wavelets” analysis where the data are described by a combination of
the CBC model with PhenomDNRT and the wavelet model. We restrict the wavelet
model to times $t_{0}$ $\in$
$\left[t_{\textrm{trig}}-0.025\text{s},t_{\textrm{trig}}+0.025\text{s}\right]$
and frequencies above 150Hz, targeting the merger portion of the signal that
is missed by PhenomDNRT.

The run labels and relevant settings are given in Table 1.

In Fig. 2 we show the signal reconstructions and parameter posteriors for
these three analyses. The left panel shows the whitened data and the 90$\%$
credible intervals for each reconstruction in LIGO Hanford (top) and LIGO
Livingston (bottom). All analyses are able to identify the GW signal and they
further lead to consistent signal reconstructions of the early portion of the
signal. However, as the signal approaches the merger phase, the CBCInsp
analysis deviates from the reference CBCIMR analysis as well as the hybrid
CBCInsp+wavelets analysis by underpredicting the strength of the signal. The
hybrid CBCInsp+wavelets agrees well with the full CBCIMR, though the
uncertainty of the former is larger in the merger phase where the signal is no
longer modeled with CBC templates but with wavelets.

The recovered source parameters are given in the right panel of Fig. 2. The
CBCIMR analysis leads to consistent results with those reported in previous
studies [131, 132, 133], while the CBCInsp analysis shows a significant bias
away from the expected posteriors. The bias is most evident in the chirp mass:
PhenomDNRT compensates for the lack of a merger and ringdown by decreasing the
chirp mass, thus leading to a longer inspiral phase in an attempt to capture
part of the missing merger cycles. Despite this bias, the CBCInsp analysis is
still not able to capture the full signal as seen on the left panel. The
CBCInsp+wavelets proof-of-concept model results in recovered parameter
posteriors that are consistent with the full CBCIMR analysis. The posteriors
are not expected to be identical as now the CBC model has access only to the
inspiral portion of the signal and therefore lacks information available to
the CBCIMR analysis. However, the wavelet model can efficiently capture the
missing portion of the signal, allowing the CBC model to recover unbiased
system parameters from the inspiral phase only.

Figure 3: Breakdown of the various components of the hybrid CBCInsp+wavelets
analysis of GW150914. We present 90% credible intervals for the whitened
signal reconstruction of the CBC model (yellow), the wavelet model (red), and
the sum of the models (black). The sum is identical to the black lines from
Fig. 2. The gray dashed line indicate the data whitened with a fair draw PSD
from the noise model posterior.

Figure 3 examines the CBCInsp+wavelets analysis in more detail and the
complementary roles of the components in the hybrid model. We plot the signal
reconstruction from each submodel as well as their sum. The sum of the
individual model components is the same as the CBCInsp+wavelets model in Fig.
2. Figure 3 shows that as expected the BBH inspiral phase is primarily
reconstructed by the CBC PhenomDNRT model. As we move from the inspiral toward
the merger time, the CBC model starts tapering off and deviating from the full
GW signal. At the same time, the wavelet model captures the part of data that
is not covered by the CBC model, such that the combination of CBC and wavelets
accurately models the full GW signal. This toy model for GW150914 demonstrates
the main concept of our hybrid analysis: the full GW signal is modeled with a
well-understood component, i.e. the CBC waveform, and a model-agnostic
component, i.e. the wavelet model, that covers any parts of the signal that
are not included in the CBC model.

##  IV Constraints on the GW170817 postmerger

We turn our attention to GW170817, the first GW detection of a BNS coalescence
[5, 30, 134]. Both the low-frequency inspiral phase and a potential short-
duration high-frequency postmerger phase of the signal have been studied
separately in detail, with the latter remaining undetected [72, 30]. Here we
apply our hybrid CBC+wavelet model to study both parts of the GW170817 signal
simultaneously and compare our results with separate analyses of the inspiral
and postmerger signals.

Label |  | Segment
---
length (s)
| Sampling
---
rate (Hz)
$f_{\textrm{low}}$(Hz) | Model | Waveform |  | Wavelet
---
window (s)
| Wavelet
---
bandwidth (Hz)
$D_{\mathrm{min}}$
CBCInsp | 64 | 4096 | 32 | CBC | PhenomDNRT | N/A | N/A | N/A
CBCInsp+wavelet | 64 | 8196 | 32 | CBC/wavelet | PhenomDNRT | (-0.125, 0.125) | (1024, 4096) | 2
wavelet-only | 64 | 8196 | 1024 | wavelet | N/A | (-0.125, 0.125) | (1024, 4096) | 2
Table 2:  Settings for the runs of Sec. IV on GW170817. All runs are relative
to GPS time 1187008882.446. From left to right, columns correspond to the run
label in Fig. 4 and Fig. 5, the segment length, the sampling rate, the low
frequency cut off, the models active, the CBC waveform, the wavelet window
around the trigger time, the wavelet bandwidth, and the minimum number of
wavelets.

We analyze 64s of data from LIGO Hanford and LIGO Livingston around GPS time
1187008882.446 [129, 130] where the prominent glitch in LIGO Livingston has
already been subtracted [135, 94]. While Virgo data are available for that
time and aided in constraining the source sky location [30], we do not
consider them due to the lower sensitivity and the fact that we fix the sky
location to the known values of $\alpha=3.446$ rad and $\delta=-0.408$ rad
[136, 137]. We also use a fixed PSD rather than marginalize over the noise
model for computational efficiency. We perform three analyses with different
models and data:

  1. 1.

A “CBCInsp” analysis with data in the frequency range (32,2048)Hz that are
described only by the CBC model with the PhenomDRT waveform thus focusing on
the inspiral signal.

  2. 2.

A “CBCInsp+wavelet” analysis with data in the wider frequency range
(32,4096)Hz that are described by a combination of the CBC model with
PhenomDNRT and the wavelet model, thus targeting the full signal. The wavelets
are restricted to (1024,4096)Hz targeting a potential high frequency
postmerger signal.

  3. 3.

A “wavelet-only” analysis where we restrict the frequency range to
(1024Hz,4096)Hz and use only the wavelet model thus focusing only on the
postmerger signal.

For both analyses with the wavelets model, we use a prior on the number of
wavelets of $D\in[2,100]$. The reason for selecting a minimum number of $2$
wavelets is different for each analysis. In the case of the hybrid analysis,
the CBC part of the model does not terminate at merger, but smoothly tapers
into the postmerger phase for a few milliseconds, see Fig 9. We therefore need
at least two wavelets, one to undo this effect of the CBC model and the other
to capture the true postmerger signal. In the wavelet-only study we again need
at least two wavelets such that one wavelet can capture the merger itself
which extends into the analysis bandwidth and the other wavelet can capture
the contribution from the postmerger signal. The full settings of our three
analyses are detailed in Table 2.

Figure 4: Results on the GW170817 data with different analyses. The left panel
presents the median, 50% and 90% credible intervals for the whitened signal
reconstructions in each detector while the right panel shows one- and two-
dimensional marginalized posterior distributions for selected CBC parameters:
mass ratio $q$, detector frame chirp mass ${\cal{M}}$, the effective spin
$\chi_{\text{eff}}$ and the tidal deformation parameter $\tilde{\Lambda}$. In
both panels blue corresponds to the hybrid CBCInsp+wavelet and magenta
corresponds to the CBCInsp analysis. In the left panel the gray dashed gives
the whitened data.

We compare the two analyses that include the inspiral of the signal in Fig. 4
which shows the signal reconstructions focused on the end of the signal (left
panel) and the marginalized posterior distribution for selected source
parameters (right) of GW170817. Despite the high SNR of the signal, the
amplitude is relatively weak compared to the noise level, however both
analyses identify the signal and lead to consistent reconstructions for the
late inspiral phase as shown on the left panel. As the binary approaches
merger the reconstruction uncertainty increases due to the decreasing detector
sensitivity at increasing frequencies. This behavior is evident in both
analyses, suggesting that the GW170817 postmerger signal remains undetected as
expected from the detector sensitivity at the time of GW170817 [30].

The right panel of Fig. 4 examines selected CBC parameters recovered from the
inspiral signal. The recovered parameters are consistent with previous results
using the same CBC model [30]. The only difference is that our
$\tilde{\Lambda}$ posterior has more support for the higher of the two modes
present in the results of [30]. This can be attributed to different uses of
prior: we use of a prior that is flat in
$\tilde{\Lambda}-\delta\tilde{\Lambda}$ with the constraint that
$\Lambda_{1}>0$ and $\Lambda_{2}>0$. This choice results in a marginalized
prior for $\tilde{\Lambda}$ that favors larger values. In contrast, [30] uses
a flat marginalized prior for $\tilde{\Lambda}$. We have verified that
reweighting our posterior to a flat marginalized prior gives consistent
results with [30]. Furthermore, the hybrid CBCInsp+wavelet analysis leads to
essentially identical results for the CBC parameters as the traditional
CBCInsp analysis. This is consistent with [74] that finds that failure to
account for the postmerger signal will not lead to biases in the CBC
parameters for typical detector sensitivities, but also serves as a sanity
check of the CBCInsp+wavelet analysis and the fact that sampling for the joint
CBC+wavelet model has converged.

After confirming that the CBCInsp+wavelet analysis matches the CBCInsp
analysis as far as the inspiral portion of the signal is concerned, we switch
to comparing the two analyses that include a wavelet model for a possible
high-frequency postmerger signal. In Fig. 5 we show the reconstructed GW
spectrum for CBCInsp+wavelet and wavelet-only. The vertical band gives the 90%
credible interval for the merger frequency estimated using Eq. 11 in [138]
given the binary component masses and tidal parameters from the CBC model
posteriors. This merger frequency is also the location at which the PhenomDNRT
template starts tapering off, so it is related to the transition between the
CBC and the wavelet models.

Figure 5: The median, 50$\%$ and 90$\%$ credible intervals of reconstructed GW
spectrum in the LIGO Hanford detector in the CBCInsp+wavelet (blue) and the
wavelet-only analysis (red). The noise amplitude spectral density is also
overplotted (black dashed line). The gray band indicates the 90% credible
interval of the merger frequency, with the dashed line indicating the peak of
the merger frequency posterior at $f_{\text{merger}}$ $\sim$ 1778Hz.

At frequencies below $\sim 1000$Hz the CBCInsp+wavelet reconstruction is
dominated by the inspiral signal which is clearly detected as the spectrum
lower limit is nonzero. With increasing frequency the detector sensitivity
decreases, resulting in larger uncertainties in the reconstructed spectrum.
Finally, at around the merger frequency, $\gtrsim 1600$Hz, the reconstruction
uncertainty is large and consistent with no detected signal. A similar picture
is drawn by the wavelet-only analysis which again places only upper limits on
the signal in frequencies above $1024$Hz. The upper limit from the
CBCInsp+wavelet and the wavelet-only analyses is comparable, though the former
is consistently lower across the frequency band, indicating more stringent
constraints on the presence of a postmerger signal. Additionally, the wavelet-
only analysis does not lead to a detection of the merger signal in the
frequency range (1000,1500)Hz, unlike the CBCInsp+wavelet case. This suggests
that this high-frequency portion of the signal is not individually detectable,
but only inferred coherently from the preceding inspiral signal.

##  V Simulated signals

Going beyond the GW170817 upper limits, in this section we study simulated
BNSs where the postmerger emission is detectable by future GW detectors. We
assume a network of two detectors at the current location of LIGO Hanford and
LIGO Livingston and a zero noise realization. Since all signals have high SNRs
we assume that an electromagnetic counterpart has been identified and the sky
location of the source is known. Following [31] we do not explicitly select
any planned GW detector such as Cosmic Explorer [139, 140, 141] or the
Einstein Telescope [142, 143] and their nominal sensitivity. We instead work
with the Advanced LIGO design sensitivity [144] and gradually lower the noise
PSD, emulating improving detector sensitivity and higher signal SNR.

Since no NR simulation of the full BNS signal as observed by ground-based
detectors exists, all simulated signals are constructed in a hybrid fashion:
the postmerger signal is obtained through NR simulations, while the premerger
signal is computed with PhenomDRT using the same parameters as the NR
simulation. We choose extrinsic parameters consistent with GW170817, namely
distance $D_{L}=40$Mpc, inclination $\iota=2.635$, polarization angle
$\psi=0$, and the known sky location. All signals are analyzed with the full
CBC+wavelets model with run settings equivalent to the CBCInsp+wavelet
analysis from Table 2.

###  V.1 Waveform from Kastaun and Ohme

The first full BNS waveform we consider was constructed and released by [145].
This hybrid waveform was constructed by combining the PhenomDNRT model for the
inspiral and the results of NR simulations for the postmerger GW signal. The
hydrodynamical simulations assume a GW170817-like initial system with the
hadronic SFHO EoS [146], resulting in a short-lived hypermassive NS. We use
the resulting hybrid waveform for a BNS system with mass ratio $q=0.9$,
detector frame component masses of $m_{1}$ = 1.438$M_{\odot}$ and $m_{2}$ =
1.294$M_{\odot}$, and dimensionless tidal parameters $\Lambda_{1}$ = 280 and
$\Lambda_{2}$ = 551, giving a binary tidal parameter $\tilde{\Lambda}$ = 396
and merger frequency $f_{\mathrm{merger}}=2110$ Hz.

We inject the signal in a detector network with sensitivities 2xDS, 4xDS, 6xDS
which respectively denote 2, 4, and 6 time improved strain sensitivity
compared to the LIGO design sensitivity (i.e., the design sensitivity divided
by 2, 4 and 6). In terms of amplitude spectral density, they correspond to
8.23$\times 10^{-24}~{}$Hz-1/2, 4.12$\times 10^{-24}~{}$Hz-1/2, and
2.72$\times 10^{-24}~{}$Hz-1/2 at 3326 Hz, which is the frequency of the main
postmerger mode of the signal, $f_{\mathrm{peak}}$. The resulting network SNRs
are 142, 384, and 425 for the premerger signal ($f<f_{\mathrm{merger}}$), and
2.8, 5.6, and 8.4 for the postmerger signal ($f\geq f_{\mathrm{merger}}$).

Figure 6: Whitened-time domain reconstruction of the simulated signal from
[145] at 6xDS. The gray dashed line gives the simulated data. The shaded
regions give the 50% and 90% credible intervals for the full CBC+wavelets
reconstruction (black), as well as the CBC (yellow) and the wavelets (red)
component of the full analysis. The postmerger signal is reconstructed by the
wavelets model. Figure 7: Median, 50%, and 90% credible intervals for the
reconstructed spectrum from the Kastaun and Ohme [145] waveform with different
detector sensitivities; from left to right we assume a strain sensitivity two,
four, and six times better than the advanced LIGO design sensitivity. The gray
line gives the simulated data, which are the same in all cases, while the
black dashed line gives the detection noise amplitude spectral density. The
dashed and dotted vertical lines denote, respectively, the dominant and
subdominant postmerger spectrum frequency.

Figure 6 shows the whitened time-domain data and reconstruction for the
injection at 6xDS, again focusing around the late stages of the signal. The
full reconstruction with the combined CBC+wavelets model accurately captures
the entire signal. We also plot the individual components of the full model,
namely the CBC and the wavelets model separately. As expected, the CBC model
captures everything up to merger and then tapers off. The wavelets overlap
with the CBC model in the taper region around merger, effectively canceling
out the ringdown-like oscillations of the CBC model that do not match the
data. At later times, the wavelet model extends to the postmerger part of the
signal, capturing its main oscillatory component.

Figure 7 shows the reconstructed spectrum from injections on different
detectors sensitivities. As the detector sensitivity increases, the signal
reconstruction becomes more accurate. The premerger signal is recovered in all
analyses given its strength, and increasing sensitivity reduces the
reconstruction uncertainty. The postmerger signal, on the other hand, is too
weak to be detected in the 2xDS case, and therefore the reconstruction is
essentially uninformative and similar to the GW170817 one. At increasing
detector sensitivity, the postmerger signal emerges from the noise. At 4xDS
the reconstructed spectrum starts identifying the main postmerger peak though
only at the $\sim 50\%$ credible level. At 6xDS the reconstructed spectrum not
only confidently identifies the main postmerger peak but also has evidence of
a subdominant peak at around 2337 Hz.

Figure 8: One- and two-dimensional marginalized posterior distributions for
selected signal parameters from the analyses of Fig. 7: detector frame chirp
mass ${\cal{M}}$, detector frame total mass $M$, tidal deformation parameter
$\tilde{\Lambda}$, and postmerger peak frequency $f_{\mathrm{peak}}$. Vertical
black lines or crosses denote the true values of each parameter.

Posteriors for select premerger and postmerger parameters are given in Fig. 8.
The peak frequency of the postmerger signal $f_{\mathrm{peak}}$ is defined in
the same way as [85, 31]: for each reconstruction posterior sample we compute
the frequency at the maximum of the spectrum after the merger. If the spectrum
possesses no maximum, then a sample is drawn from the $f_{\mathrm{peak}}$
prior. All posteriors are consistent with the injected values, and
uncertainties decrease with increasing detector sensitivity as expected. The
$f_{\mathrm{peak}}$ posterior at 2xDS is essentially the prior, consistent
with the fact that the postmerger signal was not detected. At 4xDS the
$f_{\mathrm{peak}}$ posterior starts exhibiting a peak at the injected value,
consistent with the partial identification of the postmerger signal in the
middle panel of Fig. 7. At 6xDS the peak frequency is accurately measured.

###  V.2 Hadronic EoS

We construct further simulated signals based on the postmerger NR simulations
from [31], which employ the conformal flatness approximation [147, 148].
Comparisons to fully relativistic studies show a very good agreement with
regards to postmerger GW frequencies but an underestimation of the postmerger
GW amplitude by some 10% as a result of using the quadrupole formula for GW
extraction [47]. Since our analysis does not rely on calibration to any NR
simulations, we expect our results below to be unaffected by such
uncertainties. Indeed, BayesWave has been shown to produce reliable results on
simulated postmerger signals made with different NR codes [81, 85].

All simulations in [31] were at the time constructed to be consistent with
GW170817, though in light of new data from NICER [14, 13, 15, 16] some of the
softest EoSs there are now disfavored. We work with EoS3 from [31], a hadronic
EoS that is consistent with radii values inferred in [26] and corresponds to
$f_{\mathrm{peak}}=2880$Hz and $R_{1.4}=12.6$ km. The inspiral portion of the
signal is again described by PhenomDNRT for a GW170817-like system with
detector frame component masses $m_{1}=m_{2}$ = 1.362$M_{\odot}$,
dimensionless tidal parameters $\Lambda_{1}=\Lambda_{2}=\tilde{\Lambda}=587$,
and zero spin.

The full waveform is constructed by aligning the projected premerger and
postmerger waveforms in an overlap interval in the time domain. We use a
transition window $[t_{1},t_{2}]$ where the full waveform transitions from the
PhenomDNRT template to the NR simulation data. The full waveform is then

$\displaystyle
h(t)=\left[1-x(t)\right]h_{\text{inspiral}}(t)+x(t)A_{\textrm{scale}}h_{\text{NR}}(t),$
(22)

where $h_{\text{inspiral}}$ denotes the PhenomDNRT inspiral model,
$h_{\text{NR}}$ denotes the postmerger NR simulation, $A_{\textrm{scale}}$ is
a scale parameter that can be varied to control the strength of the postmerger
signal, and $x$ ensures a smooth transition through a Planck taper function

$\displaystyle x(t)=\begin{cases}0&t\leq t_{1}\\\
\left[1+\exp{\left(\frac{t_{2}-t}{t-t_{1}}+\frac{t_{2}-t}{t-t_{2}}\right)}\right]^{-1}&t_{1}<t<t_{2}\\\
1&t\geq t_{2}.\end{cases}$ (23)

Figure 9: Comparison of the premerger (PhenomDNRT) and postmerger (NR
simulation) waveform that are merged into a hybrid full waveform. The yellow
band indicates the overlap region where the hybrid waveform transitions from
being composed of the PhenomDNRT model to the postmerger NR simulation model.
The maximum amplitude of the NR waveform used here is scaled to be 0.9 times
the premerger peak amplitude.

We fix the detector sensitivity to 2xDS and vary the amplitude of the
postmerger signal through the scale parameter $A_{\textrm{scale}}$. The main
motivation for this is that it allows us to vary the postmerger signal
strength, while keeping the premerger SNR (and thus the computation cost)
manageable. This scaling further allows us to address the underestimation of
the signal amplitude within the conformal flatness approximation mentioned
above. Larger sets of simulations do not find a very tight relation between
the postmerger amplitude and binary parameters [79], and simulations in
general may over- or under-estimate the postmerger amplitude to some extent.

In what follows, we keep the premerger SNR constant at 112, while the
postmerger SNR varies and results are presented as a function of the ratio of
the premerger to the postmerger peak amplitudes. We explore three values for
the ratio of the peak postmerger to the peak premerger amplitudes: 0.7, 0.9
and 1.1 (corresponding to $A_{\textrm{scale}}$ = 1.5, 2.0 and 2.5 and
postmerger SNRs of 6.2, 8.1 and 10.0, respectively). We show an example of the
above construction process in Fig. 9.

Figure 10: Median, 50$\%$ and 90$\%$ credible intervals of the reconstructed
spectrum for the hadronic EoS3 signal from Sec. V.2 and for pre/postmerger
peak amplitude ratios of 0.7 (left), 0.9 (middle), 1.1 (right) in the LIGO
Hanford detector. The gray line gives the simulated data. The noise amplitude
spectral density is also overplotted with a black dashed line. Figure 11: One-
and two-dimensional marginalized posterior distributions for selected pre and
postmerger parameters for the analyses of Fig. 10. The parameters shown are
the source-frame chirp mass ${\cal{M}}$ and total mass $M$, the tidal
deformation parameter $\tilde{\Lambda}$ and the postmerger peak frequency
$f_{\mathrm{peak}}$. Vertical black lines or crosses denote the true values of
each parameter. The gray region is the expected value for $f_{\mathrm{peak}}$
given the premerger inferred parameters and the fit of [85] that assumes
hadronic EoSs. The premerger and postmerger results are consistent with
expectations for hadronic EoSs.

In Figure 10 we compare the reconstructed spectra from our analyses for the
three cases of post/premerger amplitude scaling. The inspiral portion of the
signal and the corresponding reconstruction are similar in the three panels.
As the amplitude of the postmerger signal increases from left to right the
reconstructed spectrum includes more detailed features, progressing from hints
of a postmerger peak on the left panel to increasingly more confident
identification in the middle and right panel. Figure 11 again shows select
recovered premerger and postmerger parameters. The gray region overlapping
with the $f_{\mathrm{peak}}$ posterior is the expected value for the
postmerger frequency given the premerger signal. Specifically we use the fit
of [85] that holds for hadronic EoSs and the binary mass and tidal
deformability as extracted from the premerger signal to compute the expected
$f_{\mathrm{peak}}$. The expected and recovered posterior values for
$f_{\mathrm{peak}}$ agree, showing that our analysis can correctly conclude
that the premerger and postmerger signals are consistent with each other given
expectations from NR simulations of hadronic EoSs.

###  V.3 EoS with phase transitions

Figure 12: Mass-radius relation for EoS DD2-SF-4 used in our study of PTs
(black line). The green shaded region gives the 90% credible level of the
mass-radius posterior derived in Legred et al. [26] using data from heavy
pulsars, GWs, and the Miller et al. [13, 15] radius results using X-ray data.
The adopted EoS is inconsistent with this posterior at the $90\%$ level for
masses above $\sim 1.8M_{\odot}$, but it is consistent to the same level with
the less constraining Riley et al. [14, 16] results.

Since the hypermassive NS that gives rise to the postmerger signal is
characterized by higher core densities and temperatures than the premerger
NSs, the postmerger signal has the potential to reveal new high-density
physics. One possibility is a strong phase transition in the EoS toward quark
degrees of freedom333We use the term “strong” phase transition to emphasize
that the exact characteristics of the hadron-quark phase transition are
currently unclear, and only a sufficiently strong transition, for instance a
first-order phase transition with large latent heat, may impact the merger
dynamics such that the occurrence of quark matter leads to an unambiguous
signature. See, e.g., [58] for a detailed discussion. In fact, quark matter
may also resemble the behavior of purely hadronic matter (commonly referred to
as the masquerade problem [149, 150]), which may not leave a very prominent
imprint on the GW signal.. Such phase transitions could occur at lower
densities and thus be detectable with premerger data only [151, 152, 153, 154,
155, 156, 157, 158, 159, 160]. However, as already mentioned above, the tidal
deformability is a steeply decreasing function of the NS compactness, making
premerger data less constraining about high-mass NSs that could contain quark
cores. postmerger data, on the other hand, can be used to probe higher
densities and NR simulation-based studies have explored the potential
signature of a high-density phase transition on the GW signal [62, 63, 64, 65,
66, 67, 68].

Figure 13: The median, 50$\%$ and 90$\%$ credible intervals of the
reconstructed spectrum for a signal with EoS3 where the inspiral signal has
been simulated with a stiffer EoS in order to mimic the phenomenology of a
strong phase transitions from Sec. V.3. We show pre/postmerger peak amplitude
ratios of 0.7 (left), 0.9 (middle), 1.1 (right). The gray line gives the
simulated data. The noise amplitude spectral density is also overplotted with
a black dashed line.  Figure 14: One- and two-dimensional marginalized
posterior distributions for selected pre and postmerger parameters for the
analyses of Fig. 13. The parameters shown are the source-frame chirp mass
${\cal{M}}$ and total mass $M$, the tidal deformation parameter
$\tilde{\Lambda}$ and the postmerger peak frequency $f_{\mathrm{peak}}$. Black
vertical black lines or crosses denote the true values of each parameter. The
gray region is the expected value for $f_{\mathrm{peak}}$ given the premerger
inferred parameters and the fit of [85] that assumes hadronic EoSs. The
premerger and postmerger results are now inconsistent with expectations for
hadronic EoSs.

The most prominent signature of a sufficiently strong phase transition would
be an increase in the postmerger peak frequency as the softening of the EoS
would lead to a more compact merger remnant [63]. It has been proposed that
such a frequency increase could be identified if one compares the premerger
and postmerger data: for hadronic EoSs the tidal deformability and peak
frequency follow the approximately EoS-insensitive relation shown in Fig. 17.
This relation could be violated for EoSs with high-density phase transitions
as the tidal deformability is determined by the hadronic part of the EoS
alone, while the postmerger frequency is affected by the phase transition [63,
66]444In fact, a postmerger frequency deviation may also occur if quark matter
is present before the merger [66].. The degree of deviation from the EoS-
insensitive relation, and thus how detectable it is, depends on the strength
of the transition [63, 58].

As a first example of such a signal, we work again with EoS3 from [31], but
the premerger data are constructed with tidal parameters that are
systematically shifted compared to their hadronic EoS value. This effectively
results in a stiffer hadronic EoS that undergoes a phase transition which
softens the postmerger signature to the level of EoS3. All intrinsic
parameters remain to be the same as the simulations of Sec. V.2 with the
exception of $\tilde{\Lambda}=\Lambda_{1}=\Lambda_{2}$ = 800.The postmerger
SNRs remain the same as Sec. V.2.

Reconstructed spectra and parameter posteriors are shown in Figs. 13 and 14
respectively. We obtain qualitatively similar results to Figs. 10 and 11 with
the main difference being that now the premerger and postmerger are now
inconsistent as expected. Figure 14 shows the expected $f_{\mathrm{peak}}$
value given the premerger constraints on the binary properties under the
assumption of a hadronic EoS. The recovered $f_{\mathrm{peak}}$ is
inconsistent with this expectation to within its measurement uncertainty,
signaling the presence of additional high-density effects in the EoS that
affect the postmerger signal.

Figure 15: The median, 50$\%$ and 90$\%$ credible intervals of the
reconstructed spectrum for a signal with the DD2-SF-4 EoS with phase
transitions from Sec. V.3 and for pre/postmerger peak amplitude ratios of 1.1
(left), 1.3 (middle), 1.5 (right). Due to the increased $f_{\mathrm{peak}}$
value toward the less sensitive detector frequency range, such a postmerger
signal would require a larger strain amplitude for detection, we therefore
display results with larger values of the pre/postmerger peak amplitude ratio
compared to Fig. 13. The corresponding postmerger SNR is 7.2, 8.1, 9.1 from
left to right. The gray line gives the simulated data. The noise amplitude
spectral density is also overplotted with a black dashed line. Figure 16: One-
and two-dimensional marginalized posterior distributions for selected pre and
postmerger parameters for the analyses of Fig. 15. The parameters shown are
the detector frame chirp mass ${\cal{M}}$ and total mass $M$, the tidal
deformation parameter $\tilde{\Lambda}$ and the postmerger peak frequency
$f_{\mathrm{peak}}$. Vertical black lines or crosses denote the true values of
each parameter. The gray region is the expected value for $f_{\mathrm{peak}}$
given the premerger inferred parameters and the fit of [85] that assumes
hadronic EoSs. The premerger and postmerger data are now as expected
inconsistent with expectations for hadronic EoSs.

We further assess how well our hybrid analysis could detect a strong phase
transition by simulating a signal with the DD2-SF-4 EoS [161, 162] and
corresponding postmerger simulation from [63], also performed with the
conformal flatness approximation. The mass-radius relation for this EoS is
given in Fig. 12 and it exhibits the characteristic radius reduction due to
strong phase transitions starting at $\sim 1.5M_{\odot}$. Given the simulated
binary masses of $1.35M_{\odot}$ for both components, the inspiral signal is
emitted by hadronic NSs, while the postmerger signal is affected by the onset
of the phase transition [63]. For reference, we also show the 90% symmetric
credible intervals of the mass-radius posterior for the EoS derived in [26]
using heavy pulsar, GW, and X-ray data. The DD2-SF-4 EoS is inconsistent with
the posterior at the 90% level for large masses as it underpredicts the radius
of a $2M_{\odot}$ NS, though it is consistent with current data at the 95%
level. If we instead used the less constraining radius results from Riley el
at. [14, 16], the EoS would be consistent with the posterior at the 90% level.
Results are presented in Figs. 15 and 16, where again we find that the
postmerger signal can be reconstructed for sufficiently loud signals, and the
inconsistency between $\tilde{\Lambda}$ and $f_{\mathrm{peak}}$ can be
identified.

Figure 17: Comparison of the recovered $\tilde{\Lambda}$-fpeak posterior
distribution from the EoS3 hadronic and PT data as well as the data with the
DD2-SF-4 EoS. The gray bands corresponds to the EoS-independent fit relating
the expected values of $\tilde{\Lambda}$-fpeak for hadronic NSs. We show 50%
(dashed) and 90% (solid) contours, while colors correspond to those of Figs.
11, 14, 16. In all cases, we find the correct agreement or disagreement with
the expected hadronic relation for sufficiently loud signals.

Finally, we elaborate on the premerger and postmerger consistency in Fig. 17.
The gray band shows the expected relation between $\tilde{\Lambda}$ and
$f_{\mathrm{peak}}$ assuming hadronic EoSs. We compute this relation by
converting the $R-f_{\mathrm{peak}}$ EoS-independent555The fit is shown in
Fig. 4 of [85] and uses EoSs that predict radii between $11.5$km and $5$km
approximately. fit by [85] to a $\tilde{\Lambda}-f_{\mathrm{peak}}$ relation.
We convert $R$ to $\Lambda$ using the scaling relation calibrated for
GW170817-like systems for our source mass parameters $M_{\text{tot}}$ and
${\cal{M}}$ given by Eq. 17 in [118]. Overplotted are the two-dimensional
posteriors for $\tilde{\Lambda}-f_{\mathrm{peak}}$ for the case of a hadronic
EoS (EoS3), the phase-transition version of the EoS3-based signal, and the
DD2-SF-4 EoS again with a strong phase transition. The injected values for the
hadronic case is consistent with expectations and this is confirmed by the
recovered values to within the statistical error. On the other hand, a strong
phase transition leads to a $\tilde{\Lambda}-f_{\mathrm{peak}}$ combination
that is inconsistent with hadronic expectations and the extracted posteriors
are able to identify this behavior.

##  VI Conclusions

We present a hybrid approach to study the full GW signal emitted during a BNS
coalescence, including a possible postmerger component. Our method models the
inspiral part of the coalescence with waveform templates as implemented in
LALSimulation and the postmerger part of the signal with a superposition of
wavelets. We do not impose phase coherence in the transition between the
template and the wavelets, however, the full signal is smooth and coherent as
the template model used already extends coherently past merger and into the
postmerger part of the signal.

Applying our method to GW170817, we demonstrate that the inspiral parameters
are consistent with previous results. We do not detect a postmerger signal for
GW170817 but find that the high-frequency portion of the signal between
1000-1500 Hz is only detected using the full hybrid analysis in contrast with
the traditional postmerger-only analyses. We apply our method to simulated
signals with a detectable postmerger component and show that our full analysis
simultaneously reconstructs both the inspiral tidal deformation and the
postmerger dominant frequency peak when the SNR is sufficiently high.

In this work, we sample the inspiral tidal parameters and postmerger peak
frequency independently, imposing no relation between them and no assumption
on the nature of the EoS. This allows us to detect a possible signature of the
hadron-quark phase transition in NS mergers, as it lead to a characteristic
frequency shift of the postmerger signal. A possible extension of this
analysis could use information extracted from a premerger signal to guide the
analysis of the postmerger signal. If one restricts to a hadronic EoS,
relations between the premerger tidal parameters and postmerger peak frequency
could be used to predict the approximate location of the dominant postmerger
peak from the inspiral tidal parameters. This could improve the prospects of
detecting the postmerger signal in weaker signals, at the expense of assuming
the EoS is similar to a hadronic one. A further possible improvement concerns
the use of “chirplets”, sine-Gaussian wavelets with an evolving frequency
[95], which might be better suited for a time-evolving postmerger frequency
mode.

Going beyond the main peak, numerical simulations of postmerger signals show
subdominant peaks in the GW spectrum, which can also be used to characterize
the properties of the merger remnant star [163, 49, 53]. Our analysis can
identify some of these peaks for sufficiently loud signals, however they are
typically less loud than the main peak and can be lost in the noise. Numerical
simulations suggest that the main and the secondary peaks in the spectrum are
related to each other [53, 84, 82]. In future work, this feature can be
utilized in the form of a prior that links the two peaks and is parametrized
in terms of the remnant compactness. Such a prior would enhance the
sensitivity of the analysis to secondary peaks as it effectively adds modeling
information about the signal similar to the approach of analytic models [76,
77, 78, 79, 80, 81, 82].

Our analysis concerns premerger SNRs in the hundreds, where systematic biases
in the waveform models could be important and would have to be mitigated in
the lead-up to such improved detectors [74, 164, 165, 166, 167, 168]. However,
as we make minimal assumptions about the postmerger waveform and use sine-
gaussian wavelets, we expect the reconstructed signal to be less susceptible
to systematics than analytic models [76, 77, 78, 79, 80, 81, 82] that are
limited by the accuracy of the NR simulations on which they are based (see for
example Fig. 4 of [80]). Our analysis instead is not limited by NR accuracy
and has yielded unbiased results on simulated signals form different codes
[81, 85]. On the other hand, analytic models can be more straightforwardly
attached to premerger waveform templates and enforce phase coherence through
merger [80].

Next generation detectors are expected to detect thousands of BNSs [169, 170,
171], the majority of which will be weak with individually undetectable
postmerger signals. While we study single loud sources here, our method can
also be applied to the expected numerous weaker BNS sources which might
dominate the overall constraints when combined [172]. Though our analysis does
not hinge on the existence of accurate models for the postmerger signal, it is
possible that biases in the premerger waveform will lead to biased inferences
about the postmerger when making use of EoS-independent relations connecting
them. The extremely sensitive observations possible with next-generation
detectors indeed require control over a wide range of potential systematic
biases, and the flexible analysis presented here helps mitigate such biases
from the postmerger signal.

###### Acknowledgements.

We thank Will Farr and Wynn Ho for many useful discussions. We also thank
Sophie Hourihane and Tyson Littenberg for discussions and assistance about
BayesWave. This research has made use of data, software and/or web tools
obtained from the Gravitational Wave Open Science Center (https://www.gw-
openscience.org), a service of LIGO Laboratory, the LIGO Scientific
Collaboration and the Virgo Collaboration. Virgo is funded by the French
Centre National de Recherche Scientifique (CNRS), the Italian Istituto
Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with
contributions by Polish and Hungarian institutes. This material is based upon
work supported by NSF’s LIGO Laboratory which is a major facility fully funded
by the National Science Foundation. The authors are grateful for computational
resources provided by the LIGO Laboratory and supported by National Science
Foundation Grants PHY-0757058 and PHY-0823459. M.W. gratefully acknowledges
support and hospitality from the Simons Foundation through the predoctoral
program at the Center for Computational Astrophysics, Flatiron Institute. The
Flatiron Institute is supported by the Simons Foundation. K.C. was supported
by NSF Grant PHY-2110111. N.J.C. was supported by NSF Grant PHY-1912053. AB
acknowledges support by the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme under grant agreement
No. 759253, by Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) - Project-ID 279384907 - SFB 1245, by DFG - Project-ID 138713538 -
SFB 881 (“The Milky Way System”, subproject A10) and by the State of Hesse
within the Cluster Project ELEMENTS. Software: gwpy [173], matplotlib [174].

## References

  * Lattimer and Prakash [2016] J. M. Lattimer and M. Prakash, The Equation of State of Hot, Dense Matter and Neutron Stars, Phys. Rept. 621, 127 (2016), arXiv:1512.07820 [astro-ph.SR] .
  * Özel and Freire [2016] F. Özel and P. Freire, Masses, Radii, and the Equation of State of Neutron Stars, Ann. Rev. Astron. Astrophys. 54, 401 (2016), arXiv:1603.02698 [astro-ph.HE] .
  * Oertel  _et al._ [2017] M. Oertel, M. Hempel, T. Klahn, and S. Typel, Equations of state for supernovae and compact stars, Rev. Mod. Phys. 89, 015007 (2017), arXiv:1610.03361 [astro-ph.HE] .
  * Baym  _et al._ [2018] G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka, From hadrons to quarks in neutron stars: a review, Rept. Prog. Phys. 81, 056902 (2018), arXiv:1707.04966 [astro-ph.HE] .
  * Abbott  _et al._ [2017a] B. P. Abbott _et al._ (LIGO Scientific Collaboration, Virgo Collaboration), GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2020] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GW190425: Observation of a Compact Binary Coalescence with Total Mass $\sim 3.4M_{\odot}$, Astrophys. J. Lett. 892, L3 (2020), arXiv:2001.01761 [astro-ph.HE] .
  * Aasi  _et al._ [2015] J. Aasi _et al._ (LIGO Scientific), Advanced LIGO, Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ [2015] F. Acernese _et al._ (Virgo Collaboration), Advanced Virgo: a second-generation interferometric gravitational wave detector, Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Abbott  _et al._ [2013] B. P. Abbott _et al._ (LIGO Scientific Collaboration, Virgo Collaboration), Prospects for Observing and Localizing Gravitational-Wave Transients with Advanced LIGO, Advanced Virgo and KAGRA, Living Rev. Rel. 19, 1 (2013), arXiv:1304.0670 [gr-qc] .
  * Antoniadis  _et al._ [2013] J. Antoniadis, P. C. Freire, N. Wex, T. M. Tauris, R. S. Lynch, _et al._ , A Massive Pulsar in a Compact Relativistic Binary, Science 340, 1233232 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Cromartie  _et al._ [2019] H. T. Cromartie _et al._ , Relativistic Shapiro delay measurements of an extremely massive millisecond pulsar, Nature Astron. 4, 72 (2019), arXiv:1904.06759 .
  * Fonseca  _et al._ [2021] E. Fonseca _et al._ , Refined Mass and Geometric Measurements of the High-mass PSR J0740+6620, Astrophys. J. Lett. 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Miller  _et al._ [2019] M. C. Miller _et al._ , PSR J0030+0451 Mass and Radius from $NICER$ Data and Implications for the Properties of Neutron Star Matter, Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Riley  _et al._ [2019] T. E. Riley _et al._ , A $NICER$ View of PSR J0030+0451: Millisecond Pulsar Parameter Estimation, Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Miller  _et al._ [2021] M. C. Miller _et al._ , The Radius of PSR J0740+6620 from NICER and XMM-Newton Data, (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Riley  _et al._ [2021] T. E. Riley _et al._ , A NICER View of the Massive Pulsar PSR J0740+6620 Informed by Radio Timing and XMM-Newton Spectroscopy, (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Raaijmakers  _et al._ [2020] G. Raaijmakers _et al._ , Constraining the dense matter equation of state with joint analysis of NICER and LIGO/Virgo measurements, Astrophys. J. Lett. 893, L21 (2020), arXiv:1912.11031 [astro-ph.HE] .
  * Dietrich  _et al._ [2020] T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Multimessenger constraints on the neutron-star equation of state and the Hubble constant, Science 370, 1450 (2020), arXiv:2002.11355 [astro-ph.HE] .
  * Landry  _et al._ [2020] P. Landry, R. Essick, and K. Chatziioannou, Nonparametric constraints on neutron star matter with existing and upcoming gravitational wave and pulsar observations, Phys. Rev. D 101, 123007 (2020), arXiv:2003.04880 [astro-ph.HE] .
  * Al-Mamun  _et al._ [2021] M. Al-Mamun, A. W. Steiner, J. Nättilä, J. Lange, R. O’Shaughnessy, I. Tews, S. Gandolfi, C. Heinke, and S. Han, Combining Electromagnetic and Gravitational-Wave Constraints on Neutron-Star Masses and Radii, Phys. Rev. Lett. 126, 061101 (2021), arXiv:2008.12817 [astro-ph.HE] .
  * Reed  _et al._ [2021] B. T. Reed, F. J. Fattoyev, C. J. Horowitz, and J. Piekarewicz, Implications of PREX-2 on the Equation of State of Neutron-Rich Matter, Phys. Rev. Lett. 126, 172503 (2021), arXiv:2101.03193 [nucl-th] .
  * Essick  _et al._ [2020] R. Essick, I. Tews, P. Landry, S. Reddy, and D. E. Holz, Direct Astrophysical Tests of Chiral Effective Field Theory at Supranuclear Densities, Phys. Rev. C 102, 055803 (2020), arXiv:2004.07744 [astro-ph.HE] .
  * Raaijmakers  _et al._ [2021] G. Raaijmakers, S. K. Greif, K. Hebeler, T. Hinderer, S. Nissanke, A. Schwenk, T. E. Riley, A. L. Watts, J. M. Lattimer, and W. C. G. Ho, Constraints on the dense matter equation of state and neutron star properties from NICER’s mass-radius estimate of PSR J0740+6620 and multimessenger observations, (2021), arXiv:2105.06981 [astro-ph.HE] .
  * Biswas [2021] B. Biswas, Impact of PREX-II, the revised mass measurement of PSRJ0740+6620, and possible NICER observation on the dense matter equation of state, (2021), arXiv:2105.02886 [astro-ph.HE] .
  * Pang  _et al._ [2021] P. T. H. Pang, I. Tews, M. W. Coughlin, M. Bulla, C. Van Den Broeck, and T. Dietrich, Nuclear-Physics Multi-Messenger Astrophysics Constraints on the Neutron-Star Equation of State: Adding NICER’s PSR J0740+6620 Measurement, (2021), arXiv:2105.08688 [astro-ph.HE] .
  * Legred  _et al._ [2021] I. Legred, K. Chatziioannou, R. Essick, S. Han, and P. Landry, Impact of the PSR J0740+6620 radius constraint on the properties of high-density matter, Phys. Rev. D 104, 063003 (2021), arXiv:2106.05313 [astro-ph.HE] .
  * Essick  _et al._ [2021a] R. Essick, I. Tews, P. Landry, and A. Schwenk, Astrophysical Constraints on the Symmetry Energy and the Neutron Skin of 208Pb with Minimal Modeling Assumptions, (2021a), arXiv:2102.10074 [nucl-th] .
  * Essick  _et al._ [2021b] R. Essick, P. Landry, A. Schwenk, and I. Tews, A Detailed Examination of Astrophysical Constraints on the Symmetry Energy and the Neutron Skin of 208Pb with Minimal Modeling Assumptions, (2021b), arXiv:2107.05528 [nucl-th] .
  * Blanchet [2014] L. Blanchet, Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Binaries, Living Rev. Rel. 17, 2 (2014), arXiv:1310.1528 [gr-qc] .
  * Abbott  _et al._ [2019a] B. Abbott _et al._ (LIGO Scientific, Virgo), Properties of the binary neutron star merger GW170817, Phys. Rev. X 9, 011001 (2019a), arXiv:1805.11579 [gr-qc] .
  * Torres-Rivas  _et al._ [2019] A. Torres-Rivas, K. Chatziioannou, A. Bauswein, and J. A. Clark, Observing the post-merger signal of GW170817-like events with improved gravitational-wave detectors, Phys. Rev. D 99, 044014 (2019), arXiv:1811.08931 [gr-qc] .
  * Chatziioannou [2020] K. Chatziioannou, Neutron star tidal deformability and equation of state constraints, Gen. Rel. Grav. 52, 109 (2020), arXiv:2006.03168 [gr-qc] .
  * Dietrich  _et al._ [2021] T. Dietrich, T. Hinderer, and A. Samajdar, Interpreting Binary Neutron Star Mergers: Describing the Binary Neutron Star Dynamics, Modelling Gravitational Waveforms, and Analyzing Detections, Gen. Rel. Grav. 53, 27 (2021), arXiv:2004.02527 [gr-qc] .
  * Bauswein and Stergioulas [2019] A. Bauswein and N. Stergioulas, Spectral classification of gravitational-wave emission and equation of state constraints in binary neutron star mergers, Journal of Physics G Nuclear Physics 46, 113002 (2019), arXiv:1901.06969 [gr-qc] .
  * Baiotti [2019] L. Baiotti, Gravitational waves from neutron star mergers and their relation to the nuclear equation of state, Progress in Particle and Nuclear Physics 109, 103714 (2019), arXiv:1907.08534 [astro-ph.HE] .
  * Sarin and Lasky [2021] N. Sarin and P. D. Lasky, The evolution of binary neutron star post-merger remnants: a review, General Relativity and Gravitation 53, 59 (2021), arXiv:2012.08172 [astro-ph.HE] .
  * Bernuzzi [2020] S. Bernuzzi, Neutron Star Merger Remnants, Gen. Rel. Grav. 52, 108 (2020), arXiv:2004.06419 [astro-ph.HE] .
  * Baumgarte  _et al._ [2000] T. W. Baumgarte, S. L. Shapiro, and M. Shibata, On the maximum mass of differentially rotating neutron stars, Astrophys. J. Lett. 528, L29 (2000), arXiv:astro-ph/9910565 .
  * Xing  _et al._ [1994] Z.-G. Xing, J. M. Centrella, and S. L. W. McMillan, Gravitational radiation from coalescing binary neutron stars, Phys. Rev. D50, 6247 (1994), arXiv:gr-qc/9411029 [gr-qc] .
  * Ruffert  _et al._ [2001] M. Ruffert, H. T. Ruffert, and H. T. Janka, Coalescing neutron stars - a step towards physical models. 3. Improved numerics and different neutron star masses and spins, Astron. Astrophys. 380, 544 (2001), arXiv:astro-ph/0106229 .
  * Shibata  _et al._ [2005] M. Shibata, K. Taniguchi, and K. Uryu, Merger of binary neutron stars with realistic equations of state in full general relativity, Phys. Rev. D 71, 084021 (2005), arXiv:gr-qc/0503119 .
  * Shibata [2005] M. Shibata, Constraining Nuclear Equations of State Using Gravitational Waves from Hypermassive Neutron Stars, Physical Review Letters 94, 201101 (2005), arXiv:gr-qc/0504082 [gr-qc] .
  * Shibata and Taniguchi [2006] M. Shibata and K. Taniguchi, Merger of binary neutron stars to a black hole: Disk mass, short gamma-ray bursts, and quasinormal mode ringing, Phys. Rev. D 73, 064027 (2006).
  * Oechslin and Janka [2007] R. Oechslin and H.-T. Janka, Gravitational waves from relativistic neutron star mergers with nonzero-temperature equations of state, Phys. Rev. Lett. 99, 121102 (2007), arXiv:astro-ph/0702228 .
  * Hotokezaka  _et al._ [2011] K. Hotokezaka, K. Kyutoku, H. Okawa, M. Shibata, and K. Kiuchi, Binary neutron star mergers: Dependence on the nuclear equation of state, Phys. Rev. D 83, 124008 (2011), arXiv:1105.4370 [astro-ph.HE] .
  * Bauswein and Janka [2012] A. Bauswein and H. T. Janka, Measuring Neutron-Star Properties via Gravitational Waves from Neutron-Star Mergers, Physical Review Letters 108, 011101 (2012), arXiv:1106.1616 [astro-ph.SR] .
  * Bauswein  _et al._ [2012] A. Bauswein, H. Janka, K. Hebeler, and A. Schwenk, Equation-of-state dependence of the gravitational-wave signal from the ring-down phase of neutron-star mergers, Phys. Rev. D 86, 063001 (2012), arXiv:1204.1888 [astro-ph.SR] .
  * Hotokezaka  _et al._ [2013] K. Hotokezaka, K. Kiuchi, K. Kyutoku, H. Okawa, Y.-i. Sekiguchi, M. Shibata, and K. Taniguchi, Mass ejection from the merger of binary neutron stars, Phys. Rev. D 87, 024001 (2013), arXiv:1212.0905 [astro-ph.HE] .
  * Takami  _et al._ [2014] K. Takami, L. Rezzolla, and L. Baiotti, Constraining the Equation of State of Neutron Stars from Binary Mergers, Physical Review Letters 113, 091104 (2014), arXiv:1403.5672 [gr-qc] .
  * Takami  _et al._ [2015] K. Takami, L. Rezzolla, and L. Baiotti, Spectral properties of the post-merger gravitational-wave signal from binary neutron stars, Phys. Rev. D 91, 064001 (2015).
  * Kastaun and Galeazzi [2015] W. Kastaun and F. Galeazzi, Properties of hypermassive neutron stars formed in mergers of spinning binaries, Phys. Rev. D 91, 064027 (2015).
  * Bernuzzi  _et al._ [2015] S. Bernuzzi, T. Dietrich, and A. Nagar, Modeling the Complete Gravitational Wave Spectrum of Neutron Star Mergers, Physical Review Letters 115, 091101 (2015), arXiv:1504.01764 [gr-qc] .
  * Bauswein and Stergioulas [2015] A. Bauswein and N. Stergioulas, Unified picture of the post-merger dynamics and gravitational wave emission in neutron star mergers, Phys. Rev. D 91, 124056 (2015).
  * Foucart  _et al._ [2016] F. Foucart, R. Haas, M. D. Duez, E. O’Connor, C. D. Ott, L. Roberts, L. E. Kidder, J. Lippuner, H. P. Pfeiffer, and M. A. Scheel, Low mass binary neutron star mergers: Gravitational waves and neutrino emission, Phys. Rev. D 93, 044019 (2016), arXiv:1510.06398 [astro-ph.HE] .
  * Lehner  _et al._ [2016] L. Lehner, S. L. Liebling, C. Palenzuela, O. L. Caballero, E. O’Connor, M. Anderson, and D. Neilsen, Unequal mass binary neutron star mergers and multimessenger signals, Class. Quant. Grav. 33, 184002 (2016), arXiv:1603.00501 [gr-qc] .
  * East  _et al._ [2016] W. E. East, V. Paschalidis, and F. Pretorius, Equation of state effects and one-arm spiral instability in hypermassive neutron stars formed in eccentric neutron star mergers, Classical and Quantum Gravity 33, 244004 (2016), arXiv:1609.00725 [astro-ph.HE] .
  * Dietrich  _et al._ [2017] T. Dietrich, M. Ujevic, W. Tichy, S. Bernuzzi, and B. Brügmann, Gravitational waves and mass ejecta from binary neutron star mergers: Effect of the mass ratio, Phys. Rev. D 95, 024029 (2017), arXiv:1607.06636 [gr-qc] .
  * Blacker  _et al._ [2020] S. Blacker, N.-U. F. Bastian, A. Bauswein, D. B. Blaschke, T. Fischer, M. Oertel, T. Soultanis, and S. Typel, Constraining the onset density of the hadron-quark phase transition with gravitational-wave observations, Phys. Rev. D 102, 123023 (2020), arXiv:2006.03789 [astro-ph.HE] .
  * Raithel  _et al._ [2021] C. A. Raithel, V. Paschalidis, and F. Özel, Realistic Finite-Temperature Effects in Neutron Star Merger Simulations, (2021), arXiv:2104.07226 [astro-ph.HE] .
  * Hammond  _et al._ [2021] P. Hammond, I. Hawke, and N. Andersson, Thermal aspects of neutron star mergers, (2021), arXiv:2108.08649 [astro-ph.HE] .
  * Radice  _et al._ [2017] D. Radice, S. Bernuzzi, W. Del Pozzo, L. F. Roberts, and C. D. Ott, Probing Extreme-Density Matter with Gravitational Wave Observations of Binary Neutron Star Merger Remnants, Astrophys. J. Lett. 842, L10 (2017), arXiv:1612.06429 [astro-ph.HE] .
  * Most  _et al._ [2019] E. R. Most, L. J. Papenfort, V. Dexheimer, M. Hanauske, S. Schramm, H. Stöcker, and L. Rezzolla, Signatures of quark-hadron phase transitions in general-relativistic neutron-star mergers, Phys. Rev. Lett. 122, 061101 (2019), arXiv:1807.03684 [astro-ph.HE] .
  * Bauswein  _et al._ [2019a] A. Bauswein, N.-U. F. Bastian, D. B. Blaschke, K. Chatziioannou, J. A. Clark, T. Fischer, and M. Oertel, Identifying a first-order phase transition in neutron star mergers through gravitational waves, Phys. Rev. Lett. 122, 061102 (2019a), arXiv:1809.01116 [astro-ph.HE] .
  * Bauswein  _et al._ [2019b] A. Bauswein, N.-U. Friedrich Bastian, D. Blaschke, K. Chatziioannou, J. A. Clark, T. Fischer, H.-T. Janka, O. Just, M. Oertel, and N. Stergioulas, Equation-of-state Constraints and the QCD Phase Transition in the Era of Gravitational-Wave Astronomy, AIP Conf. Proc. 2127, 020013 (2019b), arXiv:1904.01306 [astro-ph.HE] .
  * Weih  _et al._ [2020] L. R. Weih, M. Hanauske, and L. Rezzolla, Postmerger Gravitational-Wave Signatures of Phase Transitions in Binary Mergers, Phys. Rev. Lett. 124, 171103 (2020), arXiv:1912.09340 [gr-qc] .
  * Bauswein and Blacker [2020] A. Bauswein and S. Blacker, Impact of quark deconfinement in neutron star mergers and hybrid star mergers, Eur. Phys. J. ST 229, 3595 (2020), arXiv:2006.16183 [astro-ph.HE] .
  * Liebling  _et al._ [2021] S. L. Liebling, C. Palenzuela, and L. Lehner, Effects of High Density Phase Transitions on Neutron Star Dynamics, Class. Quant. Grav. 38, 115007 (2021), arXiv:2010.12567 [gr-qc] .
  * Prakash  _et al._ [2021] A. Prakash, D. Radice, D. Logoteta, A. Perego, V. Nedora, I. Bombaci, R. Kashyap, S. Bernuzzi, and A. Endrizzi, Signatures of deconfined quark phases in binary neutron star mergers, (2021), arXiv:2106.07885 [astro-ph.HE] .
  * Lioutas  _et al._ [2021] G. Lioutas, A. Bauswein, and N. Stergioulas, Frequency deviations in universal relations of isolated neutron stars and postmerger remnants, Phys. Rev. D 104, 043011 (2021), arXiv:2102.12455 [astro-ph.HE] .
  * Margalit and Metzger [2019] B. Margalit and B. D. Metzger, The Multi-Messenger Matrix: the Future of Neutron Star Merger Constraints on the Nuclear Equation of State, Astrophys. J. Lett. 880, L15 (2019), arXiv:1904.11995 [astro-ph.HE] .
  * Bauswein  _et al._ [2020] A. Bauswein, S. Blacker, V. Vijayan, N. Stergioulas, K. Chatziioannou, J. A. Clark, N.-U. F. Bastian, D. B. Blaschke, M. Cierniak, and T. Fischer, Equation of state constraints from the threshold binary mass for prompt collapse of neutron star mergers, Phys. Rev. Lett. 125, 141103 (2020), arXiv:2004.00846 [astro-ph.HE] .
  * Abbott  _et al._ [2017b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Search for Post-merger Gravitational Waves from the Remnant of the Binary Neutron Star Merger GW170817, Astrophys. J. Lett. 851, L16 (2017b), arXiv:1710.09320 [astro-ph.HE] .
  * Abbott  _et al._ [2019b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Search for gravitational waves from a long-lived remnant of the binary neutron star merger GW170817, Astrophys. J. 875, 160 (2019b), arXiv:1810.02581 [gr-qc] .
  * Dudi  _et al._ [2018] R. Dudi, F. Pannarale, T. Dietrich, M. Hannam, S. Bernuzzi, F. Ohme, and B. Brügmann, Relevance of tidal effects and post-merger dynamics for binary neutron star parameter estimation, Phys. Rev. D 98, 084061 (2018), arXiv:1808.09749 [gr-qc] .
  * Vretinaris  _et al._ [2020] S. Vretinaris, N. Stergioulas, and A. Bauswein, Empirical relations for gravitational-wave asteroseismology of binary neutron star mergers, Phys. Rev. D 101, 084039 (2020), arXiv:1910.10856 [gr-qc] .
  * Hotokezaka  _et al._ [2013] K. Hotokezaka, K. Kiuchi, K. Kyutoku, T. Muranushi, Y.-i. Sekiguchi, M. Shibata, and K. Taniguchi, Remnant massive neutron stars of binary neutron star mergers: Evolution process and gravitational waveform, Phys. Rev. D 88, 044026 (2013), arXiv:1307.5888 [astro-ph.HE] .
  * Bauswein  _et al._ [2016] A. Bauswein, N. Stergioulas, and H.-T. Janka, Exploring properties of high-density matter through remnants of neutron-star mergers, Eur. Phys. J. A 52, 56 (2016), arXiv:1508.05493 [astro-ph.HE] .
  * Bose  _et al._ [2018] S. Bose, K. Chakravarti, L. Rezzolla, B. S. Sathyaprakash, and K. Takami, Neutron-star Radius from a Population of Binary Neutron Star Mergers, Phys. Rev. Lett. 120, 031102 (2018), arXiv:1705.10850 [gr-qc] .
  * Tsang  _et al._ [2019] K. W. Tsang, T. Dietrich, and C. Van Den Broeck, Modeling the postmerger gravitational wave signal and extracting binary properties from future binary neutron star detections, Phys. Rev. D 100, 044047 (2019), arXiv:1907.02424 [gr-qc] .
  * Breschi  _et al._ [2019] M. Breschi, S. Bernuzzi, F. Zappa, M. Agathos, A. Perego, D. Radice, and A. Nagar, kiloHertz gravitational waves from binary neutron star remnants: time-domain model and constraints on extreme matter, Phys. Rev. D 100, 104029 (2019), arXiv:1908.11418 [gr-qc] .
  * Easter  _et al._ [2020] P. J. Easter, S. Ghonge, P. D. Lasky, A. R. Casey, J. A. Clark, F. H. Vivanco, and K. Chatziioannou, Detection and parameter estimation of binary neutron star merger remnants, Phys. Rev. D 102, 043011 (2020), arXiv:2006.04396 [astro-ph.HE] .
  * Soultanis  _et al._ [2021] T. Soultanis, A. Bauswein, and N. Stergioulas, Analytic model of the spectral properties of gravitational waves from neutron star merger remnants, (2021), arXiv:2111.08353 [astro-ph.HE] .
  * Clark  _et al._ [2014] J. Clark, A. Bauswein, L. Cadonati, H. T. Janka, C. Pankow, and N. Stergioulas, Prospects For High Frequency Burst Searches Following Binary Neutron Star Coalescence With Advanced Gravitational Wave Detectors, Phys. Rev. D 90, 062004 (2014), arXiv:1406.5444 [astro-ph.HE] .
  * Clark  _et al._ [2016] J. A. Clark, A. Bauswein, N. Stergioulas, and D. Shoemaker, Observing Gravitational Waves From The Post-Merger Phase Of Binary Neutron Star Coalescence, Class. Quant. Grav. 33, 085003 (2016), arXiv:1509.08522 [astro-ph.HE] .
  * Chatziioannou  _et al._ [2017] K. Chatziioannou, J. A. Clark, A. Bauswein, M. Millhouse, T. B. Littenberg, and N. Cornish, Inferring the post-merger gravitational wave emission from binary neutron star coalescences, Phys. Rev. D 96, 124035 (2017), arXiv:1711.00040 [gr-qc] .
  * Cornish and Littenberg [2015] N. J. Cornish and T. B. Littenberg, BayesWave: Bayesian Inference for Gravitational Wave Bursts and Instrument Glitches, Class. Quant. Grav. 32, 135012 (2015), arXiv:1410.3835 [gr-qc] .
  * Cornish  _et al._ [2021] N. J. Cornish, T. B. Littenberg, B. Bécsy, K. Chatziioannou, J. A. Clark, S. Ghonge, and M. Millhouse, BayesWave analysis pipeline in the era of gravitational wave observations, Phys. Rev. D 103, 044006 (2021), arXiv:2011.09494 [gr-qc] .
  * LIGO Scientific Collaboration and Virgo Collaboration [2018] LIGO Scientific Collaboration and Virgo Collaboration, BayesWave, https://git.ligo.org/lscsoft/bayeswave (2018).
  * Kanner  _et al._ [2016] J. B. Kanner, T. B. Littenberg, N. Cornish, M. Millhouse, E. Xhakaj, F. Salemi, M. Drago, G. Vedovato, and S. Klimenko, Leveraging waveform complexity for confident detection of gravitational waves, Phys. Rev. D 93, 022002 (2016).
  * Littenberg  _et al._ [2016] T. B. Littenberg, J. B. Kanner, N. J. Cornish, and M. Millhouse, Enabling high confidence detections of gravitational-wave bursts, Phys. Rev. D 94, 044050 (2016).
  * Bécsy  _et al._ [2017] B. Bécsy, P. Raffai, N. J. Cornish, R. Essick, J. Kanner, E. Katsavounidis, T. B. Littenberg, M. Millhouse, and S. Vitale, Parameter Estimation for Gravitational-wave Bursts with the BayesWave Pipeline, Astrophys. J. 839, 15 (2017), arXiv:1612.02003 [astro-ph.HE] .
  * Tsang  _et al._ [2018] K. W. Tsang, M. Rollier, A. Ghosh, A. Samajdar, M. Agathos, K. Chatziioannou, V. Cardoso, G. Khanna, and C. Van Den Broeck, A morphology-independent data analysis method for detecting and characterizing gravitational wave echoes, Phys. Rev. D 98, 024023 (2018), arXiv:1804.04877 [gr-qc] .
  * Tsang  _et al._ [2020] K. W. Tsang, A. Ghosh, A. Samajdar, K. Chatziioannou, S. Mastrogiovanni, M. Agathos, and C. Van Den Broeck, A morphology-independent search for gravitational wave echoes in data from the first and second observing runs of Advanced LIGO and Advanced Virgo, Phys. Rev. D 101, 064012 (2020), arXiv:1906.11168 [gr-qc] .
  * Pankow  _et al._ [2018] C. Pankow _et al._ , Mitigation of the instrumental noise transient in gravitational-wave data surrounding GW170817, Phys. Rev. D98, 084016 (2018), arXiv:1808.03619 [gr-qc] .
  * Millhouse  _et al._ [2018] M. Millhouse, N. J. Cornish, and T. Littenberg, Bayesian reconstruction of gravitational wave bursts using chirplets, Phys. Rev. D97, 104057 (2018), arXiv:1804.03239 [gr-qc] .
  * Ghonge  _et al._ [2020] S. Ghonge, K. Chatziioannou, J. A. Clark, T. Littenberg, M. Millhouse, L. Cadonati, and N. Cornish, Reconstructing gravitational wave signals from binary black hole mergers with minimal assumptions, (2020), arXiv:2003.09456 [gr-qc] .
  * Dálya  _et al._ [2021] G. Dálya, P. Raffai, and B. Bécsy, Bayesian Reconstruction of Gravitational-wave Signals from Binary Black Holes with Nonzero Eccentricities, Class. Quant. Grav. 38, 065002 (2021), arXiv:2006.06256 [astro-ph.HE] .
  * Chatziioannou  _et al._ [2021a] K. Chatziioannou, M. Isi, C.-J. Haster, and T. B. Littenberg, Morphology-independent test of the mixed polarization content of transient gravitational wave signals, Phys. Rev. D 104, 044005 (2021a), arXiv:2105.01521 [gr-qc] .
  * Chatziioannou  _et al._ [2021b] K. Chatziioannou, N. Cornish, M. Wijngaarden, and T. B. Littenberg, Modeling compact binary signals and instrumental glitches in gravitational wave data, Phys. Rev. D 103, 044013 (2021b), arXiv:2101.01200 [gr-qc] .
  * Green [1995] P. J. Green, Reversible jump Markov chain Monte Carlo computation and Bayesian model determination, Biometrika 82, 711 (1995), http://oup.prod.sis.lan/biomet/article-pdf/82/4/711/699533/82-4-711.pdf .
  * Lartillot and Philippe [2006] N. Lartillot and H. Philippe, Computing Bayes Factors Using Thermodynamic Integration, Systematic Biology 55, 195 (2006), https://academic.oup.com/sysbio/article-pdf/55/2/195/26557316/10635150500433722.pdf .
  * LIGO Scientific Collaboration, Virgo Collaboration [2018] LIGO Scientific Collaboration, Virgo Collaboration, LALSuite, https://git.ligo.org/lscsoft/lalsuite (2018).
  * Husa  _et al._ [2016] S. Husa, S. Khan, M. Hannam, M. Pürrer, F. Ohme, X. Jiménez Forteza, and A. Bohé, Frequency-domain gravitational waves from nonprecessing black-hole binaries. I. New numerical waveforms and anatomy of the signal, Phys. Rev. D 93, 044006 (2016), arXiv:1508.07250 [gr-qc] .
  * Khan  _et al._ [2016] S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer, X. Jiménez Forteza, and A. Bohé, Frequency-domain gravitational waves from nonprecessing black-hole binaries. II. A phenomenological model for the advanced detector era, Phys. Rev. D 93, 044007 (2016), arXiv:1508.07253 [gr-qc] .
  * Dietrich  _et al._ [2017] T. Dietrich, S. Bernuzzi, and W. Tichy, Closed-form tidal approximants for binary neutron star gravitational waveforms constructed from high-resolution numerical relativity simulations, Phys. Rev. D 96, 121501 (2017), arXiv:1706.02969 [gr-qc] .
  * Dietrich  _et al._ [2019a] T. Dietrich _et al._ , Matter imprints in waveform models for neutron star binaries: Tidal and self-spin effects, Phys. Rev. D99, 024029 (2019a), arXiv:1804.02235 [gr-qc] .
  * Dietrich  _et al._ [2019b] T. Dietrich, A. Samajdar, S. Khan, N. K. Johnson-McDaniel, R. Dudi, and W. Tichy, Improving the NRTidal model for binary neutron star systems, Phys. Rev. D 100, 044003 (2019b), arXiv:1905.06011 [gr-qc] .
  * Pürrer  _et al._ [2016] M. Pürrer, M. Hannam, and F. Ohme, Can we measure individual black-hole spins from gravitational-wave observations?, Phys. Rev. D 93, 084042 (2016), arXiv:1512.04955 [gr-qc] .
  * Chatziioannou  _et al._ [2018a] K. Chatziioannou, G. Lovelace, M. Boyle, M. Giesler, D. A. Hemberger, R. Katebi, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, and B. Szilágyi, Measuring the properties of nearly extremal black holes with gravitational waves, Phys. Rev. D 98, 044028 (2018a), arXiv:1804.03704 [gr-qc] .
  * Racine [2008] E. Racine, Analysis of spin precession in binary black hole systems including quadrupole-monopole interaction, Phys. Rev. D78, 044021 (2008), arXiv:0803.1820 [gr-qc] .
  * Favata [2014] M. Favata, Systematic parameter errors in inspiraling neutron star binaries, Phys. Rev. Lett. 112, 101101 (2014), arXiv:1310.8288 [gr-qc] .
  * Wade  _et al._ [2014] L. Wade, J. D. Creighton, E. Ochsner, B. D. Lackey, B. F. Farr, _et al._ , Systematic and statistical errors in a bayesian approach to the estimation of the neutron-star equation of state using advanced gravitational wave detectors, Phys. Rev. D 89, 103012 (2014), arXiv:1402.5156 [gr-qc] .
  * Binnington and Poisson [2009] T. Binnington and E. Poisson, Relativistic theory of tidal Love numbers, Phys. Rev. D80, 084018 (2009), arXiv:0906.1366 [gr-qc] .
  * Poisson [2021] E. Poisson, Compact body in a tidal environment: New types of relativistic Love numbers, and a post-Newtonian operational definition for tidally induced multipole moments, Phys. Rev. D 103, 064023 (2021), arXiv:2012.10184 [gr-qc] .
  * Cutler and Flanagan [1994] C. Cutler and E. E. Flanagan, Gravitational waves from merging compact binaries: How accurately can one extract the binary’s parameters from the inspiral wave form?, Phys. Rev. D 49, 2658 (1994), arXiv:gr-qc/9402014 .
  * Callister [2021] T. Callister, A Thesaurus for Common Priors in Gravitational-Wave Astronomy, (2021), arXiv:2104.09508 [gr-qc] .
  * Chatziioannou  _et al._ [2018b] K. Chatziioannou, C.-J. Haster, and A. Zimmerman, Measuring the neutron star tidal deformability with equation-of-state-independent relations and gravitational waves, Phys. Rev. D97, 104036 (2018b), arXiv:1804.03221 [gr-qc] .
  * Zhao and Lattimer [2018] T. Zhao and J. M. Lattimer, Tidal Deformabilities and Neutron Star Mergers, Phys. Rev. D 98, 063020 (2018), arXiv:1808.02858 [astro-ph.HE] .
  * Carson  _et al._ [2019] Z. Carson, K. Chatziioannou, C.-J. Haster, K. Yagi, and N. Yunes, Equation-of-state insensitive relations after GW170817, Phys. Rev. D 99, 083016 (2019), arXiv:1903.03909 [gr-qc] .
  * Cornish [2021] N. J. Cornish, Rapid and Robust Parameter Inference for Binary Mergers, Phys. Rev. D 103, 104057 (2021), arXiv:2101.01188 [gr-qc] .
  * Ter Braak [2006] C. J. F. Ter Braak, A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: easy Bayesian computing for real parameter spaces, Statistics and Computing 16, 239 (2006).
  * Cornish [2010] N. J. Cornish, Fast Fisher Matrices and Lazy Likelihoods, arXiv e-prints , arXiv:1007.4820 (2010), arXiv:1007.4820 [gr-qc] .
  * Raymond and Farr [2014] V. Raymond and W. M. Farr, Physically motivated exploration of the extrinsic parameter space in ground-based gravitational-wave astronomy, (2014), arXiv:1402.0053 [gr-qc] .
  * Veitch  _et al._ [2015] J. Veitch _et al._ , Parameter estimation for compact binaries with ground-based gravitational-wave observations using the LALInference software library, Phys. Rev. D 91, 042003 (2015), arXiv:1409.7215 [gr-qc] .
  * Littenberg and Cornish [2015] T. B. Littenberg and N. J. Cornish, Bayesian inference for spectral estimation of gravitational wave detector noise, Phys. Rev. D91, 084034 (2015), arXiv:1410.3852 [gr-qc] .
  * Chatziioannou  _et al._ [2019] K. Chatziioannou, C.-J. Haster, T. B. Littenberg, W. M. Farr, S. Ghonge, M. Millhouse, J. A. Clark, and N. Cornish, Noise spectral estimation methods and their impact on gravitational wave measurement of compact binary mergers, Phys. Rev. D 100, 104004 (2019).
  * Abbott  _et al._ [2016a] B. Abbott _et al._ (LIGO Scientific, Virgo), Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ [2017c] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Effects of waveform model systematics on the interpretation of GW150914, Class. Quant. Grav. 34, 104002 (2017c), arXiv:1611.07531 [gr-qc] .
  * Gravitational Wave Open Science Center () [GWOSC] Gravitational Wave Open Science Center (GWOSC), https://www.gw-openscience.org/.
  * Abbott  _et al._ [2021] R. Abbott _et al._ (LIGO Scientific, Virgo), Open data from the first and second observing runs of Advanced LIGO and Advanced Virgo, SoftwareX 13, 100658 (2021), arXiv:1912.11716 [gr-qc] .
  * Abbott  _et al._ [2016b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Properties of the Binary Black Hole Merger GW150914, Phys. Rev. Lett. 116, 241102 (2016b), arXiv:1602.03840 [gr-qc] .
  * Abbott  _et al._ [2016c] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Binary Black Hole Mergers in the first Advanced LIGO Observing Run, Phys. Rev. X 6, 041015 (2016c), [Erratum: Phys.Rev.X 8, 039903 (2018)], arXiv:1606.04856 [gr-qc] .
  * Abbott  _et al._ [2018a] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, (2018a), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ [2018b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GW170817: Measurements of neutron star radii and equation of state, Phys. Rev. Lett. 121, 161101 (2018b), arXiv:1805.11581 [gr-qc] .
  * [135] BayesWave Glitch Subtraction for GW170817, https://dcc.ligo.org/LIGO-T1700406/public.
  * Soares-Santos  _et al._ [2017] M. Soares-Santos, et al., Dark Energy Survey, and Dark Energy Camera GW-EM Collaboration, The Electromagnetic Counterpart of the Binary Neutron Star Merger LIGO/Virgo GW170817. I. Discovery of the Optical Counterpart Using the Dark Energy Camera, Astrophysical Journal, Letters 848, L16 (2017), arXiv:1710.05459 [astro-ph.HE] .
  * Abbott  _et al._ [2017d] B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, and et al., Multi-messenger observations of a binary neutron star merger, The Astrophysical Journal 848, L12 (2017d), arXiv:1710.05833 .
  * Dietrich  _et al._ [2019] T. Dietrich, S. Khan, R. Dudi, S. J. Kapadia, P. Kumar, A. Nagar, F. Ohme, F. Pannarale, A. Samajdar, S. Bernuzzi, G. Carullo, W. Del Pozzo, M. Haney, C. Markakis, M. Pürrer, G. Riemenschneider, Y. E. Setyawati, K. W. Tsang, and C. Van Den Broeck, Matter imprints in waveform models for neutron star binaries: Tidal and self-spin effects, Phys. Rev. D 99, 024029 (2019), arXiv:1804.02235 [gr-qc] .
  * Abbott  _et al._ [2017e] B. P. Abbott _et al._ (LIGO Scientific Collaboration), Exploring the Sensitivity of Next Generation Gravitational Wave Detectors, Class. Quant. Grav. 34, 044001 (2017e), arXiv:1607.08697 [astro-ph.IM] .
  * Reitze  _et al._ [2019a] D. Reitze _et al._ , The US Program in Ground-Based Gravitational Wave Science: Contribution from the LIGO Laboratory, Bull. Am. Astron. Soc. 51, 141 (2019a), arXiv:1903.04615 [astro-ph.IM] .
  * Reitze  _et al._ [2019b] D. Reitze _et al._ , Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO, Bull. Am. Astron. Soc. 51, 035 (2019b), arXiv:1907.04833 [astro-ph.IM] .
  * Hild  _et al._ [2011] S. Hild _et al._ , Sensitivity studies for third-generation gravitational wave observatories, Classical and Quantum Gravity 28, 094013 (2011), arXiv:1012.0908 [gr-qc] .
  * Punturo  _et al._ [2010] M. Punturo _et al._ , The third generation of gravitational wave observatories and their science reach, Classical and Quantum Gravity 27, 084007 (2010).
  * Barsotti  _et al._ [2018] L. Barsotti, P. Fritschel, M. Evans, and S. Gras,  _Updated Advanced LIGO sensitivity design curve_, Tech. Rep. (2018).
  * Kastaun and Ohme [2021] W. Kastaun and F. Ohme, Numerical inside view of hypermassive remnant models for GW170817, Phys. Rev. D 104, 023001 (2021), arXiv:2103.01586 [astro-ph.HE] .
  * Hempel and Schaffner-Bielich [2010] M. Hempel and J. Schaffner-Bielich, A statistical model for a complete supernova equation of state, Nuclear Physics A 837, 210 (2010), arXiv:0911.4073 [nucl-th] .
  * Wilson  _et al._ [1996] J. R. Wilson, G. J. Mathews, and P. Marronetti, Relativistic numerical model for close neutron-star binaries, Phys. Rev. D 54, 1317 (1996).
  * Isenberg and Nester [1980] J. Isenberg and J. Nester, Canonical Gravity, in  _General Relativity and Gravitation. Vol. 1. One hundred years after the birth of Albert Einstein. Edited by A. Held. New York, NY: Plenum Press, p. 23, 1980_ , edited by A. Held (1980) p. 23.
  * Alford  _et al._ [2005] M. Alford, M. Braby, M. W. Paris, and S. Reddy, Hybrid stars that masquerade as neutron stars, Astrophys. J. 629, 969 (2005), arXiv:nucl-th/0411016 .
  * Alford  _et al._ [2015] M. G. Alford, G. F. Burgio, S. Han, G. Taranto, and D. Zappalà, Constraining and applying a generic high-density equation of state, Phys. Rev. D 92, 083002 (2015), arXiv:1501.07902 [nucl-th] .
  * Del Pozzo  _et al._ [2013] W. Del Pozzo, T. G. F. Li, M. Agathos, C. Van Den Broeck, and S. Vitale, Demonstrating the feasibility of probing the neutron star equation of state with second-generation gravitational wave detectors, Phys. Rev. Lett. 111, 071101 (2013), arXiv:1307.8338 [gr-qc] .
  * Agathos  _et al._ [2015] M. Agathos, J. Meidam, W. Del Pozzo, T. G. F. Li, M. Tompitak, J. Veitch, S. Vitale, and C. Van Den Broeck, Constraining the neutron star equation of state with gravitational wave signals from coalescing binary neutron stars, Phys. Rev. D 92, 023012 (2015), arXiv:1503.05405 [gr-qc] .
  * Chatziioannou  _et al._ [2015] K. Chatziioannou, K. Yagi, A. Klein, N. Cornish, and N. Yunes, Probing the Internal Composition of Neutron Stars with Gravitational Waves, Phys. Rev. D 92, 104008 (2015), arXiv:1508.02062 [gr-qc] .
  * Han and Steiner [2019] S. Han and A. W. Steiner, Tidal deformability with sharp phase transitions in (binary) neutron stars, Phys. Rev. D 99, 083014 (2019), arXiv:1810.10967 [nucl-th] .
  * Chen  _et al._ [2020] H.-Y. Chen, P. M. Chesler, and A. Loeb, Searching for exotic cores with binary neutron star inspirals, Astrophys. J. Lett. 893, L4 (2020), arXiv:1909.04096 [astro-ph.HE] .
  * Chatziioannou and Han [2020] K. Chatziioannou and S. Han, Studying strong phase transitions in neutron stars with gravitational waves, Phys. Rev. D 101, 044019 (2020), arXiv:1911.07091 [gr-qc] .
  * Han and Prakash [2020] S. Han and M. Prakash, On the Minimum Radius of Very Massive Neutron Stars, Astrophys. J. 899, 164 (2020), arXiv:2006.02207 [astro-ph.HE] .
  * Zhang and Li [2019] N.-B. Zhang and B.-A. Li, Implications of the Mass $M=2.17^{+0.11}_{-0.10}$M⊙ of PSR J0740+6620 on the Equation of State of Super-dense Neutron-rich Nuclear Matter, Astrophys. J. 879, 99 (2019), arXiv:1904.10998 [nucl-th] .
  * Pang  _et al._ [2020] P. T. H. Pang, T. Dietrich, I. Tews, and C. Van Den Broeck, Parameter estimation for strong phase transitions in supranuclear matter using gravitational-wave astronomy, Phys. Rev. Res. 2, 033514 (2020), arXiv:2006.14936 [astro-ph.HE] .
  * Drischler  _et al._ [2021] C. Drischler, S. Han, J. M. Lattimer, M. Prakash, S. Reddy, and T. Zhao, Limiting masses and radii of neutron stars and their implications, Phys. Rev. C 103, 045808 (2021), arXiv:2009.06441 [nucl-th] .
  * Fischer  _et al._ [2018] T. Fischer, N.-U. F. Bastian, M.-R. Wu, P. Baklanov, E. Sorokina, S. Blinnikov, S. Typel, T. Klähn, and D. B. Blaschke, Quark deconfinement as a supernova explosion engine for massive blue supergiant stars, Nature Astronomy 2, 980 (2018), arXiv:1712.08788 [astro-ph.HE] .
  * Bastian [2021] N.-U. F. Bastian, Phenomenological quark-hadron equations of state with first-order phase transitions for astrophysical applications, Phys. Rev. D 103, 023001 (2021), arXiv:2009.10846 [nucl-th] .
  * Stergioulas  _et al._ [2011] N. Stergioulas, A. Bauswein, K. Zagkouris, and H.-T. Janka, Gravitational waves and non-axisymmetric oscillation modes in mergers of compact object binaries, Monthly Notices of the Royal Astronomical Society 418, 427 (2011), arXiv:1105.0368 [gr-qc] .
  * Samajdar and Dietrich [2019] A. Samajdar and T. Dietrich, Waveform systematics for binary neutron star gravitational wave signals: Effects of spin, precession, and the observation of electromagnetic counterparts, Phys. Rev. D 100, 024046 (2019), arXiv:1905.03118 [gr-qc] .
  * Gamba  _et al._ [2021] R. Gamba, M. Breschi, S. Bernuzzi, M. Agathos, and A. Nagar, Waveform systematics in the gravitational-wave inference of tidal parameters and equation of state from binary neutron star signals, Phys. Rev. D 103, 124015 (2021), arXiv:2009.08467 [gr-qc] .
  * Chatziioannou [2021] K. Chatziioannou, Uncertainty limits on neutron star radius measurements with gravitational waves, (2021), arXiv:2108.12368 [gr-qc] .
  * Kunert  _et al._ [2021] N. Kunert, P. T. H. Pang, I. Tews, M. W. Coughlin, and T. Dietrich, Quantifying modelling uncertainties when combining multiple gravitational-wave detections from binary neutron star sources, (2021), arXiv:2110.11835 [astro-ph.HE] .
  * 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.HE] .
  * Regimbau  _et al._ [2012] T. Regimbau _et al._ , A Mock Data Challenge for the Einstein Gravitational-Wave Telescope, Phys. Rev. D 86, 122001 (2012), arXiv:1201.3563 [gr-qc] .
  * Sachdev  _et al._ [2020] S. Sachdev, T. Regimbau, and B. S. Sathyaprakash, Subtracting compact binary foreground sources to reveal primordial gravitational-wave backgrounds, Phys. Rev. D 102, 024051 (2020), arXiv:2002.05365 [gr-qc] .
  * Adhikari  _et al._ [2020] R. X. Adhikari _et al._ (LIGO), A cryogenic silicon interferometer for gravitational-wave detection, Class. Quant. Grav. 37, 165003 (2020), arXiv:2001.11173 [astro-ph.IM] .
  * Haster  _et al._ [2020] C.-J. Haster, K. Chatziioannou, A. Bauswein, and J. A. Clark, Inference of the neutron star equation of state from cosmological distances, Phys. Rev. Lett. 125, 261101 (2020), arXiv:2004.11334 [gr-qc] .
  * Macleod  _et al._ [2020] D. Macleod, A. L. Urban, S. Coughlin, T. Massinger, M. Pitkin, paulaltin, J. Areeda, E. Quintero, T. G. Badger, L. Singer, and K. Leinweber, gwpy/gwpy: 1.0.1 (2020).
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2d graphics environment, Computing In Science & Engineering 9, 90 (2007).