#  gwbench: a novel Fisher information package for gravitational-wave
benchmarking

S. Borhanian  Institute for Gravitation and the Cosmos, Department of Physics,
Pennsylvania State University, University Park, PA 16802, USA  Theoretisch-
Physikalisches Institut, Friedrich-Schiller-Universität Jena, 07743, Jena,
Germany<EMAIL_ADDRESS>
###### Abstract

We present a new Python package, gwbench, implementing the well-established
Fisher information formalism as a fast and straightforward tool for the
purpose of gravitational-wave benchmarking, i.e. the estimation of signal-to-
noise ratios and measurement errors of gravitational waves observed by a
network of detectors. Such an infrastructure is necessary due to the high
computational cost of Bayesian parameter estimation methods which renders them
less effective for the scientific assessment of gravitational waveforms,
detectors, and networks of detectors, especially when determining their
effects on large populations of gravitational-wave sources spread throughout
the universe. gwbench further gives quick access to detector locations and
sensitivities, while including the effects of Earth’s rotation on the latter,
as well as waveform models and their derivatives, while giving access to the
host of waveforms available in the LSC Algorithm Library. With the provided
functionality, gwbench is relevant for a wide variety of applications in
gravitational-wave astronomy such as waveform modeling, detector development,
cosmology, and tests of general relativity.

##  I Introduction and motivation

The initial detections of gravitational waves (GWs) emitted during the
coalescences of black holes and neutron stars GW150914 ; GW170817  have
confirmed a fundamental prediction of general relativity and opened a new
observational window into the universe in the same breath. Further detections
by the LIGO and Virgo observatories have pushed the number of observed GW
signals into the decades GWTC1 , and hence set off the era of GW astronomy.
This success has ignited the desire for detections that allow for better
estimation of the information buried in the signals that the GW observatories
record, thus allowing us to reveal new phenomena and novel physics as well as
shed light into greater depths of the cosmos.

Such desires pose challenges on many essential components in the figurative GW
detection machinery, with the development of new waveform (WF) models and the
planning of future detector designs being two important tasks to tackle. The
former is imperative to lower the systematic biases in what kind of signals
the observatories will be able to detect, to increase the information that can
be extracted from the detector output, and possibly to lower the computational
cost of WF model evaluation. The inclusion of higher-order spherical harmonic
modes beyond the quadrupole London:2017bcn ; Cotesta:2018fcv  is one such
example illustrating the impact of improved WF models as it enhances the
detectability of binaries whose orbits are highly inclined with respect to the
line of sight and allows for an improved parameter estimation by lifting the
degeneracy between luminosity distance and inclination angle.

The detector challenge is even more fundamental: new, ground-breaking analyses
require ‘richer’ data to forge through. The current, so called second-
generation (2G) of GW detectors, LIGO Hanford and Livingston aLIGO  as well as
Virgo Virgo , will be joined by two more detectors, KAGRA KAGRA  and LIGO-
India ligo_india , and further sensitivity improvements are planned to go in
effect in the years to follow aPlusWP ; Aasi:2013wya . Nevertheless the limits
of these facilities to address future, scientific goals are already apparent
LISA_science_case ; ET_science_case . Thus the GW community has been pursuing
the conceptualization of the next generation of detectors. The proposals range
from the Voyager effort, envisioning a multitude of upgrades to the LIGO
facilities Voyager , to completely new, third-generation (3G) observatories,
namely the Einstein Telescope ET_1 ; ET_2 ; ET_science_case  and the Cosmic
Explorer effort CE_proposal . In addition, efforts are under way to search for
GW from space, with the Laser Interferometer Space Antenna (LISA)
LISA_science_case .

Independently of what WF models shall be developed and adopted, or which
facilities continue or start operating, the GW community will need tools to
_benchmark_ the scientific output of any combination of these two components
in a quick and efficient way. One such tool is the Fisher information
formalism (FIF) Cutler:1994ys ; Poisson:1995ef  which has proven to be a
viable method for GW benchmarking. Its major caveat is its dependence on high
signal-to-noise ratios to give reliable results, thus making it mostly
suitable to exceptional events in current detectors. Nevertheless, the
formalism is very useful and widely applicable in studies entailing future
generations of detectors. For an extended review on the caveats of the FIF we
refer to Vallisneri:2007ev . Further, developments to improve the speed of
full Bayesian parameter estimation, for example via the application of machine
learning techniques Chua:2019wwt ; Gabbard:2019rde ; Green:2020dnx ;
Green:2020hst , could make the Bayesian framework a viable benchmarking tool
in the future.

In this work, we present a new Python package, gwbench, that implements the
FIF in an easy-to-use manner and further provides means to compute and access
a host of quantities that are necessary for benchmarking. For a selection of
GW detectors, a so-called _network_ , and given a WF model, gwbench can
compute: plus and cross polarizations of the WF, detector power spectral
densities, antenna patterns, location phase factors, detector responses,
detector and network signal to noise ratios, and measurement errors in the WF
model parameters from the FIF. Two particularly substantial features are the
inclusion of WFs from the LSC Algorithm Library (LAL) lalsuite ; wette2020
and the capability to include the effects of Earth’s rotation in the detector
antenna patterns.

The paper is structured as follows: Section II introduces the Cosmic Explorer
trade study which stimulated the development of this package. In Section III
we briefly review the FIF and aforementioned core concepts utilized by
gwbench. Next, we show a straightforward example application of gwbench to
showcase its functionality in Section IV and proceed to present a few
validation tests of the core methods of gwbench in Section V. In Section VI,
we describe the application of the code in two previous papers and conclude
this work in Section VII.

##  II Cosmic Explorer Trade Study

The Cosmic Explorer Project111https://cosmicexplorer.org/ aims to develop
among other things the science case for a next-generation GW detector, the
Cosmic Explorer (CE). CE is a promising US proposal for a ground-based
detector design which is targeting to answer a wide variety of scientific
questions Sathyaprakash:2019yqt ; Kalogera:2019sui ; Sathyaprakash:2019rom ;
Sathyaprakash:2019nnu ; Barack:2018yly  as part of a global network of GW
observatories. These questions are grouped under three overarching _science
goals_ : (i) dynamics of dense matter and extreme environments, (ii) black
holes and neutron stars through cosmic time, and (iii) extreme gravity and
fundamental physics. The questions encompass a broad range of fields relevant
to GW astronomy including but not restricted to neutron star physics,
cosmology, tests of general relativity, and the feasibility of multi-messenger
observations.

gwbench was developed to perform an exhaustive _trade study_ to gauge CE’s
performance in reaching the set science goals. The trade study examines some
457 global GW detector networks combining 50 CE configurations with plausible
background detectors such as upgraded 2G facilities or the Einstein Telescope.
Besides the number of networks, the computational requirements are also driven
by the list of science questions to be answered.

Hence, we decided to implement the FIF in gwbench, as it allows us to perform
the required tasks effectively for many events and different waveform models
in each network of interest. The measurement error estimates computed from the
FIF are processed to obtain pre-defined performance metrics which in return
allow us to evaluate the science potential of the various network
configurations. Such metrics are for example the signal-to-noise ratios,
90%-credible sky areas, expected event rates, and measurement error estimates
on WF parameters, but also on inferred quantities such as neutron star radii.

As such, gwbench has an important role to play for the development of CE’s
science case. The broad range of science goals that can be examined with the
code showcases that gwbench is a powerful tool that has a wide variety of
applications in gravitational-wave astronomy. We present two of our recent
studies that were enabled by gwbench in Section VI.

##  III Concepts for gravitational wave benchmarking

gwbench implements a variety of standard concepts and aims to make them easily
accessible to the user. We briefly review them LRR_Sathyaprakash_2009  in this
section.

#### Waveform polarizations:

Gravitational waves are emitted by a host of astrophysical and cosmological
phenomena, with the coalescence of compact binaries as the most prominent
example. The incoming GWs can be decomposed into two polarizations, _plus_
$h_{+}(t;\bm{\Theta}_{GW})$ and _cross_ $h_{\times}(t;\bm{\Theta}_{GW})$, and
carry an abundance of information about their source systems, encoded in the
parameter vector $\bm{\Theta}_{GW}$. This vector depends on the WF model. In
case of the coalescence of two black holes the standard set of such
parameters,
$\bm{\Theta}_{GW}=\\{m_{1},m_{2},\bm{\chi_{1}},\bm{\chi_{2}},D_{L},\iota,t_{c},\phi_{c}\\}$
contains the companion masses $m_{1},m_{2}$ and dimensionless spin vectors
$\bm{\chi_{1}},\bm{\chi_{2}}$, the binary’s luminosity distance $D_{L}$, the
inclination angle $\iota$ of its orbital plane with respect to the line of
sight, as well as two integration parameters: the time $t_{c}$ and phase
$\phi_{c}$ of coalescence.

In the remainder of this manuscript, we will use the Fourier transforms of the
WF and other quantities in the frequency-domain:

$\displaystyle
H_{+/\times}(f)=\int_{-\infty}^{\infty}\mathrm{d}t\,h_{+/\times}(t)\,\mathrm{e}^{\mathrm{i}\,2\pi
ft}.$ (1)

#### Antenna patterns and location phase factor:

The quadrupolar nature of GWs results in sky-position dependent responses from
detectors to incoming waves. The responses are governed by the so-called
_antenna pattern functions_ $F_{+}$ and $F_{\times}$ which are functions of
the source’s sky position—right ascension $\alpha$ and declination
$\delta$—and polarization angle $\psi$. Besides, the relative motion between
the source and detector results in a frequency-dependence in the Fourier-
domain Wen:2010cr ; Zhao:2017cbb ; albeit being negligible for signals that
are short compared to the time scale of the motion.

Another group of important quantities that encode the spatial distribution of
detectors in networks are the _location phase factors_ $F_{lp}$ which shift
the phase of the detector responses appropriately for each detector:

$\displaystyle F_{lp}=\exp\left[\mathrm{i}2\pi
f/c\,\hat{\bm{r}}(\alpha,\delta)\cdot\bm{d}(\alpha_{d},\beta_{d})\right],$ (2)

where $\hat{\bm{r}}$ and $\bm{d}$ are the unit vector pointing to the source
and the position vector of the detector, respectively, if measured from the
center of Earth. $\alpha_{d}$ and $\beta_{d}$ are the detector’s longitude and
latitude, respectively.

#### Detector response:

