# Waveform uncertainty quantification and interpretation for gravitational-
wave astronomy

Jocelyn Read  Nicholas and Lee Begovich Center for Gravitational Wave Physics
and Astronomy, California State University Fullerton, Fullerton CA 92831 USA
<EMAIL_ADDRESS>
###### Abstract

We demonstrate how to quantify the frequency-domain amplitude and phase
accuracy of waveform models, $\delta A$ and $\delta\phi$, in a form that could
be marginalized over in gravitational-wave inference using techniques
currently applied for quantifying calibration uncertainty. For concreteness,
waveform uncertainties affecting neutron-star inspiral measurements are
considered, and post-hoc error estimates from a variety of waveform models are
made by comparing time-domain and frequency-domain analytic models with
multiple-resolution numerical simulations. These waveform uncertainty
estimates can be compared to GW170817 calibration envelopes or to Advanced
LIGO and Virgo calibration goals. Signal-specific calibration and waveform
uncertainties are compared to statistical fluctuations in gravitational-wave
observatories, giving frequency-dependent modeling requirements for detectors
such as Advanced LIGO Plus, Cosmic Explorer, or Einstein Telescope. Finally,
the distribution of waveform error for the GW170817 posterior is computed from
tidal models and compared to the constraints on $\delta\phi$ or $\delta A$
from GWTC-1 by Edelman et. al. In general, $\delta\phi$ and $\delta A$ can
also be interpreted in terms of unmodeled astrophysical energy transfer within
or from the source system.

##  1 Introduction

Waveform models are critical for estimating the properties of source systems
in gravitational-wave astronomy. In particular for neutron-star mergers,
modeling uncertainties will be potentially significant for tidal inference in
coming observation runs [1, 2, 3]. In the analysis of gravitational-wave
observations, these waveform systematics have so far been assessed by
comparing or combining inference results from different waveform model
families. Examples include the first binary-black-hole merger GW150914 [4] and
the first neutron-star merger GW170817 [5]. Model selection can also be
applied to the predictions of variant waveform families, and marginalization
through a mixture model allows the gravitational-wave observation to inform
multi-waveform inference [6, 7, 8].

These methods generally take each waveform family’s prediction as an exact
model at each point in parameter space. They do not account for modeling
uncertainty within each waveform approximant, or its dependence on frequency.
Here, we focus on quantifying uncertainty from waveform model choices in the
signal domain, in a form that would allow marginalization over the full range
of predictions compatible with current theory and simulation. Observed strain
data could then inform the best-fit waveform model frequency-by-frequency.
Waveform uncertainty estimates in the signal domain can also be used to guide
analysis efforts by providing an indication of how accurately the source will
need to be modeled to interpret future signals.

In this work, we present a frequency-dependent characterization of waveform
error in amplitude and phase, $\delta A$ and $\delta\phi$, and demonstrate its
application to existing waveform models of neutron star inspirals with tidal
contributions. To demonstrate current levels of uncertainty in modeling the
physics of waveform generation, we compare multiple waveform approximants with
each other and with a high-resolution numerical simulation from the CoreDB
waveform library in Section 4. The model differences are compared to reference
noise spectral densities, anticipating the requirements of future
observations. We demonstrate a frequency-dependent condition for the effect of
waveform modeling uncertainty to be smaller than statistical fluctuation,
given a candidate signal amplitude and a noise power spectral density.

Finally, we show two implications of this result for the interpretation of the
double neutron star merger GW170817 [9]. First, we estimate the size of
current model uncertainties by evaluating the differences between waveform
families in the posterior region of parameter space. We compare these to the
envelope of calibration uncertainties and the differences within the posterior
distribution of waveforms to assess the relative impact of waveform systematic
error. Second, we use this method to interpret results from Edelman et. al.
[10], which constrain how the amplitude and phase of GWTC-1 signals depart
from the waveform models used for analysis. Observed departures from the
waveform models can be mapped to physical effects: unmodeled luminosity or
internal energy contributions for the source system. This allows the largest
median $\delta\phi$ at $\sim 60\,$Hz to be interpreted as an possible
unmodeled energy transfer of $\delta E=\Delta E/E\sim 0.001$ relative to the
orbital binding energy of $E\simeq-0.006M_{\odot}c^{2}$. In general, the
recovered bounds on $\delta A$ and $\delta\phi$ limit the amount of unmodeled
energy transfer compatible with the observed signals.

##  2 Background

To interpret the observations of gravitational-wave astronomy, interferometer
strain data $\bm{d}=\bm{h}+\bm{n}$ is assumed to be generated by an
astrophysical strain signal $\bm{h}$ and background noise fluctuations
$\bm{n}$. The likelihood of a detector measurement $\bm{d}$, given a
particular incident wave $\bm{h}$, derives from the likelihood function for
$\bm{n}=\bm{d}-\bm{h}$. In practice, it will be computed from the Fourier
transforms of the data and signal time series, evaluated a set of discrete
frequencies:

$p(\bm{d}|\bm{h})\propto\exp\left(-\sum_{i}2\Delta
f\frac{|d_{i}-h_{i}|^{2}}{S_{n}(f_{i})}\right)$ (1)

where $S_{n}(f)$ is the power spectral density of the noise [11, 12].

This motivates the definition of a noise-weighted inner product [13] in either
discrete or continuous form,

$\left<\bm{h_{1}},\bm{h_{2}}\right>=2\Delta
f\sum_{i}\frac{h_{1i}^{*}h_{2i}+h_{1i}h_{2i}^{*}}{S_{n}(f_{i})}=4\Re\int_{0}^{f_{\text{max}}}df\frac{\tilde{h}_{1}(f)^{*}\tilde{h}_{2}(f)}{S_{n}(f)}$
(2)

This inner product definition sets $\left<\bm{n},\bm{n}\right>=1$. The
expected signal-to-noise ratio $\varrho$ of an incident astrophysical signal
$\bm{h}$ is given by $\varrho^{2}=\left<\bm{h},\bm{h}\right>$. We can then re-
write the likelihood of the data $\bm{d}$ given an incident wave $\bm{h}$ as

$p(\bm{d}|\bm{h})\propto\exp\left(-\left<\bm{d}-\bm{h},\bm{d}-\bm{h}\right>\right)$
(3)

Two waveforms are considered “indistinguishable” in a given detector if the
difference between them is smaller than the noise:
$\delta\bm{h}=\bm{h}_{1}-\bm{h}_{2}$ satisfies
$\left<\delta\bm{h},\delta\bm{h}\right><1$. Often when comparing waveform
families, especially those which agree at leading post-Newtonian order, we
will have $\varrho_{1}\simeq\varrho_{2}:=\varrho$, and indistinguishability is
estimated by the “mismatch” after waveform differences are minimized over
relative shifts in time and phase:

${\rm min}_{\Delta
t_{c},\Delta\phi_{c}}\left[\left<\delta\bm{h},\delta\bm{h}\right>\right]\gtrsim
2\varrho^{2}\left(1-{\rm max}_{\Delta
t_{c},\Delta\phi_{c}}\left[\left<\bf{h_{1}},\bf{h_{2}}\right>\right]/\sqrt{\left<\bf{h_{1}},\bf{h_{1}}\right>\left<\bf{h_{2}},\bf{h_{2}}\right>}\right),$
(4)

where the mismatch is in brackets on the RHS. Gravitational waveform
uncertainty requirements have thus been presented in terms of the mismatch
between the true waveform and the model used when inferring source properties
at a given signal-to-noise ratio $\varrho$ [14, 15]. Mismatch is an integrated
quantity over all frequencies.

In contrast, calibration uncertainties that affect the inferred $\bm{h}$ are
explicitly computed as functions of frequency [16, 17, 18] and can be
marginalized over for gravitational-wave inference [19, 11, 20]. To do this,
calibration uncertainty is cast in terms of amplitude and phase errors in the
inferred detector strain compared to the true conditions of the local
spacetime

$\tilde{h}_{\text{meas}}=\tilde{h}_{\text{true}}(f)(1+\delta
A_{\text{cal}}(f))\exp(i\delta\phi_{\text{cal}}(f))$ (5)

which differentiate the measured $\tilde{h}_{\text{meas}}$ from the true
strain $\tilde{h}_{\text{true}}(f)$ generated by an incident astrophysical
wave. These corrections arise from the error budget of mapping between
astrophysical and instrumentally measured strain [21, 22].

An approach similar to the calibration framework has been developed to
constrain differences between the true signal and the waveform model used for
inference: Edelman et al [10] estimated the size of unmodeled signal
contributions for the observations of GWTC-1 [23]. To do this, coherent
deviations across all detectors are also expressed in the form

$\tilde{h}_{\text{astr}}(f)=\tilde{h}_{\text{model}}(f)(1+\delta
A(f))\exp(i\delta\phi(f))$ (6)

where we are again using $\phi$ for the frequency-domain phase. In that work,
corrections were modeled as splines for $\delta A(f)$ and $\delta\phi(f)$, and
frequency-dependent departures from the baseline $\bm{h}$ model were
constrained for signals observed by the LIGO [24] and Virgo [25] observatories
.

In this work, we demonstrate the frequency-dependent uncertainty in amplitude
and phase $\delta A(f)$ and $\delta\phi(f)$ coming from existing waveform
models. These uncertainties are explicitly connected to the underlying time-
domain model uncertainties and the implications for the underlying physics of
the source system. To do this, we follow standard procedures in the
literature, but avoid common choices that assume a post-Newtonian framework
with infinite coalescence frequency. We define all source characteristics in
terms of weak-field gravitational-wave observables. Expressing model
uncertainty in this form will allow explicit marginalization during the
inference of source properties, and also allows the interpretation of waveform
differences in terms of far-field source energetics.

##  3 Waveform assumptions and model implications

Waveform models come from time-domain physics of the emitting source and the
response of the observatory’s interferometer. Here, we review the mapping from
the source properties to the signal prediction, to demonstrate how error
intrinsic to the source physics can be disentangled from the reference time
and phase $t_{c}$ and $\phi_{c}$ of a specific signal. The fundamental
differences coming from varying model families are those that are independent
of overall shifts in signal time and phase, and lead to a residual waveform
model error that can be additionally marginalized over in inference.

We adopt a description of the time-domain waveform, similar to that in [26],
with the assumption that we can write the source emission as a multipolar
expansion of oscillatory mode functions. Specifically, we assume the $h_{+}$
and $h_{\times}$ emitted along a line of observation from a given source can
be characterized by an expansion in spin-weighted spherical harmonics [27,
28]:

$h_{+}(t)-ih_{\times}(t)=\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{\ell}h_{\ell
m}(t)Y_{-2}^{\ell m}(\iota,\varphi)$ (7)

where $\iota$ is the inclination angle and $\varphi$ is the azimuthal angle of
the line of sight from the source. Each $h_{\ell m}$ component is decomposed
into a real time-domain amplitude and phase as

$h_{\ell m}(t)=\mathcal{A}_{\ell m}(t)e^{-i\psi_{\ell m}(t)}$ (8)

in the usual form for characterizing gravitational radiation in numerical
simulations, described for example for the lalsimulation numerical relativity
injection infrastructure [29] following numerical data formats of [30]. We
assume the time-domain amplitude scales as
$\mathcal{A}(t)=\mathcal{A}_{0}(t)\left(d_{0}/d\right)$ for a source at
luminosity distance $d$ in terms of $\mathcal{A}_{0}(t)$ at a reference
distance $d_{0}$ [31].

