11institutetext:  Department of Astrophysics/IMAPP, Radboud University, P O
Box 9010, NL-6500 GL Nijmegen, The Netherlands 22institutetext: SRON,
Netherlands Institute for Space Research, Sorbonnelaan 2, NL-3584 CA Utrecht,
The Netherlands 33institutetext: Institute of Astronomy, KU Leuven,
Celestijnenlaan 200D, B-3001 Leuven, Belgium 44institutetext: Leiden
Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the
Netherlands 55institutetext: Key Laboratory of Particle Astrophysics,
Institute of High Energy Physics, Chinese Academy of Sciences, 19B Yuquan
Road, Beijing 100049, People’s Republic of China

# The Gravitational Wave Universe Toolbox:

A software package to simulate observation of the Gravitational Wave Universe
with different detectors

Shu-Xu Yi<EMAIL_ADDRESS>Gijs Nelemans 112233 Christiaan Brinkerink 11
Zuzanna Kostrzewa-Rutkowska 441122 Sjoerd T. Timmer 11 Fiorenzo Stoppa 11
Elena M. Rossi 44 Simon F. Portegies Zwart 44

###### Abstract

Context. As the importance of Gravitational Wave (GW) astrophysics increases
rapidly, astronomers interested in GW who are not experts can have the need to
get a quick idea of what GW sources can be detected by certain detectors and
the accuracy of the measured parameters.

Aims. The GW-Toolbox is a set of easy-to-use, flexible tools to simulate
observations on the GW universe with different detectors, including ground-
based interferometers (advanced LIGO, advanced VIRGO, KAGRA, Einstein
Telescope, Cosmic Explorer and also customised interferometers), space-borne
interferometers (LISA and a customised design), pulsar timing arrays mimicking
the current working ones (EPTA, PPTA, NANOGrav, IPTA) and future ones. We
include a broad range of sources such as mergers of stellar mass compact
objects, namely black holes, neutron stars and black hole-neutron stars,
supermassive black hole binaries mergers and inspirals, Galactic double white
dwarfs in ultra-compact orbit, extreme mass ratio inspirals and Stochastic GW
backgrounds.

Methods. We collect methods to simulate source populations and determine their
detectability with various detectors. The paper aims at giving a comprehensive
description of the methodology and functionality of the GW-Toolbox.