The aforementioned _detector response_ $H$ is the GW strain in the frequency-
domain measured by a detector, combining the WF polarizations with the
detector’s antenna patterns and location phase factor. It differs from the
actually recorded detector output $s(t)=h(t)+n(t)$ in the time-domain or
$S(f)=H(f)+N(f)$ in the frequency-domain, which also contains contributions
from the detector noise $n(t)$ or $N(f)$, respectively. The detector response
is given by

$\displaystyle H(f;\bm{\Theta})$ $\displaystyle=$ $\displaystyle
F_{lp}(f;\alpha,\delta)\cdot$ (3)
$\displaystyle[H_{+}(f;\bm{\Theta}_{GW})\,F_{+}(f;\alpha,\delta,\psi)+$
$\displaystyle\phantom{[}H_{\times}(f;\bm{\Theta}_{GW})\,F_{\times}(f;\alpha,\delta,\psi)],$

where $\bm{\Theta}=\\{\bm{\Theta}_{GW},\alpha,\delta,\psi\\}$.

#### Power spectral density:

The sensitivity of a detector determines its ability to measure the strain
$H$. It is governed by systematics due to detector noise $n(t)$ and
characterized by the noise auto-correlation
$\kappa=\overline{n(t_{1})\,n(t_{2})}$222The average is taken with respect to
an ensemble of noise realizations. The times $t_{1}$ and $t_{2}$ are any two
arbitrary times at which the detector noise is evaluated or measured.. In the
following, we assume the noise to be stationary and Gaussian, with mean
$\overline{n}=0$. The former assumption ensures that $\kappa$ only depends on
the time-difference $t^{*}=t_{1}-t_{2}$. Thus, the _one-sided power spectral
density_ (PSD) $S_{n}$—the Fourier transform of $\kappa$ for positive
frequencies and thus its frequency-domain counterpart—is given by

$\displaystyle
S_{n}(f)=\frac{1}{2}\int_{-\infty}^{\infty}\mathrm{d}t^{*}\kappa(t^{*})\,\mathrm{e}^{\mathrm{i}2\pi
ft^{*}},\text{ with }f>0.$ (4)

$S_{n}(f)$ indicates the sensitivity of the detector at different GW
frequencies.

#### Noise-weighted scalar product and signal-to-noise ratio:

The existence of noise in the detector output affects the detection and
measurement of the GW signal strain, obscuring the difference between various
detector responses to some level. One way to assess the _overlap_ of these
responses in a manner that takes the detector’s sensitivity into account is by
defining a _noise-weighted scalar product_ LRR_Sathyaprakash_2009 ,

$\displaystyle\langle
H,G\rangle=2\int_{0}^{\infty}\mathrm{d}f\,\frac{H(f)G^{*}(f)+H^{*}(f)G(f)}{S_{n}(f)},$
(5)

where $H,G$ are two detector responses, e.g. for two different WF models or
the same model evaluated at different parameters.

Given this scalar product, we can calculate the _signal-to-noise ratio_ (SNR)
$\rho$ of a given detector response $H$ for a given detector PSD $S_{n}$ as

$\displaystyle\rho^{2}=\langle
H,H\rangle=4\int_{0}^{\infty}\mathrm{d}f\,\frac{|H(f)|^{2}}{S_{n}(f)},$ (6)

indicating the relative loudness of the detector response compared to the
noise in the detector output LRR_Sathyaprakash_2009 .

#### Fisher information formalism:

The _FIF_ or _Fisher analysis_ is well-established in GW data analysis and
discussed in detail in Cutler:1994ys ; Poisson:1995ef . Here, we give a brief
summary of how it allows us to estimate parameter error bounds.

Since we assumed that the noise $n=s-h$ behaves Gaussian with zero mean, the
same holds for its Fourier transform $N=S-H$. Thus, the probability of the
noise can be written as

$\displaystyle
p(\bm{\Theta})=p^{0}(\bm{\Theta})\,\mathrm{e}^{-\frac{1}{2}\langle
S-H(\bm{\Theta}),S-H(\bm{\Theta})\rangle},$ (7)

where $\bm{\Theta}$ and $p^{0}$ are the parameter vector of the detector
response $H$ and the prior on these parameters, respectively. Then the value
of $\bm{\Theta}$ at peak probability is a good estimate of the true value
$\bm{\Theta^{*}}$ for the given GW signal $S,$ assuming a large SNR.

The probability is the greatest when the exponential is the
largest333Poisson:1995ef  shows that this holds even for non-flat priors
$p^{0}$.. Around the maximal value, the exponent $E=\langle S-H,S-H\rangle$
can be expanded as

$\displaystyle
E(\bm{\Theta})=E(\bm{\Theta^{*}})+\frac{1}{2}\left.\frac{\partial^{2}E(\bm{\Theta})}{\partial\Theta_{i}\partial\Theta_{j}}\right|_{\bm{\Theta}=\bm{\Theta^{*}}}\Delta\Theta_{i}\,\Delta\Theta_{j}+\dots,$
(8)

where $\Delta\Theta_{i}=\Theta_{i}-\Theta^{*}_{i}$. The Hessian

$\frac{\partial^{2}E(\bm{\Theta})}{\partial\Theta_{i}\partial\Theta_{j}}=2\,\langle\partial_{\Theta_{i}}H(\bm{\Theta}),\partial_{\Theta_{j}}H(\bm{\Theta})\rangle+\langle\partial_{\Theta_{i}}\partial_{\Theta_{j}}H(\bm{\Theta}),N\rangle$
(9)

reduces for noise with zero-mean to the first summand in which we identify the
_Fisher information matrix_ (FIM) $\Gamma$

$\Gamma_{ij}\equiv\langle\partial_{\Theta_{i}}H(\bm{\Theta}),\partial_{\Theta_{j}}H(\bm{\Theta})\rangle.$
(10)

Equations (8) and (10) inserted in the probability (7) yield

$p(\bm{\Theta})\sim\mathrm{e}^{-\frac{1}{2}\Gamma_{ij}\,\Delta\Theta_{i}\,\Delta\Theta_{j}}.$
(11)

Thus the assumption of Gaussianity motivates the association of the FIM to be
the inverse of the _covariance matrix_ $\Sigma\equiv\Gamma^{-1}$. The diagonal
and off-diagonal elements of $\Sigma$ denote the variances and covariances of
the parameters, respectively, due to the uncertainty introduced by the
detector noise and give $1\sigma$-error estimates via
$\sigma_{\Theta_{i}}=\sqrt{\Sigma_{ii}}$.

#### Detector networks and associated quantities:

The above definitions of the SNR $\rho$ and FIM $\Gamma$ depend on quantities
such as the detector response and PSD and are therefore inherently detector-
specific. In a network, detectors work in tandem to improve detectability and
parameter estimation which can be captured by defining _network versions_ of
the SNR and FIM

$\displaystyle\rho_{\text{net}}^{2}$
$\displaystyle=\sum_{n=1}^{N_{d}}\,\rho_{n}^{2},$ (12)
$\displaystyle\Gamma_{\text{net}}$
$\displaystyle=\sum_{n=1}^{N_{d}}\,\Gamma_{n},$ (13)

where $N_{d}$ is the number of detectors in the network and $\rho_{n}$ and
$\Gamma_{n}$ are the SNR and FIM of the n-th detector, respectively.
Consequently, the network covariance matrix and error estimates follow from
$\Gamma_{\text{net}}$ as defined before.

#### Condition number of the Fisher information matrix:

If the covariance matrix is obtained via numerical inversion, the FIM needs to
be well-conditioned which can be checked with the condition number
$c_{\Gamma}=e_{max}/e_{min}$. Here, $e_{max},e_{min}$ are the largest and
smallest eigenvalues of $\Gamma$. If the condition number exceeds a critical
value—a conservative choice for double-precision is $c_{crit}=10^{15}$
Rodriguez:2011aa —the inversion result cannot be trusted generically and
should thus be disregarded.

##  IV Gravitational wave benchmarking with gwbench

###  IV.1 Package overview

The installation instructions for gwbench and the source code are publicly
available on GitLab: https://gitlab.com/sborhanian/gwbench. The package
consists of 33 modules which implement the concepts discussed in Section III
as well as several WFs models. In this section, we will demonstrate with the
aid of a few examples how to use the main module of gwbench, network, to
perform GW benchmarking. These examples are available as Python scripts at
https://gitlab.com/sborhanian/gwbench/-/tree/master/example_scripts.

network contains the Network class, see Appendix A, which is designed to
handle the concepts introduced in Section III—detector initialization (antenna
patterns, location phase factors, and PSDs), WF evaluation, and both SNR and
error benchmarking—by conducting and facilitating the tasks of the other
modules in the package. We will refer to the appendices as needed to discuss
these modules and aspects of the package.

###  IV.2 Basic application using gwbench

We begin with a step-by-step showcase of the basic usage of the Network module
following the example script quick_start.py. In this example we estimate the
SNR, measurement errors, and 90%-credible sky area for a signal with the same
masses as the first detected GW event GW150914 GW150914 , if it was measured
by three detectors at the LIGO sites in Hanford, WA and Livingston, LA and the
Virgo site in Cascina, Italy. Each detector is set to the aLIGO detector
sensitivities aLIGO  and we perform the Fisher analysis with the WF model
TaylorF2 1991PhRvD..44.3819S ; Poisson:1997ha ; Mikoczi:2005dn ; Arun:2008kb ;
Buonanno:2009zt ; Bohe:2013cla ; Bohe:2015ana ; Mishra:2016whh , see Appendix
B.

We begin by importing the numpy and gwbench.network modules and specify the
network_spec which tells the Network which Detector instances to initialize:
[’aLIGO_H’,’aLIGO_L’,’aLIGO_V’] will result in three detectors with aLIGO
technology, located and oriented like the Hanford, Livingston, and Virgo
detectors (refer to the Appendices C and D for other sites and sensitivities
available in gwbench).

⬇

import numpy as np

from gwbench import network

network_spec = [’aLIGO_H’,’aLIGO_L’,’aLIGO_V’]

net = network.Network(network_spec)

While the Network class is the main interaction point for the user, single-
detector tasks are relayed to the Detector class. Appendix A compares these
two classes to show their similarities and how they synergize to conduct the
GW benchmarking together.

Next, we choose the WF model of interest, TaylorF2, and initialize the
Waveform object inside our Network. Waveform is another important class that
handles WF selection and evaluation. Appendix B contains a brief description
of this class and the full list of available WF models in gwbench, at the time
of writing this paper.

⬇

wf_model_name = ’tf2’

net.set_wf_vars(wf_model_name=wf_model_name)

Finally, we need to set a few Network variables like the frequency array and
injection parameters in the detector frame, used for WF evaluation. Further,
the code needs to know with respect to which of these parameters to take
derivatives for the FIM. Besides, we can set the Network to convert the
derivatives $\partial_{p}$ and errors $\sigma_{p}$ for certain parameters $p$
into the cosine ($\partial_{\cos(p)}$ and $\sigma_{\cos(p)}$) or logarithmic
($\partial_{\log(p)}$ and $\sigma_{\log(p)}$) versions. The final option tells
the code whether to take Earth’s rotation into account for the computation of
the antenna pattern functions, see Appendix C.

⬇

f = np.arange(5.,61.5,2**-4)

inj_params = {

’Mc’: 30.9,

’eta’: 0.247,

’chi1z’: 0,

’chi2z’: 0,

’DL’: 475,

’tc’: 0,

’phic’: 0,

’iota’: np.pi/4,

’ra’: np.pi/4,

’dec’: np.pi/4,

’psi’: np.pi/4,

’gmst0’: 0

}

deriv_symbs_string = ’Mc eta DL tc phic iota ra dec psi’

conv_cos = (’iota’,’dec’)

conv_log = (’Mc’,’DL’)

use_rot = 0

net.set_net_vars(

f=f, inj_params=inj_params,

deriv_symbs_string=deriv_symbs_string,

conv_cos=conv_cos, conv_log=conv_log,

use_rot=use_rot

)

The high frequency cutoff of $61.5\,\text{Hz}$ corresponds to the GW frequency
at the innermost stable circular orbit of a compact binary coalescence with
total detector-frame mass $M=\mathcal{M}/\eta^{3/5}=71.5\,M_{\odot}$ for the
example injections with detector-frame chirp mass
$\mathcal{M}=30.9\,M_{\odot}$ and symmetric mass ratio $\eta=0.247$. The
detector-frame masses correspond to redshifted sources-frame masses of
GW150914, i.e. $m_{\rm det}=m_{\rm source}(1+z)$ where we used a fiducial
value of $z=0.1$ and computed the luminosity distance via Planck18 cosmology
from the astropy package astropy:2013 ; astropy:2018 . gmst0 is the Greenwich
Mean Sidereal Time when the signal passes Earth’s center.

The Network is prepared for benchmarking and we calculate the WF polarizations
and their derivatives with respect to the chosen parameters with the following
function calls.

⬇

net.calc_wf_polarizations()

net.calc_wf_polarizations_derivs_num()

The first line only computes the plus and cross polarizations in the
frequency-domain, net.hfp and net.hfc, while the second obtains a dictionary
net.del_hfpc containing the derivatives of the two WF polarizations, in
addition to that. The structure of these _derivative dictionaries_ is detailed
in Appendix E.

Next, we setup antenna patterns, location phase factors, and PSDs for each
Detector in the Network .

⬇

net.setup_ant_pat_lpf_psds()

The calculated quantities can be accessed in each Detector instance via
net.detectors[i].Fp, net.detectors[i].Fc, net.detectors[i].Flp, and
net.detectors[i].psd for $\texttt{i}=0,...,N_{d}-1$. We note that the Detector
automatically truncates the frequency array, if necessary, to the range
dictated by the given detector’s PSD. This truncation is saved in
net.detectors[i].f for $\texttt{i}=0,...,N_{d}-1$ and translates to all
consecutive computations which depend on the frequency, e.g. in the detector
response or the location phase factor, and is necessary to avoid array-length
mismatches within a Detector instance. On the other hand, it can result in a
mismatch when comparing frequency-dependent quantities of two Detector
instances or between the Network and one of its Detectors.

With antenna patterns and location phase factors prepared, we can calculate
the detector responses and their derivatives, analogously to the case for the
WF polarizations.

⬇

net.calc_det_responses()

net.calc_det_responses_derivs_num()

The evaluated responses and derivative dictionaries can be accessed in each
Detector instance via net.detectors[i].hf and net.detectors[i].del_hf for
$\texttt{i}=0,...,N_{d}-1$.

With the previously computed PSDs, we obtain readily the detector and network
SNRs.

⬇

net.calc_snrs()

The function calculates both $\rho$ and $\rho^{2}$ for the total network and
each detector in the network. The values are saved in the respective Network
and Detector instances.

Finally, we compute the error estimates on the parameters specified in
deriv_symbs_string and the estimates of the 90%-credible sky area.

⬇

net.calc_errors()

net.calc_sky_area_90()

net.print_detectors()

net.print_network()

calc_errors contains several steps: It starts by calculating the FIMs for all
the Detector instances and then the Network one from those. Next, it computes
and saves the condition number $c_{\Gamma}$ for each FIM $\Gamma$, and checks
if the matrices are well-conditioned ($c_{\Gamma}<10^{15}$). If this criterion
is fulfilled, the code will proceed to invert each FIM to obtain all the
covariance matrices $\Sigma$ and consequently the error estimates. Finally, it
calculates the inversion accuracy $\epsilon=||\Gamma\cdot\Sigma-I||_{\text{\rm
max}}$, where $I$ and $||\cdot||_{\text{\rm max}}$ are the identity and
maximum matrix norm, respectively.

Analogously to the SNR computations, the Network and Detector instances save
their respective quantities. The final two function calls simply print the
contents of all the Detector instances inside the Network as well its own
contents.

#### Further remarks:

For the sake of simplicity, we omitted non-default function arguments in the
example above. The Network class is designed to perform ‘default’-benchmarking
after setting all the necessary variables via net.set_wf_vars and
net.set_net_vars without the need to pass further variables during the
consecutive function calls. Nevertheless, certain functions give the user the
option to tune the functionality of the code. In the following, we will
present these options for the main GW benchmarking functions showcased in the
example.

####  net.setup_ant_pat_lpf_psds:

This function takes the following three arguments (default values) F_lo
(-np.inf), F_hi (np.inf), and psd_file_dict (None). The first two are user-
settable frequency cutoffs which truncate the frequency array, if they lie
within the PSD’s intrinsic frequency range. The third argument allows the user
to specify PSD files that are not part of gwbench. It has to be passed as

⬇

psd_file_dict = {

det_key1:{’psd_file’:file1,’is_asd’:bool1}, ...

}

where the det_key are the exact keys that each Detector contains, the file are
the relative paths to the PSD files, and the booleans bool tell the code
whether the file contains PSDs $S_{n}$ (bool=0) or amplitude spectral
densities $\sqrt{S_{n}}$ (bool=1). The PSD files should be formatted like the
txt-files in ’gwbench/noise_curves’. A caveat of using external PSD files is
that the Network and Detector labels might not reflect this choice and will
need manual changes. Fortunately, these can be done quite easily, if needed,
in run time.

gwbench can naturally be used to just compute the antenna patterns, location
phase factors, and PSDs of detectors as shown in the example script
compute_ant_pat_lpf_psd.py.

####  net.calc_det_responses_derivs_num, net.calc_wf_polarizations_derivs_num:

These functions calculate the derivatives of the detector responses and WF
polarizations numerically via the numdifftools package numdifftools . The
numeric derivatives require four special arguments: the step size step, the
differentiation scheme method, the differentiation order order, and the
derivative order n. The respective default values in gwbench are 1e-7,
’central’, 2, and 1. We refer to the documentation of numdifftools for further
information. gwbench further allows the use of symbolic derivatives for
certain waveforms via the sympy package. Their evaluation is faster, but
requires a few extra steps compared to the numeric ones. Their usage is
showcased in Section IV.3, while Section V.1 compares the output of the
numeric and symbolic differentiation methods.

####  net.calc_snrs, net.calc_errors, and net.calc_sky_area_90:

These three functions have a common argument only_net which tells the code
whether to calculate the respective quantities just for the Network or also
for each Detector within. The default is 0, computing both network and
detector quantities.

Besides, net.calc_errors takes two more arguments (default values) cond_sup
(1e15) and by_element (0). The first determines the condition number supremum
$c^{*}$ that the code uses to assess whether a computed FIM $\Gamma$ is well-
conditioned ($c_{\Gamma}<c^{*}$). Passing None sets it to infinity,
effectively making every FIM well-conditioned. The second argument tells the
code whether to obtain a single inversion error $\epsilon$ for the entire FIM
via the maximum norm (default) or $n$ errors $\epsilon_{i}$. $n$ is the
dimension of the FIM while each $\epsilon_{i}$ is the maximum of the elements
of the $i$-th row and column combined. The single inversion error $\epsilon$
is equal to the maximum of the $n$ errors.

###  IV.3 Symbolic derivatives

#### Compute dictionaries of lambdified sympy derivatives:

While numeric derivatives are computed and evaluated in the same,
computationally expensive step, symbolic derivatives via sympy are first
computed generically and then evaluated for a given set of arguments. The
latter is faster than the numeric approach presented in Section IV.2, but
requires the generation of the so-called _lambdified_ derivative functions
ahead of time. In the following, we will showcase an example script that can
handle this step: generate_lambdified_functions.py.

The user has to specify the WF, the derivative variables, the detector
locations of interest, and whether to take Earth’s rotation into account. By
default, the script is setup to (i) use the WF model TaylorF2, (ii) take
derivatives with respect to chirpmass $\mathcal{M}$, symmetric mass ratio
$\eta$, luminosity distance $D_{L}$, time and phase of coalescence
$t_{c},\,\phi_{c}$, inclination angle $\iota$, right ascension $\alpha$,
declination $\delta$, and polarization angle $\psi$, (iii) use all locations
available in gwbench, and (iv) include Earth’s rotation.

Next, the code initializes a Waveform object and checks whether the derivative
variables actually form a subset of the WF and antenna pattern arguments. The
actual computation is then handled by

⬇

dr.generate_det_responses_sym(

wf,deriv_symbs_string,locs=locs,use_rot=use_rot)

which takes the Waveform, the derivative variables, the locations of interest,
and the rotation setting.

The function generates two sets of dictionaries containing different _generic
derivative expressions_ stored as lambdified sympy functions. The first type
contains the derivatives of the WF polarizations while the second stores the
derivatives of the detector responses for the detector locations specified in
the script. These dictionaries are then saved in a folder named
lambdified_functions using the dill package. The output files contain one
dictionary each and the file names are structured as

⬇

’par_deriv_WFM_{wf_model_name}_

VAR_{variables}_DET_{location}.dat’

where {wf_model_name} and {variables} are set according to the values passed
to generate_det_responses_sym. {location} is replaced by pl_cr in case of the
WF polarizations.

#### Symbolic derivative evaluation for Fisher analysis:

sym_gw_benchmarking.py showcases how to load and evaluate the symbolic
derivatives generated with generate_lambdified_functions.py and how to use the
injections module to generate random samples of input parameters for the
detector response evaluation (see Appendix F). The script is mimicked by
num_gw_benchmarking.py which differs by computing the derivatives numerically;
it can be seen as a more generic expansion of the example shown in Section
IV.2.

Most of the script simply prepares the benchmarking and hence we focus on the
usage of symbolic derivatives, in order to avoid repetition. As before, the
Network is the central hub and requires the same setup as shown in Section
IV.2 to prepare the Network and Detector instances. In contrast to the earlier
numeric example, the function call

⬇

net.calc_det_responses_derivs_num()

is replaced by

⬇

net.load_det_responses_derivs_sym()

net.calc_det_responses_derivs_sym()

which first loads the generated lambdified derivative functions and then
evaluates them. This is all that is needed besides ensuring that the variables
wf_model_name, wf_other_var_dic, deriv_symbs_string, and use_rot are set to
the same values as during the generation of the lambdified derivatives.

##  V Validating gwbench

Figure 1: The plots show the relative residuals $r$, see Equation (16),
between the numeric and symbolic partial derivatives of the detector responses
for 11 (12) input parameters and WF model tf2 (tf2_tidal). The $\sim 20$ very
poor residuals in the chirpmass derivatives correspond to ill-conditioned
Fisher matrices and are therefore disregarded for further analysis.

In the following, we present benchmarking results obtained with gwbench. We
compare numeric and symbolic derivatives and show their effects on the error
estimates. We further present error estimates from various WF models for
particularly interesting real events, namely GW170817 GW170817 , GW190521
GW190521 , and GW190814 GW190814 .

###  V.1 Comparison of numeric and symbolic differentiation methods

gwbench contains both numeric and symbolic differentiation methods. The former
can be used with all WF models, while the latter only works with some of the
implemented WFs, as specified in Appendix B. In order to compare the two
approaches for the purposes of GW benchmarking, we will compute the partial
derivatives of the detector responses $H$ and the error estimates for two WF
models, TaylorF2 (tf2) and TaylorF2 with tidal effects (tf2_tidal)
Lackey:2014fwa .

Both WFs share the same set of input parameters,
$(\mathcal{M},\eta,\chi_{1,z},\chi_{2,z},D_{L},t_{c},\phi_{c},\iota,\alpha,\delta,\psi)$,
with the tidal WF further depending two tidal deformability parameter
combinations $(\tilde{\Lambda},\delta\tilde{\Lambda})$ Wade:2014vqa

$\displaystyle\tilde{\Lambda}$
$\displaystyle=\frac{8}{13}\left[(1+7\eta-31\eta^{2})\,(\Lambda_{1}+\Lambda_{2})\,+\right.$
(14)
$\displaystyle\phantom{=\frac{8}{13}}\left.\sqrt{1-4\eta}\,(1+9\eta-11\eta^{2})\,(\Lambda_{1}-\Lambda_{2})\right],$
$\displaystyle\delta\tilde{\Lambda}$ $\displaystyle=\frac{1}{2\cdot
1319}\cdot$ (15)
$\displaystyle\left[\sqrt{1-4\eta}\,(1319-13272\eta+8944\eta^{2})\,(\Lambda_{1}+\Lambda_{2})\,+\right.$
$\displaystyle\left.(1319-15910\eta+32850\eta^{2}+3380\eta^{3})\,(\Lambda_{1}-\Lambda_{2})\right],$

where the component tidal deformabilities $\Lambda_{1}$ and $\Lambda_{2}$ are
defined as $\Lambda_{i}=\frac{2}{3}k_{2}C_{i}^{5}$ with Love number $k_{2}$,
compactness $C_{i}=\frac{m_{i}}{R_{i}}$, mass $m_{i}$, and radius $R_{i}$.

The partial derivatives are taken with respect to all available parameters
except $\delta\tilde{\Lambda}$ for the tidal WF. The detector responses are
evaluated for three detectors, each set to a 40km compact-binary optimized CE1
PSD, see Appendix D.

We performed this comparison for a total of 1000 random sets of input
parameters for each WF. In both cases we chose the injections to be uniformly
distributed in co-moving volume up to a redshift of $z=0.5$ as well as over
sky positions and orientation angles. The masses and spins are also sampled
uniformly for both WFs albeit the ranges are chosen to resemble binary black
hole or binary neutron star systems, respectively: $5\,M_{\odot}\leq
m_{1},m_{2}\leq 100\,M_{\odot}$ and $-0.75\leq\chi_{1,z},\chi_{2,z}\leq 0.75$
for tf2 and $1\,M_{\odot}\leq m_{1},m_{2}\leq 2\,M_{\odot}$ and
$-0.05\leq\chi_{1,z},\chi_{2,z}\leq 0.05$ for tf2_tidal.

Figure 2: The plots show the relative residuals $\tilde{r}$, see Equation
(17), between the $1\sigma$-error bounds, estimated via the FIF from numeric
or symbolic derivatives, for 11 (12) input parameters and WF model tf2
(tf2_tidal).

Figure 1 shows the relative residuals of all partial derivatives of the
detector responses for both waveforms and all three detectors. The residuals
are defined via

$r_{x}=\frac{\sqrt{\sum_{i}(\partial_{x}^{\text{num}}H(f_{i})-\partial_{x}^{\text{sym}}H(f_{i}))^{2}}}{\sqrt{\sum_{i}(\partial_{x}^{\text{sym}}H(f_{i}))^{2}}}$
(16)

where $i$ runs over the input frequency array and
$x\in\\{\ln(\mathcal{M}),\eta,\chi_{1,z},\chi_{2,z},\ln(D_{L}),t_{c},\phi_{c},\cos(\iota),\alpha,\cos(\delta),\psi,\ln(\tilde{\Lambda})\\}$.

The inclusion of tidal effects results in two trends for the relative
residuals of the partial derivatives in Figure 1: The errors of the
derivatives with respect to the so-called intrinsic parameters,
$(\mathcal{M},\eta,\chi_{1,z},\chi_{2,z},t_{c},\phi_{c},\tilde{\Lambda},\delta\tilde{\Lambda})$,
worsen but are also more tightly constrained, while the derivatives with
respect to the extrinsic parameters $(D_{L},\iota,\alpha,\delta,\psi)$ appear
to not be affected. The only exceptions are the chirpmass and symmetric mass
ratio derivatives, with the former showing improved residuals with tidal
effects and the latter showing more scatter. The reason for the discrepancy in
the behavior of the intrinsic and extrinsic derivatives is found in the
correlation of the intrinsic parameters to each other, thus worsening the
numeric derivatives when adding another intrinsic parameter, whereas the
extrinsic parameters stay unaffected.

Overall, the partial derivatives appear to be well-behaved with relative
residuals well below $10^{-4}$. The tidal parameter derivatives perform the
worst, but stay below an error of 10% for most of the samples. Finally, the
$\sim 20$ very poor residuals in chirpmass derivatives actually correspond to
ill-conditioned Fisher matrices and are therefore disregarded for further
analysis.

Figure 2 shows the relative residuals

$\tilde{r}_{x}=\left|\frac{\sigma_{x}^{\text{num}}-\sigma_{x}^{\text{sym}}}{\sigma_{x}^{\text{sym}}}\right|$
(17)

of the $1\sigma$-errors bounds estimated from the FIF using the numeric and
symbolic derivatives of the detector responses.

While the relative residuals on the parameter error estimates from Figure 2 do
not show the actual discrepancy between the two derivative methods, they do
show their effect on the final product of the Fisher analysis. It is important
to note that the number of injections shown has decreased due to some
parameter combinations yielding ill-conditioned Fisher matrices for which the
further analysis is aborted. In fact, it appears that the inclusion of tidal
effects has a stabilizing effect for the FIF, resulting in many more well-
conditioned Fisher matrices. We remark that the mentioned differences in the
behavior of the residuals between tf2 and tf2_tidal could be simply a symptom
of the different input parameter sets which we leave for a future study.

It appears that the correlations between intrinsic or extrinsic parameters
smudges the differences in the numeric and symbolic derivatives within each of
these parameter groups and the inclusion of tidal effects has little influence
on the residuals. The only exception is presented in the case of the
luminosity distance error bounds where the tidal WF resulted in improved
errors. Overall, the relative residuals between the parameter error bounds,
estimated with numeric and symbolic derivative methods, stay under 1% and
mostly under 0.1%. Thus, we can use either approach confidently to perform the
FIF.

Figure 3: The figure shows the histograms of four quantities (SNR $\rho$,
$1\sigma$-fractional error estimates of the chirpmass
$\Delta\mathcal{M}/\mathcal{M}$ and luminosity distance $\Delta D_{L}/D_{L}$,
and the 90%-credible sky area $\Omega$) computed with gwbench for three real
events (GW170817, GW190521, and GW190814). The histograms are generated from
1000 realizations of the sky positions and binary orientations while masses,
spins, and luminosity distances are set to the measured values. The networks
are set the same detector locations, LIGO Hanford, LIGO Livingston, and Virgo,
and PSDs from the actual events. The vertical lines are the measured values
from the respective events GWTC1 ; GW170817 ; GW190521 .

###  V.2 Application to real events

GW170817 and GW190814 are among the loudest events observed by the LIGO and
Virgo detectors so far with SNRs $\rho_{\mathrm{GW170817}}=32$ and
$\rho_{\mathrm{GW190814}}=25$, respectively. These are exceptional candidates
to check the FIF and GW benchmarking against. GW190521 is another interesting
event representing signals of average loudness with
$\rho_{\mathrm{GW190521}}=14$ and very short duration as well as heavy
binaries with merger- and ringdown-dominated signals. Neither of these events
have well determined sky positions and binary orientations. We applied gwbench
to 1000 realizations of the various angles while we fixed the masses, spins,
and luminosity distances of the systems to the reported median values GWTC1 ;
GW170817 ; GW190521 . All three events were detected by three detectors
located at the LIGO Hanford, LIGO Livingston, and Virgo sites. Since the
actual sensitivities of these detectors during the detections are determined
by the facilities, we used the respective PSDs for this analysis GWTC1_PSDs
and ran the FIF for 4 different WF models: tf2_tidal for GW170817, IMRPhenomD
Husa:2015iqa ; Khan:2015jqa  and IMRPhenomHM London:2017bcn  for GW190521, and
tf2, IMRPhenomD, and IMRPhenomHM for G190814.

We compare the reported values for the SNR $\rho$, the fractional error
estimates of the chirpmass $\Delta\mathcal{M}/\mathcal{M}$ and luminosity
distance $\Delta D_{L}/D_{L}$, as well as the 90%-credible sky area $\Omega$
to our benchmarking results. While this comparison contains only four
quantities, we still perform the FIF with respect to 11 (GW190521 and
GW190814) or 12 (GW170817) parameters, in order to not underestimate the error
bounds. These 11 and 12 parameters are the same sets as in Section V.1 above
for tf2 and tf2_tidal, respectively. Besides, in order to allow for an easier
comparison with the measured values we converted the fractional errors from
the standard $1\sigma$-output of the FIF to the respective 90% credible
bounds.

The benchmarking posteriors for the 1000 realizations are presented in Figure
3 indicating that they capture the measured errors well for at least one of
the WF models per event. The figure shows that tf2_tidal performs really well
for an inspiral dominated, binary neutron star signal like GW170817. This does
not hold for signals from the heavier binaries. For such systems, the error
estimates from tf2 do not capture the measured values due to missed post-
inspiral contributions (GW190814), while the FIF completely failed for a
strongly merger- and ringdown-dominated signal (GW190521). In fact, for both
of these systems it was necessary to use IMRPhenomHM, a WF with higher modes,
and not just IMRPhenomD in order to fully capture the measured errors.

These findings are in line with the WFs used for the actual parameter
estimation analyses of these events and they show that gwbench provides a
consistent tool for the purposes of GW benchmarking.

##  VI Example applications

gwbench has already been applied in two works: Borhanian:2020vyr ;
Gupta:2020lxa . In Borhanian:2020vyr  we examined the possibility to use dark
sirens—GW sources without an electromagnetic counterpart—to resolve the
tension in the measurement of the Hubble constant. In Gupta:2020lxa  we
studied how multiband observations of GW signals with the space-based LISA and
ground-based 3G observatories will open up new possibilities to conduct
multiparameter tests of general relativity in the line of Arun:2006hn .

###  VI.1 “Dark sirens to resolve the Hubble-Lemaître tension”

We examined the prospect to use dark sirens to perform precision measurements
of the Hubble constant and thus provide a GW-assisted method to resolve the
Hubble-Lemaître tension. The proposed method makes use of the potential of
future observatories to measure sky positions of close-by dark sirens well
enough to pin-point the host galaxy of the GW event. This provides two
independent measurements, the distance from the GW signal and the galaxy’s
redshift from electromagnetic follow-up observations, which compound to
estimate the Hubble constant.

We simulated three populations of dark sirens: one resembled binary black
holes from the GWTC-1 catalog GWTC1 , the second focused on the heavy sub-
population of such BBHs with component masses $m_{1},m_{2}>25\,M_{\odot}$,
while the third was GW190814-like GW190814 . The GWTC-1 catalog summarizes the
detections from the first two observing runs of the LIGO and VIRGO detectors
and GW190814 was the most asymmetric binary observed to date.

For each population we performed the Fisher analysis in four networks
containing three possible detector technologies of the future (’Advanced+’,
’Voyager’, and 3G). In order to leverage the sensitivities of such future
networks we used the lalsimulation WF IMRPhenomHM London:2017bcn  which also
models subdominant, higher-order spherical harmonic modes beyond the
quadrupole. These higher modes improve the distance measurement by breaking
the degeneracy between the luminosity distance and the inclination angle of
the binary orbit to the line of sight, thus ultimately yielding better
measurements of the Hubble constant.

Using gwbench we performed the FIF with numeric derivatives for $10^{4}$
binaries per population. We estimated the number of detections to expect from
the studied networks that yield precise measurements of the dark siren’s sky
position and propagated the error bounds on the luminosity distance
measurement to estimate the accuracy of a Hubble constant measurement. We
showed that, if redshift errors are well controlled, dark siren observations
will yield $\lesssim 2\%$-accurate Hubble constant measurements already in the
Advanced+ era and thus resolve the Hubble-Lemaître tension.

###  VI.2 “Multiparameter tests of general relativity using multiband
gravitational-wave observations”

Parameterized post-Newtonian tests Arun:2006hn  present a crucial resource to
perform theory-agonistic tests of general relativity which can probe
generically potential deviations of GW signals from Einstein’s theory. This
class of tests is especially important when verifiable GW predictions from
alternate theories of gravity are sparse. The tests are based on modifications
to the expansion coefficients in the post-Newtonian approximation to the WFs
emitted by compact binary coalescences, such that these deviations vanish if
general relativity is correct. Unfortunately, the current tests are
handicapped by the strong correlations between the introduced deviation
parameters. Therefore, parameterized post-Newtonian tests are limited to
single-coefficient modifications for now. We showed that the combination of
the information present in the low-frequency LISA and high-frequency 3G bands
breaks this degeneracy. Thus such multiband observations allow for true
multiparameter tests of general relativity.

In our study we examined stellar mass binary black hole coalescences. Their
inspiral signals will last several years in the LISA band ($\sim
0.1-100\,\mathrm{mHz})$ before they reach the 3G band ($\sim$ Hz-kHz) where
they finally merge. Unfortunately, these signals, while extremely loud in 3G
networks, will be faint to non-visible in LISA, making the ‘joint’ observation
a difficult task. The issue lies in the detection methodology for such quiet
signals. The current state-of-the-art are blind, matched filtering searches
gstlal_live ; pycbc_live  without a-priori knowledge of the GW signal. Such
blind searches require GW template banks which cover a large volume in the
parameter space. Thus, the sizes of the banks become computationally
infeasible for the signals of interest in the LISA band Moore:2019pke .

Fortunately, we can leverage the high-fidelity detections with 3G detectors to
dig out these quiet signals from LISA’s archival data. In contrast to a blind
search the 3G detection confirms the existence of the signal we are searching
for and further provides information about it, allowing for the generation of
signal-specific template banks that decrease the computational cost immensely.
This is enabled in two-ways: A 3G network will estimate all WF parameters
except chirpmass and symmetric mass ratio to better precision than LISA could,
thus decreasing the parameter dimension of the template bank to two and fixing
the other parameters to the 3G measurement. The search volume can be further
narrowed down to within the error bounds estimated from the 3G detection for
both chirpmass and symmetric mass ratio.

Using gwbench, we tackled the question of feasibility to search the archival
LISA data following a 3G detection, strengthening our proposal to perform
multiband multi-parameter tests of general relativity. For this purpose we
generated a set of 500,000 binary black holes up to a redshift of $z=10$.
Roughly 200 of these binaries should be visible in LISA with an SNR $\geq 4$
which we set as the threshold for visibility in the archival search. gwbench
allowed us to estimate the SNR and $1\sigma$ error bounds on the WF parameters
for a network consisting of one Einstein Telescope and two $40\,\text{km}$,
CE2 detectors in compact binary optimization.

Given the LISA SNR $\rho_{\text{LISA}}$ and 3G errors on chirpmass
$\sigma_{M_{c}}$ and symmetric mass ratio $\sigma_{\eta}$, we used gwbench to
calculate the detector-dependent volume $V$ via the metric $g$ Owen:1995tm ;
Owen:1998dk

$g=\rho_{\text{LISA}}^{-2}\det\left(\langle\partial_{x}H,\partial_{y}H\rangle_{x,y\in\\{M_{c},\eta\\}}\right)$
(18)

as

$V=\int\limits_{M_{c}-2\sigma_{M_{c}}}^{M_{c}+2\sigma_{M_{c}}}\mathrm{d}M_{c}\int\limits_{\eta-2\sigma_{\eta}}^{\eta+2\sigma_{\eta}}\mathrm{d}\eta\,\sqrt{g}.$
(19)

The inner product in the definition of the metric is taken with respect to the
LISA PSD. Finally, the number of templates $N$ is given by the ratio of the
total volume to the volume of the error ball around each template for a given
minimal match $mm$

$N=\frac{V}{1-mm}.$ (20)

Using $mm=0.95$, we found that such targeted searches of LISA’s archival data
following the detection of a binary black hole merger in a network of 3G
detectors should allow for template bank sizes of the order $10^{3}$ –
$10^{4}$. This is substantially less than the value of $10^{12}$ predicted
before Moore:2019pke . Ultimately, gwbench enabled us to check the feasibility
of multibanding to perform multi-parameter tests of general relativity by
providing the tools to calculate estimates of both the 3G errors bounds and
the LISA template numbers.

##  VII Summary

In this work we presented gwbench, a Python package for fast and
straightforward applications of the Fisher information formalism for the
purposes of GW benchmarking. The package is written with the usage for
detector networks in mind. Thus, it is structured around the Network class
which acts as a hub for user interaction and facilitates all the possible
computations. The package gives easy access to a variety of GW benchmarking
quantities such as the frequency-domain WF polarizations, antenna pattern
functions and location phase factors as well as noise PSDs for various
detector locations and technologies, frequency-domain detector responses,
SNRs, Fisher and covariance matrices, and the product of the FIF, the
$1\sigma$ error bounds for input parameters of choice.

There are two types of WFs available, both of which support numeric
differentiation via the numdifftools package: LAL WFs as well as those
implemented gwbench. The latter further support symbolic derivatives via sympy
which are faster and more accurate. The inclusion of LAL WF support is very
important as it gives access to a host of WF models that are well tested and
commonly used in the literature. gwbench contains nine existing, planned, or
fiducial detector locations and 23 different detector technologies ranging
from aLIGO to various CE sensitivity curves. The addition of user-defined
locations and technologies is also possible. Another specialty of this package
is the inclusion of the effects of Earth’s rotation in the antenna patterns
and detector responses, and thus the FIF. This is particularly important for
long-lived GW signals such as from binary neutron stars.

gwbench is and has been used for several research projects of which we
presented two that are in public already. In both cases, gwbench allowed us to
generate and calculate parameter errors and even template numbers for a
matched filtered search for large numbers of sources and detector networks.
gwbench represents a unique code that implements the FIF and GW benchmarking
in an easy-to-use manner for a wide range of problems, WFs, and network
options and allows for straightforward application in large-scale simulations
for improved statistics of the scientific claims.

###### Acknowledgements.

I thank Bangalore Sathyaprakash, Anuradha Gupta, Patrick Godwin, and everyone
in the Cosmic Explorer Project for helpful comments and discussions during the
development of gwbench. I further thank Kevin Kuns, Varun Srivastava, Evan
Hall, Matthew Evans, and Stefan Ballmer for providing the sensitivity curves
used in gwbench. I also thank all front line workers combating the CoVID-19
pandemic and those maintaining ordinary life without whose support this work
would not have been possible. I acknowledge the support by NSF grant
PHY-1836779. Computing resources for the development and testing of this code
were provided by the Pennsylvania State University.

##  Appendix A Network and Detector classes

GW detector networks consist, as per name, of a number of detectors that work
in tandem to observe GWs to improve signal strength, sky coverage, and
parameter estimation. Leaning on this understanding, the Network class is
designed as a GW benchmarking hub for a chosen set of detectors. The class
structure allows for straightforward handling of the relevant quantities and
network methods and thus provides the user with tools to load, calculate, and
manipulate quantities like waveforms, SNRs, or errors estimates.

The Detector class is written with similar syntax and methodology as the
Network class. This is very natural because network quantities like SNR or
error estimates are handled in the same manner as the respective single
detector ones. This choice allows for clearer codes and faster development,
but also better understanding of both classes, which helps when using Network
objects. Thus, Table 1 compares the _instance variables_ of both classes side-
by-side. The top section of the table lists the variables which are set by the
user and tell the Network

  * •

the label of the Network and what location and technology to use for a
Detector instance,

  * •

which Detector instances to initialize,

  * •

what frequency array and injection parameters to use for function evaluation,
and with respect to which parameters to take derivatives for error estimation
via FIF,

  * •

what Waveform to initialize,

  * •

which parameters $x$ to convert to $\cos(x)$ or $\ln(x)$ to estimate errors on
the latter, and whether to take Earth’s rotation into account.

The bottom section shows the variables which can be loaded and computed with
instance methods based on the user-set variables:

  * •

Detector-specific quantities like PSDs, antenna patterns, and location phase
factors are loaded and stored in each Detector instance.

  * •

The WF polarizations and the respective derivative dictionaries are detector-
agnostic and saved in the Network.

  * •

Each Detector handles the respective detector response and derivative
dictionaries.

  * •

The detector and network SNRs (and their squares) are stored in the respective
instances, where each Detector also saves the integrand of the square SNR.

  * •

Each detector and the network contain their respective copy of the FIM, its
condition number, a well-conditioned flag, the covariance matrix, the
inversion error, and an error dictionary containing the $1\sigma$ errors
extracted from the covariance matrix.

##  Appendix B Implemented waveform models and the Waveform class

Table 1: Network vs Detector variables. |  Network |  Detector
---|---|---
_set by user:_
labels |  .label (label to initialize network) |  .det_key, .tec, .loc (full detector key, technology and location labels)
detectors |  .detectors, .det_keys (list of Detector objects and detector keys) |  —
global quantities |  .f, .inj_params, .deriv_symbs_string, .deriv_variables (frequency array, injection parameters, derivative parameters string and list) |  .f (truncated frequency array)
waveform |  .wf (Waveform object) |  —
analysis settings |  .conv_cos, .conv_log, .use_rot (lists of parameters $x$ to convert to $\cos(x),\ln(x)$, consider Earth’s rotation) |  —
_loaded or calculated by network methods:_
detector specific quantities |  — |  .psd, .Fp, .Fc, .Flp (detector PSD, antenna patterns, location phase factor)
waveform polarizations |  .hfp, .hfc, .del_hfpc, .del_hfpc_expr (plus and cross polarizations, derivatives dictionary, sympy expressions dict) |  —
detector responses |  — |  .hf, .del_hf, .del_hf_expr (detector response, derivatives dictionary, sympy expressions dict)
SNR |  .snr, .snr_sq (SNR, SNR square) |  .snr, .snr_sq, .d_snr_sq (SNR, SNR square, SNR square integrand)
error quantities |  .fisher, .cond_num, .wc_fisher, .cov, .inv_err, .errs (Fisher information matrix, Fisher condition number, well conditioned boolean, covariance matrix, inversion error, the parameter errors) |  .fisher, .cond_num, .wc_fisher, .cov, .inv_err, .errs (Fisher information matrix, Fisher condition number, well conditioned boolean, covariance matrix, inversion error, the parameter errors)
Table 2: The available WFs are presented here, together with their labels, their input parameters, and whether they allow for numeric or symbolic differentiation. gwbench uses the same labeling convention for ’approximant’ in the LAL WFs as specified in lalsimulation. Check the _Approximant_ -list for frequency-domain approximants, implemented in lalsimulation, under _Enumerations_ at https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_inspiral__h.html.  Waveforms |  Waveform label |  Input parameters |  Other parameter dictionary
---|---|---|---
_LAL WFs — only numeric differentiation_
BBH |  ’lal_bbh’ |  ’f Mc eta chi1x chi1y chi1z chi2x chi2y chi2z DL tc phic iota’ |  {’approximant’:name_string}
BNS |  ’lal_bns’ |  ’f Mc eta chi1x chi1y chi1z chi2x chi2y chi2z DL tc phic iota lam_t delta_lam_t’ |  {’approximant’:name_string}
_Coded WFs — numeric and symbolic differentiation_
TaylorF2 |  ’tf2’ |  ’f Mc eta chi1z chi2z DL tc phic iota’ |  None
TaylorF2 + tidal |  ’tf2_tidal’ |  ’f Mc eta chi1z chi2z DL tc phic iota lam_t delta_lam_t’ |  None

gwbench handles WF selection and evaluation within the Waveform class. The
available WFs are listed in Table 2 together with their labels, their input
parameters, and whether they allow for numeric or symbolic differentiation.
One of the strong suites of gwbench is the implementation of WF models from
LAL in the Waveform class and the capability to use them for the FIF with
numeric derivatives. The class contains five instance variables

  * •

.wf_model_name: internal name of the WF model,

  * •

.wf_other_var_dic: extra input parameters,

  * •

.wf_symbs_string: string containing the basic parameters (space-separated),

  * •

.hfpc_np: numpy function,

  * •

.hfpc_sp: sympy function,

and few methods of which two are of particular interest to the user:

  * •

.get_sp_expr(): load the sympy function to obtain a symbolic expression which
can be manipulated,

  * •

.eval_np_func(f,inj_params): evaluates the numpy function for the given
frequency array f and injection parameters inj_params.
(The injection parameters inj_params can be either provided as a dictionary
with keys equal to the parameters defined in .wf_symbs_string or as a list
with the same ordering given by .wf_symbs_string.)

Table 3: The available detector locations are presented together with the respective labels and detector angles (longitude $\alpha_{d}$, latitude $\beta_{d}$, and orientation of the y-arm with respect to due East $\gamma_{d}$). The five locations H, L, V, K, I are chosen in accordance with the implementation in LAL (https://lscsoft.docs.ligo.org/lalsuite/lal/_l_a_l_detectors_8h_source.html). The four locations E, C, N, and S are fiducial example sites and do not represent candidates under active consideration. Label | Location | angles (rad)
---|---|---
$\alpha_{d}$ | $\beta_{d}$ | $\gamma_{d}$
H | LIGO Hanford (Washington, USA) | $-2.084060$ | $0.810795$ | $-2.513290$
L | LIGO Livingston (Louisiana, USA) | $-1.584310$ | $0.533423$ | $-1.261590$
V | Virgo (Cascina, Italy) | $0.183338$ | $0.761512$ | $2.802430$
K | KAGRA (Kamioka, Japan) | $2.394200$ | $0.632682$ | $2.087480$
I | LIGO India (Hingoli, India) | $1.334013$ | $0.248418$ | $1.570800$
E | Einstein Telescope (Cascina, Italy) | $0.183338$ | $0.761512$ | $2.802430$
C | Cosmic Explorer (main) (Idaho, USA) | $-1.969170$ | $0.764918$ | $0.000000$
N | Cosmic Explorer North (New Mexico, USA) | $-1.858430$ | $0.578751$ | $-1.047200$
S | Cosmic Explorer South (New South Wales, Australia) | $2.530730$ | $-0.593412$ | $0.785398$

Only the wf_model_name and wf_other_var_dic are needed to fully initialize a
Waveform instance:

⬇

from gwbench import wf_class

wf = wf_class.Waveform(

wf_model_name = wf_model_name,

wf_other_var_dic = wf_other_var_dic)

The LAL WF models do not provide any of the aforementioned symbolic
functionality, enabled by sympy, as this requires a direct implementation of
the WF model to ensure proper sympy function calls.

##  Appendix C Implemented detector locations - antenna patterns and location
phase factor

Table 3 summarizes the nine available detector locations in gwbench together
with the respective labels used in the code, while Figure 4 presents their
spread across the world map. There are five 2G detector sites at LIGO Hanford
(Washington, USA), LIGO Livingston (Louisiana, USA), VIRGO (Cascina, Italy),
KAGRA (Kamioka, Japan), and LIGO India (Hingoli, India) as well as four
fiducial sites for the 3G detectors. The Einstein Telescope is set to the same
coordinates as VIRGO, while the three Cosmic Explorer locations Main, North,
and South are set to sites in Idaho (USA), New Mexico (USA), and New South
Wales (Australia)444The four fiducial 3G detector locations are example sites
and do not represent candidates under active consideration..

Figure 4: The available detector locations shown on the world map.

Figure 5: The _top_ shows the amplitude spectral densities (ASD)
$\sqrt{S_{n}(h)}$ available in gwbench, see Figure 6, as functions of the
frequency between 1 and 4096 Hz. The _bottom_ presents the observable
redshifts for equal-mass binaries at SNR 100 (solid) and 12 (dotted) as
functions of the total binary mass between 1 and 1000 $M_{\odot}$ (in the
source frame). The systems are assumed to be optimally orientated for an
L-shaped detector with the respective ASD. The SNRs were obtained for the LAL
WF IMRPhenomHM. The shaded regions show the respective ranges for events with
SNRs between 12 and 100.

The computation of the antenna patterns $F_{+},F_{\times}$ and location phase
factors $F_{lp}$ for these locations is implemented in the two modules
antenna_pattern_np.py and antenna_pattern_sp.py which differ in their intended
usage: the former is used to calculate the respective quantities directly
whereas the latter is required for the generation of the symbolic derivatives
via sympy.

The core concepts for a single detector are demonstrated in
LRR_Sathyaprakash_2009 , whereas Schutz:2011tw  outlines the more intricate
methods required for the analysis of a network of detectors where calculations
are carried out in a universal frame centered at Earth’s center. In either
case, $F_{+}$, $F_{\times}$, and $F_{lp}$ depend on a few signal-specific
quantities: right ascension $\alpha$, declination $\delta$, and polarization
angle $\psi$. In a network, these functions further depend on the arrival time
$t_{E}$ of the signal at Earth’s center as well as three detector- and
location-specific angles: the longitude $\alpha_{d}$ and latitude $\beta_{d}$
of the detector site as well as the angle $\gamma_{d}$ of the detector’s y-arm
with respect to due East555Since the Einstein Telescope consists of three
V-shaped detectors ET1, ET2, and ET3, that form a equilateral triangle,
$\gamma_{d}$ is fiducially set to correspond to ET1. The fixed geometry allows
the calculation of the detector quantities for all three detectors..

All these considerations omit the effects of Earth’s rotation on the antenna
patterns and location phases factors. These effects are negligible for short-
lived, transient signals lasting only a few minutes or less. They cannot be
neglected for 3G detectors which will see signals that can be up to a few days
long. Hence, gwbench contains the functionality to take Earth’s rotation,
simplified as a constant circular rotation with a period of one day, into
account which follows the calculations detailed in Wen:2010cr ; Zhao:2017cbb .

##  Appendix D Implemented detector technologies - detector PSDs

Figures 6 and 5 show the 23 available PSDs trade_study_psd  in gwbench. The
former is a graphical summary presenting the hierarchy of these PSDs in terms
of detector generation (2G and 3G) and design concepts (Advanced+ vs Voyager
vs ET vs CE and more) while the latter shows the amplitude spectral densities
$\sqrt{S_{n}(h)}$ as functions of frequencies. Figure 5 further presents the
observable redshift ranges for equal-mass binaries with total mass between 1
and 1000 $M_{\odot}$ (source frame), if observed with SNRs between 12 and 100.
Table 4 presents the respective label and key options for these PSDs.

##  Appendix E Derivative dictionaries

gwbench obtains the partial derivatives of the WF polarizations and detector
responses either numerically or symbolically using the numdifftools or sympy
packages, respectively. The former calculates the evaluated derivatives
directly which is computationally expensive but applicable to most WFs. sympy
on the other hand works with symbolic expressions that require modifications
to the WF source code. In return they enable the analytic computation of
generic derivative expressions which can be stored in the form of _lambdified
functions_ for fast evaluation in run time.

Table 4: The available detector PSDs are presented together with the
respective labels and internal keys as used by gwbench. CBO stands for
compact-binary optimized and PMO for post-merger optimized. The former PSD
will be tuned towards increased sensitivity at low frequencies, while the
latter is geared towards high frequencies. There are two examples illustrating
the nomenclature of the CE configurations.

Label | Internal keys | PSD
---|---|---
- | aLIGO | aLIGO
+ | A+,V+,K+ | A+ (adV+/KAGRA+ for V/K location)
vc | Voyager-CBO | Voyager, CBO
vp | Voyager-PMO | Voyager, PMO
E | ET | ET-D (same as location)
sLo | <s>-<L>-<o> | CE configurations (see below)

CE configurations - <s>-<L>-<o>
---
| stage s | length L | optimizations o
Label option | i / a | 1 / 2 / 3 / 4 | c / p
Key option |  CE1 / CE2 |  10 / 20 / 30 / 40 |  CBO / PMO
PSD setting | CE1 / CE2 |  $(10$ / $20$ / $30$ / $40)\,\text{km}$ | CBO / PMO

CE configurations - Examples
---
Label | Internal key | PSD
i2p | CE1-20-PMO | CE1, 20 km, PMO
a4c | CE2-40-CBO | CE2, 40 km, CBO

The evaluated partial derivatives are stored in dictionaries which are
structured as follows

⬇

{’del_variable1_hf’:<numpy array>,

’del_variable2_hf’:<numpy array>,

... }

in case of the detector responses and

⬇

{’del_variable1_hfp’:<numpy array>,

’del_variable1_hfc’:<numpy array>,

’del_variable2_hfp’:<numpy array>,

’del_variable2_hfc’:<numpy array>,

... }

for the WF polarizations. These dictionaries have the same structure in case
of the symbolic derivative expressions, but store lambdified sympy functions
instead of numpy arrays. Further, they contain two extra keys, ’variables’ and
’deriv_variables’, that store the input parameters and the differentiation
variables.

##  Appendix F injections module

Figure 6: Detector technologies/PSDs available in gwbench. The included 2G
PSDs encompass aLIGO as well as the planned Advanced+ and proposed Voyager
upgrades. The 3G PSDs encompass the ET-D curve for the Einstein Telescope and
various Cosmic Explorer configurations, subdivided in the CE1 and CE2 stages.
CBO and PMO stand for compact-binary or post-merger optimized, i.e. focus on
low- or high-frequency sensitivity, respectively.

The module injections contains functions to randomly sample input parameters
appropriate for GW problems: chirpmass and symmetric mass ratio, binary
component spins, sky location and binary orientation angles, as well as
redshifts and luminosity distances. The inputs of the sampler functions are
the same across the board: dictionaries specifying details of the parameter
space to sample, the number of injections num_injs to sample, and the seed
value seed for the random number generator. In the following we show the usage
of these functions and explain what the respective dictionaries contain.

#### Mass sampler:

The chirpmass and symmetric mass ratio sampling is handled by mass_sampler:

⬇

from gwbench import injections

Mcs, etas = injections.mass_sampler(

mass_dict,num_injs,seed)

The mass_dict contains three definite keys:

  * •

’dist’: specifies the mass distribution (values: ’gaussian’, ’power’,
’power_uniform’, and ’uniform’).

  * •

’mmin’: specifies the minimum mass.

  * •

’mmax’: specifies the maximum mass.

The Gaussian and power-law distributions need further inputs, namely ’mean’
and ’sigma’ or ’alpha’, respectively. The first two specify the mean and
standard deviation for the Gaussian distribution while the ’alpha’ is the
power-law index.

#### Spin sampler:

The binary spin components sampling is handled by spin_sampler:

⬇

chi1x, chi1y, chi1z, chi2x, chi2y, chi2z =

injections.spin_sampler(spin_dict,num_injs,seed)

The spin_dict contains four keys.

  * •

’geom’: specifies the geometry (values: ’cartesian’ and ’spherical’).

  * •

’dim’: specifies spin dimension (values: 1 and 3; one-dimensional spin
sampling disregards ’geom’ and samples using Cartesian geometry; further 1
returns zeroes for the x- and y-components to sample aligned injections).

  * •

’chi_lo’: specifies the minimum spin value (smallest spin value in each
dimension for ’cartesian’ geometry or smallest spin magnitude for ’spherical’
geometry).

  * •

’chi_hi’: specifies the maximum spin value (largest spin value in each
dimension for ’cartesian’ geometry or largest spin magnitude for ’spherical’
geometry).

#### Angle sampler:

The sky location and orientation angles sampling is handled by angle_sampler:

⬇

iotas, ras, decs, psis = injections.angle_sampler(

num_injs,seed)

#### Redshift and luminosity distance sampler:

The redshift and luminosity distance sampling is handled by
redshift_lum_distance_sampler:

⬇

zs, DLs = injections.redshift_lum_distance_sampler(

cosmo_dict,num_injs,seed)

The cosmo_dict contains three definite keys

  * •

’sampler’: specifies the redshift sampling method (values: ’uniform’,
’uniform_comoving_volume_inversion’ and ’uniform_comoving_volume_rejection’;
the inversion method is faster, while the rejection method is more accurate,
especially for very jagged redshift distributions).

  * •

’zmin’: specifies the minimum redshift.

  * •

’zmax’: specifies the maximum redshift.

By default the code samples redshifts from the redshift distribution
astropy.cosmology.Planck18 provided in the astropy package. The user can
change this behavior by adding three further keys to the cosmo_dict whose
values are then passed to astropy.cosmology.LambdaCDM:

  * •

’H0’: specifies Hubble constant.

  * •

’Om0’: specifies the mass density $\Omega_{m,0}$ at the current time.

  * •

’Ode0’: specifies dark energy density $\Omega_{\Lambda,0}$ at the current
time.

#### Injection parameter sampler:

injections further contains a function that combines the four other samplers
into one function call: injections_CBC_params_redshift. Hence it depends on
the same input as the previous four methods:

⬇

data = injections.injections_CBC_params_redshift(

cosmo_dict,mass_dict,spin_dict,redshifted,

num_injs,seed,file_path)

The two additional inputs, redshifted and file_path, control whether the
sampled masses are already redshifted (1) according to the sampled redshifts
or not (0) as well as if the sampled data should be stored in a file at the
specified path (default is None), respectively.

## References

  * (1) LIGO Scientific, Virgo collaboration, _Observation of Gravitational Waves from a Binary Black Hole Merger_ , _Phys. Rev. Lett._ 116 (2016) 061102 [1602.03837].
  * (2) LIGO Scientific, Virgo collaboration, _GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral_ , _Phys. Rev. Lett._ 119 (2017) 161101 [1710.05832].
  * (3) LIGO Scientific, 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 (2019) 031040 [1811.12907].
  * (4) L. London, S. Khan, E. Fauchon-Jones, C. García, M. Hannam, S. Husa et al., _First higher-multipole model of gravitational waves from spinning and coalescing black-hole binaries_ , _Phys. Rev. Lett._ 120 (2018) 161102 [1708.00404].
  * (5) R. Cotesta, A. Buonanno, A. Bohé, A. Taracchini, I. Hinder and S. Ossokine, _Enriching the Symphony of Gravitational Waves from Binary Black Holes by Tuning Higher Harmonics_ , _Phys. Rev. D_ 98 (2018) 084028 [1803.10701].
  * (6) LIGO Scientific collaboration, _Advanced LIGO_ , _Class. Quant. Grav._ 32 (2015) 074001 [1411.4547].
  * (7) VIRGO collaboration, _Advanced Virgo: a second-generation interferometric gravitational wave detector_ , _Class. Quant. Grav._ 32 (2015) 024001 [1408.3978].
  * (8) KAGRA collaboration, _Interferometer design of the KAGRA gravitational wave detector_ , _Phys. Rev. D_ 88 (2013) 043007 [1306.6747].
  * (9) C. for Indian Initiative in Gravitational-wave Observations (IndIGO), _Ligo-india, proposal_ , (2011), https://dcc.ligo.org/LIGO-M1100296/public.
  * (10) LSC, _Instrument science white paper_ , (2016), https://dcc.ligo.org/LIGO-T1600119/public.
  * (11) KAGRA, LIGO Scientific, VIRGO collaboration, _Prospects for Observing and Localizing Gravitational-Wave Transients with Advanced LIGO, Advanced Virgo and KAGRA_ , _Living Rev. Rel._ 21 (2018) 3 [1304.0670].
  * (12) eLISA collaboration, _The Gravitational Universe_ , [1305.5720].
  * (13) M. Maggiore et al., _Science Case for the Einstein Telescope_ , _JCAP_ 03 (2020) 050 [1912.02622].
  * (14) LIGO collaboration, _A Cryogenic Silicon Interferometer for Gravitational-wave Detection_ , _Class. Quant. Grav._ 37 (2020) 165003 [2001.11173].
  * (15) M. Punturo et al., _The Einstein Telescope: A third-generation gravitational wave observatory_ , _Class. Quant. Grav._ 27 (2010) 194002.
  * (16) M. Punturo et al., _The third generation of gravitational wave observatories and their science reach_ , _Class. Quant. Grav._ 27 (2010) 084007.
  * (17) D. Reitze et al., _Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO_ , _Bull. Am. Astron. Soc._ 51 (2019) 035 [1907.04833].
  * (18) 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 (1994) 2658 [gr-qc/9402014].
  * (19) E. Poisson and C. M. Will, _Gravitational waves from inspiraling compact binaries: Parameter estimation using second postNewtonian wave forms_ , _Phys. Rev. D_ 52 (1995) 848 [gr-qc/9502040].
  * (20) M. Vallisneri, _Use and abuse of the Fisher information matrix in the assessment of gravitational-wave parameter-estimation prospects_ , _Phys. Rev. D_ 77 (2008) 042001 [gr-qc/0703086].
  * (21) A. J. K. Chua and M. Vallisneri, _Learning Bayesian posteriors with neural networks for gravitational-wave inference_ , _Phys. Rev. Lett._ 124 (2020) 041102 [1909.05966].
  * (22) H. Gabbard, C. Messenger, I. S. Heng, F. Tonolini and R. Murray-Smith, _Bayesian parameter estimation using conditional variational autoencoders for gravitational-wave astronomy_ , [1909.06296].
  * (23) S. R. Green and J. Gair, _Complete parameter inference for GW150914 using deep learning_ , _Mach. Learn. Sci. Tech._ 2 (2021) 03LT01 [2008.03312].
  * (24) S. R. Green, C. Simpson and J. Gair, _Gravitational-wave parameter estimation with autoregressive neural network flows_ , _Phys. Rev. D_ 102 (2020) 104057 [2002.07656].
  * (25) L. S. Collaboration et al., _Lsc algorithm library suite_ , (2020), https://doi.org/10.7935/GT1W-FZ16.
  * (26) K. Wette, _SWIGLAL: Python and octave interfaces to the lalsuite gravitational-wave data analysis libraries_ , _Accepted in SoftwareX_ (2020) .
  * (27) B. Sathyaprakash et al., _Extreme Gravity and Fundamental Physics_ , [1903.09221].
  * (28) V. Kalogera et al., _Deeper, Wider, Sharper: Next-Generation Ground-Based Gravitational-Wave Observations of Binary Black Holes_ , [1903.09220].
  * (29) B. Sathyaprakash et al., _Multimessenger Universe with Gravitational Waves from Binaries_ , [1903.09277].
  * (30) B. Sathyaprakash et al., _Cosmology and the Early Universe_ , [1903.09260].
  * (31) L. Barack et al., _Black holes, gravitational waves and fundamental physics: a roadmap_ , _Class. Quant. Grav._ 36 (2019) 143001 [1806.05195].
  * (32) B. Sathyaprakash and B. Schutz, _Physics, Astrophysics and Cosmology with Gravitational Waves_ , _Living Rev. Rel._ 12 (2009) 2 [0903.0338].
  * (33) L. Wen and Y. Chen, _Geometrical Expression for the Angular Resolution of a Network of Gravitational-Wave Detectors_ , _Phys. Rev. D_ 81 (2010) 082001 [1003.2504].
  * (34) W. Zhao and L. Wen, _Localization accuracy of compact binary coalescences detected by the third-generation gravitational-wave detectors and implication for cosmology_ , _Phys. Rev. D_ 97 (2018) 064031 [1710.05325].
  * (35) C. L. Rodriguez, I. Mandel and J. R. Gair, _Verifying the no-hair property of massive compact objects with intermediate-mass-ratio inspirals in advanced gravitational-wave detectors_ , _Phys. Rev. D_ 85 (2012) 062002 [1112.1404].
  * (36) B. S. Sathyaprakash and S. V. Dhurandhar, _Choice of filters for the detection of gravitational waves from coalescing binaries_ , _Phys. Rev. D_ 44 (1991) 3819.
  * (37) E. Poisson, _Gravitational waves from inspiraling compact binaries: The Quadrupole moment term_ , _Phys. Rev. D_ 57 (1998) 5287 [gr-qc/9709032].
  * (38) B. Mikoczi, M. Vasuth and L. A. Gergely, _Self-interaction spin effects in inspiralling compact binaries_ , _Phys. Rev. D_ 71 (2005) 124043 [astro-ph/0504538].
  * (39) K. Arun, A. Buonanno, G. Faye and E. Ochsner, _Higher-order spin effects in the amplitude and phase of gravitational waveforms emitted by inspiraling compact binaries: Ready-to-use gravitational waveforms_ , _Phys. Rev. D_ 79 (2009) 104023 [0810.5336].
  * (40) A. Buonanno, B. Iyer, E. Ochsner, Y. Pan and B. Sathyaprakash, _Comparison of post-Newtonian templates for compact binary inspiral signals in gravitational-wave detectors_ , _Phys. Rev. D_ 80 (2009) 084043 [0907.0700].
  * (41) A. Bohé, S. Marsat and L. Blanchet, _Next-to-next-to-leading order spin–orbit effects in the gravitational wave flux and orbital phasing of compact binaries_ , _Class. Quant. Grav._ 30 (2013) 135009 [1303.7412].
  * (42) A. Bohé, G. Faye, S. Marsat and E. K. Porter, _Quadratic-in-spin effects in the orbital dynamics and gravitational-wave energy flux of compact binaries at the 3PN order_ , _Class. Quant. Grav._ 32 (2015) 195010 [1501.01529].
  * (43) C. K. Mishra, A. Kela, K. Arun and G. Faye, _Ready-to-use post-Newtonian gravitational waveforms for binary black holes with nonprecessing spins: An update_ , _Phys. Rev. D_ 93 (2016) 084054 [1601.05588].
  * (44) Astropy collaboration, _Astropy: A Community Python Package for Astronomy_ , _Astron. Astrophys._ 558 (2013) A33 [1307.6212].
  * (45) Astropy collaboration, _The Astropy Project: Building an Open-science Project and Status of the v2.0 Core Package_ , _Astron. J._ 156 (2018) 123 [1801.02634].
  * (46) P. A. Brodtkorb and J. D’Errico, _Documentation: numdifftools_ , (2019), https://numdifftools.readthedocs.io/en/latest/index.html.
  * (47) LIGO Scientific, Virgo collaboration, _GW190521: A Binary Black Hole Merger with a Total Mass of $150~{}M_{\odot}$_, _Phys. Rev. Lett._ 125 (2020) 101102 [2009.01075].
  * (48) LIGO Scientific, Virgo collaboration, _GW190814: Gravitational Waves from the Coalescence of a 23 Solar Mass Black Hole with a 2.6 Solar Mass Compact Object_ , _Astrophys. J. Lett._ 896 (2020) L44 [2006.12611].
  * (49) B. D. Lackey and L. Wade, _Reconstructing the neutron-star equation of state with gravitational-wave detectors from a realistic population of inspiralling binary neutron stars_ , _Phys. Rev. D_ 91 (2015) 043002 [1410.8866].
  * (50) L. Wade, J. D. Creighton, E. Ochsner, B. D. Lackey, B. F. Farr, T. B. Littenberg 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 (2014) 103012 [1402.5156].
  * (51) LSC, _Power spectral densities (psd) release for gwtc-1_ , (2019), https://dcc.ligo.org/LIGO-P1900011/public.
  * (52) S. Husa, S. Khan, M. Hannam, M. Pürrer, F. Ohme, X. Jiménez Forteza et al., _Frequency-domain gravitational waves from nonprecessing black-hole binaries. I. New numerical waveforms and anatomy of the signal_ , _Phys. Rev. D_ 93 (2016) 044006 [1508.07250].
  * (53) S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer, X. Jiménez Forteza et al., _Frequency-domain gravitational waves from nonprecessing black-hole binaries. II. A phenomenological model for the advanced detector era_ , _Phys. Rev. D_ 93 (2016) 044007 [1508.07253].
  * (54) S. Borhanian, A. Dhani, A. Gupta, K. G. Arun and B. S. Sathyaprakash, _Dark Sirens to Resolve the Hubble–Lemaître Tension_ , _Astrophys. J. Lett._ 905 (2020) L28 [2007.02883].
  * (55) A. Gupta, S. Datta, S. Kastha, S. Borhanian, K. Arun and B. Sathyaprakash, _Multiparameter tests of general relativity using multiband gravitational-wave observations_ , _Phys. Rev. Lett._ 125 (2020) 201101 [2005.09607].
  * (56) K. Arun, B. R. Iyer, M. Qusailah and B. Sathyaprakash, _Probing the non-linear structure of general relativity with black hole binaries_ , _Phys. Rev. D_ 74 (2006) 024006 [gr-qc/0604067].
  * (57) C. Messick et al., _Analysis Framework for the Prompt Discovery of Compact Binary Mergers in Gravitational-wave Data_ , _Phys. Rev. D_ 95 (2017) 042001 [1604.04324].
  * (58) A. H. Nitz, T. Dal Canton, D. Davis and S. Reyes, _Rapid detection of gravitational waves from compact binary mergers with PyCBC Live_ , _Phys. Rev. D_ 98 (2018) 024050 [1805.11174].
  * (59) C. J. Moore, D. Gerosa and A. Klein, _Are stellar-mass black-hole binaries too quiet for LISA?_ , _Mon. Not. Roy. Astron. Soc._ 488 (2019) L94 [1905.11998].
  * (60) B. J. Owen, _Search templates for gravitational waves from inspiraling binaries: Choice of template spacing_ , _Phys. Rev. D_ 53 (1996) 6749 [gr-qc/9511032].
  * (61) B. J. Owen and B. Sathyaprakash, _Matched filtering of gravitational waves from inspiraling compact binaries: Computational cost and template placement_ , _Phys. Rev. D_ 60 (1999) 022002 [gr-qc/9808076].
  * (62) B. F. Schutz, _Networks of gravitational wave detectors and three figures of merit_ , _Class. Quant. Grav._ 28 (2011) 125023 [1102.5421].
  * (63) K. Kuns, V. Srivastava, E. Hall, M. Evans and S. Ballmer (2020), https://dcc.cosmicexplorer.org/CE-T2000007.