The astrophysical strain $\bf{h}$ measured by a single detector is a
projection on to the detector frame of the incident time domain polarizations
$h_{+}(t)$ and $h_{\times}(t)$,

$h(t)=F_{+}(\alpha,\delta,\psi_{p})h_{+}(t)+F_{\times}(\alpha,\delta,\psi_{p})h_{\times}(t)$
(9)

with the specific detector’s antenna response functions $F_{+,\times}$ that
depend on the source’s sky location (right ascension $\alpha$ and declination
$\delta$) and a polarization orientation angle $\psi_{p}$ of the source
relative to the interferometer. In a single detector observation, the
polarization angle is degenerate with phase. A sky location directly above or
below an interferometer with orthogonal arms gives
$F_{+}^{2}+F_{\times}^{2}=1$.

The incident $\bm{h}$ can be written as the real part of a sum over the
spherical harmonic mode amplitudes.

$h(t)=\sum_{\ell m}Q_{\ell m}h_{\ell m}$ (10)

where $Q_{\ell m}$ captures the sky-location-dependent detector response to
each mode given the source orientation. $Q_{\ell m}$ will be constant for
short-duration transients; for long-duration signals a transfer function
encoding the response of the detector would need to be included in the signal
model [32].

The leading-order quadrupole modes $\ell,m=2,\pm 2$ are dominant for many
gravitational-wave sources. We consider circular, non-precessing, near-equal-
mass binaries as the sources in this demonstration, so will restrict to this
case — neutron-star binary signals have masses within a range of roughly
1–2.5$\,M_{\odot}$ [33] and are well-described by the $2,\pm 2$ mode. For an
optimal sky location directly above or below the detector, with detector arms
aligned with the plus polarization, the sum of the $\ell,m=2,\pm 2$ modes
gives

$h_{+}(t)-ih_{\times}(t)=\sqrt{5/16\pi}\left(\left(1+\cos^{2}\iota\right)\cos
2\varphi+2i\cos\iota\sin 2\varphi\right)h(t)$ (11)

which yields $h(t)$ from $h_{22}(t)$ with the conversion factor
$Q_{22}=\sqrt{5/4\pi}$ for an optimally oriented, face-on source. If a general
sky location and inclination are considered,
$Q_{22}=\sqrt{5/4\pi}\left(F_{+}^{2}\left(1+\cos^{2}\iota\right)^{2}/4+F_{\times}\cos^{2}\iota\right)^{1/2}$.
Overall, measured signal amplitude for quadrupole sources will be scaled
relative to the optimal face-on and overhead configuration by the effective
distance
$d_{\text{eff}}=d\left(F_{+}^{2}\left(1+\cos^{2}\iota\right)^{2}/4+F_{\times}\cos^{2}\iota\right)^{-1/2}$
which combines the effects of luminosity distance, sky location, and the
source’s orientation angles [34].

We focus now on signals from inspiraling compact binaries, where the frequency
of the signal sweeps slowly upward as the binary evolves over multiple cycles.
We are interested in modeling the uncertainty given known physical properties
of the source — for example, those characterizing the masses, spins,
eccentricity, and tides — which we call intrinsic properties. These properties
will determine the amplitude of gravitational waves as function of the
emission frequency, and the gravitational-wave emission also drives a change
in the emission frequency. We show how to compute the Fourier-domain signal
$\tilde{h}(f)$, with $f$ the Fourier frequency, from integrations of the
instantaneous gravitational-wave frequency $F=\dot{\psi}/2\pi$ of the time-
domain phase. This results in the integration constants $t_{c}$ and $\phi_{c}$
that fix the arrival time and phase of a specific source. The intrinsic
properties of the source model generate the characteristic $\dot{F}(F)$.

A relatively slowly-varying amplitude allows the use of the stationary phase
approximation (SPA) in determining the frequency-domain Fourier transform for
each $m>0$ and $f>0$ [13, 35, 26] when the two conditions

$\displaystyle\left|\frac{d}{dt}\ln\mathcal{A}_{\ell m}(t)\right|$
$\displaystyle\ll\left|\dot{\psi}(t)\right|$ $\displaystyle|\ddot{\psi}|$
$\displaystyle\ll\dot{\psi}^{2}$ (12)

are satisfied. Overdots denote a time derivative. Error from the SPA is
expected to be smaller than windowing artifacts for gravitational-wave
inspirals [35]. For the reference neutron-star binary system used later in
this work, all waveforms including those from numerical simulation will
satisfy LHS$<0.1\times$RHS for both conditions, only reaching 0.1 as the stars
collide. This supports the use of instantaneous frequency as a characteristic
of the time-domain system. The SPA has been shown to give faithful frequency-
domain representations for state-of-the-art time-domain waveforms up to merger
[36]. Applications to LISA signals demonstrate higher order contributions that
could be used for extension through the merger phase [37].

To apply the SPA, we expand the transform integral around time where the
instantaneous frequency $F$ of the incident wave matches the Fourier frequency
$f$ of interest, specifically $T_{\ell m}$ defined by

$\dot{\psi}_{lm}\left(T_{\ell m}\right)=2\pi f$ (13)

The resulting frequency domain waveform is $\tilde{h}(f)=A(f)\exp(-i\phi(f))$,
with

$\displaystyle A(f)$ $\displaystyle=\sum_{\ell m}Q_{\ell
m}\left(\frac{2\pi}{\ddot{\psi}_{\ell m}(T_{\ell
m})}\right)^{1/2}\mathcal{A}_{\ell m}(T_{\ell m})$ (14a)
$\displaystyle\phi(f)$ $\displaystyle=\frac{\pi}{4}+\psi_{lm}(T_{\ell m})-2\pi
fT_{\ell m}$ (14b)

The SPA allows direct use of model predictions for the time-domain amplitude
and phase $\mathcal{A}_{\ell m}(t)$ and $\psi_{\ell m}(t)$ to calculate the
corresponding frequency-domain waveform.

The functions that enter into $\tilde{h}$ are relative to an explicit
coalescence time $t_{c}$, which is traditionally defined in the $f\to\infty$
limit following Cutler and Flanagan [13]. Since $f\to\infty$ does not happen
in all models, we here define this as the time at which the chosen waveform
model reaches a specific reference coalescence frequency $f_{c}$. The
LIGO/Virgo software lalsimulation instead uses a convention where
$t_{c}=t_{\text{peak}}$ is defined by the maximum amplitude of the waveform
[29]. This choice will imply a reference $f_{c}=F(t_{\text{peak}})$ which
varies between waveform models. For fixed $t_{c}$ and $\phi_{c}$, differences
between waveform models will depend on the $f_{c}$ assumed in the parameter
estimation, so will be sensitive to the definition of $t_{c}$. However, the
residual phase error as defined below will be independent of $f_{c}$.

For each time-domain waveform model, the following time and phase observables
are generated from the wave’s instantaneous frequency $F$ and its time
derivative $\dot{F}$,

$\displaystyle T(f)$ $\displaystyle=t_{c}-\int_{f}^{f_{c}}dF\,T^{\prime}(F)$
(15a) $\displaystyle\psi(f)$
$\displaystyle=\psi_{c}-2\pi\int_{f}^{f_{c}}dF\,F\,T^{\prime}(F)$ (15b)
$\displaystyle=\psi_{c}-2\pi\left(f_{c}t_{c}-fT(f)-\int_{f}^{f_{c}}T(F)dF\right)$
(15c)

defining the function $T^{\prime}(F):=\left(\dot{F}(F)\right)^{-1}=d\,T(F)/dF$
to streamline notation. These observables are entirely defined from the far-
field wave emission, and can be interpreted as $t_{c}-$(time to $f_{c}$) and
$\psi_{c}-$(time-domain phase accumulation remaining before $f_{c}$). After
the instantaneous frequency is calculated, $(T-t_{c})$ and $(\psi-\psi_{c})$
can be read in directly from numerical simulation data.

Substituting $T(f)$ and $\psi(f)$ into Eqs. 14 yields the frequency-domain
form for the quadrupole waveform in terms of the characteristic functions, as
we begin to write $\mathcal{A}(f)$ for $\mathcal{A}(T(f))$:

$\displaystyle A(f)$
$\displaystyle=Q(\bm{\theta}_{\text{ext}})\left(T^{\prime}(f)\right)^{1/2}\mathcal{A}(f)/2$
(16a) $\displaystyle\phi(f)$ $\displaystyle=\frac{\pi}{4}+\psi(f)-2\pi fT(f)$
(16b)
$\displaystyle=\frac{\pi}{4}+\psi_{c}-2\pi\left(f_{c}t_{c}-\int_{f}^{f_{c}}dF\,T(F)\right)$
(16c) $\displaystyle=\phi_{c}-2\pi
ft_{c}+2\pi\int_{f}^{f_{c}}d\tilde{f}\int_{\tilde{f}}^{f_{c}}dF\,T^{\prime}(F)\
$ (16d)

Setting $f=f_{c}$ in Eq. 16c defines the signal coalescence phase
$\phi_{c}=\frac{\pi}{4}+\psi_{c}-2\pi f_{c}t_{c}$. Note that $\phi_{c}=0$ does
not correspond to $\psi_{c}=0$. The explicit second integration over frequency
$\tilde{f}$ demonstrates how $\phi(f)$ depends on the underlying model’s
$T^{\prime}(F)$ and the two integration constants $t_{c}$ and $\phi_{c}$ in
Eq. 16c. The $t_{c}$ and $\phi_{c}$ that characterize an observed signal
emerge as integration constants in Eq. 16d, a form that allows marginalization
over reference phase $\phi_{c}$ and time $t_{c}$ in gravitational-wave
inference [38, 11]. The remaining intrinsic parameter contribution is entirely
from $T^{\prime}(F)$.

Following the form used in calibration marginalization, we write frequency-
domain signal uncertainty terms as $\delta A(f)$ and $\delta\phi(f)$

$\tilde{h}_{\text{true}}(f)=\tilde{h}_{\text{model}}(f)(1+\delta
A(f))\exp(i\delta\phi(f))$ (17)

The intrinsic model dependence of $\tilde{h}$ will be derived from the
characteristic $\mathcal{A}$ and $T^{\prime}=1/\dot{F}$ functions of the
waveform family in question. We can therefore express the impact of modeling
differences on the signal in terms of these functions, specifically in the
form

$\displaystyle\mathcal{A}_{\text{true}}(f)$
$\displaystyle=\mathcal{A}_{\text{model}}(f)\left(1+\delta\mathcal{A}(f)\right),$
(18a) $\displaystyle T^{\prime}_{\text{true}}(f)$
$\displaystyle=T^{\prime}_{\text{model}}(f)(1+\delta T^{\prime}(f))$ (18b)

It is computationally useful to track differences between models through the
intermediate functions $T(f)=T_{\text{model}}(f)+\delta T(f)$ and
$\psi(f)=\psi_{\text{model}}(f)+\delta\psi(f)$. As the contribution to $T(f)$
and $\psi(f)$ are defined relative to $f_{c}$ in Eqs. 15, the differences
$\delta T(f)$ and $\delta\psi(f)$ will also be relative to $f_{c}$. The signal
uncertainty at fixed $t_{c}$ and $\phi_{c}$ is