Results. The GW-Toolbox produces results that are consistent with previous
work in literature, and the tools can be accessed with a website interface
(gw-universe.org) or as a python package
(https://bitbucket.org/radboudradiolab/gwtoolbox). In the future, it will be
upgraded with more functions.

###### Key Words.:

Gravitational waves, stars:neutron, stars:black holes, stars:white dwarfs,
methods:numerical, galaxies:evolution, large-scale structure of the Universe

##  1 Introduction

Gravitational Wave (GW) astrophysics rises quickly into a pivotal branch of
astronomy. A growing number of researchers from different fields find GW
interesting and relevant to their various scientific goals (e.g. Petiteau et
al., 2013; Lasky et al., 2016; Caprini & Figueroa, 2018; Nelemans, 2018;
McWilliams et al., 2019; Barausse et al., 2020). The first opened frequency
window covers $\sim 10-1000$ Hz (a.k.a., the kHz band), which corresponds to
the working frequency range of ground-based interferometers, e.g., the
operating 2nd generation detectors advanced Laser Interferometer
Gravitational-Wave Observatory (LIGO)/Virgo interferometer (Virgo)/Kamioka
Gravitational Wave Detector (KAGRA) (Harry & LIGO Scientific Collaboration,
2010; Acernese et al., 2015; Kagra Collaboration et al., 2019) and the planned
3rd generation detectors e.g., Einstein Telescope (ET, Punturo et al. 2010)
and Cosmic Explorer (CE, Reitze et al. 2019). Space-borne interferometers,
e.g., Laser Interferometer Space Antenna (LISA, Amaro-Seoane et al. 2017) and
other similar projects, e.g., DECIGO (Kawamura et al., 2006), Taiji (Mei et
al., 2020) and Tianqin (Luo et al., 2016), will cover the GW spectrum in the
frequency range $\sim 10^{-3}-1$ Hz (a.k.a., mHz band). At even lower
frequencies, pulsar timing arrays (PTAs, Hobbs & Dai 2017, e.g.; Dahal 2020,
e.g.) are used to probe GWs with frequencies around $10^{-8}-10^{-5}$ Hz
((a.k.a., nHz band). There are also attempts to search for GW at frequencies
higher than the kHz regime, i.e., MHz/GHz (Aggarwal et al., 2020). Although
this frequency range is another important part of the GW Universe, which is
full of opportunities to discovery new physics and phenomena, we do not
include detectors and sources of this frequency range in the GW-Toolbox at
this stage, due to its larger uncertainties.

Together, these detectors are sensitive to a very broad range of GW sources,
where higher frequency detectors are sensitive to smaller objects. In the
frequency range $\sim 10-1000$ Hz, there are inspiral and mergers of stellar
mass compact object binaries, spinning neutron stars and supernovae explosions
(see e.g., Clark et al. 1979; Thorne 1987; Schutz 1989; Lipunov et al. 1997;
Portegies Zwart & McMillan 2000; Belczynski et al. 2002; Ott 2009; Goetz &
Riles 2011; Postnov & Yungelson 2014). On the other hand, the early inspiral
phase of stellar mass compact object binaries occupy the frequency range $\sim
10^{-3}-1$ Hz (Sesana, 2016). Also white dwarfs binaries, which are not
compact enough to be detected by kHz detectors, are prominent sources for mHz
detectors (e.g., Nelemans et al., 2001; Ruiter et al., 2010; Sesana, 2016;
Lamberts et al., 2018; Tauris, 2018; Sesana et al., 2020; Lau et al., 2020).
Since BH sizes increase with mass, the mHZ detectors are sensitive to mergers
of supermassive BHs (SMBHs $10^{3}-10^{8}\,M_{\odot}$ (e.g., Sesana et al.,
2005; Berti et al., 2006) and extreme mass ratio inspirals (EMRIs, see Gair et
al. 2004; Amaro-Seoane et al. 2007). Again, in the early inspiral phase, SMBH
binaries occupy the low frequency end of the spectrum $\sim 10^{-8}-10^{-5}$
Hz, either as individual sources or as a stochastic background (Sazhin, 1978;
Detweiler, 1979; Hellings & Downs, 1983; Jenet et al., 2005; Sesana et al.,
2008).

The first detection of GW was made by the LIGO/Virgo Collaboration (LVC) in
2015 (Abbott et al., 2016). The event GW150914 originates from the merger of a
binary black hole (BBH) and already had interesting astrophysical implications
(see Abbott et al., 2016; Belczynski et al., 2016). Since then, there have
been 90 GW events detected (as reported by LVC Abbott et al. 2019b, 2020c,
2021b, and there are more candidates reported by other groups from the public
strain data). These events include mergers of BBH, double neutron stars (DNS,
e.g., GW170817 and GW190425) and Black Hole-Neutron stars (BHNS, e.g.
GW190426, GW200105, GW200115 (Abbott et al., 2021a)). The population of
detected BBH mergers provides important clues on stellar formation and
evolution history (e.g., The LIGO Scientific Collaboration et al., 2020).
Among the detected sources, there are several unique systems which provide
input to and even challenge current stellar evolution theory (e.g., GW190412,
GW190814) and provide insight into the nature of neutron star matter
(GW170817, GW190425). The multi-messenger observations of the DNS merger event
GW170817/GRB 170817A/AT 2017gfo brought huge progress on astrophysics,
fundamental physics and cosmology (e.g. Abbott et al., 2017a, b, c, d, e, f;
Coulter et al., 2017; Pian et al., 2017; Abbott et al., 2018b, a, 2019a,
2019c, 2019d, 2020a). Future detectors like ET will vastly increase the
detection reach and thus allow even broader science investigations (e.g.
Maggiore et al., 2020).

Also PTA observations are routinely happening, and although there has no
conclusive evidence on GW detection with PTA111NanoGrav found evidence for a
common stochastic signal across pulsars, but there is no significant evidence
of that being a GW (Arzoumanian et al., 2020)., the communities are putting
more and more stringent upper limits on the nHz GW signals (van Haasteren et
al., 2011; Demorest et al., 2013; Lentati et al., 2015; Arzoumanian et al.,
2018; Verbiest et al., 2016; Aggarwal et al., 2019, 2020) and already ruled
out some theories of galaxy evolution (Shannon et al., 2013, 2015). Future
studies, in particular including SKA will make significant steps in the
science that can be done (e.g. Janssen et al., 2015).

In the following sections we describe the GW-Toolbox, while in section 5 we
show the results and compare these to the literature. In sections 6 and 7 we
discuss the caveats, further plans and summarize the paper.

##  2 The Gravitational Wave Universe Toolbox

The GW-Toolbox (website: gw-universe.org)222Python package repository:
bitbucket.org/radboudradiolab/gwtoolbox. is a set of tools for a broad
audience to quickly simulate the observation on the GW universe with different
detectors. The results include the expected number of detections, synthetic
catalogues and parameter uncertainties. We include three classes of GW
detectors, namely Earth-based interferometers, space-borne interferometers,
and PTAs. In each of these classes, the GW-Toolbox has the following detectors
with default and customised settings:

  * •

Earth-based interferometers:

    * –

Advanced LIGO in O3

    * –

Advanced LIGO at final design sensitivity

    * –

Advanced Virgo at final design sensitivity

    * –

KAGRA

    * –

Einstein telescope

    * –

Cosmic Explorer

    * –

A LIGO/Virgo-Like interferometer with customised parameters

  * •

Space-borne interferometers:

    * –

default LISA

    * –

LISA-like spacecraft with customised parameters

  * •

Pulsar Timing Arrays:

    * –

Existing EPTA

    * –

Existing PPTA

    * –

Existing NANOGrav

    * –

Existing IPTA

    * –

any of the existing PTAs plus simulated new pulsars to mimick future PTAs

The three classes of detectors correspond to the three main frequency regimes,
namely kHz, mHz and nHz. In each of these regimes, we include the following
sources.

  * •

kHz GW:

    * –

Binary BH (BBH) mergers

    * –

Double Neutron Star (DNS) mergers

    * –

BH-Neutron Star (BHNS) mergers

  * •

mHz GW:

    * –

Supermassive BH Binary Mergers (SMBBH)

    * –

Close Galactic White Dwarf binaries insprials (GWD)

    * –

Extreme Mass Ratio Inspirals (EMRIs, Stellar mass BHs inspiraling into
supermassive BHs)

  * •

nHz GW:

    * –

Individual SMBBH inspiral

    * –

A stochastic GW background

See figure 1 for a summary of the detectors and sources.

Figure 1: A summary of detectors and sources included in the GW-Toolbox

Figure 2: Screen shot of the GWT start page (left) and ground-based
observatories page with results for advanced LIGO (right).

In practice, the logical flow of the GW-Toolbox is to first select the
detector class and the detector parameters, then choose the source class,
select its parameters and run. For some of the sources, the underlying
cosmological model is also relevant to the simulation, where users are able to
select a certain build-in cosmology or to input parameters for a self-defined
$\Lambda$-CDM cosmological model. Examples in this paper are simulated with
the cosmology model that corresponds to Planck DR15 (Planck Collaboration et
al., 2016) We use the astropy package (Astropy Collaboration et al., 2013,
2018) with the cosmology model “Planck15”. Changing to more up-to-date
“Planck18” cosmology model does not significantly change the simulation
catalogues.

The results in general are the expected number of detections of the selected
sources, and a synthetic catalogue with/without uncertainties on parameters.
Plots and histograms of parameters of the catalogue can be made in-browser,
and the figures and catalogue can also be exported. In figure 2, a screen shot
of the GW-Toolbox website with the top-level selection is shown (left panel),
in addition to an example of the kHz detector selection page (right panel).

This paper aims at giving a comprehensive description of the methodology and
the functionality of the GW-Toolbox. The paper is organised with a similar
structure as the set-up of the GW-Toolbox. For each class of sources, a model
for the Universe is made in which the sources are distributed in space, time
and other relevant parameters (such as mass, spins etc) according to a pre-
defined population model. The user is allowed to change the population model,
according to a parametrised formalism. In section 3, we give description on
each of the population models.

The resulting GW source population is then simulated to be observed against a
selected GW detector from the above list. After source and detector are
determined, the _detected_ sources are selected using a signal-to-noise ratio
(SNR) criterion. Using the Fisher Matrix formalism, the uncertainties of
parameters of the detected sources are estimated. In section 4, we describe
how the response of the detector is represented, and the algorithm to generate
the catalogue of detected sources. For PTAs, we describe our simplified
representation of the pulsar array, the properties of the timing noises and
the observation campaign.

##  3 Implementation 1: The Universe model

We first turn to the implementation of the Universe model, in which we select
source population models from the literature and in some cases, allow users to
change or adapt the source population. There are many reviews of GW sources
that discuss these in more detail, such as Sathyaprakash & Schutz (2009);
González et al. (2013) and below we discuss the relevant ones. We start with
sources for ground-based detectors, where we concentrate on compact binary
mergers detectable with ground-based detectors, and then move to space-borne
detectors, with supermassive binary black hole mergers, white dwarf compact
binaries and extreme mass-ratio inspirals (EMRIs) as sources, and finish with
pulsar timing arrays for which we discuss individual supermassive black hole
binaries as well as stochastic background.

###  3.1 Sources for Ground Based detectors

There have been many studies of the formation and population of sources for
ground based detectors (e.g. Phinney, 1991; Schneider et al., 2001; Fryer &
Kalogera, 2001; Belczynski et al., 2002; Dominik et al., 2015; de Mink &
Belczynski, 2015) and the population depends in a complicated way on many
uncertain aspects of binary evolution, the formation of binary stars and even
the metallicity evolution of the cosmic star formation history. We take a
simple, yet flexible approach by using a parametrised description of the
source populations (but see e.g. Moe & Di Stefano 2017; Chruślińska et al.
2019, 2020 for discussions of some of the complications).

####  3.1.1 Binary Black Hole Mergers

The population of BBH mergers is the most prominent in current GW detectors
(Abbott et al., 2020c). The source population is characterised by the merger
rate as a function of redshift and the distribution of masses and spins (e.g.
Mapelli et al., 2017; Farr et al., 2011; Kovetz et al., 2017; Talbot & Thrane,
2018; Postnov & Kuranov, 2019; Abbott et al., 2019e).

In the population model for BBH, the merger rate density is expressed as:

$\dot{n}(z,m_{1},q,\chi)=\mathcal{R}(z)f(m_{1})\pi(q)P(\chi),$ (1)

where $f(m_{1})$ is the mass function of the primary (heavier) black hole,
$\pi(q)$ and $P(\chi)$ are the probability distributions of the mass ratio
$q\equiv m_{2}/m_{1}$ and the effective spin $\chi$ respectively.
$\mathcal{R}$ as function of $z$ is often refer to as the cosmic merger rate
density. We take the parameterisation as in (Vitale et al., 2019):

$\mathcal{R}(z_{m})=\mathcal{R}_{n}\int^{\infty}_{z_{m}}\psi(z_{f})P(z_{m}|z_{f})dz_{f},$
(2)

where $\psi(z)$ is the Madau-Dickinson star formation rate:

$\psi(z)=\frac{(1+z)^{\alpha}}{1+(\frac{1+z}{C})^{\beta}},$ (3)

with $\alpha=2.7$, $\beta=5.6$, $C=2.9$ (Madau & Dickinson, 2014), and
$P(z_{m}|z_{f},\tau)$ is the probability that the BBH merger at $z_{m}$ if the
binary is formed at $z_{f}$, which we refer to as the distribution of delay
time with the form (Vitale et al., 2019):

$P(z_{m}|z_{f},\tau)=\frac{1}{\tau}\exp{\left[-\frac{t_{f}(z_{f})-t_{m}(z_{m})}{\tau}\right]}\frac{dt}{dz}.$
(4)

In the above equation, $t_{f}$ and $t_{m}$ are the look back time
corresponding to $z_{f}$ and $z_{m}$ respectively.

We give plots of $\mathcal{R}(z)$ with different $\mathcal{R}_{n}$ and $\tau$
in figure 3. The default parameters are set to $\mathcal{R}_{n}=13\,$Gpc-3
yr-1 and $\tau=3$ Gyrs, which are compatible with the local merger rate of BBH
found in O3a (The LIGO Scientific Collaboration et al., 2020).

Figure 3: $\mathcal{R}(z)$, [eq. (2)] with different $\mathcal{R}_{0}$ and
$\tau$.

Although some information is known about the masses of the observed BBH (The
LIGO Scientific Collaboration et al., 2020), we assume a generic mass function
$p(m_{1})$:

$\displaystyle p(m_{1})\propto$
$\displaystyle\begin{cases}\exp\left(-\frac{c}{m_{1}-\mu}\right)(m_{1}-\mu)^{-\gamma},&m_{1}\leq
m_{\text{cut}}\\\ 0,&m_{1}>m_{\text{cut}}\\\ \end{cases}$ (5)

The distribution of $p(m_{1})$ is defined for $m_{1}>\mu$, which has a power
law tail of index $-\gamma$ and a cut-off above $m_{\rm cut}$. When
$\gamma=3/2$, the distribution becomes a shifted Lévy distribution. $p(m_{1})$
peaks at $m_{1}=c/\gamma+\mu$. We set $\mu=3$, $\gamma=2.5$, $c=6$,
$m_{cut}=95\,M_{\odot}$ as default, which result in simulated catalogue that
fits with the observed one (see section 5.1). The normalization of $p(m_{1})$
is

$c^{1-\gamma}\Gamma(\gamma-1,\frac{c}{m_{\text{cut}}-\mu}),$

where $\Gamma(a,b)$ is the upper incomplete gamma function;

In order to provide more flexibility, we also provide an alternative
$p(m_{1})$, which has an extra Gaussian peak component $p_{\rm{peak}}(m_{1})$
on top of that in equation (5):

$p_{\rm{peak}}(m_{1})=A_{\rm{peak}}\exp\left[-\left(\frac{m_{1}-m_{\rm{peak}}}{\sigma_{\rm{peak}}}\right)^{2}\right],$
(6)

the normalization of the alternative $p(m_{1})$ is:

$c^{1-\gamma}\Gamma(\gamma-1,\frac{c}{m_{\text{cut}}-\mu})+\sqrt{2\pi}\sigma_{\rm{peak}}A_{\rm{peak}}.$

We denote the population model without/with the peak component in the mass
function as BBH-PopA/B. Our default parameters for the peak component are
$A_{\rm{peak}}=0.002$, $m_{\rm{peak}}=40\,M_{\odot}$,
$\sigma_{\rm{peak}}=1\,M_{\odot}$, which are compatible with that implied from
GWTC-3. In figure 4, we plot the mass distributions for BBH-PopA/B.

Figure 4: The primary BH mass distribution in the default models of BBH-
PopA/B, see eqns (5), (6). Here we mark $\mu$, $\mu+c/\gamma$ and $m_{\rm
peak}$ in the plot, to give an intuition of these quantities.

For $\pi(q)$ we use a uniform distribution between [$q_{\rm cut}$, 1] and
assume $\chi$ follows a truncated Gaussian distribution centered at zero with
standard deviation $\sigma_{\chi}$, which is limited between -1 and 1 as
demanded by general relativity. The actual mass ratio distribution is still
poorly constrained from observation. The discovery of the GW190412 with
asymmetric masses implies that the mass ratio distribution can be quite board
(Abbott et al., 2020b). The Gaussian spin model is consistent with the finding
in The LIGO Scientific Collaboration et al. (2020). The default parameters we
use are $q_{\rm cut}=0.4$ and $\sigma_{\chi}=0.1$. Those parameters in the
population model can all be reset by users in both the web interface or in the
python package.

####  3.1.2 Double neutron star mergers

For the population model of DNS mergers, the merger rate density is similarly
expressed as:

$\dot{n}(z,m_{a},m_{b},\chi)=\mathcal{R}(z)p(m_{a})p(m_{b})P(\chi),$ (7)

where $\mathcal{R}(z)$ is taking the same form as in BBH population model, but
with a different default setting: $\mathcal{R}_{n}=300\,$Gpc-3 yr-1 and
$\tau=3$ Gyrs, which are compatible with the local merger rate of DNS found
with GWTC-2;

$p(m_{a,b})$ is the mass function of the neutron stars. Note that we use
$m_{a,b}$ instead of $m_{1,2}$. The latter are the primary and secondary stars
based on their masses, while the former we do not distinguish between the
primary and the secondary. We assume both $m_{a}$ and $m_{b}$ are following
the same mass function. We use a truncated Gaussian with upper and lower cuts
as the mass function. The default parameters are the mean
$\overline{m}=1.4\,M_{\odot}$, the mass dispersion
$\sigma_{m}=0.5\,M_{\odot}$, upper cut $m_{\rm cut,low}=1.1\,M_{\odot}$,
$m_{\rm cut,high}=2.5\,M_{\odot}$; we apply also a truncated Gaussian spin
model, with a smaller dispersion $\sigma_{\chi}=0.05$. These simplified
choices are roughly agree with observations (Valentim et al., 2011; Özel et
al., 2012; Kiziltan et al., 2013; Miller & Miller, 2015; Farrow et al., 2019;
Zhang et al., 2019), while they can be substitute with more realistic models.

####  3.1.3 Neutron star/black hole mergers

For the population model of BHNS, we again assume the merger rate density to
be:

$\dot{n}(z,m_{1},m_{2},\chi)=\mathcal{R}(z)f(m_{\bullet})p(m_{\rm n})P(\chi),$
(8)

where $\mathcal{R}(z)$ is taking the same form as for BBH and DNS, with a
different default setting: $\mathcal{R}_{n}=45\,$Gpc-3 yr-1 and $\tau=3$ Gyrs;
which is broadly consistent with the merger rate, which is implied by the
number of BHNS detection in LVK O3a+b (Abbott et al., 2021a). $f(m_{\bullet})$
is the mass function of the BH, which has the same function forms same as BBH.
We denote the population model without/with the peak component in
$f(m_{\bullet})$ as BHNS-PopA/B; The default parameters for $f(m_{\bullet})$
are $\mu=3$, $\gamma=2.5$, $c=6$, $m_{cut}=95\,M_{\odot}$,
$A_{\rm{peak}}=0.002$, $m_{\rm{peak}}=40\,M_{\odot}$,
$\sigma_{\rm{peak}}=1\,M_{\odot}$. $p(m_{n})$ is the mass function of the NS,
which is the same as in the DNS case.

###  3.2 Sources for space-borne detectors

####  3.2.1 Suppermassive Black Hole Binaries

In the last two decades, it has been established that in the center of most
galaxies there is a supermassive black hole (SMBH, with mass from
$10^{4}\,M_{\odot}$ to $10^{10}\,M_{\odot}$). Since the mergers of galaxies
are thought to be ubiquitous under the hierarchical clustering process, there
are expected to be close binaries of supermassive binary black holes (SMBBH)
in the merger galaxies, which emit GW during their inspiral and merger phase
(e.g. Colpi, 2014). Such SMBBH insprials are the main targets of LISA (Amaro-
Seoane et al., 2017), TianQin (Luo et al., 2016) and PTA (Jenet et al., 2004,
2005), since the frequency of their GW falls in the $10^{-8}-1$ Hz range. We
use the SMBBH merger catalogues from Klein et al. (2016) (Klein16 hereafter),
which are based on Barausse (2012). There are three population models being
considered, namely pop3, Q3_nodelays and Q3_delays. They mainly differ in the
origin of their SMBH seeds, and whether the delays between SMBH mergers and
galaxy mergers are included (see Klein16 for a detailed description; see also
a review on supermassive BH formation and evolution by Inayoshi et al. 2020).
For pop3, the SMBH seeds are from PopIII stars (light seeds), and the delay
between SMBH and galaxy mergers is accounted; for Q3_nodelays and Q3_delays,
the SMBH seeds are assumed to be from the collapse of protogalactic disks
(heavy seeds). The former accounts for the delay between SMBH and galaxy
mergers, while the latter does not. For each population model, there are ten
catalogues, each corresponding to a realisation of all sources in the Universe
within five years. The number of total events for each population is $\sim
890$ for pop3, $\sim 630$ for Q3_nodelays and $\sim 40$ for Q3_delays. They
are compatible with the reported averaged merger rates in Klein16 (pop3:
175.36 yr-1, Q3_nodelays: 121.8 yr-1 and Q3_delays: 8.18 yr-1). In figure 5 we
plot $M_{z,\rm tot}$ (red-shifted total mass) and $z$ of SMBBH mergers that
occur in the Universe over a timescale of five years for three population
models as a direct demonstration of the population models. The distribution
agrees with those shown in Fig 3 of Klein16.

Figure 5: $M_{z,\rm tot}$ vs. $z$ of SMBBH mergers that occur in the Universe
over a timescale of five years for three population models from Klein16. The
data corresponds to one realisation.

####  3.2.2 Double White Dwarfs

Another important population of LISA sources are Ultra Compact Galactic
Binaries. Among those Galactic binaries, close double white dwarfs (DWDs) are
the dominant, and are long expected to be promising targets for LISA and other
space-borne GW detectors (e.g. Nelemans et al., 2001; Yu & Jeffery, 2010;
Nissanke et al., 2012; Toonen et al., 2012; Korol et al., 2017; Lamberts et
al., 2019; Huang et al., 2020).

We use a synthetic catalogue of close DWD in the whole Galaxy (Nelemans et
al., 2001). In figure 6 we plot the joint distribution density of the
frequencies $f_{\rm s}=2/P$ (where $P$ is the orbital period of the binaries)
and the intrinsic amplitudes $A=2\left(G\mathcal{M}\right)^{5/3}\left(\pi
f\right)^{2/3}/(c^{4}d)$ of GW emitted from binaries in the catalogue. The
contours mark levels which correspond to uniform iso-proportions of the
density. The total number of sources in the catalogue is $\sim 2.6\times
10^{7}$.

Beside the synthetic catalogue, we also include a catalogue of 81 known DWDs
(Huang et al., 2020) a.k.a. verification binaries (VBs, see also Kupfer et al.
2018). Those VBs are plotted along in figure 6 with star markers. Note that
the distribution of the verification binaries do not follow that of the
simulated DWD in the whole Galaxy, due to very strong observational biases
which favour close by sources to be detected as VBs.

Figure 6: The distribution density of GW properties of DWD binaries in the
simulated catalogue (black contours), as function of GW frequency $f_{\rm s}$
and and GW amplitude $A$. The blue stars are known DWDs as verification
binaries.

####  3.2.3 Extreme Mass-Ratio Inspirals

In the nuclei region of galaxies, surrounding the supermassive black holes
(SMBH), there are abundant stellar populations. Among them, compact objects,
including stellar mass BHs, NS and white dwarfs, can inspiral into the central
SMBH, radiating a large amount of energy in GW. These systems are referred to
as extreme mass ratio insprials (EMRIs). EMRIs are very interesting targets
for space-borne detectors such as LISA (see Amaro-Seoane et al., 2007).

In the GW-Toolbox, we use the simulated catalogues from Babak et al. (2017)
(Babak17 hereafter) and use their population models M1-M11 (skipping M7,
explained later). These populations differ in several ingredients: the mass
function and spin distribution of SMBHs, the $M-\sigma$ relation, the ratio of
plunges to EMRIs and the characteristic mass of the compact objects (see
Babak17) for a detailed description of models). For each population model,
there are ten realizations of the catalogues, which contain detectable EMRIs
within one year with the assumption of standard LISA noise properties. The
distributions of $\mu$ (mass of the stellar BH), $M$ (mass of the massive BH)
and $D$ (luminosity distance) in the catalogues are plotted in figure 7. The
histogram for each population is averaged among the ten realizations. Since we
are using an SNR limited sample, instead of a complete one of the whole
Universe (which is $\sim 10$ times larger), the GW-Toolbox will give an
underestimated detectable number and an incomplete catalogue of detections,
especially when using a lower SNR cutoff or a more sensitive LISA
configuration. For now, we exclude M7 and M12 from the Toolbox, because in
their population models the direct plunges are ignored, therefore the total
number of EMRIs are about one order of magnitude larger than others, which
will make the computation take too long.

Figure 7: The averaged histograms of $\mu$, $M$ and $D$ of EMRIs happen in one
year, assuming different population models from Babak17

###  3.3 Sources for Pulsar Timing Arrays

####  3.3.1 Individual Massive Black Hole

Long before (millions or tens of millions of years, depends on their chirp
mass; see Sesana & Vecchio 2010; Petiteau et al. 2013; Burke-Spolaor et al.
2019) the GW from supermassive black hole binaries entering the band of space-
borne GW detectors, it lies in the PTA frequency range. If there would be such
MBH binaries sufficiently close to the Earth, PTAs could detect their signals.
Since no sources are known yet, we incorporate this source class in the GW-
Toolbox in the form of a free form in which the user can fill in the frequency
and GW amplitude, and the GW-Toolbox will determine if this binary, as a
monochromatic GW source, can be detected by the selected PTA.

####  3.3.2 Stochastic background

The second class of target for PTA are Stochastic GW background (SGWB). It can
originate from the incoherent overlapping of many unresolvable SMBBH (Phinney,
2001; Sesana et al., 2008; McWilliams et al., 2014), the relic of primordial
GW (Grishchuk, 1976, 1977; Linde, 1982; Starobinsky, 1980), or the collision
of cosmic-strings (Damour & Vilenkin, 2001, 2005; Siemens et al., 2006, 2007;
Ölmez et al., 2010). Each of these gives rise to a power-law GW signal

$h_{c}^{2}(f)=Cf^{\gamma},$ (9)

where the index $\gamma$ corresponds to the origin of SGWB. For incoherent
overlapping of SMBBH, $\gamma=-2/3$; for relic GW, $\gamma=-1$ and for cosmic-
strings, $\gamma=-7/6$. We also enable users to customize $\gamma$.

##  4 Implementation 2: Gravitational Wave detectors

###  4.1 Ground-based interferometers

####  4.1.1 Noise model of interferometers

For ground-based interferometers, the GW-Toolbox integrates the design
performance of advanced LIGO (aLIGO), Advanced Virgo (AdV), KAGRA, CE and ET
instruments. The noise models for the above-mentioned interferometers are
taken from the following resources:

  * •

aLIGO:
https://dcc.ligo.org/LIGO-T1800044/public, see also The LIGO Scientific
Collaboration et al. (2015);

  * •

adV and KAGRA:
https://dcc.ligo.org/LIGO-T1500293/public;

  * •

CE:
https://dcc.ligo.org/LIGO-P1600143/public;

  * •

ET:
http://www.et-gw.eu/index.php/etsensitivities, see also Hild et al. (2008,
2010, 2011);

In the upper panel of figure 8, we plot the noise curves that are used for the
default detectors. aLIGO-O3 corresponds to the noise performance in O3 period,
while aLIGO-design is the designed sensitivity of the aLIGO; CE1 and CE2
correspond to the expected first and second stages of CE respectively; for ET
we use the ET-D-sum curve.

Figure 8: Upper Panel: Noise Curves of the default detectors: aLIGO-O3
corresponds to the noise performance in O3 period, while aLIGO-design is the
designed sensitivity of the aLIGO; CE1 and CE2 correspond to the expected
first and second stages of CE respectively; ET-D-sum curve is applied for ET.
Lower Panel: aLIGO in design vs. Customised settings (arm length=4 km, laser
power=125 W)

In addition, the GW-Toolbox employs the package FINESSE to calculate
$S_{\rm{n}}$ of a LIGO-like or a ET-like interferometer with customised
settings (Brown et al., 2020). Users can define the following parameters of
the detector:

  * •

Arm Length

  * •

Laser power

  * •

Arm mirror mass

  * •

Arm mirror transmission coefficient

  * •

Signal recycling mirror transmission coefficient

  * •

Power recycling mirror transmission coefficient

  * •

Signal recycling phase factor

  * •

Power recycling cavity length

  * •

Signal recycling cavity length

The most important parameters that affect the sensitivity are the arm length
and laser power. In the lower panel of figure 8, we show the effects of
varying the arm length and the laser power starting from the design aLIGO
sensitivity. One of the other influential parameters is the mass of the arm
mirror. Heavier arm mirrors will decrease the noise in the low frequency ends
slope and leave the high frequency end unaffected.

####  4.1.2 SNR of GW from compact binary merger

The core of the method with which the GW-Toolbox determines detectability of
sources is to compare the SNR threshold $\rho_{\star}$ with that of a source,
which can be calculated with (Maggiore, 2008):

$\rho^{2}=4\int^{f_{\rm{high}}}_{f_{\rm{low}}}\frac{|h^{2}(f)|}{S_{\rm{n}}(f)}df,$
(10)

where $h(f)$ is the frequency domain response of the interferometer to the GW
signal, and $S_{\rm n}$ is the noises power density. For a binary system, the
detector response can be expressed as333in general, the detector response to
the ($\ell$, $m$) multipole mode is  $h^{\ell m}=F_{+}h^{\ell
m}_{+}+F_{\times}h^{\ell m}_{\times}.$ When we limit the waveform to the
dominant (quadrupole) mode, the response takes the form in equation 10, with
the terms $(1+\cos^{2}\iota)/2$, $\cos\iota$ from the $h^{22}_{+}$ and
$h^{22}_{\times}$ respectively. For expressions of higher order mode $h^{\ell
m}_{+,\times}$, see Mehta et al. (2017); Mills & Fairhurst (2021). We discuss
the influence of including the higher order modes in Section 6.:

$h(f)=\mathcal{C}\sqrt{\left(\frac{1+\cos^{2}\iota}{2}\right)^{2}F_{+}^{2}+\cos\iota^{2}F_{\times}^{2}}A(f)e^{-i(\Psi(f)+\phi_{\rm{p}})}.$
(11)

In the above equation, the constant

$\mathcal{C}=\frac{1}{2}\sqrt{\frac{5}{6}}\frac{\left(G\mathcal{M}\right)^{5/6}}{c^{3/2}\pi^{2/3}D_{\rm{L}}},$
(12)

where $\mathcal{M}$ is the red-shifted chirp mass of the binary, $D_{\rm{L}}$
is the luminosity distance, $\iota$ is the inclination angle between the
orbital angular momentum and the line of sight; $A(f)$ is the frequency
dependence of the GW amplitude. As a proof of principle of the GW-Toolbox, we
apply the waveform approximant, which has hybrid degrees of simplification.
For the amplitude $A(f)$, we use a broken power-law:

$A(f)=\begin{cases}f^{-7/6},&f\leq f_{\rm trans},\\\ f^{-2/3},&w_{m}f_{\rm
trans}<f\leq f_{\rm cut};\\\ 0,&f\geq f_{\rm cut},\\\ \end{cases}$ (13)

where $w_{m}$ is the scaling factor that to make $A(f)$ continuous. The above
formula represents the inspiral, merger and ring-down phases of a circular,
non-spinning point mass binary. It corresponds to the $A(f)$ in the
approximant IMRPhenomD (Ajith et al., 2011), when the effective spin $\chi=0$.
The transition frequency $f_{\rm trans}$ and the cut-off frequency $f_{\rm
cut}$ corresponds to $f_{1}$ and $f_{2}$ of Ajith et al. (2011). $A(f)$ and
$m_{1},m_{2}$ and $D$ via $\mathcal{C}$ essentially determine the SNR of
sources, because the effective spin has a much weaker effect on the SNR.
Therefore the dependence on $\chi$ can be separated from that on $m_{1},m_{2}$
and $D$. It makes the catalogue sampling process easier and more effective
(the Markov Chain mixes faster with less dimension, see below).

$\Psi(f)$ is the phase of the waveform. It is essential for evaluating the
uncertainties on the intrinsic parameters. Therefore we use $\Psi(f)$ from the
frequency domain waveform IMRPhenomD (Ajith et al., 2011) and restore its
dependency on $\chi$. It corresponds to a hybrid waveform of a circular binary
with parallel spin, that matches post Newtonian and numerical relativity
waveforms.

$F_{+,\times}$ are the antenna patterns of the interferometer, which are
function of position angles of the source $(\theta,\varphi)$ and the
polarization angle of the GW $\psi$, For LIGO/Virgo/KAGRA-like
interferometers, which have two perpendicular arms, the antenna patterns are:

$\displaystyle F_{+,90^{\circ}}$
$\displaystyle=\frac{1}{2}(1+\cos^{2}\theta)\cos 2\varphi\cos
2\psi+\cos\theta\sin 2\varphi\sin 2\psi$ $\displaystyle F_{\times,90^{\circ}}$
$\displaystyle=\frac{1}{2}(1+\cos^{2}\theta)\cos 2\varphi\sin
2\psi+\cos\theta\sin 2\varphi\cos 2\psi,$ (14)

and for ET-like interferometers with $60^{\circ}$ angles between the arms
(Regimbau et al., 2012):,

$\displaystyle F_{+,60^{\circ}}$
$\displaystyle=-\frac{\sqrt{3}}{2}F_{+,90^{\circ}}$ $\displaystyle
F_{\times,60^{\circ}}$
$\displaystyle=\frac{\sqrt{3}}{2}F_{\times,90^{\circ}}.$ (15)

ET will have three nested interferometers, $60^{\circ}$ rotated with respect
to each other. The antenna pattern for each interferometers are
$F_{i,+,\times}(\theta,\varphi,\psi)=F_{0,+,\times}(\theta,\varphi+2/3i\pi,\psi)$,
where $i=0,1,2$ is the index of the interferometers, and $F_{0,+,\times}$ are
those in equation (15). The joint response can be calculated with equation
(11), where the antenna pattern squared should be substituted with:

$F^{2}_{+,\times}=\sum^{2}_{i=0}F^{2}_{i,+,\times}.$ (16)

In our treatment, the modulation on the antenna patterns due to the Earth
rotation is neglected, for the duration of the GW signal is typically much
less than the period of Earth rotation. The situation would change for 3G
detectors, which has much broader frequency window and the ability to better
resolve the waveform.

####  4.1.3 Determining the sample of detected sources

Given the differential cosmic merger rate density for compact binary mergers
$\dot{n}$, the theoretical number distribution for each source class in the
catalogue is:

$N_{\rm{D}}(\mathbf{\Theta},\theta,\varphi,\psi,\iota)=\frac{\Delta
T}{1+z}\frac{dV_{\rm{c}}}{4\pi
dz}\dot{n}(\mathbf{\Theta},\theta,\varphi,\iota,\psi)\mathcal{H}(\rho^{2}-\rho^{2}_{\star}).$
(17)

where $\Delta T$ is the time span of observation, $dV_{\rm{c}}/dz$ is the
differential cosmic comoving volume (volume per redshift), $\mathcal{H}$ is
the Heaviside step function and $\rho_{\star}$ is the SNR threshold,
$\mathbf{\Theta}$ denotes the intrinsic parameters and the luminosity distance
of the source. Marginalising over the directional parameters and assuming that
$\dot{n}$ is isotropic yields

$N_{\rm{D}}(\mathbf{\Theta})=\frac{T}{1+z}\frac{dV_{\rm{c}}}{dz}\dot{n}(\mathbf{\Theta})\mathcal{D}(\mathbf{\Theta}),$
(18)

where

$\mathcal{D}(\mathbf{\Theta})=\oiint d\Omega
d\Omega^{\prime}\mathcal{H}(\rho^{2}-\rho^{2}_{\star})/(4\pi)^{2},$ (19)

is the detectability of the source, which is determined by the detector
properties and the waveform of the source. Since we use the same waveform for
BBH, DNS and BH-NS, the difference among these three populations are only in
the cosmic merger rate $\dot{n}$ discussed above in Sections 3.1.1-3.1.3.

The total number of expected events catalogue is:

$N_{\rm{tot}}=\int d\mathbf{\Theta}N_{D}(\mathbf{\Theta}),$ (20)

and the number of detections thus is Poisson realisation of the expectation
value $N_{\rm{D}}(\mathbf{\Theta})$. The synthetic catalogue is then obtained
by a Markov Chain-Monte Carlo sampling from $N_{\rm{D}}(\mathbf{\Theta})$. We
apply the elliptical slice sampling algorithm (Murray et al., 2010), which
converges faster to the target distribution comparing with the traditional
Metropolis-Hasting algorithm (Neal et al., 1999) and requires less tuning on
the initial parameters.

We also give the estimated uncertainties on the parameters using the Fisher
Information Matrix (FIM): the covariance matrix is related to the Fisher
matrix with:

$\left<\delta\mathbf{\Theta}_{i}\delta\mathbf{\Theta}_{j}\right>=\mathcal{F}^{-1}_{ij},$
(21)

where the Fisher matrix is defined as:

$\mathcal{F}_{ij}=(\partial h/\partial\mathbf{\Theta}_{i}|\partial
h/\partial\mathbf{\Theta}_{j}).$ (22)

The partial derivatives in the above equation are calculated numerically:

$\frac{\partial
h}{\partial\mathbf{\Theta}_{i}}=\frac{h(\mathbf{\Theta}_{i}+\Delta\mathbf{\Theta}_{i})-h(\mathbf{\Theta}_{i}-\Delta\mathbf{\Theta}_{i})}{2\Delta\mathbf{\Theta}_{i}}.$
(23)

In the GW-Toolbox, we use $\Delta\Theta_{i}=10^{-8}\Theta_{i}$, as it is small
enough to give stable results. We only calculate the FIM for intrinsic
parameters ($m_{1,z}$, $m_{2,z}$ and $\chi$) and an overall scaling factor
that represents the effect of all extrinsic parameters. As a result, we cannot
give an estimate of the uncertainties of extrinsic parameters, i.e., the
distance, sky locations, inclination and polarization angle, for these
parameters contribute to the overall scaling factor in a highly correlated
way. In reality, the sky locations are determined largely by triangulation
with a network of detectors, and the precision of other extrinsic parameters
are also depended on triangulation. In the current version of the GW-Toolbox
we do not include a detector network. The calculation of FIM is thus
simplified. However, in order to obtained an estimation on the uncertainties
on the source frame masses, one still needs to propagate the uncertainties on
the redshift. The uncertainties on the redshift itself is also an interesting
quantity that users want to obtain from the synthetic catalogue of GW-Toolbox.
We give evaluation on $\delta z$ from empirical relations. For 2G detectors,
we use $\delta z=0.5z$ roughly represents the trend of $\delta z-z$ in GWTC-3;
for 3G detectors, we use $\delta z=0.017z+0.012$, which is a fit from the
simulation results of Vitale & Evans (2017). $\delta z$ is propagated to the
uncertainties on the source frame masses with the following equation:

$\delta m_{i}=m_{i}\sqrt{\big{(}\frac{\delta
m_{i,z}}{m_{i,z}}\big{)}^{2}+\big{(}\frac{\delta z}{1+z}\big{)}^{2}}.$ (24)

It is worth mentioning that although FIM is a quick, simple and widely applied
method of uncertainties, it something gives overestimated uncertainties
comparing with those from full Bayesian inference (e.g. Veitch et al., 2015),
especially for events with low SNR (e.g. Vallisneri, 2008; Rodriguez et al.,
2013).

The returned catalogue is composed by the mean values of the parameters and
their estimated uncertainties. A more realistic simulation of the observed
catalogue can be obtained by shifting the mean values according to the
corresponding uncertainties, which is straightforward and easy. Since the
uncertainties given here are conservative, the GW-Toolbox return the un-
shifted catalogue, and let the user to decide whether or how to further shift
the catalogue.

###  4.2 Space-borne interferometers

The Space-borne interferometers module of the Toolbox enables users to
simulate observations with LISA-like space-borne GW observatories (see Barke
et al., 2015). We work with the codes of the LISA Data Challenge (LDC,
https://lisa-ldc.lal.in2p3.fr, a successor program of the earlier Mock LISA
Data Challenge (Babak et al., 2010)), and make it possible for users to
customize the arm length, the laser power and the telescope diameter of LISA.
In the ground-based interferometers section the theoretical probability
distribution of parameters of the detectable sources are first calculated,
afterwards samples are drawn from such distribution as synthetic catalogues of
observations. The procedure for LISA-like detectors is different: we go
through pre-generated synthetic catalogues of different source populations in
the whole Universe/Galaxy and calculate the SNR of each source to be detected
by LISA. The SNR is still calculated with:

$\rho^{2}=4\int\frac{|h^{2}(f)|}{S_{\rm{n}}(f)}df,$ (25)

where $h(f)$ is the LISA response to a waveform of a source, and $S_{n}(f)$ is
the noise power spectrum density (PSD). The time-delay interference (TDI)
channels are combinations of data streams such that the noises arising from
the fluctuation of the laser frequency can be exactly cancelled while the
signal in GW can be preserved (Tinto & Armstrong, 1999; Armstrong et al.,
1999; Estabrook et al., 2000). In the GW-Toolbox, we consider the LISA
responses and the noise spectral density in the first generation TDI-$X$
channel (Armstrong et al., 1999). The noise spectrum will be introduced in the
next section. Three classes of sources are included in the GW-Toolbox for
LISA, namely: mergers of SMBBH, resolved DWDs in the Galaxy and EMRIs.
Waveforms and the corresponding LISA responses will be introduced in the
following subsections. Examples of synthetic observations are given and
compared with the literature in Section 5. Uncertainties of the parameters are
again estimated with the FIM method.

####  4.2.1 Noise TDI

The PSD of the noise TDI-$X$ response is formulated as (Armstrong et al.,
1999):

$S_{X}(f)=[4\sin^{2}(2x)+32\sin^{2}x]S^{\rm{accel}}_{y}+16\sin^{2}xS^{\rm{optical}}_{y},$
(26)

where $x=2\pi fL/c$, and $L$ is the arm length of LISA and $c$ is the speed of
light, $S^{\rm{accel}}_{y}$ and $S^{\rm{optical}}_{y}$ are the fractional
frequency fluctuations due to the acceleration noise of the spacecrafts and
the optical meteorology system noise respectively. For the acceleration noise,
we use (Amaro-Seoane et al., 2017):

$S^{\rm{accel}}_{a}=9\times
10^{-30}\frac{[\text{m}\,\text{s}^{-2}]^{2}}{[\rm{Hz}]}\left(1+\left(\frac{[0.4\text{mHz}]}{f}\right)^{2}\right)\left(1+\left(\frac{f}{[8\,\text{mHz}]}\right)^{4}\right).$
(27)

Note that the above noise is in the form of acceleration. To convert it into a
fractional frequency fluctuation, one needs to divide it by a factor
$4\pi^{2}f^{2}c^{2}$ (Armstrong et al., 1999), resulting in

$\displaystyle S^{\rm{accel}}_{y}$ (28) $\displaystyle=\frac{3.9\times
10^{-44}}{[\text{Hz}]}\left(1+\left(\frac{[0.4\text{mHz}]}{f}\right)^{2}\right)\left(\left(\frac{[8\,\text{mHz}]}{f}\right)^{2}+\left(\frac{f}{[8\,\text{mHz}]}\right)^{2}\right).$
(29)

The noise of the optical metrology system can be split

$S^{\rm{optical}}_{y}=S^{\rm{ops}}+S^{\rm{opo}},$ (30)

where $S^{\rm{ops}}$ is the laser shot noise, which scales with the arm length
$L$, the laser power $P$ and the diameter of the telescope $D$ as (Amaro-
Seoane et al., 2017):

$S^{\rm{ops}}=5.3\times
10^{-38}\times\left(\frac{f}{[\rm{Hz}]}\right)^{2}\frac{[2W]}{P}\left(\frac{L}{2.5[\rm{Gm}]}\right)^{2}\left(\frac{[0.3\rm{m}]}{D}\right)^{2}\,\rm{Hz}^{-1};$
(31)

and

$S^{\rm{opo}}=2.81\times 10^{-38}(f/[\rm{Hz}])^{2}\,\rm{Hz}^{-1}$ (32)

denotes the contribution from other noise in the optical meteorology system.

We also include the TDI-$X$ noise PSD originating from the foreground GW
emission from Galactic DWDs (GWD). In practice this confusion noise will be
modulated with the orbital phase of the spacecraft. For simplicity, we adopt
an analytic approximation for the averaged equal-arm Michelson PSD of GWD as
(Robson et al., 2019):

$S_{\rm GWD}(f,T_{\rm{obs}})=Af^{-7/3}e^{-f^{\alpha}+\beta f\sin(\kappa
f)}\left(1+\tanh(\gamma(f_{k}-f))\right),$ (33)

Note that the noise depends on the observation duration, because for longer
observations, more and more individual DWDs can be resolved and removed from
the confusion noise foreground. This time dependent is represented by using
different parameters with different $T_{\rm{obs}}$ as shown in table 1.

|  $\leq$6 mo |  $\leq$1 yr |  $\leq$2 yr |  $\leq 4$ yr
---|---|---|---|---
$\alpha$ | 0.133 | 0.171 | 0.165 | 0.138
$\beta$ | 243 | 292 | 299 | -221
$\kappa$ | 482 | 1020 | 611 | 521
$\gamma$ | 917 | 1680 | 1340 | 1680
$f_{k}$ | 0.00258 | 0.00215 | 0.00173 | 0.00113
Table 1: The parameters of the confusion noise of the unresolved GWD
background.

The amplitude $A$ is fixed at $9\times 10^{-45}$ for $T_{\rm{obs}}\leq 4$
years, and is set to zero for larger $T_{\rm{obs}}$ since the foreground then
disappears.

The equal-arm Michelson response (fractional displacement) PSD term $S_{\rm
GWD}$ can be converted to the fractional frequency by timing a factor $x^{2}$
and then added to $S^{\rm{ops}}$ calculated with equation (26). The upper
panel of figure 9 shows the square root of the noise PSD with various LISA
parameters.

Figure 9: Upper Panel: Square root of the PSD noise TDI-$X$ corresponding to
different LISA configurations. The solid curves are the total PSD, while the
dashed curves are the contribution from the confusion GWDs. The bundle of
curves in the same colour correspond to $T_{\rm{obs}}=1,2,4,5$ years from top
to bottom; Lower Panel: Sensitivity curves corresponding to different arm
length and $T_{\rm{obs}}$. The solid curves are the total curve, while the
dashed curves are the contribution from the confusion GWDs. The bundle of
curves in the same colour correspond to $T_{\rm{obs}}=1,2,4,5$ years from top
to bottom.

Note that our $S_{X}$ should not be confused with the PSD in the Michelson
response. The latter is more commonly applied and sometimes referred as the
sensitivity curve. The GW-Toolbox also provide the latter with the following
analytic model (Robson et al., 2019):

$S_{\rm{n}}=\frac{10}{3L^{2}}\left(S^{\rm{op}}_{\rm{dis}}+2(1+\cos^{2}(x))\frac{S^{\rm{acc}}_{\rm{a}}}{(2\pi
f)^{4}}\right)\left(1+\frac{3}{5}x^{2}\right)+S_{\rm{GWD}},$ (34)

where $S^{\rm{op}}_{\rm{dis}}$ is the noise in the optical system in term of
the displacement, which can be converted into the previous Doppler
$S^{\rm{optical}}_{y}$ by multiplication by a factor $2\pi f/c$. We plot the
sensitivity curves corresponding to different arm length and $T_{\rm{obs}}$ in
figure 9. In the appendix, we give a summary plot of the conversion among
different detector responses.

####  4.2.2 TDI response to the SMBBH waveform

The TDI-$X$ response of LISA due to a GW from a SMBBH merger is calculated
using the LDC code (Babak et al., 2010), where the IMRPhenomD waveform is
adopted (Ajith et al., 2011). Figure 10 shows the modulus of the TDI-$X$
responses in the frequency domain, for three different sources. The parameters
of the example sources are listed in Table 2. The low-frequency limit
corresponds to the time to coalescence at the beginning of the observation,
and the dips at the high frequency end are due to the term $\sin(x)$ when
converting to TDI. The sample cadence is fixed to 5 s, which corresponds to a
high frequency cut-off at $0.1$ Hz. For systems with heavy BHs, like #2 in the
example (masses $>10^{5}M_{\odot}$), the frequency at coalescence is lower
than the cut-off frequency, therefore the current cadence will not lose any
power from the signal; On the other hand, for systems with light BHs, like #1
in the example (masses $<50,000M_{\odot}$), the high frequency part ($>0.1$
Hz) of the waveform will be lost. However, the decrease of SNR is less than 1%
comparing to that using a cadence of 1 s. Therefore it is acceptable to fix
the cadence to 5 s for all sources, in order to perform the simulation fast.

Figure 10: Modulus of frequency domain TDI-$X$ responses to GW from different SMBBHs (whose parameters are listed in Table 2). The solid curves correspond to LISA with 5 Gm laser arms, and dashed curves correspond to a 2.5 Gm arms configuration. # |  $\beta$ (rad) |  $\lambda$ (rad) |  $\theta_{\chi_{1}}$ (rad) |  $\theta_{\chi_{2}}$ (rad) |  $\varphi_{\chi_{1}}$ (rad) |  $\varphi_{\chi_{2}}$ (rad) | $\chi_{1}$ | $\chi_{2}$ | $m_{1}\,(M_{\odot})$ | $m_{2}\,(M_{\odot})$ |  $\theta_{L}$ (rad) |  $\varphi_{L}$ (rad) | $z$ |  $t_{\rm c}$ (yr)
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
1 | -1.3 | 0.44 | 0.8 | 2.6 | 4.5 | 5.98 | 0.8 | 0.2 | 37695 | 4582 | 2.18 | 1.3 | 0.069 | 0.76
2 | -0.44 | 4.7 | 0.08 | 0.037 | 4.26 | 5.48 | 0.04 | 0.2 | 420555 | 298237 | 1.14 | 3.26 | 5.1 | 0.0038
3 | -0.01 | 2.7 | 0.12 | 0.12 | 0.72 | 6.08 | 0.6 | 0.13 | 76476 | 28854 | 1.673 | 1.47 | 2.8 | 0.5
Table 2: Parameters of example sources corresponding the figure 10. The
meaning of the parameters are: $\beta$-Ecliptic Latitude; $\lambda$-Ecliptic
Longitude; $\theta_{{\chi_{1}}/{\chi_{1}}}$-Polar angle of spin 1/2;
$\chi_{1,2}$-Spin 1/2; $m_{1,2}$-(Intrinsic) masses of primary/secondary BH;
$\theta_{L}$-Initial polar Angle of the orbital angular momentum;
$\varphi_{L}$ Initial azimuthal Angle of the orbital angular momentum; $z$
red-shift; $t_{\rm{c}}$ time to coalescence at the beginning of observation;

####  4.2.3 TDI Waveform of DWD

We first derive the frequency domain TDI-$X$ waveform in a monochromatic plane
wave approximation. The equal-arm Michelson response of a plane GW in the
long-wavelength region can be approximated as a sine wave with a orbit
averaged amplitude $\langle\mathcal{A}\rangle$ and frequency $f_{\rm s}$. The
relation between $\langle\mathcal{A}\rangle$ and the intrinsic amplitude $A$
of the source binary can be found in equations (A12,A13) of Korol et al.
(2017).

The Fourier transform of such a signal with the duration $T_{\rm obs}$ is:

$\tilde{h}_{\rm Mich}(f)=\frac{1}{2}\langle\mathcal{A}\rangle T_{\rm obs}{\rm
sinc}((f-f_{\rm s})T_{\rm obs}).$ (35)

Note that here we use the convention that ${\rm sinc}(x)=\sin(\pi x)/(\pi x)$,
such that the integration of $|\tilde{h}(f)|^{2}$ equals $T_{\rm obs}A^{2}$.

To convert the equal-arm Michelson into TDI-$X$, we multiply by a factor
$4x\sin x$, where $x=f(2\pi L/c)$.

$X(f)=2x\sin x\langle\mathcal{A}\rangle T_{\rm obs}{\rm sinc}\left[(f-f_{\rm
s})T_{\rm obs}\right]$ (36)

In figure 11 we show the waveforms of a DWD with $A=10^{-20}$, $f_{\rm
s}=10^{-3}$ Hz, calculated analytically from equation (36) and compare them
with those calculated numerically with the LDC code (which is based on Cornish
& Littenberg 2007). Although simplified, the overall agreement between the two
is good.

From equations (25,36), we obtain an approximated squared SNR expression:

$\rho^{2}_{\rm approx.}=\frac{16x^{2}\sin^{2}x<\mathcal{A}^{2}>T_{\rm
obs}}{S_{X}(f_{\rm s})}.$ (37)

When we replace the TDI-$X$ noise PSD with the Michelson noise PSD, and drop
the $16x^{2}\sin^{2}x$ term, the above equation becomes equation (10) of Korol
et al. (2017).

Figure 11: Frequency domain LISA responses to GW from a DWD: The blue and
oranges lines correspond to two different LISA configurations; solid and
dashed lines correspond to responses calculated with numerical and analytical
methods respectively.

####  4.2.4 TDI Waveform for EMRIs

The analytic kludge (AK) waveforms (Barack & Cutler, 2004) for EMRIs are
applied and the corresponding TDI-$X$ responses are calculated with the code
package EMRI_Kludge_Suite444https://github.com/alvincjk/EMRI_Kludge_Suite
(Chua & Gair, 2015; Chua et al., 2017). As examples, in figure 12 we plot the
frequency domain TDI-$X$ responses to the AK waveforms, which correspond to
three EMRIs systems and two $L$ designs of LISA. The first system (blue
curves) has a supermassive BH with mass $M=10^{6}\,M_{\odot}$ and stellar mass
BH with mass $m=20\,M_{\odot}$. Here the masses are all measured in the
observer’s frame, i.e., red-shifted. The frequency domain response corresponds
to a time-domain waveform simulated from the semi-latus rectum $p=8GM/c^{2}$
to the final plunge. The time resolution is $dt=25$ s, which is set to ensure
that the highest frequency cut-off set by $1/(2dt)$ is larger than the Kepler
frequency around the innermost stable circular orbit (ISCO) of the
supermassive BH. The lower frequency cut corresponds to initial orbital
inspiral, and the higher frequency cut corresponds to the orbital frequency at
the plunge, which approximates to the Kepler frequency at ISCO; The second
system (orange curves) has a supermassive BH with mass $M=10^{5}\,M_{\odot}$
and stellar mass BH with mass $m=20\,M_{\odot}$. The initial semi-latus rectum
is also $p=8GM/c^{2}$. Since the Kepler frequency at ISCO is inversely
proportional to the mass of the supermassive BH, we use $dt=2.5$ s. The third
system (green curves) is identical with the second one. The difference is the
initial semi-latus rectum of the third one is $p=20GM/c^{2}$. To track the
evolution from this larger initial semi-latus rectum to the final plunge, the
simulation includes $\sim 10$ times longer time steps. As a result, more low-
frequency components are included in the third waveforms than the second.
Other physical parameters are identical for the three system and are (values
used in parenthesis) (Barack & Cutler, 2004):

  * •

$s$: dimension-less spin of the massive BH ($s=0.5$);

  * •

$e$: the initial eccentricity ($e=0$);

  * •

$\iota$: the initial inclination ($\iota=0.524$ rad);

  * •

$\gamma$: the pericenter angle in AK Waveform ($\gamma=0$);

  * •

$\psi$: the initial phase ($\psi=0.785$ rad);

  * •

$\theta_{\rm{S}}$: the sky position polar angle of source in an ecliptic-based
coordinate system, equals to $\pi/2$ minus the ecliptic latitude
($\theta_{\rm{S}}=0.785$ rad);

  * •

$\phi_{\rm{S}}$: the sky position azimuth angle of source in an ecliptic-based
coordinate system, equals to the ecliptic longitude ($\phi_{\rm{S}}=0.785$
rad);

  * •

$\theta_{\rm{K}}$: the polar angle of the massive BH spin
($\theta_{\rm{K}}=1.05$ rad);

  * •

$\phi_{\rm{K}}$: the azimuth angle of the massive BH spin
($\phi_{\rm{K}}=1.05$ rad);

  * •

$\alpha$: the azimuthal direction of the orbital angular momentum
($\alpha=0$);

  * •

$D$: luminosity distance ($D=1$ Gpc).

As mentioned above, if the frequency domain waveform were to be simulated in
real time for every candidate event, it is difficult to reconcile both the
speed of simulation and to include the full GW signal from the beginning of
the observation to the plunge. As a solution, we generate the frequency domain
TDI waveform corresponding to each candidate event in the catalogue in
advance, and store their modulus in files. The pre-generated TDI waveform
corresponds to signal from the beginning of the observation to the final
plunge. The initial semi-latus is calculated with a Newtonian formula equation
4.136 of Maggiore (2008) according to its masses, eccentricity and time to
plunge at the beginning of observation. The pre-calculated TDI corresponds to
a LISA arm length 2.5 Gm. The conversion to a different LISA arm length can be
done by rescaling with $x_{1,i}\sin x_{1,i}/(x_{0,i}\sin x_{0,i})$, where:

$x_{0,i}=2\pi f_{i}L_{\rm{default}}/c,$ (38)

and

$x_{1,i}=2\pi f_{i}L_{\rm{new}}/c.$ (39)

Figure 12: Frequency domain TDI-$X$ responses to EMRI AK waveform, which
correspond to three EMRIs systems and two $L$ designs of LISA

###  4.3 Pulsar Timing Arrays

Pulsars are rotating neutron stars. Some of the known pulsars which are very
stable, i.e., their spin period only changes a tiny fraction in a very long
epoch. Therefore, the arrival time of each pulse from such a pulsar can be
modeled with high precision. The passing of a series of GW will cause
additional changes to the expected time of arrival of the pulses (TOAs), and
thus provide a way to detect GW, at frequencies between $10^{-8}$ and
$10^{-5}$ Hz, where the lower frequency limit corresponds to the observation
span of years, and the high frequency limit corresponds to the average cadence
of a couple of days (see Sazhin, 1978; Detweiler, 1979; Hellings & Downs,
1983; Jenet et al., 2005). A pulsar timing array is a group of pulsars, which
are stable and have been monitored with radio telescopes for a long time. The
existing PTA consortia are EPTA (Kramer & Champion, 2013), PPTA (Hobbs, 2013),
NANOGrav (McLaughlin, 2013) and IPTA (Hobbs et al., 2010). The standard
procedure of pulsar timing is to first fit a timing model to the TOAs of
individual pulsars, that includes inaccuracies in the pulsars’ astrometric
parameters, a model for spin evolution, refractive effects of interstellar
medium and solar wind, the orbital and spin motion of the Earth, delays due to
general relativity etc. (see Hobbs et al., 2006). The difference between the
timing model and the observed TOAs, the timing residuals, are used to extract
information about possible GW signals with frequentist (Jenet et al., 2005;
Babak & Sesana, 2012; Ellis et al., 2012) or Bayesian methods (van Haasteren
et al., 2009; Ellis, 2013). Here we want to use a simplified way to represent
the properties of PTAs, without the need to make use of the full time series
of the timing residuals, and obtain results which agree to an order of
magnitude with the published results. We base our method on measuring the
excess power from GW over analytic timing noise power spectra. Such a practice
was also used by some early work (Jenet et al., 2004; Yardley et al., 2010; Yi
et al., 2014).

####  4.3.1 Representing the timing noises

Suppose that we have already removed every known effect that contributed to
differences in the TOAs, the timing residuals that remain are purely intrinsic
to the pulsars due to their spin irregularity. Previous studies found that
such timing noise can be decomposed into a red noise component and a white
noise component (Hobbs et al., 2010). The red noise component can be
represented with a power-law spectrum, with increasing power towards the lower
frequencies, while the white noise component has a frequency independent power
level. In the GW-Toolbox, we use the following equation to represent the noise
spectrum density of the timing residuals of an individual pulsar:

$S_{total}(f)=\sigma^{2}_{\rm w}/(f_{\rm high}-f_{\rm low})+S_{\rm n,red}(f),$
(40)

where $\sigma_{\rm w}$ is the level of the white noise, $f_{\rm high}=N/(2T)$
is the high frequency cut-off defined by the observation cadence and $f_{\rm
low}=1/T$ is the low frequency cut-off defined by the inverse of the duration,
$S_{\rm n,red}(f)$ is the red noise component, which has a power-law form:

$S_{\rm
n,red}(f)=\frac{A^{2}_{\rm{red}}}{12\pi^{2}}\left(\frac{f}{\text{yr}^{-1}}\right)^{-\alpha}.$
(41)

Therefore, we define the noise spectrum of a pulsar with five parameters,
namely: $N$ the number of observations, $T$ the duration of observation,
$A_{\rm red}$ the normalization of the red noise, $\alpha$ the power index of
the red noise, and $\sigma_{\rm w}$ the level of the white noise. The last
three are intrinsic properties of the pulsar. These parameters for the pulsars
in the above mentioned PTAs are fitted and published (Desvignes et al., 2016;
Porayko et al., 2019; Alam et al., 2020). The GW-Toolbox includes 42 pulsars
in EPTA, 26 pulsars in PPTA, 47 pulsars in NANOGrav and 87 pulsars in IPTA.
Besides the pulsars in the current PTAs, the GW-Toolbox also includes
simulated future observations, with customised observation cadence and
duration, and an increasing number of newly discovered pulsars during the
observation period. The parameters (sky coordinates RA, DEC and noise
parameters $A_{\rm{red}}$, $\alpha$, $\sigma_{\rm{w}}$) of the simulated
pulsars are assigned in the following way: randomly select two pulsars from
the current PTA with replacement, and draw a uniformly random number between
the parameters of the selected pair of pulsars, and assign the random variable
as the corresponding parameter of the new pulsar. In this method, the noise
properties and sky distribution of the new pulsars reflect those of the known
pulsars.

In figure 14, we plot the noise spectra density of pulsars in the PTAs used by
the GW-Toolbox. Blue curves correspond to known pulsars in existing PTAs, and
orange curves are simulated new pulsars. In figure 13, we plot the sky
coordinates of pulsars in the PTA. The blue dots indicate the pulsars from
IPTA, and orange dots indicate simulated new pulsars. The size of the markers
is proportional to the number of TOAs of the corresponding pulsar.

Figure 13: Galactic distribution of the PTAs. The blue dots indicate the
pulsars from current IPTA, and orange dots indicate simulated new pulsars. The
size of the markers is proportional to the number of TOAs of the corresponding
pulsar. Figure 14: The noise spectra of the pulsars. Blue curves correspond to
known pulsars in existing PTAs, and orange curves are simulated new pulsars.

Figure 15: aLIGO-design-1 year: The simulated catalogue from one year
observation by aLIGO with designed noise performance on BBH, marked with blue
crosses; aLIGO-O3-1 year: The simulated catalogue from one year observation by
aLIGO with O3 noise performance on BBH, marked with orange dots; GWTC-3 (BBH):
BBH events in GWTC-3, marked with green star symbols.

###  4.4 PTA detections

The GW from an individual supermassive BH can be approximated with a
monochromatic wave. The timing residuals induced in the $i$-th pulsar is (Yi
et al., 2014):

$\mathcal{A}_{i}=\frac{h_{\rm{s}}}{\omega}\left(1+\cos\theta\right)\sqrt{\cos^{2}2\psi\left(\frac{1+\cos^{2}\iota}{2}\right)^{2}+\sin^{2}2\psi\cos^{2}\iota}$
(42)

where $\theta$, $\iota$, $\Psi$ are the angle between the direction towards
the pulsar and the GW source, the inclination of the source binary plan and
the polarization angle respectively. The SNR of the GW in the $i$-th pulsar
is:

$\rho^{2}_{i}=\mathcal{A}_{i}^{2}/S_{n,i}(f)\times T_{i},$ (43)

where $S_{n,i}(f)$ is the noise spectrum density of the pulsar, and $T_{i}$ is
the observation duration. The total SNR squared of a PTA is:

$\rho^{2}=\sum\rho_{i}^{2}.$ (44)

The effects of a stochastic GW background (SGWB) to the timing residual are an
additional red noise:

$h_{c}^{2}(f)=Cf^{\gamma},$ (45)

where $h_{\rm c}$ is the characteristic GW strain at the frequency 1 yr-1, and
the index $\gamma$ depends on the origin of the SGWB. For incoherent
overlapping of MBH, $\gamma=-2/3$; for relic GW, $\gamma=-1$ and for cosmic
strings, $\gamma=-7/6$. Besides the additional red noise, it is also expected
that the timing residuals due to the SGWB are correlated between pairs of
pulsars. The correlation as a function of the angular separation between the
pair is:

$\Gamma_{0}=3\Big{\\{}\frac{1}{3}+\frac{1-\cos\xi}{2}[\ln(\frac{1-\cos\xi}{2}-\frac{1}{6}]\Big{\\}},$
(46)

which is referred as the Hellings and Downs Curve (Hellings & Downs, 1983).
The SNR squared in a pair of pulsars is (Anholm et al., 2009):

$\rho_{1,2}=\frac{H_{0}^{2}}{4\pi^{2}}\sqrt{2T\int^{\infty}_{0}df\frac{\Omega^{2}_{\rm{gw}}(f)\Gamma^{2}_{0}}{f^{6}P_{1}(f)P_{2}(f)}},$
(47)

where $H_{0}$ is the Hubble constant, $\Omega_{\rm gw}(f)$ is the energy
density of the SGWB relative to the critical density of the Universe. The
relation between $h_{c}(f)$ and $\Omega_{\rm gw}(f)$ is:

$h^{2}_{c}(f)=\frac{3H_{0}^{2}}{2\pi^{2}}\frac{1}{f^{2}}\Omega_{\rm{gw}}(f).$
(48)

In equation (47) $P_{1,2}(f)=S_{\rm n1,2}(f)f^{2}$. The total SNR squared is
the summation of SNR squared over all pairs in the PTA.

##  5 Results and examples

In order to test and validate the GW-Toolbox, we discuss the outcome of the
simulations for the different source populations and different detectors, and
compare these with the literature where possible.

###  5.1 Examples for stellar mass BBH detected with Earth-based detectors

For BBHs, we generate catalogues for a one year observation of aLIGO, both
with the noise spectrum of O3 (aLIGO-O3) and that of the design (aLIGO-
design). The underlying population model is BBH-PopB with the default
parameters (delay time of 3 Gyr and a mass function with an extra Gaussian
peak at 40$M_{\odot}$, see section 3.1.1). The simulated number of detection
of aLIGO-O3 is 99 (with a Poisson expectation 87.4). It is compatible with the
rate of detection in O3a ($\sim 36$ BBH in six month of observation with duty
cycle $\sim 80\%$, see Abbott et al. (2020c)); the simulated number of
detection for the aLIGO-design is 379 (expectation 376.0). In the panels of
figure 15, we plot the simulated catalogues in the $z-m_{1}$ and $m_{1}-m_{2}$
planes. To compare, we also plot the BBH events in GWTC-3(Abbott et al.,
2021b), which is comprised by GWTC-1 (Abbott et al., 2019b) and events
detected in O3a+b. The simulated sets agree well with the observed ones,
except that the uncertainties of parameters are in general larger in our
simulation comparing with those in the real observation in GWTC-3 (up to $\sim
50\%$), especially in the low SNR region. The origin of this overestimation is
double folded: 1, the intrinsic difference between FIM and Bayesian methods,
as mentioned above; 2, the localisation information from triangulation is not
included in our method as were in real catalogue. Therefore, the estimated
uncertainties can serve as a conservative upper limits.

We also simulate a catalogue of one month observation with ET. In this case,
the number of detection is 1923 (expectation $1.9\times 10^{3}$). In the
panels of figure 16, we plot the distribution of the catalogue in $z$ and
$m_{1}$. We also plot the distribution of BBH mergers in the whole Universe as
defined by the population model for comparison. From this example, it is clear
that ET will probe the distribution of sources throughout the Universe very
well, as was shown before (e.g. Vitale & Evans, 2017). We will study this in
more detail in a forthcoming paper (Yi et al. in preparation)

Figure 16: $m_{1}$ distribution (upper panel) and the number density as
function of redshift (lower panel) in the catalogue for one month of
observation of BBH mergers by ET (solid blue). The integral of the latter is
the total number of detection. As a comparison, we plot that of the whole
Universe within one month in dashed orange.

Figure 17: The same as figure 15, but for ten years observation of DNS with
eLIGO (design and O3 sensitivity).

Figure 18: The total mass distribution (upper panel) and the number density as
function of red-shift (right panel) in the catalogue of one year observation
of DNS mergers by ET (solid blue). The integral of the latter is the total
number of detection. The dashed orange curve is that in the whole universe
within one year.

Figure 19: aLIGO-design-10 year: The simulated catalogue from ten years’
observation by aLIGO with designed noise performance on BHNS, marked with blue
crosses; aLIGO-O3-10 year: The simulated catalogue from ten year observation
by aLIGO with O3 noise performance on BHNS, marked with orange dots.

###  5.2 Examples for double neutron stars detected with Earth-based detectors

The number of DNS mergers detected so far is small (2-3). Therefore, we
generate catalogues for ten years of detection of DNS mergers by aLIGO, both
with the noise spectrum in O3 (aLIGO-O3) and that at design (aLIGO-design).
The number of detection of aLIGO-O3 is 56 (expectation 48.7) in ten years,
which is in accordance with the real detection rate in O3a (1-2 in six month);
the number for aLIGO-design is 238 (expectation 236.3). In the panels of
figure 17, we plot the simulated catalogues in the $z-m_{1}$ and $m_{1}-m_{2}$
planes. We also plot the DNS events in GWTC-3 (GW170817, GW190425) to compare.
It is difficult to draw strong conclusions with so few detection, but broadly
the results of the GW-Toolbox agree with the observations so far. To look
further in the future, we also simulate the catalogue of one year observation
by ET. The expected number of detection is 168455 (expectation $1.68\times
10^{5}$). In figure 18, we plot the distribution of the catalogue in redshifts
and the total masses. We plot together the distribution of DNS mergers in the
whole Universe as defined by the population model for comparison. As we can
see from the upper panel of figure 18, the detected mass distribution is
shifted to the high mass side, due to higher detectability; and in the lower
panel of figure 18, we see the portion of detectable DNS merger decreases
towards higher redshift, as expected.

###  5.3 Example for neutron star/black hole mergers detected with Earth-based
detectors

In a similar way as above, we generate catalogues for ten years aLIGO
observation for BHNS mergers, both with the noise spectrum in O3 (aLIGO-O3)
and that at design (aLIGO-design). The underlying population model is BHNS-
PopB with the default parameters. The event number in the simulated catalogue
of aLIGO-O3 is 29 (a Poisson random with expectation value 32.0), which is
compatible with the rate found in the O3 period (2-3 in one year); the number
for aLIGO-design is 553 (expectation value 517.2). In the top and bottom
panels of figure 19, we plot the simulated catalogues in $z-m_{\bullet}$ and
$m_{\bullet}-m_{\rm{n}}$ planes. A more sensitive detector would help to
properly characterize this population.

We also simulate a catalogue for one year observation by ET. The number of
events is 6251 (expectation $6.32\times 10^{4}$). In the panels of figure 20,
we plot the distribution of the catalogue in $z$, $m_{\bullet}$ and
$m_{\rm{n}}$. We also plot the distribution of BHNS mergers in the whole
Universe as defined by the population model for comparison. We see the the
effect that the fraction of detectable BHNS mergers increases towards higher
masses, and decreases towards higher redshift, although ET probes essentially
the whole distribution.

Figure 20: The Probability density distribution as function of $m_{\bullet}$
(upper panel), $m_{\rm{n}}$ (middle panel) and the number density as function
of red-shift (lower panel) for 10 years observation of BHNS mergers by ET. The
integral of the latter is the total number of detection. The orange lines show
the intrinsic distributions.

###  5.4 Examples for SMBBH mergers detected with LISA

Figure 21: The red-shifted chirp masses and the redshift of the total events
(markers) and detected ones (squares on top of markers) by the standard LISA
configuration in five years. Orange, blue and green markers correspond to
Q3_nodelay, pop3 and Q3_delays population models respectively. The underlying
SNR cutoff is 8.

We now turn to the space-based detectors. The GW-Toolbox simulates the
observed catalogue of SMBBH mergers by going through the catalogue in the
Universe for a given observation duration, calculate the SNR for each SMBBH
and select against the SNR cut-off, in this case SNR=8 as default. In figure
21, we plot the red-shifted chirp masses and the redshift of the total events
and detected ones (squares on top of markers) by the standard LISA
configuration in five years, corresponding to different population models. As
we can see from figure 21, the detection horizon of our default LISA passes
though the Pop3 population and below the Q3_delays and Q3_nodelays population.
As a result, almost all events in Q3_delays and Q3_nodelays population are
detectable, while the detectable fraction of Pop3 changes significantly with
different LISA noise settings.

Figure 22: Upper: The distribution of relative uncertainties on the total
masses in the catalogue of MBH mergers detected by LISA in five year. Solid
lines are for default LISA and dashed lines are for 5Gm arm length; Middle
Same as the upper panel, but on the relative uncertainties on the luminosity
distances; Lower: Same as the upper two panels, but on uncertainties of the
sky location (deg2).

For these sources, we also calculate the uncertainties in the parameters. In
panels of figure 22, we plot the distribution of the relative uncertainties of
total masses, distances and sky localization $d\Omega$ in units of square
degrees. The uncertainties are all estimated with a FIM method. The
uncertainties span a wide range, but a significant fraction has quite well
determined masses while only a small fraction has well determined distance and
sky position. Our findings are in general in agree with previous results in
Klein16.

Figure 23: The detected number of DWDs with the default LISA and a larger LISA with 5 million km arms as function of the observation duration. We use a threshold SNR $\rho_{\star}=7$ Name |  $f$ (Hz) |  $\beta$ (deg) |  $\lambda$ (deg) | $A$ |  $d\Omega$ (deg2) | SNR
---|---|---|---|---|---|---
J0806 |  $6.2\times 10^{-3}$$\pm$$3.1\times 10^{-8}$ | -4.704 | 120.442 |  $1.2\times 10^{-22}$$\pm$$6.5\times 10^{-23}$ | $6.9\times 10^{-3}$ | 91.851
V407 Vul |  $3.5\times 10^{-3}$$\pm$$3.1\times 10^{-8}$ | 46.783 | 294.995 |  $5.9\times 10^{-23}$$\pm$$1.1\times 10^{-22}$ | $1.2\times 10^{-1}$ | 65.930
ES Cet |  $3.2\times 10^{-3}$$\pm$$3.1\times 10^{-8}$ | -20.334 | 24.612 |  $4.7\times 10^{-22}$$\pm$$1.0\times 10^{-22}$ | $1.7\times 10^{-1}$ | 46.347
ZTF J153932.16+502738.8 |  $4.8\times 10^{-3}$$\pm$$3.1\times 10^{-8}$ | 66.162 | 205.031 |  $3.1\times 10^{-24}$$\pm$$1.8\times 10^{-22}$ | $1.3\times 10^{-2}$ | 188.051
SDSS J065133.34+284423.4 |  $2.6\times 10^{-3}$$\pm$$3.1\times 10^{-8}$ | 5.805 | 101.340 |  $1.9\times 10^{-23}$$\pm$$1.6\times 10^{-22}$ | 2.306 | 15.613
SDSS J093506.92+441107.0 |  $1.6\times 10^{-3}$$\pm$$3.1\times 10^{-8}$ | 28.091 | 130.980 |  $7.6\times 10^{-22}$$\pm$$3.0\times 10^{-22}$ | 19.776 | 10.653
Table 3: Detection of the verification binaries with the LISA detector.

###  5.5 Examples for DWDs detected with LISA

In order to simulate the DWDs observed by LISA, we go through a simulated
catalogue of all Galactic DWDs and calculate their SNR with the analytic
approximation in equation (37). Then we select the sources with SNR larger
than the SNR threshold of 10 as the detected sources. We do the same for the
verification binaries (VBs). In Table 3, we list the detected VBs and the
estimated uncertainties of parameters. When calculating the uncertainties with
FIM, the numerical waveform calculated with LDC are used. The results are
compatible with earlier work, the SNR a bit lower than Kupfer et al. (2018),
due the use of a slightly different LISA sensitivity. Since the size of the
DWD catalogue is huge ($\sim 2.6\times 10^{7}$), we go through a smaller sub-
catalogue that are randomly drawn from the whole synthetic DWD catalogue
instead, and scale the number of detection from this small sub-sample to the
whole GWD catalogue to obtain the total expected number of detection. In
figure 23, we plot the number of detection as function of the observation
duration. Our results are about a factor of 0.5 lower than the earlier results
e.g. in Nelemans et al. (2001); Nissanke et al. (2012); Korol et al. (2017).
We attribute this to a slight difference in the LISA noise model or the DWD
catalogue.

###  5.6 Examples for EMRIs detected with LISA

The signal from an EMRI can fall into the detectable frequency range of LISA
from an early phase untill the final plunge, which can span quite a long
duration from months to years. In order to guarantee the speed of simulation,
we use a frequency domain waveform (modulus) generated in advance and stored
in files. In each simulation, we go through all the candidate events in the
catalogue, read in the corresponding TDI waveform modulus and calculate the
SNR against the noise curve. The detected events are then selected against an
SNR-threshold. In the EMRIs catalogues that we are using, there are also pre-
calculated SNR for each system, which corresponds to a slightly different LISA
setup and waveform (see Babak17). In figure 24, we compare their pre-
calculated SNRs with ours of the same catalogue for one year observation of
the M1 population with our default LISA. Our calculated SNR values disperse
within a factor two around the values of Babak17. We attribute this dispersion
to the slightly different LISA noise and waveform (AK Schwarzchild versus AK
Kerr, see Babak17).

Figure 24: The pre-calculated SNR by Babak17 comparing with SNR of this work
for the same catalogue, for one years’ observation on M1 population.

In figure 25, we plot the histogram of the SNR of the bright EMRIs in the
catalogue of Babak17 and that calculated from this work. The underlying
population model is M1.

Figure 25: Histogram of the SNR of the bright EMRIs in catalogue of Babak17.
The blue histogram indicate the SNR calculated in this work, and the orange
histogram is that pre-calculated by Babak17 using a slightly different LISA
noise and waveform.

In panels of figure 26, we plot histograms of relative uncertainties of the
masses $\mu$, $M$, distance $D$ and sky location $\Omega$ (in units of square
degrees) in a catalogue detected by the default LISA. The observation duration
is two years, and the SNR cut-off is set to 20 and the population is M1. The
uncertainties are estimated with FIM (section 4.1.2). In the calculation of
FIM, derivatives of the complex waveform relative to all the relevant
parameters are needed (equation 23). Therefore, if we were to use the pre-
calculated waveform for the uncertainty estimation, the storing files would be
$\sim 20$ times larger in size than those used for the SNR calculation. On
that account, we calculate the late stage waveform in real time and use them
for the uncertainty evaluation. In general the parameters of EMRIs are very
well determined, expect in some cases the sky position. The estimated level of
uncertainties are in agreement with Babak17.

Figure 26: Histogram of relative uncertainties on the mass of the lighter
component $\mu$, the mass of the heavier component $M$, the distance $D$ and
the uncertainties of the sky location $d\Omega$ in a catalogue of EMRIs
detected by the default LISA for a two year observation. The SNR cut-off is
set to 20.

####  5.6.1 Results for PTAs

For PTAs, we calculate the detection limits for individual SMBBHs and a
stochastic background based on the different PTA configurations. Given a
certain PTA, a SNR cut-off and the coordinates of the source, we can give a
sensitivity curve for GW emitted by an individual SMBBH, as a function of
frequency. In figure 27, we plot the sky-averaged sensitivity curve to
individual sources of EPTA, PPTA, NANOGrav, IPTA and a simulated future PTA
(labeled “future”). For the future PTA, we assume daily observation on the
IPTA pulsars for 10 more years, with two more new pulsars being added to the
PTA per year (the setting of the future PTA can be customized by users). The
corresponding $\rho_{\rm cri}=10$. The results are in agreement with published
ones (Babak et al., 2016; Schutz & Ma, 2016; Aggarwal et al., 2019).

Figure 27: Sensitivity plot for an individual SMBBH as function of frequency,
for $\rho_{\rm cri}=10$ and averaged over the celestial sphere for different
PTAs.

In Table 4, we list the upper limits on the SGWB of different origins, the
corresponding $\rho_{\rm cri}=100$. They are broadly in agreement with
published results (to within a factor $\sim 2$), which are also listed in the
table in the parentheses.

|  SBHBH ($A_{\rm yr^{-1}}$) |  Cosmic Strings ($\Omega_{\rm gw,yr^{-1}}h^{2}$) |  Relic ($\Omega_{\rm gw,yr^{-1}}h^{2}$)
---|---|---|---
EPTA (Lentati et al., 2015) |  $1.4\times 10^{-15}$($3\times 10^{-15}$) | $2.5\times 10^{-10}$ |  $4.7\times 10^{-10}$($1.2\times 10^{-9}$)
PPTA (Shannon et al., 2015) |  $2.45\times 10^{-15}$($1\times 10^{-15}$) | $6.4\times 10^{-10}$ | $1.2\times 10^{-9}$
NANOGrav (Arzoumanian et al., 2018) |  $1.6\times 10^{-15}$($1.45\times 10^{-15}$) | $3.6\times 10^{-10}$ |  $6.3\times 10^{-10}$($3.4\times 10^{-10}$)
IPTA | $9.9\times 10^{-16}$ | $1.0\times 10^{-10}$ | $2.0\times 10^{-10}$
Table 4: The upper limits set with different PTAs to SGWB of different
origins, $\rho_{\rm cri}=100$. Numbers in parentheses are values in
literature. The most recent reported upper limit on the GW originates from
Cosmic String are in terms of the cosmic-string tension ($G\mu$). The
conversion from $G\mu$ to $\Omega_{\rm gw,yr^{-1}}h^{2}$ depends on models and
the reconnection probability $p$, so we cannot give literature values.

##  6 Caveats & Discussion

We have implemented and described a first version of the GW-Toolbox, which
still has a few caveats and is missing some ingredients that we plan to
implement in the (near) future. The most important caveats of the current
version are:

  * •

Higher order multipole modes of the waveform: Theoretically, GW radiation are
treated with spin-weighed spherical decomposition (Thorne, 1980). The leading
term is the $\ell=2$, $|m|=2$ (quadrupole) term. In the current GW-Toolbox, we
apply a waveform that neglects the higher order modes beyond quadrupole. High
order modes will extend the spectrum of a GW to high frequencies (London et
al., 2018), and their significance increase with higher mass inequality,
higher total mass and higher orbital inclination angle (Mills & Fairhurst,
2021). Including higher order modes in the waveform can break the degeneracy
between the inclination angle and the distance. thus results in more accurate
estimate of the parameters of the source (London et al., 2018). In the O3
observation, there are already some GW events with evidence of higher order
modes (GW170729: Chatziioannou et al. 2019; GW190412: Abbott et al. 2020b;
GW190814: Abbott et al. 2020d). Mills & Fairhurst (2021) found that a few
percent of binaries will be detected with higher order modes by aLIGO with
designed sensitivity. For 3G GW detectors, Divyajyoti et al. (2021) found that
the SNR contribution from higher order modes in some GWTC-3 systems would
surpass 10% of the total SNR. For larger total mass and mass ratio binaries,
the higher order modes are expecting to contribute a significant fraction of
the SNR. Therefore, neglecting higher order modes in the waveform will results
in an underestimate of the detection rate of sources, especially for those
with large total mas and mass ratio. In an updated GW-Toolbox, the waveform
will be replaced with one that includes the higher order modes, e.g., PhenomHM
(London et al., 2018).

  * •

Uncertainties estimation for Earth-based detectors: the uncertainties are
estimated in a hybrid way: those of the intrinsic parameters ($m_{i,z}$ and
$\chi$) are calculated with FIM, with all of the amplitude parameters treated
as one. In addition, the uncertainty on the red-shift ($\delta z$) is
estimated with empirical relations. Then the source-frame masses uncertainties
are calculated with error propagation equation (Equation 24). Equation (24) is
valid only when the covariance between $dm_{i,z}$ and $dz$ is vanishing. Since
$m_{i,z}$ are largely determined from matching the phases of the waveform,
while the constraint on $z$ comes from the fitting amplitude of the waveform,
we expect that the dependence between $dm_{i,z}$ and $dz$ are weak. A more
self-consistent and accurate non-Bayesian treatment would be: to treat all
extrinsic parameters separately in the FIM (the dimensional of FIM will thus
be enlarged), and calculate the covaraince among all parameters explicitly.
Then impose the constraints from the triangulation localization and obtain the
uncertainties on all parameters.

  * •

Populations and Waveforms of EMRIs: In order to return the simulated catalogue
of detection in a tolerable time for a website user, we use catalogues of
EMRIs in which only bright ones are included (pre-calculated
SNR${}_{\rm{tot}}>20$). Therefore, the user should not set an SNR cut-off
lower than $\sim 15$, otherwise the returned synthetic catalogue is
incomplete. In order to compare with previous results of Babak17, the waveform
we employed is analytic kludge, which is known to be fast but less realistic.
In the GW-Toolbox, it can be replaced with a more accurate waveform augmented
analytic kludge (AAK) easily, since the latter can also be simulated with the
same package EMRI_Kludge_Suite. For the sake of the speed of simulation, we
use the pre-calculated frequency domain TDI waveform for the SNR calculation.
When estimating the uncertainties using FIM methods of the detected sources,
we only include the late stage of their waveform. As shown in the above
sections, we will miss some of the low frequency section corresponding to the
early inspiral stage, which is in fact detectable by LISA, and results in
underestimation of the accuracy of parameters inference. As shown in examples
in the above section, such underestimation is not severe.

  * •

PTAs: Our sensitivity curves and upper limits are given according to an SNR
threshold and the SNR of GW are calculated based on a simplified parameterised
noise spectra. On the other hand, upper limits are reported in literature with
confidence levels. Due to the very different nature between the methods, the
correspondence between the confidence level in literature and our SNR
threshold is difficult to explore. In our examples, we set the
$\rho_{\rm{cri}}=10$ for continuous GW in plotting the sensitivity curves, and
$\rho_{\rm{cri}}=100$ for SGWB, in order to obtain results which are in order
of magnitude in accordance with literature. For continuous GW, the sensitivity
scales with the $\rho_{\rm{cri}}$; while for SGWB, upper limits on the
characteristic strain scale with the square root of $\rho_{\rm{cri}}$.

We plan to include a number of additions to the GW-Toolbox in the future. The
first ones involve several proposed space-borne detectors, in particular
DECIGO, Taiji and Tianqin. In the ground-based and space-borne modules, we
will integrate SNR calculators for individual sources that are specified by
the user. On a longer time scale, we plan to include more GW sources, e.g.,
supernovae explosions, single spinning and recycling neutron stars, multiple
black holes encounters and catalogs of SMBBH for PTAs. We are also working to
extend the GW-Toolbox with electro-magnetic counterparts, e.g., to return the
fluence of short GRBs and peak fluxes of kilonovae. Triangulation of a network
of ground-based detectors will also be developed.

In the next step, the GW-Toolbox will have the ability to simulate
observations of different evolutionary phases of the same population in
different GW frequency ranges. For instance, each population of compact object
mergers corresponds to a population of persistent GW source from the earlier
phase. The former are targets of ground-based interferometers, while the
latter are targets of space-borne interferometers. Another instance is the
close orbit and inspiral-merger phases of SMBBH, which can be observed with
PTA and LISA respectively. We will also include simulated observation on SGWB
with ground-based and LISA-like detectors.

##  7 summary

In this paper, we introduce the GW-Toolbox (www.gw-universe.org), a set of
tools that quickly simulates a Universe with kHz/mHz/nHz GW source
populations, and observations with different GW detectors, i.e., ground-based
interferometers, space-borne interferometers and pulsar timing arrays. We
hereby summarize the functionalities and methodologies of the GW-Toolbox for
each module:

  * •

The module for ground-based interferometers simulates observation of mergers
of compact objects, including binary black holes (BBH), double neutron stars
(DNS) and black hole-neutron stars (BHNS). The detectors include default
aLIGO, Virgo, KAGRA, Einstein Telescope and Cosmic Explorer, and a user
customised LIGO/Virgo-like and ET-like detectors. The noise curves of the
default detectors are taken from the literature, and those of user customised
detectors are simulated with the FINESSE software. After the noise curve and
antenna patterns are determined, we calculate the optimal SNR for all sources
in the selected source class. A simplified IMRPhenomD waveform that assumes
zero effective spin is employed in this step. With a certain SNR-threshold of
detection, we marginalize the geometrical parameters and obtain the
detectability $\mathcal{D}(m_{1},m_{2},z,\chi)$ as function of the source’s
masses, redshift (luminosity distance) and effective spin. The product of
$\mathcal{D}(m_{1},m_{2},z,\chi)$ and the user-selected probability density
function (p.d.f) of the source population defines the p.d.f of the detectable
sources, $N_{d}(m_{1},m_{2},z,\chi)$. A synthetic catalogue of observations is
obtained with a MCMC sampling from the $N_{d}(m_{1},m_{2},z,\chi)$. We use
Fisher Information Matrix (FIM) method to estimate the uncertainties of the
parameters of events. In the process of calculating the FIM, we apply the
complete IMRPhenomD waveform phases (for aligned spins).

  * •

The module for space-borne interferometers simulates observations with LISA or
a customised LISA-like configuration. The noise power density in the TDI-$X$
response channel is calculated with an analytical formula, which includes
acceleration noise, laser shot noise, other optical Meteorology noises and
confusion noise due to Galactic double white dwarf (DWD) foreground. The
targets we include are the inspiral of Supermassive Binary Black Holes
(SMBBH), individual resolvable Galactic DWD and Extreme Mass Ratio Insprials
(EMRIs). For SMBBH, we calculate the TDI-$X$ LISA responses of a GW source
with LDC codes. The optimal SNR is subsequently calculated. There are three
population models being considered, namely Pop3, Q3_nodelays and Q3_delays,
taken from Babak17. There are ten realizations of the simulated catalogues of
SMBBH mergers in the Universe within five years for each population model. The
GW-Toolbox will re-sample from the catalogue according to user specified
observation duration, and calculate the SNR for each source in the sample. A
synthetic detection catalogue is thus returned based on a SNR threshold of
detection set by the user. The uncertainties are estimated with FIM. For DWDs,
we use an analytic equation to calculate the modulus TDI-$X$ LISA response to
a series of sinusoidal GWs, and therefore the SNR. We consider two samples of
DWDs, namely the verification DWDs and a simulated entire Galactic population.
For the former sample, the SNR is calculated individually in the catalogue,
and a detected catalogue is returned according to a SNR threshold. For the
latter sample, due to its huge number, we randomly draw a sub-sample from it,
and find the catalogue of detections in the sub-sample. The total expected
number of detection is obtained by rescaling the number of events in the
returned catalogue. The uncertainties are also estimated with FIM, where we
use LDC codes for the complex TDI LISA response, instead of using the
analytical equation as in the SNR computation. For EMRIs, we make use of the
EMRI_Kludge_Suite for the TDI LISA response, and therefore the SNR. We
calculate the SNR for each source in pre-simulated catalogues of EMRIs of
different populations, and select those with SNR surpassing the SNR threshold.
The uncertainties are again computed with FIM.

  * •

In the PTA module, we include four currently operating PTAs: EPTA, PPTA,
NANOGrav and IPTA. For the pulsars in these PTAs, we use the following
parameters to represent their noises properties and observation campaigns: the
levels of white noise, the level of red noise, the red noise spectrum index,
total observation duration and averaged interval between observations. We
allow users to include new pulsars which will be discovered in the course of
future observations. The sky-locations of the new pulsars and their noise
properties are randomly assigned according to the distributions of those of
the known pulsars. In this module, the GW-Toolbox computes the SNR of a series
of monochromatic GWs with given frequency and amplitude, which corresponds to
a GW from the orbital motion of SMBBH. Another function of this module is to
evaluate the upper limit that a PTA can set to Stochastic GW Background (SGWB)
from different origins.

In the (near) future, the GW-Toolbox will be extended with new standard
detectors, triangulation of a network of ground-based detectors, new source
classes and electro-magnetic counterparts and the ability to ”observe” the
same source model with different detectors. In this way, the GW-Toolbox will
provide even more functionality to give users quick idea of the power of
different GW detectors for their favourite source population.

## Acknowledgements

We thank the many researchers that provided results for the different source
populations and detectors that made it possible to collect all these together
in the GW-Toolbox. We also thank the referee for comments and detailed testing
that greatly improved the paper and the GW-Toolbox. We thank in particular the
(Mock) LISA Data Challenge team and the participating groups that made it
possible to include so many results in the space module. Special thanks to
Stas Babak and Antoine Petiteau for help and support. We thank our colleagues
Paul Groot, Samaya Nissanke, Sarah Caudill, Chris van den Broeck, Gemma
Janssen, Antonia Rowlinson, Peter Jonker, Selma de Mink, Marc Klein-Wolt and
Jess Broderick for the initial discussions that led to this project. This
research was made possible by support from the Dutch National Science Agenda,
NWA Startimpuls – 400.17.608

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102. doi:10.1103/PhysRevLett.116.061102
  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, ApJ, 818, L22. doi:10.3847/2041-8205/818/2/L22
  * Abbott et al. (2017a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Phys. Rev. Lett., 119, 161101. doi:10.1103/PhysRevLett.119.161101
  * Abbott et al. (2017b) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 848, L12. doi:10.3847/2041-8213/aa91c9
  * Abbott et al. (2017c) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 848, L13. doi:10.3847/2041-8213/aa920c
  * Abbott et al. (2017d) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Nature, 551, 85. doi:10.1038/nature24471
  * Abbott et al. (2017e) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 850, L39. doi:10.3847/2041-8213/aa9478
  * Abbott et al. (2017f) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 850, L40. doi:10.3847/2041-8213/aa93fc
  * Abbott et al. (2018a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Phys. Rev. Lett., 121, 161101. doi:10.1103/PhysRevLett.121.161101
  * Abbott et al. (2018b) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Phys. Rev. Lett., 120, 091101. doi:10.1103/PhysRevLett.120.091101
  * Abbott et al. (2019a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Physical Review X, 9, 011001. doi:10.1103/PhysRevX.9.011001
  * Abbott et al. (2019b) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Physical Review X, 9, 031040. doi:10.1103/PhysRevX.9.031040
  * Abbott et al. (2019c) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Phys. Rev. Lett., 122, 061104. doi:10.1103/PhysRevLett.122.061104
  * Abbott et al. (2019d) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Phys. Rev. Lett., 123, 011102. doi:10.1103/PhysRevLett.123.011102
  * Abbott et al. (2019e) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, ApJ, 882, L24. doi:10.3847/2041-8213/ab3800
  * Abbott et al. (2020a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020, Classical and Quantum Gravity, 37, 045006. doi:10.1088/1361-6382/ab5f7c
  * Abbott et al. (2020b) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, Phys. Rev. D, 102, 043015. doi:10.1103/PhysRevD.102.043015
  * Abbott et al. (2020c) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, arXiv:2010.14527
  * Abbott et al. (2020d) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, ApJ, 896, L44. doi:10.3847/2041-8213/ab960f
  * Abbott et al. (2021a) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021, ApJ, 915, L5. doi:10.3847/2041-8213/ac082e
  * Abbott et al. (2021b) Abbott, R., Abbott, T. D., Acernese, F., et al. 2021, link: https://dcc.ligo.org/public/0170/P2000318/007/o3b_catalog.pdf
  * Acernese et al. (2015) Acernese, F., Agathos, M., Agatsuma, K., et al. 2015, Classical and Quantum Gravity, 32, 024001. doi:10.1088/0264-9381/32/2/024001
  * Aggarwal et al. (2020) Aggarwal, N., Aguiar, O. D., Bauswein, A., et al. 2020, arXiv:2011.12414
  * Aggarwal et al. (2019) Aggarwal, K., Arzoumanian, Z., Baker, P. T., et al. 2019, ApJ, 880, 116. doi:10.3847/1538-4357/ab2236
  * Aggarwal et al. (2020) Aggarwal, K., Arzoumanian, Z., Baker, P. T., et al. 2020, ApJ, 889, 38. doi:10.3847/1538-4357/ab6083
  * Alam et al. (2020) Alam, M. F., Arzoumanian, Z., Baker, P. T., et al. 2020, arXiv:2005.06495
  * Armstrong et al. (1999) Armstrong, J. W., Estabrook, F. B., & Tinto, M. 1999, ApJ, 527, 814
  * Amaro-Seoane et al. (2007) Amaro-Seoane, P., Gair, J. R., Freitag, M., et al. 2007, Classical and Quantum Gravity, 24, R113. doi:10.1088/0264-9381/24/17/R01
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv:1702.00786
  * Anholm et al. (2009) Anholm, M., Ballmer, S., Creighton, J. D. E., et al. 2009, Phys. Rev. D, 79, 084030. doi:10.1103/PhysRevD.79.084030
  * Arzoumanian et al. (2018) Arzoumanian, Z., Baker, P. T., Brazier, A., et al. 2018, ApJ, 859, 47. doi:10.3847/1538-4357/aabd3b
  * Ajith et al. (2011) Ajith, P., Hannam, M., Husa, S., et al. 2011, Phys. Rev. Lett., 106, 241101
  * Arzoumanian et al. (2020) Arzoumanian, Z., Baker, P. T., Blumer, H., et al. 2020, ApJ, 905, L34. doi:10.3847/2041-8213/abd401
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33. doi:10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123. doi:10.3847/1538-3881/aabc4f
  * Barack & Cutler (2004) Barack, L. & Cutler, C. 2004, PRD, 69, 082005
  * Barke et al. (2015) Barke, S., Wang, Y., Esteban Delgado, J. J., et al. 2015, Classical and Quantum Gravity, 32, 095004. doi:10.1088/0264-9381/32/9/095004
  * Babak et al. (2010) Babak, S., Baker, J. G., Benacquista, M. J., et al. 2010, Classical and Quantum Gravity, 27, 084009. doi:10.1088/0264-9381/27/8/084009
  * Babak & Sesana (2012) Babak, S. & Sesana, A. 2012, Phys. Rev. D, 85, 044034. doi:10.1103/PhysRevD.85.044034
  * Babak et al. (2016) Babak, S., Petiteau, A., Sesana, A., et al. 2016, MNRAS, 455, 1665. doi:10.1093/mnras/stv2092
  * Babak et al. (2017) Babak, S., Gair, J., Sesana, A., et al. 2017, Phys. Rev. D, 95, 103012. doi:10.1103/PhysRevD.95.103012
  * Barausse (2012) Barausse, E. 2012, MNRAS, 423, 2533. doi:10.1111/j.1365-2966.2012.21057.x
  * Barausse et al. (2020) Barausse, E., Berti, E., Hertog, T., et al. 2020, General Relativity and Gravitation, 52, 81. doi:10.1007/s10714-020-02691-1
  * Belczynski et al. (2002) Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407. doi:10.1086/340304
  * Belczynski et al. (2016) Belczynski, K., Holz, D. E., Bulik, T., et al. 2016, Nature, 534, 512. doi:10.1038/nature18322
  * Berti et al. (2006) Berti, E., Cardoso, V., & Will, C. M. 2006, Phys. Rev. D, 73, 064030. doi:10.1103/PhysRevD.73.064030
  * Brown et al. (2020) Brown, D. D., Jones, P., Rowlinson, S., et al. 2020, SoftwareX, 12, 100613. doi:10.1016/j.softx.2020.100613
  * Burke-Spolaor et al. (2019) Burke-Spolaor, S., Taylor, S. R., Charisi, M., et al. 2019, A&A Rev., 27, 5. doi:10.1007/s00159-019-0115-7
  * Caprini & Figueroa (2018) Caprini, C. & Figueroa, D. G. 2018, Classical and Quantum Gravity, 35, 163001. doi:10.1088/1361-6382/aac608
  * Chatziioannou et al. (2019) Chatziioannou, K., Cotesta, R., Ghonge, S., et al. 2019, Phys. Rev. D, 100, 104015. doi:10.1103/PhysRevD.100.104015
  * Chruślińska et al. (2019) Chruślińska, M., Nelemans, G., & Belczynski, K. 2019, MNRAS, 482, 5012. doi:10.1093/mnras/sty3087
  * Chruślińska et al. (2020) Chruślińska, M., Jeřábková, T., Nelemans, G., et al. 2020, A&A, 636, A10. doi:10.1051/0004-6361/202037688
  * Chua & Gair (2015) Chua, A. J. K. & Gair, J. R. 2015, Classical and Quantum Gravity, 32, 232002. doi:10.1088/0264-9381/32/23/232002
  * Chua et al. (2017) Chua, A. J. K., Moore, C. J., & Gair, J. R. 2017, Phys. Rev. D, 96, 044005. doi:10.1103/PhysRevD.96.044005
  * Clark et al. (1979) Clark, J. P. A., van den Heuvel, E. P. J., & Sutantyo, W. 1979, A&A, 72, 120
  * Colpi (2014) Colpi, M. 2014, Space Sci. Rev., 183, 189. doi:10.1007/s11214-014-0067-1
  * Cornish & Littenberg (2007) Cornish, N. J. & Littenberg, T. B. 2007, Phys. Rev. D, 76, 083006. doi:10.1103/PhysRevD.76.083006
  * Coulter et al. (2017) Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017, Science, 358, 1556. doi:10.1126/science.aap9811
  * Dahal (2020) Dahal, P. K. 2020, Journal of Astrophysics and Astronomy, 41, 8. doi:10.1007/s12036-020-9625-y
  * Damour & Vilenkin (2001) Damour, T. & Vilenkin, A. 2001, Phys. Rev. D, 64, 064008. doi:10.1103/PhysRevD.64.064008
  * Damour & Vilenkin (2005) Damour, T. & Vilenkin, A. 2005, Phys. Rev. D, 71, 063510. doi:10.1103/PhysRevD.71.063510
  * Demorest et al. (2013) Demorest, P. B., Ferdman, R. D., Gonzalez, M. E., et al. 2013, ApJ, 762, 94. doi:10.1088/0004-637X/762/2/94
  * Desvignes et al. (2016) Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, MNRAS, 458, 3341. doi:10.1093/mnras/stw483
  * Detweiler (1979) Detweiler, S. 1979, ApJ, 234, 1100. doi:10.1086/157593
  * de Mink & Belczynski (2015) de Mink, S. E. & Belczynski, K. 2015, ApJ, 814, 58. doi:10.1088/0004-637X/814/1/58
  * Divyajyoti et al. (2021) Divyajyoti, Baxi, P., Mishra, C. K., et al. 2021, Phys. Rev. D, 104, 084080. doi:10.1103/PhysRevD.104.084080
  * Dominik et al. (2015) Dominik, M., Berti, E., O’Shaughnessy, R., et al. 2015, ApJ, 806, 263. doi:10.1088/0004-637X/806/2/263
  * Ellis et al. (2012) Ellis, J. A., Siemens, X., & Creighton, J. D. E. 2012, ApJ, 756, 175. doi:10.1088/0004-637X/756/2/175
  * Ellis (2013) Ellis, J. A. 2013, Classical and Quantum Gravity, 30, 224004. doi:10.1088/0264-9381/30/22/224004
  * Estabrook et al. (2000) Estabrook, F. B., Tinto, M., & Armstrong, J. W. 2000, Phys. Rev. D, 62, 042002
  * Farr et al. (2011) Farr, W. M., Sravan, N., Cantrell, A., et al. 2011, ApJ, 741, 103. doi:10.1088/0004-637X/741/2/103
  * Farrow et al. (2019) Farrow, N., Zhu, X.-J., & Thrane, E. 2019, ApJ, 876, 18. doi:10.3847/1538-4357/ab12e3
  * Ferrarese & Merritt (2000) Ferrarese, L. & Merritt, D. 2000, ApJ, 539, L9. doi:10.1086/312838
  * Fryer & Kalogera (2001) Fryer, C. L. & Kalogera, V. 2001, ApJ, 554, 548. doi:10.1086/321359
  * Gair et al. (2004) Gair, J. R., Barack, L., Creighton, T., et al. 2004, Classical and Quantum Gravity, 21, S1595. doi:10.1088/0264-9381/21/20/003
  * Goetz & Riles (2011) Goetz, E. & Riles, K. 2011, Classical and Quantum Gravity, 28, 215006. doi:10.1088/0264-9381/28/21/215006
  * González et al. (2013) González, G., Viceré, A., & Wen, L. 2013, Frontiers of Physics, 8, 771. doi:10.1007/s11467-013-0329-5
  * Grishchuk (1976) Grishchuk, L. P. 1976, Soviet Journal of Experimental and Theoretical Physics Letters, 23, 293
  * Grishchuk (1977) Grishchuk, L. P. 1977, Eighth Texas Symposium on Relativistic Astrophysics, 302, 439. doi:10.1111/j.1749-6632.1977.tb37064.x
  * Harry & LIGO Scientific Collaboration (2010) Harry, G. M. & LIGO Scientific Collaboration 2010, Classical and Quantum Gravity, 27, 084006. doi:10.1088/0264-9381/27/8/084006
  * Hellings & Downs (1983) Hellings, R. W. & Downs, G. S. 1983, ApJ, 265, L39. doi:10.1086/183954
  * Hild et al. (2008) Hild, S., Chelkowski, S., & Freise, A. 2008, arXiv:0810.0604
  * Hild et al. (2010) Hild, S., Chelkowski, S., Freise, A., et al. 2010, Classical and Quantum Gravity, 27, 015003. doi:10.1088/0264-9381/27/1/015003
  * Hild et al. (2011) Hild, S., Abernathy, M., Acernese, F., et al. 2011, Classical and Quantum Gravity, 28, 094013. doi:10.1088/0264-9381/28/9/094013
  * Hobbs et al. (2006) Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655. doi:10.1111/j.1365-2966.2006.10302.x
  * Hobbs et al. (2010) Hobbs, G., Archibald, A., Arzoumanian, Z., et al. 2010, Classical and Quantum Gravity, 27, 084013. doi:10.1088/0264-9381/27/8/084013
  * Hobbs et al. (2010) Hobbs, G., Lyne, A. G., & Kramer, M. 2010, MNRAS, 402, 1027. doi:10.1111/j.1365-2966.2009.15938.x
  * Hobbs (2013) Hobbs, G. 2013, Classical and Quantum Gravity, 30, 224007. doi:10.1088/0264-9381/30/22/224007
  * Hobbs & Dai (2017) Hobbs, G. & Dai, S. 2017, arXiv:1707.01615
  * Huang et al. (2020) Huang, S.-J., Hu, Y.-M., Korol, V., et al. 2020, arXiv:2005.07889
  * Hughes et al. (2021) Hughes, S. A., Warburton, N., Khanna, G., et al. 2021, arXiv:2102.02713
  * Inayoshi et al. (2020) Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27. doi:10.1146/annurev-astro-120419-014455
  * Janssen et al. (2015) Janssen, G., Hobbs, G., McLaughlin, M., et al. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 37
  * Jenet et al. (2004) Jenet, F. A., Lommen, A., Larson, S. L., et al. 2004, ApJ, 606, 799. doi:10.1086/383020
  * Jenet et al. (2005) Jenet, F. A., Hobbs, G. B., Lee, K. J., et al. 2005, ApJ, 625, L123. doi:10.1086/431220
  * Kagra Collaboration et al. (2019) Kagra Collaboration, Akutsu, T., Ando, M., et al. 2019, Nature Astronomy, 3, 35. doi:10.1038/s41550-018-0658-y
  * Kawamura et al. (2006) Kawamura, S., Nakamura, T., Ando, M., et al. 2006, Classical and Quantum Gravity, 23, S125. doi:10.1088/0264-9381/23/8/S17
  * Kiziltan et al. (2013) Kiziltan, B., Kottas, A., De Yoreo, M., et al. 2013, ApJ, 778, 66. doi:10.1088/0004-637X/778/1/66
  * Klein et al. (2016) Klein, A., Barausse, E., Sesana, A., et al. 2016, Phys. Rev. D, 93, 024003
  * Kormendy & Ho (2013) Kormendy, J. & Ho, L. C. 2013, ARA&A, 51, 511. doi:10.1146/annurev-astro-082708-101811
  * Korol et al. (2017) Korol, V., Rossi, E. M., Groot, P. J., et al. 2017, MNRAS, 470, 1894
  * Kovetz et al. (2017) Kovetz, E. D., Cholis, I., Breysse, P. C., et al. 2017, Phys. Rev. D, 95, 103010. doi:10.1103/PhysRevD.95.103010
  * Kramer & Champion (2013) Kramer, M. & Champion, D. J. 2013, Classical and Quantum Gravity, 30, 224009. doi:10.1088/0264-9381/30/22/224009
  * Kupfer et al. (2018) Kupfer, T., Korol, V., Shah, S., et al. 2018, MNRAS, 480, 302. doi:10.1093/mnras/sty1545
  * Lamberts et al. (2018) Lamberts, A., Garrison-Kimmel, S., Hopkins, P. F., et al. 2018, MNRAS, 480, 2704. doi:10.1093/mnras/sty2035
  * Lamberts et al. (2019) Lamberts, A., Blunt, S., Littenberg, T. B., et al. 2019, MNRAS, 490, 5888. doi:10.1093/mnras/stz2834
  * Lasky et al. (2016) Lasky, P. D., Mingarelli, C. M. F., Smith, T. L., et al. 2016, Physical Review X, 6, 011035. doi:10.1103/PhysRevX.6.011035
  * Lau et al. (2020) Lau, M. Y. M., Mandel, I., Vigna-Gómez, A., et al. 2020, MNRAS, 492, 3061. doi:10.1093/mnras/staa002
  * Lentati et al. (2015) Lentati, L., Taylor, S. R., Mingarelli, C. M. F., et al. 2015, MNRAS, 453, 2576. doi:10.1093/mnras/stv1538
  * Linde (1982) Linde, A. D. 1982, Physics Letters B, 108, 389. doi:10.1016/0370-2693(82)91219-9
  * Lipunov et al. (1997) Lipunov, V. M., Postnov, K. A., & Prokhorov, M. E. 1997, New A, 2, 43. doi:10.1016/S1384-1076(97)00007-9
  * London et al. (2018) London, L., Khan, S., Fauchon-Jones, E., et al. 2018, Phys. Rev. Lett., 120, 161102. doi:10.1103/PhysRevLett.120.161102
  * Luo et al. (2016) Luo, J., Chen, L.-S., Duan, H.-Z., et al. 2016, Classical and Quantum Gravity, 33, 035010. doi:10.1088/0264-9381/33/3/035010
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415
  * Maggiore (2008) Maggiore, M. Gravitational Waves: Volume 1: Theory and Experiments. (OUP Oxfor,2008)
  * Maggiore et al. (2020) Maggiore, M., Van Den Broeck, C., Bartolo, N., et al. 2020, J. Cosmology Astropart. Phys., 2020, 050. doi:10.1088/1475-7516/2020/03/050
  * Magorrian et al. (1998) Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285. doi:10.1086/300353
  * Mapelli et al. (2017) Mapelli, M., Giacobbo, N., Ripamonti, E., et al. 2017, MNRAS, 472, 2422. doi:10.1093/mnras/stx2123
  * McLaughlin (2013) McLaughlin, M. A. 2013, Classical and Quantum Gravity, 30, 224008. doi:10.1088/0264-9381/30/22/224008
  * McWilliams et al. (2014) McWilliams, S. T., Ostriker, J. P., & Pretorius, F. 2014, ApJ, 789, 156. doi:10.1088/0004-637X/789/2/156
  * McWilliams et al. (2019) McWilliams, S. T., Caldwell, R., Holley-Bockelmann, K., et al. 2019, arXiv:1903.04592
  * Mei et al. (2020) Mei, J., Bai, Y.-Z., Bao, J., et al. 2020, arXiv:2008.10332
  * Mehta et al. (2017) Mehta, A. K., Mishra, C. K., Varma, V., et al. 2017, Phys. Rev. D, 96, 124010. doi:10.1103/PhysRevD.96.124010
  * Miller & Miller (2015) Miller, M. C. & Miller, J. M. 2015, Phys. Rep, 548, 1. doi:10.1016/j.physrep.2014.09.003
  * Mills & Fairhurst (2021) Mills, C. & Fairhurst, S. 2021, Phys. Rev. D, 103, 024042. doi:10.1103/PhysRevD.103.024042
  * Moe & Di Stefano (2017) Moe, M. & Di Stefano, R. 2017, ApJS, 230, 15. doi:10.3847/1538-4365/aa6fb6
  * Murray et al. (2010) Murray, I., Adams, R. P., MacKay, D. J. C.. , Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, PMLR 9:541-548
  * Neal et al. (1999) Neal, R.M. 1999, Bayesian Statistics, 6, 475
  * Nelemans (2018) Nelemans, G. 2018, arXiv:1807.01060
  * Nelemans et al. (2001) Nelemans, G., Yungelson, L. R., & Portegies Zwart, S. F. 2001, A&A, 375, 890. doi:10.1051/0004-6361:20010683
  * Nissanke et al. (2012) Nissanke, S., Vallisneri, M., Nelemans, G., et al. 2012, ApJ, 758, 131. doi:10.1088/0004-637X/758/2/131
  * Ölmez et al. (2010) Ölmez, S., Mandic, V., & Siemens, X. 2010, Phys. Rev. D, 81, 104028. doi:10.1103/PhysRevD.81.104028
  * Ott (2009) Ott, C. D. 2009, Classical and Quantum Gravity, 26, 204015. doi:10.1088/0264-9381/26/20/204015
  * Özel et al. (2012) Özel, F., Psaltis, D., Narayan, R., et al. 2012, ApJ, 757, 55. doi:10.1088/0004-637X/757/1/55
  * Petiteau et al. (2013) Petiteau, A., Babak, S., Sesana, A., et al. 2013, Phys. Rev. D, 87, 064036. doi:10.1103/PhysRevD.87.064036
  * Phinney (1991) Phinney, E. S. 1991, ApJ, 380, L17. doi:10.1086/186163
  * Phinney (2001) Phinney, E. S. 2001, astro-ph/0108028
  * Pian et al. (2017) Pian, E., D’Avanzo, P., Benetti, S., et al. 2017, Nature, 551, 67. doi:10.1038/nature24298
  * Planck Collaboration et al. (2016) Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016, A&A, 594, A1. doi:10.1051/0004-6361/201527101
  * Porayko et al. (2019) Porayko, N. K., Noutsos, A., Tiburzi, C., et al. 2019, MNRAS, 483, 4100. doi:10.1093/mnras/sty3324
  * Portegies Zwart & McMillan (2000) Portegies Zwart, S. F. & McMillan, S. L. W. 2000, ApJ, 528, L17. doi:10.1086/312422
  * Postnov & Yungelson (2014) Postnov, K. A. & Yungelson, L. R. 2014, Living Reviews in Relativity, 17, 3. doi:10.12942/lrr-2014-3
  * Postnov & Kuranov (2019) Postnov, K. A. & Kuranov, A. G. 2019, MNRAS, 483, 3288. doi:10.1093/mnras/sty3313
  * Punturo et al. (2010) Punturo, M., Abernathy, M., Acernese, F., et al. 2010, Classical and Quantum Gravity, 27, 194002. doi:10.1088/0264-9381/27/19/194002
  * Regimbau et al. (2012) Regimbau, T., Dent, T., Del Pozzo, W., et al. 2012, Phys. Rev. D, 86, 122001. doi:10.1103/PhysRevD.86.122001
  * Reitze et al. (2019) Reitze, D., Adhikari, R. X., Ballmer, S., et al. 2019, BAAS
  * Robson et al. (2019) Robson, T., Cornish, N. J., & Liu, C. 2019, Classical and Quantum Gravity, 36, 105011. doi:10.1088/1361-6382/ab1101
  * Rodriguez et al. (2013) Rodriguez, C. L., Farr, B., Farr, W. M., et al. 2013, Phys. Rev. D, 88, 084013. doi:10.1103/PhysRevD.88.084013
  * Ruiter et al. (2010) Ruiter, A. J., Belczynski, K., Benacquista, M., et al. 2010, ApJ, 717, 1006. doi:10.1088/0004-637X/717/2/1006
  * Sathyaprakash & Schutz (2009) Sathyaprakash, B. S. & Schutz, B. F. 2009, Living Reviews in Relativity, 12, 2. doi:10.12942/lrr-2009-2
  * Sazhin (1978) Sazhin, M. V. 1978, Sov. Ast., 22, 36
  * Sesana (2016) Sesana, A. 2016, Phys. Rev. Lett., 116, 231102. doi:10.1103/PhysRevLett.116.231102
  * Sesana et al. (2005) Sesana, A., Haardt, F., Madau, P., et al. 2005, ApJ, 623, 23. doi:10.1086/428492
  * Sesana et al. (2008) Sesana, A., Vecchio, A., & Colacino, C. N. 2008, MNRAS, 390, 192. doi:10.1111/j.1365-2966.2008.13682.x
  * Sesana & Vecchio (2010) Sesana, A. & Vecchio, A. 2010, Phys. Rev. D, 81, 104008. doi:10.1103/PhysRevD.81.104008
  * Sesana et al. (2020) Sesana, A., Lamberts, A., & Petiteau, A. 2020, MNRAS, 494, L75. doi:10.1093/mnrasl/slaa039
  * Schneider et al. (2001) Schneider, R., Ferrari, V., Matarrese, S., et al. 2001, MNRAS, 324, 797. doi:10.1046/j.1365-8711.2001.04217.x
  * Schutz (1989) Schutz, B. F. 1989, Classical and Quantum Gravity, 6, 1761. doi:10.1088/0264-9381/6/12/006
  * Schutz & Ma (2016) Schutz, K. & Ma, C.-P. 2016, MNRAS, 459, 1737. doi:10.1093/mnras/stw768
  * Shannon et al. (2013) Shannon, R. M., Ravi, V., Coles, W. A., et al. 2013, Science, 342, 334. doi:10.1126/science.1238012
  * Shannon et al. (2015) Shannon, R. M., Ravi, V., Lentati, L. T., et al. 2015, Science, 349, 1522. doi:10.1126/science.aab1910
  * Siemens et al. (2006) Siemens, X., Creighton, J., Maor, I., et al. 2006, Phys. Rev. D, 73, 105001. doi:10.1103/PhysRevD.73.105001
  * Siemens et al. (2007) Siemens, X., Mandic, V., & Creighton, J. 2007, Phys. Rev. Lett., 98, 111101. doi:10.1103/PhysRevLett.98.111101
  * Starobinsky (1980) Starobinsky, A. A. 1980, Physics Letters B, 91, 99. doi:10.1016/0370-2693(80)90670-X
  * Talbot & Thrane (2018) Talbot, C. & Thrane, E. 2018, ApJ, 856, 173. doi:10.3847/1538-4357/aab34c
  * Tauris (2018) Tauris, T. M. 2018, Phys. Rev. Lett., 121, 131105. doi:10.1103/PhysRevLett.121.131105
  * Timpano et al. (2006) Timpano, S. E., Rubbo, L. J., & Cornish, N. J. 2006, Phys. Rev. D, 73, 122001. doi:10.1103/PhysRevD.73.122001
  * Tinto & Armstrong (1999) Tinto, M. & Armstrong, J. W. 1999, Phys. Rev. D, 59, 102003
  * The LIGO Scientific Collaboration et al. (2015) LIGO Scientific Collaboration, Aasi, J., Abbott, B. P., et al. 2015, Classical and Quantum Gravity, 32, 074001. doi:10.1088/0264-9381/32/7/074001
  * The LIGO Scientific Collaboration et al. (2020) The LIGO Scientific Collaboration, the Virgo Collaboration, Abbott, R., et al. 2020, arXiv:2010.14533
  * Thorne (1980) Thorne, K. S. 1980, Reviews of Modern Physics, 52, 299. doi:10.1103/RevModPhys.52.299
  * Thorne (1987) Thorne, K. S. 1987, Three Hundred Years of Gravitation, 330
  * Toonen et al. (2012) Toonen, S., Nelemans, G., & Portegies Zwart, S. 2012, A&A, 546, A70. doi:10.1051/0004-6361/201218966
  * van Haasteren et al. (2009) van Haasteren, R., Levin, Y., McDonald, P., et al. 2009, MNRAS, 395, 1005. doi:10.1111/j.1365-2966.2009.14590.x
  * van Haasteren et al. (2011) van Haasteren, R., Levin, Y., Janssen, G. H., et al. 2011, MNRAS, 414, 3117. doi:10.1111/j.1365-2966.2011.18613.x
  * Vitale & Evans (2017) Vitale, S. & Evans, M. 2017, Phys. Rev. D, 95, 064052. doi:10.1103/PhysRevD.95.064052
  * Vitale et al. (2019) Vitale, S., Farr, W. M., Ng, K. K. Y., et al. 2019, ApJ, 886, L1
  * Yardley et al. (2010) Yardley, D. R. B., Hobbs, G. B., Jenet, F. A., et al. 2010, MNRAS, 407, 669. doi:10.1111/j.1365-2966.2010.16949.x
  * Valentim et al. (2011) Valentim, R., Rangel, E., & Horvath, J. E. 2011, MNRAS, 414, 1427. doi:10.1111/j.1365-2966.2011.18477.x
  * Vallisneri (2008) Vallisneri, M. 2008, Phys. Rev. D, 77, 042001. doi:10.1103/PhysRevD.77.042001
  * Veitch et al. (2015) Veitch, J., Raymond, V., Farr, B., et al. 2015, Phys. Rev. D, 91, 042003. doi:10.1103/PhysRevD.91.042003
  * Verbiest et al. (2016) Verbiest, J. P. W., Lentati, L., Hobbs, G., et al. 2016, MNRAS, 458, 1267. doi:10.1093/mnras/stw347
  * Yi et al. (2014) Yi, S., Stappers, B. W., Sanidas, S. A., et al. 2014, MNRAS, 445, 1245. doi:10.1093/mnras/stu1826
  * Yu & Jeffery (2010) Yu, S. & Jeffery, C. S. 2010, A&A, 521, A85. doi:10.1051/0004-6361/201014827
  * Zhang et al. (2019) Zhang, J., Yang, Y., Zhang, C., et al. 2019, MNRAS, 488, 5020. doi:10.1093/mnras/stz2020

##  Appendix A Conversion among different LISA responses

In above sections, when working with LISA responses to gravitational wave
signal and noises, we often need to convert among different kinds of LISA
responses. We summarise the conversion relationship in figure 28.

Figure 28: Conversion between different kinds of LISA responses