$\displaystyle 1+\delta A(f)$ $\displaystyle=\left(1+\delta
T^{\prime}(f)\right)^{1/2}\left(1+\delta\mathcal{A}(f)\right)$ (19a)
$\displaystyle\delta\phi(f)$ $\displaystyle=\delta\psi(f)-2\pi f\delta
T(f)=2\pi\int_{f}^{f_{c}}d\tilde{f}\int_{\tilde{f}}^{f_{c}}dF\,T^{\prime}(F)\delta
T^{\prime}(F)$ (19b)

The historic choice of $f_{c}\to\infty$ results in better-behaved modeling
uncertainties at fixed $t_{c}$ and $\phi_{c}$ for post-Newtonian expansions
than lower reference frequencies.

To determine the impact of model uncertainties on inference, it’s important to
consider that the integration constants $\phi_{c}$ and $t_{c}$ will be
searched over when identifying incident waves and marginalized over when
determining the intrinsic source properties [11, 16, 38]. The impact of
waveform uncertainty on measurements will be reduced by this marginalization.
For example, a toy modified-GR waveform model that produced only an overall
constant shift $\delta\phi(f)=\delta\phi_{0}$ across all frequencies would not
change any recovered parameters other than $\phi_{c}$. Similarly, an overall
shift of the form $\delta\phi(f)=2\pi ft_{0}$ would be entirely absorbed by
marginalization over $t_{c}$.

To find the residual phase error for a measurement scenario, consider a model
$\tilde{h}(f)$ with errors $\delta A(f)$ and $\delta\phi(f)$. The amplitude
$\delta A(f)$ in the Fourier domain is independent of shifts in time and
phase. However, when inferring the properties of this model from a signal that
differs by $\delta\phi(f)$, any model uncertainty described by $\phi_{0}+2\pi
ft_{0}$ will be absorbed by marginalization over shifts in time and phase
during inference. In the absence of noise, the maximum likelihood value of
$\phi_{0}$ and $t_{0}$ will be at the minimum of

$\displaystyle\left<h-h_{\text{model}}|h-h_{\text{model}}\right>$
$\displaystyle\propto\sum_{i}\frac{\left|A(f_{i})\right|^{2}\left|1-\exp{i\left(\delta\phi(f_{i})-\phi_{0}-2\pi
f_{i}t_{0}\right)}\right|^{2}}{S_{n}(f_{i})}$ (20a)
$\displaystyle\simeq\sum_{i}\frac{\left|A(f_{i})\right|^{2}}{S_{n}(f_{i})}\left|\delta\phi(f_{i})-\phi_{0}-2\pi
f_{i}t_{0}\right|^{2}$ (20b)

where the approximation applies when residual model error is small enough that
that a small angle approximation can apply across relevant frequencies. The
maximum likelihood condition becomes equivalent to a weighted least squares
fit of $\delta\phi(f)$ to the linear function $\phi_{0}+2\pi ft_{0}$, with
weights $|A(f)|^{2}/S_{n}(f)$ corresponding to the expected variance at each
frequency. Subtracting off the fit leaves a residual model error
$\delta\phi_{\text{res}}$. This characterizes the frequency-dependent waveform
phase uncertainty that should be marginalized over in addition the existing
time and phase marginalization.

##  4 Uncertainty estimates for current models

A proper error estimate should rise from careful analysis of the range of
viable analytic choices, universal relation uncertainties, numerical
resolution error, neglected effects in modeling, or other choices made by the
waveform modelers. The effect of a particular range of “reasonable” choices
will determine the frequency-dependent error distribution of the model-
characteristic functions $T^{\prime}(F)=1/\dot{F}(F)$ and $\mathcal{A}(F)$, as
in Eqs. 18, where $F$ is defined as instantaneous frequency from the phase
evolution of the waveform model. These distributions will best be estimated
directly by those creating state-of-the-art models, and could be released in
concert with a waveform implementation.

To illustrate the size of current modeling uncertainties, however, we show an
example comparing waveform models from lalsimulation [39] and simulations from
the CoRe library [28]. This illustrates both the current model differences as
well as a possible application of this error budget calculation $A(f)$ and
$\phi(f)$.

We start with a high-resolution CoreDB numerical simulation waveform from
Dietrich et. al., BAM:0095 [40], which has masses $m_{1}=m_{2}=1.349998$, zero
spin, and the SLy equation of state. This equation of state has
$\Lambda_{1}=\Lambda_{2}=390.1104$ at this mass, and we generate tidal
waveform models with matching $\bm{\theta}_{\text{int}}$ parameters using
TEOBResumS [41, 42] and SEOBNRv4T [43, 44]. We calculate instantaneous
frequency with second order accurate central differences using numpy.gradient
for each model. For the numerical simulation data, we fit $f(t)$ with a
B-spline using scipy.interpolate with smoothing to reduce sub-orbital
oscillations in the derivative quantity. Another gradient gives $\dot{F}$.

We show also the equivalent model characteristics of the IMRPhenomP model [45]
with numerical tidal contributions NRTidal [46, 47] as well as the TaylorF2
[48] waveforms with post-Newtonian tides [49] using the analytic stationary
phase approximation inversion

$\displaystyle T(f)$ $\displaystyle=\phi^{\prime}(f)/(2\pi)-t_{c}^{\prime}$
(21a) $\displaystyle\psi(f)$ $\displaystyle=\phi(f)+2\pi
f\,T(f)-\phi_{c}^{\prime}$ (21b) $\displaystyle\dot{F}(F)$
$\displaystyle=\frac{1}{T^{\prime}(F)}=\frac{2\pi}{\phi^{\prime\prime}(F)}$
(21c) $\displaystyle\mathcal{A}(F)$
$\displaystyle=\frac{1}{Q}\sqrt{\dot{F}(F)}\left|\tilde{h}(F)\right|$ (21d)

where we choose $\phi_{c}^{\prime}$ and $t_{c}^{\prime}$ to set the $T$ and
$\psi$ functions to zero at our chosen $f_{c}$.

Figure 1: The characteristic functions $\dot{F}(F)$ and the time-domain strain
amplitude $\mathcal{A}(F)$ as a function of instantaneous frequency. These
functions characterize the waveform predictions of each model independent of
overall shifts in time and phase. Results are shown for a set of time-domain
models, frequency-domain models, and two resolutions of numerical simulation
data. Some nonphysical oscillation is visible for the numerical data from the
start of the simulations. All simulate a zero-spin double neutron star system
with $m_{1}=m_{2}=1.35M_{\odot}$ and the SLy equation of state, except for the
binary black hole (BBH) which uses SEOBNRv4 [50] for the same masses. The
difference in time-domain amplitude $\mathcal{A}$ for TaylorF2 at these
frequencies is because it includes only the leading-order amplitude term. The
two light-shaded NR simulation resolutions demonstrate consistent frequency-
domain amplitude and phase predictions from early tidal departure through to
merger frequencies.

Fig. 1 shows the resulting $\dot{F}$ and $\mathcal{A}$ for all models. As will
be discussed in Section 7, differences in $\dot{F}$ can be interpreted in
terms of additional energy losses as a function of frequency, for example
those that rise from tidal contributions relative to the BBH model. The
amplitude $\mathcal{A}$ is connected to the gravitational-wave luminosity at a
give frequency, with differences arising for example from changes in the
source quadrupole moment.

We note that this comparison of numerical data and semi-analytic models is
independent of any waveform alignment choices such as shifts in time and
phase; it is closely related to comparisons of orbital frequency derivative
used in numerical simulation analyses (e.g. [51]) cast in terms of waveform
characteristics. Fitting these functions directly from a combination of semi-
analytic and numerical-simulation information could therefore inform a
numerical-relativity calibrated waveform model which is independent of
hybridization choices and extends through to the post-merger phase.

For this set of models, we note that TEOBResumS terminates at its $t_{c}$ with
a lower $f_{c}$ than other approximants. For the first phase comparison, we
therefore choose TEOBResumS as our baseline and set
$f_{c}=f_{\text{peak,TEOB}}\simeq 1787$ Hz as a reference frequency which can
be chosen to be common for all models.

Figure 2: Phase differences of the predicted waveform $\tilde{h}$, relative to
the TEOBResumS model. Phase is measured relative to $f_{c}=1783$ as defined
for TEOBResumS, so all models have $T(f_{c})=0$ and $\phi(f_{c})=0$ at this
frequency. All simulations and models have the same mass and tide parameters,
except for BBH which sets the tides to 0.

We compute $T(f)$ and $\psi(f)$ relative to the common reference frequency and
regenerate the Fourier-domain phase. Although smoothing was require to reduce
differentiation noise seen in Fig. 1, the integrated $\psi(f)$ relative to the
reference frequency can be read directly from numerical data. This gives from
each waveform the amplitude and phase of a Fourier-domain signal
$\tilde{h}(f)$ with a consistent $t_{c}$ and $\phi_{c}$, using the stationary
phase approximation of Sec. 3. Once the signal-domain functions are defined,
we can compute relative amplitude error $\delta A$ and the phase shift
$\delta\psi$ relative to the baseline TEOBResumS model that set $f_{c}$. The
phase error resulting from waveform differences is relative to the reference
frequency, and naturally goes to zero at that frequency. Reults are shown in
Fig. 2.

To determine measurement impact here, we calculate the residual phase
differences that must remain after marginalization over time and phase. We
calculate residual phase differences by removing a weighted least squares fit
of $\delta\phi(f)$ to $\phi_{0}+2\pi ft_{0}$ with variance
$S_{n}(f)/A^{2}(f)$, and show the resulting differences with the “APLUS”
$S_{n}$ [52] in Fig. 3. The residual phase error is not very sensitive to the
specific ground-based detector noise spectrum used to fit $\phi_{0}$ and
$t_{0}$; very similar residual phase is found using $S_{n}(f)$ from other
ground-based spectra. As a cross-check, recalling that most of the signal-to-
noise accumulates through the early inspiral [53] and would “anchor” the
$t_{c}$ and $\phi_{c}$ used to describe the models, we verified that similar
residual error can also be found by subtracting a linear fit to $\delta\phi$
in the “bucket” of ground-based detector sensitivities, between 50 and 150 Hz.

Figure 3: Relative signal amplitude differences $\delta A(f)$ and residual
signal phase differences $\delta\phi_{\text{res}}(f)$ of waveform models
relative to the reference model TEOBResumS. Residual phase is that remaining
after removing a least-squares fit of $\delta\phi(f)=\phi_{0}+2\pi ft_{0}$
weighted by the expected variance $\propto A(f)^{2}/S_{n}(f)$ in A+; the
choice of ground-based detector has minimal impact on the result. Numerical
simulations are too short to compute residual phase; phase differences for
both resolutions from TEOBResumS with $f_{c}=700$ Hz are shown for reference.
The estimated equivalent contribution of noise fluctuations in current and
next-generation detectors for a signal at $d_{\text{eff}}=100$ Mpc is shown
following discussion of Section 5.

Calibration envelopes from O4 [21] give $\delta A$ better than 11.29% and
$\delta\phi$ better than 9.18 deg or 0.160 radians from 20-2000 Hz. Comparing
the results of Fig. 3, we see that waveform uncertainties for these neutron-
star models will be smaller than calibration uncertainties at low frequency.
Current waveform differences for $\delta\phi$ between TEOBResumS and SEOBNRv4
may exceed calibration error above 500 Hz and between IMRPhenomD_NRTidalv2 and
TEOBResumS above 1 kHz. Waveform differences from $\delta A$ remain below
calibration uncertainty past 1 kHz for tidal models that go beyond leading
order.

##  5 Detector-dependent accuracy requirements for calibration and modeling

For our reference $\delta A(f)$ and $\delta\phi(f)$, we can compute the
resulting waveform differences and estimate whether they have the potential to
affect measurements in a detector with a given noise power spectral density
$S_{n}(f)$ and signal amplitude $A(f)$. This motivates frequency-dependent
modeling requirements to bring systematic uncertainties below the level of
frequency-dependent statistical fluctuations implied by a future detector’s
power spectral density $S_{n}(f)$ for a given signal. These criteria will
apply equally well to the $\delta A(f)$ and $\delta\phi(f)$ of calibration
envelopes, similarly applied to the amplitude of a specific signal.

We choose three reference noise PSDs in this work. First, we denote as
“17.08.17” the LIGO Livingston power spectral density during the observation
of GW170817 [5], as documented in [54], to set the scale of background noise
in 2017, and compare with assessments of waveform error in determining
GW170817’s source properties. Second, we show as “APLUS” the A+ detector
design of [52], as a reference for the goal sensitivities of
LIGO/Virgo/Kagra’s O5 observing run in the later 2020s [55, 24, 25]. Finally,
as a reference for what might be anticipated for next-generation
gravitational-wave observatories taking data in the 2030s, we show results
using a wide-band Cosmic Explorer noise configuration “CE” [56, 57, 58], which
has comparable sensitivity to the Einstein Telescope [59].

Figure 4: Waveform differences of semi-analytic waveform models from phase
differences relative to TEOBResumS for a signal at $d_{\text{eff}}=100\,$Mpc.
This uses the same $\delta A(f)$ and $\delta\phi(f)$ as Fig. 3, with its each
effect propagated through to $|\delta h(f)|$, including for $\delta\phi(f)$
the motion in and out of phase in $A(f)\sqrt{2(1-\cos\delta\phi(f))}$ which
gives the bumpy structure. A similar plot style was used in [60] to compare
waveform differences to detector noise. A ratio of the noise and signal
quantities squared is integrated to determine whether waveforms are
‘distinguishable’ by the criteria of Ref. [14].

Recall that the signal coming from the difference between two waveforms is
$\delta\bm{h}=\bm{h_{1}}-\bm{h_{2}}$. In general, we are interested in the
magnitude of the waveform difference
$\sqrt{\left<\delta\bm{h},\delta\bm{h}\right>}$, which is the signal-to-noise
ratio of the waveform differences. For a relative difference in amplitude
$\delta A(f)$, this is given by

$\left<\delta\bm{h},\delta\bm{h}\right>=4\int_{0}^{f_{\text{max}}}df\,\frac{|\delta\tilde{h}(f)|^{2}}{S_{n}(f)}=4\int_{0}^{f_{\text{max}}}df\,\frac{A^{2}(f)\delta
A^{2}(f)}{S_{n}(f)}$ (22)

and for a frequency-dependent phase difference $\delta\phi(f)$ we have

$\displaystyle\left<\delta\bm{h},\delta\bm{h}\right>$
$\displaystyle=4\int_{0}^{f_{\text{max}}}df\,\frac{A^{2}(f)\left|1-\exp{\left(i\delta\phi(f)\right)}\right|^{2}}{S_{n}(f)}$
(23a)
$\displaystyle=\int_{0}^{f_{\text{max}}}df\,\frac{2A^{2}(f)\left(1-\cos\delta\phi\right)}{S_{n}(f)}$
(23b)

Requirements for $\delta\phi(f)$ and $\delta A(f)$ are set from a particular
level of background noise fluctuations $S_{n}(f)$. The frequency-bin variance
$fS_{n}(f)$ comes from the power spectral density $S_{n}(f)$. The
characteristic strain squared of the signal in the frequency bin is
$4f^{2}|\delta\tilde{h}(f)|^{2}$, and is defined so that the signal-to-noise
of the waveform difference $\left<\delta\bm{h}|\delta\bm{h}\right>$ will be
less than one, and the waveforms are indistinguishable as defined in Section
2, if the characteristic strain $4f^{2}|\delta\tilde{h}(f)|^{2}<fS_{n}(f)$ at
all frequencies - see [61, 62] for explanations of characteristic strain.

We show this comparison recast in terms of the more commonly plotted amplitude
spectral density $\sqrt{S_{n}(f)}$ and the effective strain amplitude
$2\sqrt{f}|\delta\tilde{h}|$ in Figure 4 for a TEOBResumS reference signal at
$d_{\text{eff}}=100\,$Mpc. All waveform parameters are set following the
reference numerical simulation of Section 4. While the IMRPhenomPv4_NRT
waveform model remains nearly indistinguishable from TEOBResumS even for A+
sensitivites, other waveform models have differences with a potential impact
on the likelihood. However, as explored in [15], the actual impact on the
inferred parameter posterior distributions depends additionally on whether or
not the model differences (or calibration errors) are orthogonal to the
waveform differences introduced by a parameter variation. The result is
similar to the waveform difference estimates of [60], but uses the maximum-
likelihood phase and amplitude method of Section 3. In this work, we can
recover a similar estimate of systematic uncertainty for state-of-the-art
analytic model variants that do not share a common $f_{\text{peak}}$.

For $\delta A(f)$, and for sufficiently small $\delta\phi(f)$ that
$\sqrt{2(1-\cos\delta\phi)}\simeq\delta\phi$, this comparison motivates a goal
for the allowable waveform and calibration errors by again rearranging the
frequency-by-frequency indistinguishability criterion to

$\displaystyle|\delta A(f)|$
$\displaystyle<\frac{\sqrt{S_{n}(f)}}{2A(f)\sqrt{f}}$ and
$\displaystyle|\delta\phi(f)|$
$\displaystyle<\frac{\sqrt{S_{n}(f)}}{2A(f)\sqrt{f}}$ (24)

For both waveform and calibration uncertainties, then we have a goal
frequency- dependent tolerance that is set by the noise spectral density and
waveform amplitude. We show the frequency-dependent indistinguishable error
level for $\delta\phi(f)$ or $\delta A(f)$ for the three reference detector
$S_{n}(f)$ and our reference neutron-star $A(f)$ at $d_{\text{eff}}=100\,$Mpc
in Fig. 5. These are translated into the shaded detector bounds in Figs. 3.

Figure 5: The size of of $\delta\phi$ (in radians) or $\delta A$ (as a
fraction) that gives signal modifications with the same characteristic strain
as that from the power spectral density for the labeled detectors, for an
incident BNS signal from $d_{\text{eff}}=100$ Mpc. The model uncertainty and
calibration requirements set by the power spectral density are specific to the
signal under consideration; the rapid rise at 2 kHz for CE reflects the
termination frequency of the inspiral model being applied here. The important
region for measurement and calibration will be set by the frequency range of
the desired signal information, as could be determined for example in Figure 2
of Ref. [53]. For signals of this amplitude, systematic uncertainties from
calibration and waveform modeling seem likely to dominate the error budget at
$\sim 10-100$ Hz in next-generation observatories.

Figure 5 shows that waveform and calibration error of $\sim 3\%$ in amplitude
(2 Deg of phase) through the tens of Hz range would be required for
calibration to remain subdominant to stochastic noise for a single-detector
$d_{\text{eff}}=100$ Mpc BNS with signal-to-noise ratio approximately 60 in
the A+ detector.

In Cosmic Explorer, the same system would have a signal-to-noise ratio of
$\simeq$1158 above 10 Hz and require waveform and calibration errors as small
as $10^{-3}$ for the impact of those uncertainties to be smaller than the
statistical uncertainty.

##  6 Range of waveform uncertainties for GW170817

Recent discussions have raised the question of whether the analysis of
gravitational-wave signals should be truncated at a specific frequency to
limit systematic uncertainty in the inference of tidal properties, as those
increase due to modeling differences at high frequencies (see for example [63,
64]). We propose instead that a more appropriate treatment of increasing
waveform uncertainty is to explicitly marginalize over the high-frequency
$\delta\phi(f)$ range implied by those model differences. Analytic
marginalization procedures or likelihood reweighting [20, 65] could make this
more tractable.

As a first estimate, we characterize the behavior $\delta A(f)$ and
$\delta\phi(f)$ for GW170817. Waveform modeling uncertainty is illustrated by
generating an equivalent of the results of Fig. 3 for $\delta A(f)$ and
$\delta\phi(f)$ between two waveform families (here,
IMRPhenomPv2_NRTidalv2[66, 67, 40] and SEOBNRv4T_surrogate[68]) at parameter
values drawn from the the posterior distributions of the bilby samples for
GW170817 [11]. Waveform amplitudes are calculated using the sampled distance,
sky position and inclination, and projected on to the LIGO Livingston (L1)
detector antenna pattern at the signal GPS time using a slightly modified
version of pesummary [69]. Residual $\delta\phi$ is calculated as described in
Section 3, by subtracting of the phase contribution of the maximum likelihood
value for an overall shift in time and phase describing the waveform
differences. To compare the effects of spin modeling and tidal modeling, we
first calculate waveform differences with spin set to zero, and then with
tides set to zero. The median $\delta A(f)$ and $\delta\phi(f)$ in both cases,
and the distribution over the tidal posterior range, are shown in Figure 6.

Figure 6:  Left: waveform amplitude and phase differences over the posterior
region of parameters in the bilby samples for GW170817 [11]. Differences are
calculated between IMRPhenomPv2_NRTidalv2[66, 67, 40] and
SEOBNRv4T_surrogate[68] and projected onto LIGO Livingston (L1). The solid red
line is the median distribution without spin terms, and 1-$\sigma$, and
2-$\sigma$ ranges around that median are shaded in grey. The dashed red lines
show median waveform differences with spin but without tidal terms. The
relatively small GWTC-1 envelopes for L1 calibration uncertainty during
GW170817 [70] are shown in orange. Right: to compare with realistic waveform
uncertainties, the distribution of amplitude and phase differences at L1
between two random draws from the posterior distribution is shown. The PSD-
based estimate of 1-$\sigma$ error from L1 noise fluctuations at the time of
GW170817 [71], as seen in the shaded region of Figure 3, is outlined in blue
in all panels.

Calibration envelopes for GW170817 from GWTC-1 [70] are also shown on the same
scale in Figure 6. Calibration uncertainties are larger than waveform
uncertainties at low frequencies, but both are small (sub 2 $\deg$ or 3%) at
low frequencies. In this comparison, waveform phase uncertainties become
larger than calibration uncertainty starting between 100 and 200 Hz, while
waveform amplitude uncertainties remain lower than calibration uncertainty up
to 700 Hz.

To compare the waveform systematics to realistic uncertainties, we also show
the distribution of $\delta A(f)$ and $\delta\phi(f)$ within the inferred
posterior distribution of waveforms from the same bilby samples. Two waveforms
are generated from two random parameter draws, both with the
IMRPhenomPv2_NRTidalv2 model that was used for inference. The waveforms are
again projected on to the L1 observatory. Here, the $t_{c}$ and $\phi_{c}$
parameters fix the full phase difference between the waveforms; each is
generated with its sample’s $\phi_{c}$, and the effect of $t_{c}$ on the
Fourier-domain phase is incorporated by subtracting $2\pi f\Delta t_{c}$ using
the difference in coalescence time between the samples.

As a comparison, a PSD-based estimate of the 1-$\sigma$ uncertainty from
statistical fluctuations is also shown, using the L1 noise curve [54] and
following the method of Section 5. We note that the PSD-based estimates of
1-$\sigma$ waveform errors are an imperfect estimate of realistic
uncertainties. The amplitude prediction shows reasonable agreement with the
most constrained frequencies in the posterior. The waveform model has limited
flexibility in amplitude before the merger, so it can not saturate the
statistical range at all frequencies. However, phase differences within the
posterior distribution have a 1-$\sigma$ level that is somewhat larger than
the PSD-based estimate.

The distribution of errors within the posterior distribution can also be
compared to the credible intervals for $\delta A$ and $\delta\phi$ found in
Edelman et al [10]. Their recovered-spline distribution for $\delta\phi$ is
flatter over frequencies at the 1-$\sigma$ level, but at the 2-$\sigma$ level
the credible intervals increase with increasing frequency (as seen in their
Figures 5 and 12) in a similar way as the posterior- and PSD-based estimates
shown here.

The waveform differences shown here come from a combination of different spin
models, different tidal models, and differences rising from a re-summed
effective-one-body waveform model vs. a frequency-domain phenomenological
waveform generation. The range of tides and the presence of spin each affect
the relative differences and their frequency dependence. The results confirm
that systematic error for neutron-star inference can be reduced by restricting
analysis to frequencies below $\sim 1$kHz, as $\delta\phi$ are smaller than
the statistical background in that range. They also confirm that amplitude
uncertainty is not significant for GW170817.

However, this type of waveform comparison could be used to generate a
distribution model for $\delta\phi$ and $\delta A$ that could be applied for
coherent frequency-dependent waveform marginalization in proposed future event
analyses. With this in mind, we computed the distribution of $\delta\phi$ and
$\delta A$ over the full range of prior samples with the
IMRPhenomPv2NRT_lowSpin_prior in the GW170817 PE release of GWTC-1 [23], using
optimized versions of the waveform models IMRPhenomPv2_NRTidal_INTERP [72] and
SEOBNRv4T_surrogate [68]. Even this low-spin prior range of parameters gives
gives very wide distribution of both $\delta\phi(f)$ and $\delta A(f)$ at
several hundred Hz; the 2-$\sigma$ range exceeds 50 Deg/% between 200 Hz and 1
kHz. Most systems in the prior have poorly-defined contributions above $\sim
1$ kHz as the stars have merged at low frequencies due to large size and tides
at larger values of the component $\Lambda$s.

A realistic implementation of marginalization would therefore require a
prescription for the distribution of $\delta\phi(f)$ that varies over the
intrinsic parameters. Such a marginalization could also reflect the increase
in modeling uncertainties away from well-calibrated regions of mass ratio.
While such a distribution could first be estimated by comparing multiple
waveform models as is done here, it may also be better derived from
theoretical estimates of intrinsic model uncertainties as discussed in Section
3.

##  7 Astrophysical interpretation of $\delta\phi$ and $\delta A$

Considering the size of the energy transfers required to modify the
characteristic functions $\dot{F}(F)$ and $\mathcal{A}(F)$ will allow us to
generally interpret constraints on $\delta\phi(f)$ and $\delta A(f)$, as might
be recovered from a coherent analysis like that of Edelman et. al. [10].

Consider the usual energy balance equation for the evolution of a source
system, with a total system energy $E(F)$ (for orbits: negative relative to
$F\to 0$) that can be written as a function of the instantaneous
gravitational-wave frequency $F$. For example, two masses $m_{1}$ and $m_{2}$
in circular orbits emitting gravitational waves at $F=\Omega_{\text{orb}}/\pi$
for orbital angular frequency $\Omega_{\text{orb}}$ have a leading order
orbital binding energy of

$E_{\text{orb}}=-\frac{m_{1}m_{2}}{m_{1}+m_{2}}\frac{c^{2}}{2}\left(\frac{G(m_{1}+m_{2})(\pi
F)}{c^{3}}\right)^{2/3}$ (25)

where the masses combine to a reduced mass and the term to the power of 2/3 is
a dimensionless frequency parameter often denoted in post-Newtonian expansions
as $x$ [73] . The gravitational-wave energy loss or luminosity
$\mathcal{L}_{\text{GW}}$ from the same orbit is

$\mathcal{L}_{\text{GW}}=\frac{32}{5}\frac{c^{5}}{G}\frac{m_{1}^{2}m_{2}^{2}}{\left(m_{1}+m_{2}\right)^{4}}\left(\frac{G(m_{1}+m_{2})(\pi
F)}{c^{3}}\right)^{10/3}$ (26)

using the same dimensionless frequency term.

Orbital evolution is driven by the energy balance at each characteristic
emission frequency, which results in the characteristic $\dot{F}$ function in
terms of the energy gradient $E^{\prime}(F)=dE(F)/dF$ and the total rate of
energy loss or luminosity $\mathcal{L}(F)$:

$\displaystyle\frac{dE(F)}{dt}$ $\displaystyle=-\mathcal{L}(F)$ and
$\displaystyle\dot{F}$ $\displaystyle=-\frac{\mathcal{L}(F)}{E^{\prime}(F)}$
(27a)

The gravitational-wave luminosity is determined by the multipole decomposition
of strain emitted by the source [74, 75]:

$\displaystyle\mathcal{L}_{\text{GW}}$ $\displaystyle=\frac{1}{4\pi}\int
d\Omega\sum_{\ell m}Y_{-2}^{\ell m}(\Omega)\left|d\dot{h}_{\ell
m}(f)\right|^{2}$ (28) $\displaystyle=\frac{1}{16\pi}\sum_{\ell
m}d^{2}\left(\dot{\mathcal{A}}_{\ell m}^{2}+(2\pi)^{2}\mathcal{A}_{\ell
m}^{2}F^{2}\right)$ (29)

where we use again $F=\dot{\psi}/2\pi$ and where $d$ the luminosity distance
at which $h$ is measured.

When the SPA applies (following Eq. 12) we can neglect the
$\dot{\mathcal{A}^{2}}$ term, so gravitational-wave luminosity in each mode
determines the waveform amplitude through

$\mathcal{A}^{2}=\frac{4}{\pi}\frac{\mathcal{L}_{\text{GW}}}{d^{2}F^{2}}$ (30)

The characteristic time-domain amplitude is directly connected to the total
luminosity of the sytem.

The energy balance implied by a waveform model could be modified by any
unmodeled physical effect which increases gravitational-wave luminosity
$\mathcal{L}_{\text{GW}}(1+\delta\mathcal{L}_{\text{GW}})$, such as the
coherent excess quadrupole moment induced by tides on each component neutron
star. Unmodeled non-gravitational-wave energy losses (including neutrino,
electromagnetic,…) could also contribute to the system’s evolution through an
additional $\mathcal{L}_{\text{tot}}(1+\delta\mathcal{L}_{\text{MM}})$ that
drives the energy balance without affecting the gravitational-wave amplitude.
Finally, models derived from or calibrated to numerical relativity simulations
might include unphysical excess energy loss due to numerical dissipation. The
true waveform would then have
$\mathcal{L}_{\text{tot}}(1-\delta\mathcal{L}_{\text{num}})$; this is
equivalent to a negative $\delta\mathcal{L}_{\text{MM}}$ needed to generate
the true model.

###  7.1 Internal energy transfers: adiabatic vs. dynamic

The waveform model may also neglect effects that change the system energy
associated with a characteristic gravitational-wave emission frequency, such
as high-order corrections to the orbital binding energy or a missing energy
reservoir. For example, there may be changes in the energy stored in a
component neutron star’s internal modes.

Consider some effect that changes total energy $E$ as a function of $F$, for
example changing the leading order energy of Eq. 25 to $E(1+\delta E)$. A
total derivative with respect to $F$ gives the correction to the energy
gradient term in the energy balance equation, $dE/dF=E^{\prime}$, as

$E^{\prime}(1+\delta E^{\prime})=E^{\prime}\left(1+\delta
E+\frac{E}{E^{\prime}}\left(\delta E\right)^{\prime}\right)$ (31)

Drawing from the extensive discussion of potential energy transfer impacts in
the literature (such as in [76, 77, 78, 44]), we interpret the first term as a
correction from adiabatic effects: those where the change in $\delta E$ with
frequency is slow. In this case, noting each $\delta$ term here represents the
size of _relative_ corrections to the original system $E^{\prime}$ and $E$,
the correction $\delta E^{\prime}_{A}=\delta E_{A}$. For adiabatic tides, the
correction $\delta E_{A}$ is negative [77]: energy going into the stellar
deformation implies a less negative total energy at a given frequency.

In contrast, the second term characterizes dynamic effects, coming from a
rapid energy change. For example, consider a resonant transfer of energy
$\Delta E_{D}$ from the orbit to a stellar mode, starting at a characteristic
gravitational-wave frequency $F_{0}$ over a bandwidth $\Delta F$. Around
$F=F_{0}$, we have a change in the gradient

$\displaystyle E^{\prime}(1+\delta E^{\prime})$ $\displaystyle\simeq
E^{\prime}-\frac{\Delta E_{D}}{\Delta F}$ (32)
$\displaystyle=E^{\prime}\left(1-\frac{E}{E^{\prime}}\frac{\delta
E_{D}}{\Delta F}\right)$ (33)

where we introduce the relative energy transfer $\delta E_{D}=\Delta
E/E_{\text{orb}}$.

Define the energy transfer timescale $t_{D}=\Delta F/\dot{F}$ and the
adiabatic decay timescale or gravitational-wave emission timescale
$t_{\text{A}}=|E/\mathcal{L}|\simeq 2F/3\dot{F}$. Using the energy balance
relationship $E^{\prime}=-\mathcal{L}/\dot{F}$, we find as has previously been
shown [78] that the impact of a dynamic energy transfer compared to adiabatic
energy transfer is amplified:

$\delta E^{\prime}_{D}=-\frac{E}{E^{\prime}}\frac{\delta E_{D}}{\Delta
F}=\frac{E}{\mathcal{L}}\frac{\dot{F}}{\Delta F}\delta
E_{D}\simeq\frac{3}{2}\frac{t_{A}}{t_{D}}\delta{E_{D}}$ (34)

Generally, then, internal energy transfers lead to a gradient term

$\delta E^{\prime}\sim\delta{E_{A}}+\frac{t_{A}}{t_{D}}\delta{E_{D}}$ (35)

where $\delta{E_{A}}E_{\text{orb}}$ and $\delta{E_{D}}E_{\text{orb}}$
characterize the amounts of energy transferred adiabatically or dynamically.

###  7.2 Signal implications

Assume that all unmodeled effects are very small corrections to the baseline
energy and luminosity, so that their effects can be linearized. We find the
impact on the characteristic functions

$\displaystyle\mathcal{A}_{\text{true}}$
$\displaystyle=\mathcal{A}_{\text{model}}\left(1+\delta\mathcal{L}_{\text{GW}}\right)$
(36a) $\displaystyle\dot{F}_{\text{true}}$
$\displaystyle=\dot{F}_{\text{model}}\left(1-\delta
E^{\prime}+\delta\mathcal{L}_{\text{GW}}+\delta\mathcal{L}_{\text{MM}}\right)$
(36b) $\displaystyle T^{\prime}_{\text{true}}$
$\displaystyle=T^{\prime}+T^{\prime}\left(\delta
E^{\prime}-\delta\mathcal{L}_{\text{GW}}-\delta\mathcal{L}_{\text{MM}}\right)$
(36c)

This is consistent with tidally driven increases in the $\dot{F}$ function
that can be seen at moderate frequencies for all binary neutron star models
relative to the binary black hole model in Fig. 1.

To confirm the interpretation of these energy transfers, we work out the
impact on the time-domain waveform following Eqs. 15 and using the energy
balance relation $T^{\prime}=-E^{\prime}/\mathcal{L}$,

$\displaystyle\delta T(f)$
$\displaystyle=\int_{f}^{f_{c}}dF\frac{E^{\prime}}{\mathcal{L}}\left(-\delta
E^{\prime}+\delta\mathcal{L}_{\text{GW}}+\delta\mathcal{L}_{\text{MM}}\right)$
(37a) $\displaystyle\delta\psi(f)$
$\displaystyle=2\pi\int_{f}^{f_{c}}dFF\frac{E^{\prime}}{\mathcal{L}}\left(-\delta
E^{\prime}+\delta\mathcal{L}_{\text{GW}}+\delta\mathcal{L}_{\text{MM}}\right)$
(37b)

To interpret them physically: since orbital $E^{\prime}<0$, this shows that
both the time to coalescence and the number of cycles before coalescence can
be shortened by additional luminosity or by a less rapid decrease of total
energy of the system with increasing $F$.

After propagation through to the signal domain, again linearizing the
different corrections, we will find from Eqs. 19 the energetics implications
for the frequency-domain signal:

$\displaystyle\delta A(f)$ $\displaystyle=\frac{1}{2}\left(\delta
E^{\prime}-\delta\mathcal{L}_{\text{GW}}-\delta\mathcal{L}_{\text{MM}}\right)+\delta\mathcal{L}_{\text{GW}}$
(38a) $\displaystyle=\frac{1}{2}\left(\delta
E^{\prime}+\delta\mathcal{L}_{\text{GW}}-\delta\mathcal{L}_{\text{MM}}\right)$
(38b) $\displaystyle\delta\phi(f)$ $\displaystyle=2\pi f\delta
T(f)-\delta\psi(f)$ (38c)
$\displaystyle=-2\pi\int_{f}^{f_{c}}d\tilde{f}\int_{\tilde{f}}^{f_{c}}dF\,\frac{E^{\prime}}{\mathcal{L}}\left(\delta
E^{\prime}-\delta\mathcal{L}_{\text{GW}}-\delta\mathcal{L}_{\text{MM}}\right)$
(38d)

Note that any energy losses that are not in gravitational waves lead to a
decrease in gravitational-wave amplitude at the corresponding frequency - the
orbit sweeps through more quickly and fewer cycles contribute. Phase
accumulates more rapidly with any additional energy losses from the system.
The relative change in $\delta A$ vs $\delta\phi$ depends on the type of
energy transfer that is added to the source model.

###  7.3 Example: dynamic energy transfer and GW170817

Applying the energetics framework above allows the interpretation of Edelman
et al [10] results for GW170817 in terms of a possible augmentation of the
inspiral waveform due to neutron-star energy transfers. While no confident
identification of a departure from the signal model was identified, the
posterior distribution of $\delta\phi(f)$ in Figure 12 in [10] shows a
1-$\sigma$ excess at around $60$ Hz of perhaps 5 Deg or order 0.1 radians. At
the same frequency, no significant $\delta A$ change is observed.

One candidate for a short-duration $\delta\phi(f)$ at a specific inspiral
frequency is dynamic energy transfer into a neutron-star mode, for example as
explored in [79, 80]. As the net impact on the Fourier-domain phase comes from
a difference between $2\pi f\delta T$ and $\delta\phi$ in Eq. 38c, a
completely instantaneous energy loss only affects amplitude at the specific
emission frequency. Consider instead the transfer of $\Delta E$ at $F_{0}$
over a small range $\Delta F$. If the excitation does not induce any coherent
quadrupole, we have a corresponding short-duration relative decrease in the
signal amplitude from Eq. 34 of

$\delta A(F_{0})=\frac{1}{2}\delta E^{\prime}=\frac{1}{2}\frac{\Delta
E}{t_{D}}\frac{1}{\mathcal{L}}=\frac{1}{2}\frac{t_{A}}{t_{D}}\delta E_{D}$
(39)

using the timescale notation of Sec. 7.1. If amplitude variation is not
observed, this limits the total unmodeled energy transfer rate to be small
compared to the luminosity at 60 Hz, which is $\sim 2\times
10^{-4}M_{\odot}c^{2}$/s or $-3\times 10^{43}$ joules/s at 60 Hz. At the same
time, there is an increase in accumulated Fourier-domain phase $\delta\phi$
which will be nonzero for frequencies between $F_{0}$ and $F_{0}+\Delta F$. We
estimate the nonzero values as

$\displaystyle\delta\phi$ $\displaystyle\simeq 2\pi\delta
E^{\prime}\int_{F_{0}}^{F_{0}+\Delta
F}d\tilde{F}\int_{\tilde{F}}^{F_{0}+\Delta F}dF\,T^{\prime}(F)$ (40)
$\displaystyle\simeq\pi\delta E^{\prime}\,\frac{\left(\Delta
F\right)^{2}}{\dot{F}}\simeq\pi\Delta
E\frac{\dot{F}}{\mathcal{L}}t_{D}=\frac{2\pi}{3}t_{D}F\delta E_{D}$ (41)

where $t_{D}F$ is the number of cycles over which the energy is transferred.

A resonant energy transfer has a dynamic timescale of
$t_{D}\sim\sqrt{t_{A}/F}$ [78], so that both $t_{A}/t_{D}$ and
$t_{D}F\sim\sqrt{t_{A}F}$. At 60 Hz, $\sqrt{t_{A}F}\simeq 40$, and the orbital
binding energy $E\simeq-0.006M_{\odot}c^{2}$ or $-1.1\times 10^{45}$ joules,
and a energy transfer of $\delta E=\Delta E/E\sim 0.001$ at 60 Hz corresponds
to roughly $\delta A=2\%$ and $\delta\phi=5$ Deg, compatible with the credible
intervals recovered for GW170817 [10]. We leave more realistic explorations of
astrophysical interpretations to future work.

##  8 Conclusion

An amplitude-phase decomposition of time-domain gravitational-wave emission
enables the use of instantaneous frequency in characterizing the underlying
model physics in terms of characteristic functions for the signal evolution
$\dot{F}$ and the time-domain amplitude $\mathcal{A}$. Resulting model
uncertainties can be propagated from the time-domain wave to the frequency-
domain signal for comparison with gravitational-wave observations. The
differences in model waveform predictions that propagate from the underlying
physics can be decoupled from the time and phase constants $t_{c}$ and
$\phi_{c}$ that characterize a specific observation by finding maximum-
likelihood values. Writing waveform uncertainties in the same form as
calibration uncertainties — namely, as frequency-dependent signal amplitude
and phase error terms $\delta A(f)$ and $\delta\phi(f)$ — would then
characterize model-dependent waveform uncertainty suitable for marginalization
in gravitational-wave inference.

Especially for observations of early inspiral, energy-balance arguments can
aid in the interpretation of coherent waveform deviations for binary signals.
We have illustrated how a common astrophysical energy transfer scenarios can
be used to interpret recovered values of $\delta\phi(f)$ and $\delta A(f)$
compared to a waveform model used for analysis; recovered bounds on these
functions limit the size of unmodeled energy transfers in the source system.

##  9 Acknowledgements

Many thanks to the IGWN Conda Distribution, numpy [81], pandas [82], watpy
[28], lalsimulation [39], pycbc [72], bilby [83], pesummary [69], and the work
of everyone that supports computational infrastructure for the LIGO-Virgo-
Kagra collaborations and the gravitational-wave community. I thank David
Radice, Wynn Ho, Sanjay Reddy, Josh Smith, Nils Andersson and the Caltech LIGO
group for helpful discussion. Thanks also to Aaron Zimmerman for useful and
detailed feedback and corrections to Section 7, and the anonymous reviewers
for especially helpful suggestions for Section 6. Read was supported during
the development of these ideas by funding from NSF PHY-1307545 and
PHY-1806962, the LIGO Lab, the Carnegie Observatories, and the Nicholas and
Lee Begovich Center for Gravitational-Wave Astronomy. This material is based
upon work supported by NSF’s LIGO Laboratory which is a major facility fully
funded by the National Science Foundation.

## References

  * [1] Nina Kunert, Peter T. H. Pang, Ingo Tews, Michael W. Coughlin, and Tim Dietrich.  Quantifying modeling uncertainties when combining multiple gravitational-wave detections from binary neutron star sources.  Physical Review D, 105(6), mar 2022.
  * [2] Rossella Gamba, Matteo Breschi, Sebastiano Bernuzzi, Michalis Agathos, and Alessandro Nagar.  Waveform systematics in the gravitational-wave inference of tidal parameters and equation of state from binary neutron-star signals.  Physical Review D, 103(12), jun 2021.
  * [3] Yiwen Huang, Carl-Johan Haster, Salvatore Vitale, Vijay Varma, Francois Foucart, and Sylvia Biscoveanu.  Statistical and systematic uncertainties in extracting the source properties of neutron star-black hole binaries with gravitational waves.  Phys. Rev. D, 103(8):083001, April 2021.
  * [4] B. P. Abbott, LIGO Scientific Collaboration, and Virgo Collaboration.  Effects of waveform model systematics on the interpretation of GW150914.  Classical and Quantum Gravity, 34(10):104002, May 2017.
  * [5] B. P. Abbott, LIGO Scientific Collaboration, and Virgo Collaboration.  Properties of the Binary Neutron Star Merger GW170817.  Physical Review X, 9(1):011001, January 2019.
  * [6] Gregory Ashton and Tim Dietrich.  The use of hypermodels to understand binary neutron star collisions.  Nature Astronomy, 6:961–967, July 2022.
  * [7] Gregory Ashton and Sebastian Khan.  Multiwaveform inference of gravitational waves.  Phys. Rev. D, 101(6):064037, March 2020.
  * [8] A. Z. Jan, A. B. Yelikar, J. Lange, and R. O’Shaughnessy.  Assessing and marginalizing over compact binary coalescence waveform systematics with RIFT.  Phys. Rev. D, 102(12):124069, December 2020.
  * [9] B. P. Abbott, LIGO Scientific Collaboration, and Virgo Collaboration.  GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral.  Phys. Rev. Lett., 119(16):161101, October 2017.
  * [10] Bruce Edelman, F. J. Rivera-Paleo, J. D. Merritt, Ben Farr, Zoheyr Doctor, Jeandrew Brink, Will M. Farr, Jonathan Gair, Joey Shapiro Key, Jess McIver, and Alex B. Nielsen.  Constraining unmodeled physics with compact binary mergers from GWTC-1.  Phys. Rev. D, 103(4):042004, February 2021.
  * [11] I. M. Romero-Shaw, C. Talbot, S. Biscoveanu, V. D’Emilio, G. Ashton, C. P. L. Berry, S. Coughlin, S. Galaudage, C. Hoy, M. Hübner, K. S. Phukon, M. Pitkin, M. Rizzo, N. Sarin, R. Smith, S. Stevenson, A. Vajpeyi, M. Arène, K. Athar, S. Banagiri, N. Bose, M. Carney, K. Chatziioannou, J. A. Clark, M. Colleoni, R. Cotesta, B. Edelman, H. Estellés, C. García-Quirós, Abhirup Ghosh, R. Green, C. J. Haster, S. Husa, D. Keitel, A. X. Kim, F. Hernandez-Vivanco, I. Magaña Hernandez, C. Karathanasis, P. D. Lasky, N. De Lillo, M. E. Lower, D. Macleod, M. Mateu-Lucena, A. Miller, M. Millhouse, S. Morisaki, S. H. Oh, S. Ossokine, E. Payne, J. Powell, G. Pratten, M. Pürrer, A. Ramos-Buades, V. Raymond, E. Thrane, J. Veitch, D. Williams, M. J. Williams, and L. Xiao.  Bayesian inference for compact binary coalescences with BILBY: validation and application to the first LIGO-Virgo gravitational-wave transient catalogue.  MNRAS, 499(3):3295–3319, December 2020.
  * [12] J. Veitch, V. Raymond, B. Farr, W. Farr, P. Graff, S. Vitale, B. Aylott, K. Blackburn, N. Christensen, M. Coughlin, W. Del Pozzo, F. Feroz, J. Gair, C. J. Haster, V. Kalogera, T. Littenberg, I. Mandel, R. O’Shaughnessy, M. Pitkin, C. Rodriguez, C. Röver, T. Sidery, R. Smith, M. Van Der Sluys, A. Vecchio, W. Vousden, and L. Wade.  Parameter estimation for compact binaries with ground-based gravitational-wave observations using the LALInference software library.  Phys. Rev. D, 91(4):042003, February 2015.
  * [13] Curt Cutler and Éanna E. Flanagan.  Gravitational waves from merging compact binaries: How accurately can one extract the binary’s parameters from the inspiral waveform\?  Phys. Rev. D, 49(6):2658–2697, March 1994.
  * [14] Lee Lindblom, Benjamin J. Owen, and Duncan A. Brown.  Model waveform accuracy standards for gravitational wave data analysis.  Phys. Rev. D, 78(12):124020, December 2008.
  * [15] Michael Pürrer and Carl-Johan Haster.  Gravitational waveform accuracy requirements for future ground-based detectors.  Physical Review Research, 2(2):023151, May 2020.
  * [16] Will M. Farr, Benjamin Farr, and Tyson M. Littenberg.  Modelling calibration errors in cbc waveforms.  Technical Report T1400682-v1, LIGO, 2014.
  * [17] Ethan Payne, Colm Talbot, Paul D. Lasky, Eric Thrane, and Jeffrey S. Kissel.  Gravitational-wave astronomy with a physical calibration model.  Phys. Rev. D, 102:122004, Dec 2020.
  * [18] Salvatore Vitale, Carl-Johan Haster, Ling Sun, Ben Farr, Evan Goetz, Jeff Kissel, and Craig Cahillane.  Physical approach to the marginalization of LIGO calibration uncertainties.  Phys. Rev. D, 103(6):063016, March 2021.
  * [19] Salvatore Vitale, Walter Del Pozzo, Tjonnie G. F. Li, Chris Van Den Broeck, Ilya Mandel, Ben Aylott, and John Veitch.  Effect of calibration errors on Bayesian parameter estimation for gravitational wave signals from inspiral binary systems in the advanced detectors era.  Phys. Rev. D, 85(6):064034, March 2012.
  * [20] Reed Essick.  Calibration uncertainty’s impact on gravitational-wave observations.  Phys. Rev. D, 105(8):082002, April 2022.
  * [21] Ling Sun, Evan Goetz, Jeffrey S Kissel, Joseph Betzwieser, Sudarshan Karki, Aaron Viets, Madeline Wade, Dripta Bhattacharjee, Vladimir Bossilkov, Pep B Covas, Laurence E H Datrier, Rachel Gray, Shivaraj Kandhasamy, Yannick K Lecoeuche, Gregory Mendell, Timesh Mistry, Ethan Payne, Richard L Savage, Alan J Weinstein, Stuart Aston, Aaron Buikema, Craig Cahillane, Jenne C Driggers, Sheila E Dwyer, Rahul Kumar, and Alexander Urban.  Characterization of systematic error in advanced ligo calibration.  Classical and Quantum Gravity, 37(22):225008, oct 2020.
  * [22] F. Acernese, et. al., and Virgo Collaboration.  Calibration of advanced Virgo and reconstruction of the detector strain h(t) during the observing run O3.  Classical and Quantum Gravity, 39(4):045006, February 2022.
  * [23] LIGO Scientific Collaboration and Virgo Collaboration.  Gwtc-1: A gravitational-wave transient catalog of compact binary mergers observed by ligo and virgo during the first and second observing runs.  Phys. Rev. X, 9:031040, Sep 2019.
  * [24] The LIGO Scientific Collaboration et al.  Advanced ligo.  Classical and Quantum Gravity, 32(7):074001, mar 2015.
  * [25] F. Acernese et al.  Advanced Virgo: a second-generation interferometric gravitational wave detector.  Class. Quant. Grav., 32(2):024001, 2015.
  * [26] Simone Mezzasoma and Nicolás Yunes.  Theory-agnostic framework for inspiral tests of general relativity with higher-harmonic gravitational waves.  Phys. Rev. D, 106(2):024026, July 2022.
  * [27] Lawrence E. Kidder.  Using full information when computing modes of post-Newtonian waveforms from inspiralling compact binaries in circular orbit.  Phys. Rev. D, 77(4):044016, February 2008.
  * [28] Alejandra Gonzalez, Francesco Zappa, Matteo Breschi, Sebastiano Bernuzzi, David Radice, Ananya Adhikari, Alessandro Camilletti, Swami Vivekanandji Chaurasia, Georgios Doulis, Surendra Padamata, Alireza Rashti, Maximiliano Ujevic, Bernd Brügmann, William Cook, Tim Dietrich, Albino Perego, Amit Poudel, and Wolfgang Tichy.  Second release of the CoRe database of binary neutron star merger waveforms.  arXiv e-prints, page arXiv:2210.16366, October 2022.
  * [29] Patricia Schmidt, Ian W. Harry, and Harald P. Pfeiffer.  Numerical relativity injection infrastructure.  Technical Report LIGO-T1500606, LIGO, 2017.
  * [30] P. Ajith, M. Boyle, D. A. Brown, S. Fairhurst, M. Hannam, I. Hinder, S. Husa, B. Krishnan, R. A. Mercer, F. Ohme, C. D. Ott, J. S. Read, L. Santamaria, and J. T. Whelan.  Data formats for numerical relativity waves.  arXiv e-prints, page arXiv:0709.0093, September 2007.
  * [31] Hsin-Yu Chen, Daniel E. Holz, John Miller, Matthew Evans, Salvatore Vitale, and Jolien Creighton.  Distance measures in gravitational-wave astrophysics and cosmology.  Classical and Quantum Gravity, 38(5):055010, March 2021.
  * [32] Sylvain Marsat and John G. Baker.  Fourier-domain modulations and delays of gravitational-wave signals.  arXiv e-prints, page arXiv:1806.10734, June 2018.
  * [33] R. Abbott et al.  Population of Merging Compact Binaries Inferred Using Gravitational Waves through GWTC-3.  Phys. Rev. X, 13(1):011048, 2023.
  * [34] Bruce Allen, Warren G. Anderson, Patrick R. Brady, Duncan A. Brown, and Jolien D. E. Creighton.  FINDCHIRP: An algorithm for detection of gravitational waves from inspiraling compact binaries.  Phys. Rev. D, 85(12):122006, June 2012.
  * [35] Serge Droz, Daniel J. Knapp, Eric Poisson, and Benjamin J. Owen.  Gravitational waves from inspiraling compact binaries: Validity of the stationary-phase approximation to the Fourier transform.  Phys. Rev. D, 59(12):124016, June 1999.
  * [36] Rossella Gamba, Sebastiano Bernuzzi, and Alessandro Nagar.  Fast, faithful, frequency-domain effective-one-body waveforms for compact binary coalescences.  Phys. Rev. D, 104(8):084058, October 2021.
  * [37] Scott A. Hughes, Niels Warburton, Gaurav Khanna, Alvin J. K. Chua, and Michael L. Katz.  Adiabatic waveforms for extreme mass-ratio inspirals via multivoice decomposition in time and frequency.  Phys. Rev. D, 103(10):104014, 2021.
  * [38] John Veitch and Walter del Pozzo.  Analytic marginalisation of phase parameter.  Technical Report T1300326-v1, LIGO, 2013.
  * [39] LIGO Scientific Collaboration.  LIGO Algorithm Library - LALSuite.  free software (GPL), 2018.
  * [40] Tim Dietrich, Sebastiano Bernuzzi, and Wolfgang Tichy.  Closed-form tidal approximants for binary neutron star gravitational waveforms constructed from high-resolution numerical relativity simulations.  Phys. Rev. D, 96(12):121501, 2017.
  * [41] Alessandro Nagar, Francesco Messina, Piero Rettegno, Donato Bini, Thibault Damour, Andrea Geralico, Sarp Akcay, and Sebastiano Bernuzzi.  Nonlinear-in-spin effects in effective-one-body waveform models of spin-aligned, inspiralling, neutron star binaries.  Phys. Rev. D, 99(4):044007, 2019.
  * [42] Sarp Akcay, Sebastiano Bernuzzi, Francesco Messina, Alessandro Nagar, Néstor Ortiz, and Piero Rettegno.  Effective-one-body multipolar waveform for tidally interacting binary neutron stars up to merger.  Phys. Rev. D, 99(4):044051, 2019.
  * [43] Jan Steinhoff, Tanja Hinderer, Tim Dietrich, and Francois Foucart.  Spin effects on neutron star fundamental-mode dynamical tides: Phenomenology and comparison to numerical simulations.  Physical Review Research, 3(3):033129, August 2021.
  * [44] Tanja Hinderer, Andrea Taracchini, Francois Foucart, Alessandra Buonanno, Jan Steinhoff, Matthew Duez, Lawrence E. Kidder, Harald P. Pfeiffer, Mark A. Scheel, Bela Szilagyi, Kenta Hotokezaka, Koutarou Kyutoku, Masaru Shibata, and Cory W. Carpenter.  Effects of Neutron-Star Dynamic Tides on Gravitational Waveforms within the Effective-One-Body Approach.  Phys. Rev. Lett., 116(18):181101, May 2016.
  * [45] Mark Hannam, Patricia Schmidt, Alejandro Bohé, Leïla Haegel, Sascha Husa, Frank Ohme, Geraint Pratten, and Michael Pürrer.  Simple model of complete precessing black-hole-binary gravitational waveforms.  Phys. Rev. Lett., 113:151101, Oct 2014.
  * [46] Tim Dietrich, Sebastian Khan, Reetika Dudi, Shasvath J. Kapadia, Prayush Kumar, Alessandro Nagar, Frank Ohme, Francesco Pannarale, Anuradha Samajdar, Sebastiano Bernuzzi, Gregorio Carullo, Walter Del Pozzo, Maria Haney, Charalampos Markakis, Michael Pürrer, Gunnar Riemenschneider, Yoshinta Eka Setyawati, Ka Wa Tsang, and Chris Van Den Broeck.  Matter imprints in waveform models for neutron star binaries: Tidal and self-spin effects.  Phys. Rev. D, 99:024029, Jan 2019.
  * [47] Tim Dietrich, Anuradha Samajdar, Sebastian Khan, Nathan K. Johnson-McDaniel, Reetika Dudi, and Wolfgang Tichy.  Improving the nrtidal model for binary neutron star systems.  Phys. Rev. D, 100:044003, Aug 2019.
  * [48] Alessandra Buonanno, Bala R. Iyer, Evan Ochsner, Yi Pan, and B. S. Sathyaprakash.  Comparison of post-newtonian templates for compact binary inspiral signals in gravitational-wave detectors.  Phys. Rev. D, 80:084043, Oct 2009.
  * [49] Justin Vines, Éanna É. Flanagan, and Tanja Hinderer.  Post-1-newtonian tidal effects in the gravitational waveform from binary inspirals.  Phys. Rev. D, 83:084051, Apr 2011.
  * [50] Alejandro Bohé et al.  Improved effective-one-body model of spinning, nonprecessing binary black holes for the era of gravitational-wave astrophysics with advanced detectors.  Phys. Rev. D, 95(4):044028, 2017.
  * [51] Sebastiano Bernuzzi, Alessandro Nagar, Marcus Thierfelder, and Bernd Brügmann.  Tidal effects in binary neutron star coalescence.  Phys. Rev. D, 86(4):044030, August 2012.
  * [52] Lisa Barsotti, Lee McCuller, Matthew Evans, and Peter Fritschel.  The a+ design curve.  Technical Report T1800042-v5, LIGO, 2018.
  * [53] Tim Dietrich, Tanja Hinderer, and Anuradha Samajdar.  Interpreting binary neutron star mergers: describing the binary neutron star dynamics, modelling gravitational waveforms, and analyzing detections.  General Relativity and Gravitation, 53(3):27, March 2021.
  * [54] LIGO Scientific Collaboration and Virgo Scientific Collaboration.  Properties of the binary neutron star merger gw170817.  Technical Report P1800061-v11, LIGO, 2018.
  * [55] B. P. Abbott, Ligo Scientific Collaboration Kagra Collaboration, and VIRGO Collaboration.  Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA.  Living Reviews in Relativity, 21(1):3, April 2018.
  * [56] Matthew Evans, Jan Harms, and Salvatore Vitale.  Exploring the sensitivity of next generation gravitational wave detectors.  Technical Report P1600143-v18, LIGO, 2018.
  * [57] David Reitze, Rana X. Adhikari, Stefan Ballmer, Barry Barish, Lisa Barsotti, GariLynn Billingsley, Duncan A. Brown, Yanbei Chen, Dennis Coyne, Robert Eisenstein, Matthew Evans, Peter Fritschel, Evan D. Hall, Albert Lazzarini, Geoffrey Lovelace, Jocelyn Read, B. S. Sathyaprakash, David Shoemaker, Joshua Smith, Calum Torrie, Salvatore Vitale, Rainer Weiss, Christopher Wipf, and Michael Zucker.  Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO.  In Bulletin of the American Astronomical Society, volume 51, page 35, September 2019.
  * [58] Matthew Evans, Rana X Adhikari, Chaitanya Afle, Stefan W. Ballmer, Sylvia Biscoveanu, Ssohrab Borhanian, Duncan A. Brown, Yanbei Chen, Robert Eisenstein, Alexandra Gruson, Anuradha Gupta, Evan D. Hall, Rachael Huxford, Brittany Kamai, Rahul Kashyap, Jeff S. Kissel, Kevin Kuns, Philippe Landry, Amber Lenon, Geoffrey Lovelace, Lee McCuller, Ken K. Y. Ng, Alexander H. Nitz, Jocelyn Read, B. S. Sathyaprakash, David H. Shoemaker, Bram J. J. Slagmolen, Joshua R. Smith, Varun Srivastava, Ling Sun, Salvatore Vitale, and Rainer Weiss.  A horizon study for cosmic explorer: Science, observatories, and community, 2021.
  * [59] Michele Maggiore, Chris Van Den Broeck, Nicola Bartolo, Enis Belgacem, Daniele Bertacca, Marie Anne Bizouard, Marica Branchesi, Sebastien Clesse, Stefano Foffa, Juan Garcí a-Bellido, Stefan Grimm, Jan Harms, Tanja Hinderer, Sabino Matarrese, Cristiano Palomba, Marco Peloso, Angelo Ricciardone, and Mairi Sakellariadou.  Science case for the einstein telescope.  Journal of Cosmology and Astroparticle Physics, 2020(03):050–050, mar 2020.
  * [60] Jocelyn S. Read, Luca Baiotti, Jolien D. E. Creighton, John L. Friedman, Bruno Giacomazzo, Koutarou Kyutoku, Charalampos Markakis, Luciano Rezzolla, Masaru Shibata, and Keisuke Taniguchi.  Matter effects on binary neutron star waveforms.  Phys. Rev. D, 88(4):044042, August 2013.
  * [61] C. J. Moore, R. H. Cole, and C. P. L. Berry.  Gravitational-wave sensitivity curves.  Classical and Quantum Gravity, 32(1):015014, January 2015.
  * [62] Jordan B. Camp and Neil J. Cornish.  Gravitational Wave Astronomy.  Annual Review of Nuclear and Particle Science, 54:525–577, December 2004.
  * [63] Tatsuya Narikawa, Nami Uchikata, Kyohei Kawaguchi, Kenta Kiuchi, Koutarou Kyutoku, Masaru Shibata, and Hideyuki Tagoshi.  Reanalysis of the binary neutron star mergers GW170817 and GW190425 using numerical-relativity calibrated waveform models.  Phys. Rev. Res., 2(4):043039, 2020.
  * [64] Rossella Gamba, Matteo Breschi, Sebastiano Bernuzzi, Michalis Agathos, and Alessandro Nagar.  Waveform systematics in the gravitational-wave inference of tidal parameters and equation of state from binary neutron-star signals.  Phys. Rev. D, 103(12):124015, June 2021.
  * [65] Ethan Payne, Colm Talbot, and Eric Thrane.  Higher order gravitational-wave modes with likelihood reweighting.  Physical Review D, 100(12), dec 2019.
  * [66] Sascha Husa, Sebastian Khan, Mark Hannam, Michael Pürrer, Frank Ohme, Xisco Jiménez Forteza, and Alejandro Bohé.  Frequency-domain gravitational waves from nonprecessing black-hole binaries. I. New numerical waveforms and anatomy of the signal.  Phys. Rev. D, 93(4):044006, 2016.
  * [67] Sebastian Khan, Sascha Husa, Mark Hannam, Frank Ohme, Michael Pürrer, Xisco Jiménez Forteza, and Alejandro Bohé.  Frequency-domain gravitational waves from nonprecessing black-hole binaries. II. A phenomenological model for the advanced detector era.  Phys. Rev. D, 93(4):044007, February 2016.
  * [68] Benjamin D. Lackey, Michael Pürrer, Andrea Taracchini, and Sylvain Marsat.  Surrogate model for an aligned-spin effective one body waveform model of binary neutron star inspirals using Gaussian process regression.  Phys. Rev. D, 100(2):024002, 2019.
  * [69] Charlie Hoy and Vivien Raymond.  PESUMMARY: The code agnostic Parameter Estimation Summary page builder.  SoftwareX, 15:100765, July 2021.
  * [70] LIGO and Virgo Scientific Collaborations.  Calibration uncertainty envelope release for gwtc-1.  Technical Report P1900040-v1, LIGO, 2019.
  * [71] LIGO and Virgo Scientific Collaborations.  Power spectral densities (psd) release for gwtc-1.  Technical Report P1900011-v1, LIGO, 2019.
  * [72] Alex Nitz, Ian Harry, Duncan Brown, Christopher M. Biwer, Josh Willis, Tito Dal Canton, Collin Capano, Thomas Dent, Larne Pekowsky, Andrew R. Williamson, Soumi De, Miriam Cabero, Bernd Machenschalk, Duncan Macleod, Prayush Kumar, Francesco Pannarale, Steven Reyes, Gareth S Cabourn Davies, dfinstad, Sumit Kumar, Márton Tápai, Leo Singer, Sebastian Khan, Stephen Fairhurst, Alex Nielsen, Shashwat Singh, Thomas Massinger, Koustav Chandra, shasvath, and veronica villa.  gwastro/pycbc: v2.0.5 release of pycbc, July 2022.
  * [73] Luc Blanchet.  Gravitational Radiation from Post-Newtonian Sources and Inspiralling Compact Binaries.  Living Reviews in Relativity, 17(1):2, December 2014.
  * [74] J. Baker, M. Campanelli, C. O. Lousto, and R. Takahashi.  Modeling gravitational radiation from coalescing binary black holes.  Phys. Rev. D, 65(12):124012, June 2002.
  * [75] Afura Taylor and Vijay Varma.  Gravitational wave peak luminosity model for precessing binary black holes.  Phys. Rev. D, 102(10):104047, November 2020.
  * [76] Wynn C. G. Ho and Dong Lai.  Resonant tidal excitations of rotating neutron stars in coalescing binaries.  MNRAS, 308(1):153–166, September 1999.
  * [77] Éanna É. Flanagan and Tanja Hinderer.  Constraining neutron-star tidal Love numbers with gravitational-wave detectors.  Phys. Rev. D, 77(2):021502, January 2008.
  * [78] Eanna E. Flanagan and Etienne Racine.  Gravitomagnetic resonant excitation of Rossby modes in coalescing neutron star binaries.  Phys. Rev. D, 75:044001, 2007.
  * [79] Nils Andersson and Wynn C. G. Ho.  Using gravitational-wave data to constrain dynamical tides in neutron star binaries.  Phys. Rev. D, 97(2):023016, January 2018.
  * [80] David Tsang, Jocelyn S. Read, Tanja Hinderer, Anthony L. Piro, and Ruxandra Bondarescu.  Resonant Shattering of Neutron Star Crusts.  Phys. Rev. Lett., 108(1):011102, January 2012.
  * [81] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  Nature, 585(7825):357–362, September 2020.
  * [82] The pandas development team.  pandas-dev/pandas: Pandas, February 2020.
  * [83] Gregory Ashton, Moritz Hübner, Paul D. Lasky, Colm Talbot, Kendall Ackley, Sylvia Biscoveanu, Qi Chu, Atul Divakarla, Paul J. Easter, Boris Goncharov, Francisco Hernandez Vivanco, Jan Harms, Marcus E. Lower, Grant D. Meadors, Denyz Melchor, Ethan Payne, Matthew D. Pitkin, Jade Powell, Nikhil Sarin, Rory J. E. Smith, and Eric Thrane.  BILBY: A User-friendly Bayesian Inference Library for Gravitational-wave Astronomy.  ApJS, 241(2):27, April 2019.