literature, one must first of all check the value of ${\cal R}_{0,{\rm BNS}}$
used.

Finally, similarly to the plots shown for BBH signals, in the left and central
panels of Fig. 17 we show scatter plots of the joint estimation errors for the
sky location and luminosity distance, and for detector–frame chirp mass and
luminosity distance, with a color code showing the redshift dependence. These
are again the basis for the statistical methods based on correlation with
galaxy catalogs and the mass function, the latter being particularly promising
given the expected narrowness of the BNS mass distribution (Taylor et al.,
2012; Finke et al., 2022). Exploiting statistical techniques is of fundamental
importance also for BNS systems even though they can have an associated
electromagnetic counterpart, since the detection of such EM emissions is
expected only for a fraction of the sources (Belgacem et al., 2019a; Ronchini
et al., 2022). In the right panels of Fig. 17 we instead report the forecast
of the relative errors attainable on the source–frame chirp mass and the
combination of the adimensional tidal deformability parameters
$\tilde{\Lambda}$ defined in Eq. (14) (with a color code giving information on
the total source–frame mass). These are the basic quantities needed to study
the equation of state of dense matter above the nuclear density. We find a
large amount of detections with percent–level accuracy attainable both at ET
alone and at ET+ 2CE detectors. This would lead to dramatic improvements in
our understanding of the equation of state of NSs [see e.g. Gupta et al.
(2022) for a recent discussion]. We again show a selection of some of the main
results presented in this section in Tab. 4.

BNS
---
Network | Detected | Analysed | $\rm SNR\geq 30$ | $\Delta d_{L}/d_{L}\leq 10\%$ |  $\Delta\Omega_{90\%}\leq$100\text{\,}{\mathrm{deg}}^{2}$$
LVK–O4 | 4 | 4 | 0 | 1 | 2
ET | 9577 | 9153 | 770 | 18 | 51
ET+2CE | 40107 | 39584 | 5289 | 2225 | 27331
Table 4: A selection of results from the analysis of the
$1\text{\times}{10}^{5}$ BNSs (corresponding to the full population in about
$1\text{\,}\mathrm{yr}$ with our choices for the parameters) at the considered
networks.

###  5.3 Neutron star–black hole binaries

Finally, we present our results for the population of NSBH systems, for which,
currently, there are only two detections, which took place in the second part
of the O3 run. In Abbott et al. (2021g), the local rate has been estimated to
be ${\cal R}_{0,{\rm NSBH}}=45^{+75}_{-33}\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ (assuming that the observed
NSBH are representative of the underlying population), or ${\cal R}_{0,{\rm
NSBH}}=130^{+112}_{-69}\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ (assuming a broad NSBH
population), while the re–analysis in Abbott et al. (2021c), using different
modelizations, produce broadly consistent rates, between $7.8$ and
$140\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$. In the
following, as a reference value, we will assume for definiteness the value
${\cal R}_{0,{\rm
NSBH}}=$45\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$, keeping
in mind that the current uncertainty on the local rate can result in an
increase or a decrease in the number of detections per year, by a factor $\sim
3$ in either direction.

For these systems, we then simulate a population of $N_{\rm
NSBH}=$4.5\text{\times}{10}^{4}$$ sources out to $z=20$ that, with our
reference value of ${\cal R}_{0,{\rm NSBH}}$ and the redshift dependence of
the merger rate that we have assumed, described in Sect. 4.1, would correspond
to the number of coalescences happening in about one year. The cumulative
distributions for the SNRs and for the parameter errors are shown in Fig. 18
and 19, for the relative fractions and the absolute number of detection,
respectively, while the redshift distributions are shown in Fig. 20 for ET and
in Fig. 21 for ET+2CE. We find that, at ${\rm SNR}\geq 12$, ET alone can
detect up to $24\%$ of the sources (corresponding, with our assumptions on the
population, to about $1.1\text{\times}{10}^{4}$ detections per year), while a
network consisting of ET + 2CE could be able to detect 66% of the full
population of sources (corresponding to about $3.0\text{\times}{10}^{4}$
detections per year). In contrast, for the current network of detectors during
O4, accounting for the duty cycle, we get 3 events, again in agreement with
the broad expectations of Abbott et al. (2020b). Assuming the current
detectors to be operational 100% of the time during O4 the number of
detections only raises to 4, thus again we do not plot the corresponding
results for the parameter estimation. Varying ${\cal R}_{0,\rm NSBH}$ in the
currently allowed range $(7.8-140)\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ (while keeping fixed the
redshift dependence of the merger rate and the other features of our
population model), the number of detection per year at LVK–O4 ranges between
one and 10, for ET varies in the range
$[$1.9\text{\times}{10}^{3}$,$3.4\text{\times}{10}^{4}$]$ and, for ET+2CE,
varies in the range $[$5.2\text{\times}{10}^{3}$,$9.2\text{\times}{10}^{4}$]$.
For parameter estimation, using ${\cal R}_{0,{\rm
NSBH}}=$45\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$, the
corresponding values of $\tilde{N}_{\rm det}$ are $9.4\text{\times}{10}^{3}$
for ET ($87\%$ of the detected events) and $2.9\text{\times}{10}^{4}$ for
ET+2CE ($99\%$).

In Fig. 22 we show scatter plots for some combinations of parameters. In
particular, in the left and central panels we show the scatter plots in the
$(\Delta\Omega_{90\%},\Delta d_{L}/d_{L})$ plane and in the $(\Delta{\cal
M}_{c}/{\cal M}_{c},\Delta d_{L}/d_{L})$ plane, respectively, with a color
code carrying information on the redshift. Also in the case of NSBHs, these
quantities are relevant when exploiting statistical techniques to extract
cosmological and astrophysical information from GWs. NSBH systems could be
particularly relevant for probing the so–called lower mass gap, whose
identification would have significant implications (Ezquiaga & Holz, 2022).
For this reason, in the right panels we show the relative errors in the
reconstruction of the source–frame chirp mass and symmetric mass ratio as a
function of the total source–frame mass. These parameters are in fact
fundamental to identify the edges of the lower mass gap, as well as
identifying the system.

The knowledge that we currently have on the distribution of NSBH systems is
quite limited. The large uncertainties characterising the astrophysical
properties of the NSBH population are apparent looking, e.g., at Fig. 14 in
Broekgaarden et al. (2021), where the authors compare the results for 420
different models, and already the shapes of the resulting distributions can be
significantly different. We thus find useful to study how our results for this
class of sources change with different assumptions. In our fiducial NSBH
model, discussed in Sect. 4.1, the BH mass distribution is relatively narrow,
and we assumed a time delay distribution $P(t_{d})\propto 1/t_{d}$, which
favours small time delays. We compare it with two of the models presented in
Broekgaarden et al. (2021): their fiducial model, characterised by a broader
distribution for the BH mass and a time delay distribution favouring higher
values, and the model denoted by O_111, which predicts smaller masses and low
time delays. These can be seen, respectively, as an “optimistic” and a
“pessimistic” scenario, since systems with higher masses and closer to the
observer are louder and easier to detect. Observe that even more optimistic or
pessimistic scenarios are possible, see e.g. Santoliquido et al. (2021b).

Figure 18: Cumulative distributions of the SNRs (for LVK–O4, ET and ET+2CE) and parameter errors (for ET and ET+2CE) for NSBH signals, using the waveform model IMRPhenomNSBH. Figure 19: As in Fig. 18, in terms of the total number of events rather than detection fraction. |  |
---|---|---
Figure 20: Redshift distributions of the NSBHs observed at ET alone, selected on the basis of different thresholds for the SNR (left panel), or setting ${\rm SNR}\geq 12$ and applying further cuts on $\Delta d_{L}/d_{L}$ (central panel), or on $\Delta\Omega_{90\%}$ (right panel). The black solid line corresponds to the total NSBH population in the astrophysical model that we have assumed. We set ${\cal R}_{0,{\rm NSBH}}=$45\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$, and show the results for the detections in one year (taking into account our assumptions of the duty cycle). Given the current uncertainly on ${\cal R}_{0,{\rm NSBH}}$, one should keep in mind that the absolute number can still change by a factor $\order{3}$. In each column, the upper panel shows the number of events per redshift bin, while the lower panel shows the corresponding cumulative distributions, normalized to the number of BNS events in our sample, $N_{\rm NSBH}=$4.5\text{\times}{10}^{4}$$. |  |
---|---|---
Figure 21: As in Fig. 20, for ET+2CE.
---
Figure 22: Scatter plots of the observed NSBH population at 3G detectors. On
the top row we report the results for ET alone and in the bottom for the ET +
2CE network. In each row, the left panel shows the distribution of errors on
the luminosity distance and sky location as a function of redhsift, the
central panel the distribution of errors on the luminosity distance and
detector–frame chirp mass as a function of redhsift, and the right panel the
distribution of errors on the source–frame chirp mass and symmetric mass ratio
as a function of the total source–frame mass. The numbers reported on the
dashed lines refer to the fraction of observed events lying inside the
corresponding region. Figure 23: Histograms of the parameter errors for NSBHs
with the ET+2CE network configuration, using the waveform model IMRPhenomNSBH,
and sampling the BH mass and system redshift from our baseline distribution,
and from two distributions taken from (Broekgaarden et al., 2021). The tails
of the distributions featuring extremely large values of the relative errors,
as in some of the panels, simply represents the contribution of events for
which the Fisher matrix approximation cannot be trusted.

To perform the comparison, we fitted the quantiles of the chosen distributions
provided in Broekgaarden et al. (2021),292929All data is publicly available at
https://zenodo.org/record/4574727#.Ymv8jy8QN70, in the distributionQuantiles
folder. with suitable metalog distributions (Keelin, 2016), capable of
reproducing PDFs with irregular features. The comparison is only made using
the ET + 2CE detector network, which has the highest number of detections, in
order to avoid a proliferation of plots and for ease of readability. We show
the results in Fig. 23, from which it is apparent that the uncertainties
characterising the distribution of these events can strongly impact our
results. In particular, as expected, the parameter errors for the events
characterised by the narrower rate distribution can be much lower, even by
more than one order of magnitude, as compared to the results obtained with our
baseline choices, which, in contrast, are in good agreement with the broader
distribution. Furthermore, also the SNRs obtained for the “optimistic”
population are on average much higher, resulting in an increased percentage of
detections, from $66\%$ to $89\%$, while the percentage of detections in the
“pessimistic” case remains consistent with our fiducial choices. The results
on NSBHs shown in this work could thus change significantly in the future,
when the population of NSBH systems will be better explored, and should only
be considered as a first step toward understanding the capabilities of GW
detectors for NSBH systems, although our estimates appear to be more on the
conservative side. Also in this case we summarise some of the main results
presented in this section in Tab. 5.

NSBH
---
Network | Detected | Analysed | $\rm SNR\geq 50$ | $\Delta d_{L}/d_{L}\leq 10\%$ |  $\Delta\Omega_{90\%}\leq$100\text{\,}{\mathrm{deg}}^{2}$$
LVK–O4 | 3 | 3 | 0 | 0 | 1
ET | 10792 | 9384 | 223 | 2 | 44
ET+2CE | 29707 | 29462 | 1785 | 2687 | 22821
Table 5: A selection of results from the analysis of the
$4.5\text{\times}{10}^{4}$ NSBHs (corresponding to the full population in
about $1\text{\,}\mathrm{yr}$ with our choices for the parameters) at the
considered networks.

###  5.4 ‘Golden binaries’ at 3G detectors

The detection rate of compact binaries at 3G detectors is huge. As we have
seen, for BBHs (keeping fixed the redshift dependence of the merger rate and
varying only the local rate) the number of detections per year should be in
the range $[$3.1\text{\times}{10}^{4}$,$8.0\text{\times}{10}^{4}$]$ at ET,
raising to $[$4.2\text{\times}{10}^{4}$,$1.1\text{\times}{10}^{5}$]$ for
ET+2CE. Similarly, for BNSs, again keeping fixed the redshift dependence of
the merger rate and varying only the local rate, the number of detections per
year are expected to be in the range
$[$9.1\text{\times}{10}^{2}$,$1.6\text{\times}{10}^{5}$]$ for ET, and
$[$3.8\text{\times}{10}^{3}$,$6.5\text{\times}{10}^{5}$]$ for ET+2CE, and even
larger numbers can be obtained with different models for the redshift
dependence of the merger rate, as we discuss in App. B. For NSBHs the
detections should be of order
$[$1.9\text{\times}{10}^{3}$,$3.4\text{\times}{10}^{4}$]$ for ET and
$[$5.2\text{\times}{10}^{3}$,$9.2\text{\times}{10}^{4}$]$ for ET+2CE.

These numbers are remarkable, and constitute one of the strength of 3G
detectors. However, such large rates pose new challenges, compared to the
situations at 2G detectors. For instance, a detection rate of
$$5\text{\times}{10}^{3}$\,{\rm BNS/yr}$ corresponds to about a BNS every
$2\text{\,}\mathrm{hr}$, while a detection rate of $8\text{\times}{10}^{4}$
BNS/yr means about one every $7\text{\,}\mathrm{min}$. As we already
mentioned, BNSs can stay in the bandwidth of 3G detectors for hours, or up to
order one day. This means that we will have to deal with overlapping signals,
as has been recognized already since several years (Regimbau & Hughes, 2009;
Regimbau et al., 2012; Meacher et al., 2016; Regimbau et al., 2017). In
particular, in Samajdar et al. (2021), assuming values of the BBH local rate
in the range $[15,37.9]\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ and BNS local rates in the
range $[80,810]\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$, consistent with the range
of values that we have used, it is found that, at a network ET+2CE, a BNS
signal will typically have tens of overlapping BBH and BNS signals, and it
will happen up to $\order{10^{4}}$ times per year that two signals enters
their final merger phase within seconds of each other (see also Himemoto et
al. (2021) for similar estimates).303030Overlapping signals also takes places
for galactic white dwarf binaries in LISA. For 3G detectors, however, the BBH
and BNS systems are observed in their full inspiral–merger–ringdown phases,
while white dwarf binaries in LISA are deep into their inspiral phase, so
several aspects of the problem are different. The bias that the overlap of
signals can introduce in the parameter reconstruction of each separate event
has been studied recently in Samajdar et al. (2021); Pizzati et al. (2022);
Himemoto et al. (2021); Antonelli et al. (2021); these works show that the
overlap of signals can indeed introduce a bias in parameter reconstruction,
when two signals with comparable SNR have a difference in merger times which
is no longer large compared to the accuracy to which merger times themselves
can be measured.

Beside problems related to the overlap of signals, another practical problem
related to the huge number of detections expected at 3G detectors is that the
computation of the full multi–dimensional posteriors of such a large number of
events might be computationally too expensive (see Smith et al. (2021) for
recent progress in this direction).

For some scientific questions, a possible solution to these problems, which
might be sub–optimal but has the advantage of simplicity, is to restrict to
events with a very large SNR, say ${\rm SNR}\geq 100$; we will refer to them
as ‘golden events’, or ‘golden binaries’. For some aspects of the science at
3G detectors such a strategy might have limited value. For instance, for
populations studies, the completeness of the sample of detections is a key
element. However, for precision studies, such as tests of General Relativity
or cosmological observations, the result can be largely dominated by the very
best events. For instance, in order to test gravity near the BH horizons and
to discriminate the BH solution of General Relativity from other possible
types of compact objects (see Cardoso et al. (2016) for general discussion)
one needs events with a sufficiently large SNR in the ringdown phase, which
would allow the extraction of the frequency and damping time of at least two
different quasi–normal modes, allowing us to perform ‘BH spectroscopy’ (Berti
et al., 2016). For this kind of questions, only events with a sufficiently
large SNR are relevant. The same can happen in cosmological studies that make
use of ‘dark sirens’, i.e. coalescing binaries without an observed
electromagnetic counterpart, by correlating them with galaxy catalogs. In this
case, the results that can be obtained for $H_{0}$, or for the parameter
$\Xi_{0}$ that characterizes modified GW propagation (Belgacem et al., 2018a,
b), is largely dominated by the events with the smallest localization volume,
that fall in a region where a galaxy catalog is sufficiently complete, see
Finke et al. (2021). Another example is provided by multi–messenger
observations; given a BNS detection rate of order one event per hour or more,
only for a small fraction of detections it will be conceivable to perform a
dedicated follow–up by electromagnetic observations.

In all these cases, for each redshift range of interest, it will be natural to
focus on the best characterized GW events. Generically, these are the ones
with the highest SNR although, depending on the application that one has in
mind, one might wish to perform the cut by selecting directly the events with
$\Delta d_{L}/d_{L}$ below a given value, or with an angular resolution below
a given value (of course, these choices are correlated, as we also see from
the scatter plots in Fig. 12). It is therefore interesting to discuss in more
detail the ensemble of detections with large SNR, or with small error on
distance, or with good angular resolution, that we generically call ‘golden
binaries’. In the left panel of Fig. 10 we showed the distribution in redshift
of BBHs at ET alone, for different cuts on the SNR, including very high
thresholds (${\rm SNR}\geq 100$, ${\rm SNR}\geq 150$ and ${\rm SNR}\geq 200$),
while in the central panel we considered the events detected, according to the
criterion ${\rm SNR}>12$, and we imposed further cuts based on $\Delta
d_{L}/d_{L}$, or (in the right panel), based on $\Delta\Omega_{90\%}$; Fig. 11
showed the same result for BBH at ET+2CE. These plots show the number of
events obtained in one year, for our reference value ${\cal R}_{0,{\rm
BBH}}=$17\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ of the
local BBH rate, taking into account the duty cycle, computed according to the
assumptions in Sect. 4.2. As we already mentioned, it is important to keep in
mind that, since we are looking for events in the tails of the distributions
(whether with respect to SNR, $\Delta d_{L}/d_{L}$, or $\Delta\Omega_{90\%}$),
changes in the rate, or in the observation time, or in the duty cycle, can
have significant effects. In particular, increasing one or the other, will in
general lead to the appearance of rarer events, i.e. events with larger SNR,
better distance measurements, and/or better angular localization.

From Fig. 10, or from Tab. 3, we can appreciate that, already at ET alone
(with our fiducial choices for the local rate, the redshift dependence of the
merger rate, and the BBH mass function) there will be $\order{800}$ BBH
detections per year with ${\rm SNR}\geq 100$, of which a significant fraction
is at $z\gtrsim 1$, and a few of them are even at redshifts as large as $z\sim
3-4$; in our sample realization, the farthest event with ${\rm SNR}\geq 100$
has $z\simeq 3.8$. In terms of accuracy on distance, there are
$\order{${10}^{3}$}$ events/yr with $\Delta d_{L}/d_{L}<0.1$, of which a large
fraction is at $z>1$; furthermore, there are $\order{200}$ events/yr with
$d_{L}$ measured to better than $5\%$, again with a large fraction of them at
$z>1$.

For ET+2CE, we see from Fig. 11 or Tab. 3 that these numbers raise by almost
one order of magnitude. There are now $\order{5000}$ events/yr with ${\rm
SNR}\geq 100$, of which $\order{150}$ events/yr have $z\gtrsim 3$. For the
luminosity distance, we find $\order{1000}$ events/yr with $d_{L}$ measured
better than $1\%$; again a significant fraction of them is at $z>1$, and, out
of them, $\order{40}$ events/yr are at redshifts as large as $z\sim 2-3$.
Thanks to triangulation, now also for the angular resolution we have
remarkable results: there are $\order{1200}$ events/yr at $z\gtrsim 1$
localized to better than $1\text{\,}{\mathrm{deg}}^{2}$, and tens of events,
still at cosmologically significant redshifts, localized to better than
$0.05\text{\,}{\mathrm{deg}}^{2}$. An especially remarkable result is that
about 4 BBH systems per year at ET+2CE could be localised in a region as small
as $0.01\text{\,}{\mathrm{deg}}^{2}$ or lower, meaning that the true host
galaxy could be identified even in absence of an EM counterpart. In fact,
considering the COSMOS2020 survey (Weaver et al., 2022), among the deepest and
most complete galaxy surveys to date, we find that, on average, only a single
galaxy is present in a region of $0.01\text{\,}{\mathrm{deg}}^{2}$ and a
redshift extension of 0.05 [in agreement with Eq. (7) of Singer et al. (2016),
see also the discussion in Borhanian et al. (2020)].313131We thank Nicola
Borghi for providing this estimate. Note that the similarity of the shapes of
the distributions for the SNR, for $\Delta d_{L}/d_{L}$ and for
$\Delta\Omega_{90\%}$ confirms that the events with high SNR, good accuracy on
$d_{L}$ and good angular localization are essentially the same, as was also
visible, on the whole ensemble of detections, from the scatter plots in Fig.
12.

For testing the near–horizon physics and the nature of BHs, as discussed
above, it is crucial to have an accurate reconstruction of the ringdown phase.
In the left panel of Fig. 24 we show the distribution of ${\rm SNR}_{\rm pm}$
for BBH systems, where ${\rm SNR}_{\rm pm}$ is defined as the portion of the
SNR obtained just from the post–merger phase, i.e. considering only the part
of the signal above the peak frequency of the waveform.323232This is the
frequency at which the amplitude of the signal, normalised to the overall
$f^{-\nicefrac{{7}}{{6}}}$ factor, reaches its maximum, determined internally
by IMRPhenomHM according to Eq. (20) of Khan et al. (2016b). We see that,
already at ET alone, there will be $\order{100}$ BBHs per year for which the
SNR of the post–merger phase is larger than 30, with a few events per year
reaching ${\rm SNR}_{\rm pm}\sim 80-100$. For ET+2CE these numbers further
rise by a factor of a few and, for the best events, ${\rm SNR}_{\rm pm}$ can
reach values of order $200$. In contrast, at LVK–O4, one can expect at most
$\order{1}$ events with ${\rm SNR}_{\rm pm}\gtrsim 10$.

(a)  (b)

Figure 24: Histograms of the SNR in the post–merger phase of BBH binaries
(left panel) and NSBH binaries (right panel), for 1 year of data. For BBHs the
computation is performed using the IMRPhenomHM waveform, while for NSBH
systems we adopted IMRPhenomNSBH.

Some properties of ‘golden events’ for BNS can be appreciated from Fig. 15 for
ET, and from Fig. 16 for ET+2CE. Given that BNS are intrinsically less loud,
now meaningful cuts are at lower values of the SNR. With our fiducial choice
for the population, at ET there are $\order{800}$ events per year at with
${\rm SNR}\geq 30$, $\order{150}$ with ${\rm SNR}\geq 50$, and $\order{20}$
with ${\rm SNR}\geq 100$, of which, in our sample, the farthest is at $z\simeq
0.19$. It is important to notice that these results have been obtained setting
${\cal R}_{0,{\rm
BNS}}=$105.5\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ and, as
we repeatedly stressed, there is an uncertainty of one order of magnitude, in
both directions, on this number. Similarly to the case of BBHs, for ET+2CE
there is a further improvement by typically one order of magnitude, both for
the number of sources with a given high value of the SNR and for those with a
small value of $\Delta d_{L}/d_{L}$. For instance, the number of BNS/yr with
${\rm SNR}\geq 50$ is $\order{1200}$. As usual, the improvement is even more
important on the localization, thanks to triangulation. While, for ET alone,
there are just about $\order{10}$ BNS/yr localized to better than
$50\text{\,}{\mathrm{deg}}^{2}$, with ET+2CE there are $\order{2400}$ BNS/yr
localized to better than $10\text{\,}{\mathrm{deg}}^{2}$, including events up
to $z\sim 1.5$, and $\order{100}$ events/yr localized to better than
$1\text{\,}{\mathrm{deg}}^{2}$, including events up to $z\sim 0.44$. To put
this result in the correct perspective, however, it is useful to recall, from
the discussion in Sect. 5.2, that, for multi–messenger observations,
localization regions of the order of hundreds of square degree, as can be
obtained by ET alone, are already sufficient to trigger the search of a
kilonova with large FOV telescopes, which are then needed to localize the
source to the arcsec precision necessary to point instruments such as the ELT
and perform spectroscopic studies. However, the remarkable improvement in
localization between ET and ET+2CE will be very important in some cosmological
applications, such as those making use of dark sirens by correlating them with
galaxy catalogs.

Finally, Fig. 20 and 21 show the analogous results for NSBHs. ‘Golden events’
for NSBHs could be particularly interesting for cosmological applications,
since these events, compared to BNSs, are louder because of the larger mass,
while they could still have an electromagnetic counterpart, if the NS is
tidally disrupted during the merger; furthermore, while the spin of NSs in
binaries is expected to be negligible, the spin of the BH in a NSBH binary in
general might not be small (see, however, Mandel & Fragos (2020); Broekgaarden
& Berger (2021); Mandel & Smith (2021)) and, if it induces a precession of the
orbital plane, it can lift the degeneracy between $d_{L}$ and $\cos\iota$,
leading to an improved distance measurement (Vitale & Chen, 2018). The plots
in Fig. 20 and 21 show a pattern intermediate between that for BBHs and that
for BNSs. This is expected since, as discussed in Sect. 4.1, for the models
that we have used, the typical BH mass in a NSBH binary is of order
$8\text{\,}\mathrm{M_{\odot}}$, so these systems are heavier than BNSs, but
lighter than typical BBHs and, in the bandwidth of 3G detectors, they are
louder than typical BNS, but weaker than typical BBHs. Finally, we observe
that NSBH binaries look less promising than BBHs for testing the post–merger
and ringdown phase of the final BH. In the right panel of Fig. 24 we show the
distribution of ${\rm SNR}_{\rm pm}$ for these systems, which shows that the
signal–to–noise ratio in the post–merger phase never exceed 10, even for
ET+2CE.333333In this case the computation of the peak frequency is less
straightforward, due to the various possible outcomes of the merger, depending
on the possible disruption of the NS and the formation of an accretion torus
(Lattimer & Schramm, 1976). We thus proceeded according to the prescriptions
in Sect. IV of Pannarale et al. (2015), and checked with an actual numerical
computation of the full waveform.

##  6 Conclusions

In this paper we have performed a comprehensive study of the capabilities of
the Einstein Telescope, alone and in a network with two Cosmic Explorer
detectors, and we have also provided forecasts for the forthcoming LVK–O4 run.
The results have been obtained using a new parameter estimation code, GWFAST,
conceived especially for application to third–generation GW detectors. Our
results broadly confirm the overall picture obtained in other recent studies
of 3G detectors (see in particular Maggiore et al. (2020); Borhanian &
Sathyaprakash (2022); Harms et al. (2022); Ronchini et al. (2022), and
references therein), on the extraordinary scientific potential of 3G
detectors. We have applied our code to study populations of BBHs, BNSs and
NSBH binaries. Compared to some of the other works devoted to forecast for 3G
detectors, our results are more tuned toward ET, including detailed studies of
what ET can do as a single detector; we used updated information on the
astrophysical population from the LVK population studies based on the GWTC–3
catalog of detections; we used state–of–the–art waveforms for all these
sources, and included the Earth’s motion in the detectors response, which is
especially relevant for BNSs; we expanded the range of parameters on which
inference is performed, including the spin of the BHs and the tidal
deformability of NSs; we also presented the first systematic study of
parameter estimation for NSBH binaries at 3G detectors.

From the technical point of view our code, which is presented in more detail
in the companion paper Iacovelli et al. (2022), has been designed so to
optimally exploit the parallel nature of the problem of computing a large
number of Fisher matrices for independent events, by vectorizing the
evaluation of the Fisher matrices even on a single CPU. This relies on the
implementation of the waveforms in Python, resulting in a gain in
computational speed, which motivates the name GWFAST. Furthermore, this allows
the use of automatic differentiation with the library JAX for the accurate
computation of derivatives. Together with this paper, we publicly release the
code GWFAST , which is available at https://github.com/CosmoStatGW/gwfast.
This paper is associated to version v1.0.1 which is archived on Zenodo
(Iacovelli & Mancarella, 2022). We also release the library WF4Py , available
at https://github.com/CosmoStatGW/WF4Py, which contains state–of–the–art
gravitational–wave waveforms in pure Python. This paper is associated to
version v1.0.0 which is archived on Zenodo (Iacovelli, 2022).

Acknowledgments. We are grateful to Marica Branchesi, Jan Harms, Ulyana
Dupletsa, Samuele Ronchini and the GSSI group for extremely useful discussions
and comparisons of our respective codes and results, all along this project.
We also thank Ssohrab Borhanian and B. S. Sathyaprakash for comparisons with
their code and Andrea Maselli, Ik Siong Heng, Mauro Pieroni and the members of
the ET OSB div. 9 for hosting discussions about comparisons of Fisher codes.
We thank Ed Porter for useful advice on the implementation of the code,
Christopher Finlay for helpful discussion about JAX usage, and Floor
Broekgaarden for a careful reading of the manuscript and extremely useful
comments. The research leading to these results has been conceived and
developed within the ET Observational Science Board (OSB). Our research is
supported by the Swiss National Science Foundation, grant 200020$\\_$191957,
and by the SwissMap National Center for Competence in Research. Computations
made use of the Yggdrasil cluster at the University of Geneva. , HEALPix
(Zonca et al., 2019).

##  Appendix A Details of the adopted population distributions

Here we report the shape and parameters of the population distributions
mentioned in the main text. All these are available in the public code
MGCosmoPop343434https://github.com/CosmoStatGW/MGCosmoPop. (Mancarella et al.,
2022), together with sampling routines, which we used to obtain the catalogs.

###  A.1 Redshift distribution

For the BBH merger distribution with redshift we adopt a model that follows
the star formation rate. As a caveat, we mention that this is just a first
approximation, and other effects could induce departures from the star
formation rate. In particular, metallicity–dependent effects might
significantly boost the merger rate at low metallicities, see Santoliquido et
al. (2021a) (in particular their Fig. 2) and Chruślińska (2022). Given the
current uncertainties, a simple model of the merger rate following the star
formation rate can be considered a meaningful starting point that, however,
will have to be updated as our understanding of the merger rate of compact
objects sharpens.

The Madau–Dickinson profile (Madau & Dickinson, 2014) for the star formation
rate density has the functional form (Madau & Fragos, 2017)

$\psi(z|\alpha_{z},\beta_{z},z_{p})\propto\frac{(1+z)^{\alpha_{z}}}{1+\left(\frac{1+z}{1+z_{p}}\right)^{\alpha_{z}+\beta_{z}}}\,.$
(A1)

Observationally, the merger rate evolution of BBHs is currently constrained
only at small redshift; assuming a profile
$\psi(z|\kappa)\propto(1+z)^{\kappa}$, the analysis of the first three LVK
observing runs gives the median value $\kappa=2.7$.353535All the values of the
population parameters given in this Appendix are taken from the results of the
combined analysis of the events coming from the O1, O2 and O3 runs, available
at https://zenodo.org/record/5655785#.YnUnPS8QN70. All values are the medians
of the estimations, contained in the file PowerLawPeakObsOneTwoThree.json.
Here we use the Madau–Dickinson profile, with the low–redshift slope fixed to
the value of $\alpha_{z}=2.7$. This is chosen to match the value of $\kappa$
since at low redshift the Madau–Dickinson profile reduces to
$\sim(1+z)^{\alpha_{z}}$. For the other parameters, we fix them to
$\beta_{z}=3$ and $z_{p}=2$, which, within current uncertainties, are typical
values used in the literature (Madau & Dickinson, 2014).

For BNS systems, instead, we convolve the Madau–Dickinson profile with a time
delay distribution $P(t_{d})\propto 1/t_{d}$, using a minimum time delay of
$20\text{\,}\mathrm{Myr}$, as suggested in Regimbau et al. (2012). To get a
faster software implementation, we further re–fit the curve obtained after the
convolution to another Madau–Dickinson profile, with new parameters
$\alpha_{z}$, $\beta_{z}$ and $z_{p}$ now determined from the fit; we obtain
the values $\alpha_{z}=1.42$, $\beta_{z}=4.62$, $z_{p}=1.84$, and this is what
we use to generate the population (the overall rescaling to the local rate is
discarded, since we always simulate a fixed amount of sources). A population
with the same parameters is also adopted for NSBH systems.

The distribution in redshift of the merger rate of compact binaries is finally

$p_{z}(z)\propto\psi(z|\alpha_{z},\beta_{z},z_{p})\frac{1}{1+z}\dfrac{{\rm
d}V_{c}}{{\rm d}z}(z)\,,$ (A2)

where ${\rm d}V_{c}/{\rm d}z$ is the comoving volume element and the factor of
$1+z$ relates the source times and the observer time. In a flat $\Lambda$CDM
Universe, the redshift is related to the luminosity distance of the source by

$d_{L}(z)=(1+z)\dfrac{c}{H_{0}}\int_{0}^{z}\dfrac{{\rm
d}\tilde{z}}{\sqrt{\Omega_{m}(1+\tilde{z})^{3}+\Omega_{\Lambda}}}\,,$ (A3)

(the contribution of radiation is negligible even at the redshifts explored by
3G detectors), and we use Planck18 parameters (Planck Collaboration et al.,
2020), i.e.
$H_{0}=$67.66\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$$,
$\Omega_{m}=0.3097$ and $\Omega_{\Lambda}=1-\Omega_{m}=0.6903$.

###  A.2 Mass distribution

The Power Law + Peak distribution (Abbott et al., 2021c) used for the BBH mass
function has the functional form, for the mass of the primary component
(always taken as the heaviest),

$p(m_{1}|\lambda_{\rm peak},\alpha,m_{\rm min},\delta_{m},m_{\rm
max},\mu_{m},\sigma_{m})=\left[(1-\lambda_{\rm
peak})\mathcal{P}(m_{1}|-\alpha,m_{\rm max})+\lambda_{\rm
peak}\mathcal{N}(m_{1}|\mu_{m},\sigma_{m})\right]S(m_{1}|m_{\rm
min},\delta_{m})$ (A4)

where $\mathcal{P}(m_{1}|-\alpha,m_{\rm max})$ is a normalized power–law
distribution with spectral index $-\alpha$ and high–mass cut–off $m_{\rm
max}$, $\mathcal{N}(m_{1}|\mu_{m},\sigma_{m})$ is a normalised Gaussian
distribution with mean $\mu_{m}$ and standard deviation $\sigma_{m}$, while
$S(m_{1}|m_{\rm min},\delta_{m})$ is a smoothing function of the form

$S(m_{1}|m_{\rm min},\delta_{m})=\begin{dcases}0&(m<m_{\rm min})\\\
[f(m-m_{\rm min},\delta_{m})+1]^{-1}&(m_{\rm min}\leq m<m_{\rm
min}+\delta_{m})\\\ 1&(m\geq m_{\rm min}+\delta_{m})\end{dcases}\,,$ (A5)

with

$f(m^{\prime},\delta_{m})={\rm
exp}\left\\{\dfrac{\delta_{m}}{m^{\prime}}+\dfrac{\delta_{m}}{m^{\prime}-\delta_{m}}\right\\}.$
(A6)

The conditional mass ratio distribution of this model is

$p(q|\beta_{q},m_{1},m_{\rm min},\delta_{m})\propto
q^{\beta_{q}}S(qm_{1}|m_{\rm min},\delta_{m}).$ (A7)

The numerical values of the parameters we used for this distribution are
$\lambda_{\rm peak}=0.039,\ \alpha=3.4,\ m_{\rm min}=5.1,\ \delta_{m}=4.8,\
m_{\rm max}=87,\ \mu_{m}=34,\ \sigma_{m}=3.6,\ \beta_{q}=1.1$.35

###  A.3 Spin distribution

The Default distribution used to simulate the BBH spins [introduced in Abbott
et al. (2019c)] assumes a Beta distribution for the spin magnitudes (taken to
be independent)

$p(\chi_{1,2}|\alpha_{\chi},\beta_{\chi})={\rm
Beta}(\alpha_{\chi},\beta_{\chi})\,,$ (A8)

while the orientations are drawn from

$p({\bf z}|\zeta,\sigma_{t})=\zeta{\cal N}_{t}({\bf
z}|0,\sigma_{t})+(1-\zeta){\rm I}({\bf z})\,,$ (A9)

where $z_{i}={\rm cos}\theta_{i}$, $\theta_{i}$ is the tilt angle between the
component spin and the binary’s orbital angular momentum, ${\cal N}_{t}({\bf
z}|0,\sigma_{t})$ is a truncated Gaussian distribution centered at 0 and with
standard deviation $\sigma_{t}$, while ${\rm I}({\bf z})$ denotes an isotropic
distribution. Though this model could be used to sample all the 6 components
of the spins of the two objects, we will always extract only the distribution
of the spin component aligned with the orbital angular momentum, setting the
others to zero. The numerical values of the parameters we used for this
distribution are $\zeta=0.66,\ \sigma_{t}=1.5,\ \mu_{\chi}=0.28,\
\sigma_{\chi}=0.03$,35 from which it follows $\alpha_{\chi}=1.6,\
\beta_{\chi}=4.12$.

##  Appendix B Comparison with other parameter estimation codes

With the recent boost in activities for Einstein Telescope and Cosmic
Explorer, parameter estimation for 3G detectors has become a particularly
important topic, both as a tool to develop the science case for these
detectors, and to inform the choices for optimal detector’s design and
detector network configuration, so several parameter estimation codes have
been developed recently. In this appendix we compare our code and our results
to those obtained with GWBENCH (Borhanian, 2021; Borhanian & Sathyaprakash,
2022), developed at Penn State, and with GWFISH, developed by the GSSI group
(Harms et al., 2022) which, together with GWFAST, are among the most complete
and advanced codes available, for parameter estimation at 3G detectors [see
also Chan et al. (2018); Grimm & Harms (2020); Nitz & Dal Canton (2021); Li et
al. (2022); Pieroni et al. (2022)].

We begin by observing that, in the context of the activities of the ET
Observational Science Board (OSB), our group, the GSSI group and the Penn
State group have analyzed a small sample of common injections, using all the
same waveform model [the TaylorF2 waveform, developed in Buonanno et al.
(2009a)], finding excellent agreement, at the level of one part in
$10^{4}-10^{5}$ for the SNR, and one part in $10^{3}$ on the parameter
estimation of individual signals. This is an important consistency check for
these rather complex codes, that had been developed totally independently.
Further activity on this is in progress.

We next compare the results presented in our paper with those presented in
Borhanian & Sathyaprakash (2022), and with those in Ronchini et al. (2022)
(that provides an extensive set of results for BNSs obtained with GWFISH). The
comparison must take into account the fact that, currently, there are large
uncertainties on several aspects, ranging from the detector configurations
considered, the assumptions on the astrophysical models, as well as several
differences in technical details, such as the waveform models used, the
assumed threshold for the SNR, down to different implementations of the Fisher
matrix inversion. To orient the reader, we therefore find useful first of all
to give a summary of several different choices adopted in these works. We
stress that all these different choices are fully legitimate, and reflect
current uncertainties. It is therefore instructive to see how the results
change when using different assumptions. We also stress that, while the
results, in our paper as well as in Borhanian & Sathyaprakash (2022) and in
Harms et al. (2022); Ronchini et al. (2022), are necessarily presented just
for a selection of few different choices, the codes themselves can be adapted
to a larger variety of choices. The main differences in the assumptions used
in these works can be summarized as follows.

Detector networks. Our work is more focused on ET. To avoid a proliferation of
plots, or of lines in each plot, we only presented the results for two 3G
configurations: ET alone, and ET+2CE. Studying ET alone allows us to assess
the strength of the ET Science Case, independently of decisions that will be
taken by different funding agencies on CE. On the other hand, combining ET
with two CE detectors allows us to examine the full strength of a 3G detector
network. The two CE detectors are assumed to be one with
$40\text{\,}\mathrm{km}$ arms, and one with $20\text{\,}\mathrm{km}$ arms,
which is the reference CE configuration (Evans et al., 2021). We set both of
them in the U.S. (the other options that is under consideration, with one of
them in Australia, allows better triangulation but, otherwise, is not very
different). In contrast, in Ronchini et al. (2022) are studied ET alone, ET+CE
($40\text{\,}\mathrm{km}$) and ET+2CE, with both CE detectors with
$40\text{\,}\mathrm{km}$ arms, and one of them located in Australia; compared
to the ET+2CE network that we have studied, the latter network therefore has a
greater horizon to coalescing binaries, because of the 40+40 km arms, and
better sky localization, because one of the CE detectors is set in Australia.
The work in Borhanian & Sathyaprakash (2022) studies a larger set of
combinations of 2G+ and 3G configurations, somewhat more tuned toward CE. ET
alone is not considered, but is studied in a so–called ‘ECS’ network,
consisting of ET and 2CE, both having $40\text{\,}\mathrm{km}$ arms, with a CE
located in Australia.363636For performing detailed comparison among our
results, it should also be observed that that the PSDs of CE that we use,
available at https://dcc.cosmicexplorer.org/CE-T2000017/public are an update
of the sensitivity curves available at
https://dcc.cosmicexplorer.org/public/0163/T2000007/005/, that have been used
in Borhanian & Sathyaprakash (2022).

Astrophysical populations. In Borhanian & Sathyaprakash (2022) are given
predictions for BBHs and BNSs, and Ronchini et al. (2022) focus on BNSs, while
we provide predictions for BBHs, BNSs and NSBHs. As we repeatedly stressed in
the text, there are currently large uncertainties on the BBH, BNS and NSBH
populations, and the fiducial choices made in different works can be sensibly
different, resulting potentially in differences as large as one order of
magnitude.

For BNSs, we assume a redshift distribution of the merger rate given by a
rather standard Madau–Dickinson distribution with a delay of 20 Myr between
formation and merger, as described in App. A; the same choice, with the same
time delay, is made in Borhanian & Sathyaprakash (2022), while Ronchini et al.
(2022) use a different model, based on a catalog generated according to the
recent study in Santoliquido et al. (2021a). Furthermore, the overall
normalizations are different. Our fiducial value for the local rate is ${\cal
R}_{0,\rm
BNS}=$105.5\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$, which is
the median value for ${\cal R}_{0,\rm BNS}$ inferred from the GWTC–3 catalog,
using the same flat mass distribution that we are using (Abbott et al.,
2021c). In contrast, Borhanian & Sathyaprakash (2022) use ${\cal R}_{0,\rm
BNS}=$320\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$, which is
the median value inferred from the GWTC–2 catalog (actually, with a flat mass
distribution, while Borhanian & Sathyaprakash (2022) use a truncated Gaussian
distribution), and Ronchini et al. (2022) use ${\cal R}_{0,\rm
BNS}=$365\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ (with the
same flat mass distribution between 1 and $2.5\text{\,}\mathrm{M_{\odot}}$
that we also use). As a result of these differences, the rate of BNS mergers
in our population model is $1\text{\times}{10}^{5}$ events/yr, in Borhanian &
Sathyaprakash (2022) is $4.7\text{\times}{10}^{5}$ events/yr, and in Ronchini
et al. (2022) is $9\text{\times}{10}^{5}$. These numbers therefore span one
full order of magnitude, and are the main reason for the differences among our
results for BNSs. Our estimates are the most conservative among these works,
but we stress again that all these assumptions are consistent with current
knowledge, and, in fact, none of them are on the extreme side. As we
mentioned, values of ${\cal R}_{0,\rm BNS}$ in the range
$(10-1700)\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ are still consistent with
current observations, so one could consider even more conservative, or even
more optimistic models.373737Earlier works, based on the GWTC–1 catalog
(Abbott et al., 2019b), used even larger values of ${\cal R}_{0,\rm BNS}$,
close to ${\cal R}_{0,\rm
BNS}=$1000\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$. For
instance, in Belgacem et al. (2019a) were used the two values ${\cal R}_{0,\rm
BNS}=$662\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ and ${\cal
R}_{0,\rm BNS}=$920\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$
given in Abbott et al. (2019b) for a flat or a Gaussian mass distribution,
respectively, while in Sathyaprakash et al. (2019) the value ${\cal R}_{0,\rm
BNS}=$1000\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ was used.

For BBHs, both us and Borhanian & Sathyaprakash (2022) use a Madau–Dickinson
distribution (with negligible delay between formation and merger for us, and a
delay of 10 Myr for Borhanian & Sathyaprakash (2022)). For the local rate,
current constraints are more stringent than for BNSs. Still, we use ${\cal
R}_{0,{\rm BBH}}=$17\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$,
which is the median value from GWTC–3, while Borhanian & Sathyaprakash (2022)
use ${\cal R}_{0,{\rm
BBH}}=$24\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ from
GWTC–2. As a result, our BBH merger rate is $7.5\text{\times}{10}^{4}$
events/yr, while in Borhanian & Sathyaprakash (2022) is
$1.2\text{\times}{10}^{5}$ events/yr.

Waveforms and detection model. A number of other differences are related to
the waveforms and the detection model assumed. For BNS, we use the
IMRPhenomD_NRTidalv2 waveform, as in Borhanian & Sathyaprakash (2022), which
includes inspiral, merger, ringdown, as well as tidal effects, while Ronchini
et al. (2022) uses TaylorF2, a simpler inspiral–only waveform, and truncates
it at a frequency of the order of $4f_{\rm ISCO}$, where $f_{\rm ISCO}$ is the
frequency at the innermost stable orbit [however, GWFISH can use all waveforms
models of the LALSimulation (LIGO Scientific Collaboration, 2018)]. We all
include the effect of the Earth’s rotation on the localization of the BNSs.
For BBHs both us and Borhanian & Sathyaprakash (2022) use IMRPhenomHM, which
includes higher modes. There is, however, a difference on the number of
parameters on which the inference is performed, both for BBHs and BNSs. In
Borhanian & Sathyaprakash (2022) the inference is performed on the standard
set of 9 parameters $\\{{\cal M}_{c},\ \eta,\ d_{L},\ \theta,\ \phi,\ \iota,\
\psi,\ t_{c},\ \Phi_{c}\\}$, while we add to this list also the two parallel
spin components $\chi_{1,z}$ and $\chi_{2,z}$ and, for BNS and NSBH systems,
we also perform the inference on the tidal deformability parameters, see Sect.
5.

Another difference concerns the threshold for the SNR. A lower threshold on
the SNR increases the number of detections, but (beside raising the
false–alarm rate in an actual run) at the level of our analysis it makes less
reliable the parameter estimation made with the Fisher matrix. In this work,
we present the bulk of our results with a detection threshold ${\rm SNR}\geq
12$, while Borhanian & Sathyaprakash (2022) use ${\rm SNR}\geq 10$; in
Ronchini et al. (2022), which is concerned with multi–messenger observation of
BNS, ${\rm SNR}\geq 8$ is used for the detection in survey mode (which can be
more appropriate to multi–messenger observations, given the temporal
coincidence with the electromagnetic signal), but it has been checked that in
pointing mode, where the sky localization is used, most of the events have
${\rm SNR}\geq 12$.

Another significant point, especially for detector networks, is the assumption
on the duty cycle. In this work, as in Ronchini et al. (2022), for 3G
detectors we have assumed an uncorrelated duty cycle of $85\%$ for each arm,
while in Borhanian & Sathyaprakash (2022) a duty cycle of $100\%$ is assumed.
As we showed in Fig. 7, including a realistic duty cycle results in long tails
of events with poor accuracy on $\Delta d_{L}/d_{L}$ and, to lesser extent, of
$\Delta\Omega_{90\%}$, which are missed assuming a $100\%$ duty cycle. The
same remark applies to the results presented for BBHs in Pieroni et al.
(2022), which agree with ours once we also assume a $100\%$ duty cycle.

As an example of the role of these different assumptions, we perform a more
detailed comparison with the results reported in Ronchini et al. (2022) for
the angular resolution of BNS at ET alone, shown in their Tab. 6, and in the
cumulative distribution functions of their Figs. A1 and A2. In particular,
they find that, at ET alone, one would detect $370\,{\rm events/yr}$ with
$\Delta\Omega_{90\%}<$100\text{\,}{\mathrm{deg}}^{2}$$, and 2 events/yr with
$\Delta\Omega_{90\%}<$1\text{\,}{\mathrm{deg}}^{2}$$, while, from the results
in Sect. 5.2, we find $51\,{\rm events/yr}$ with
$\Delta\Omega_{90\%}<$100\text{\,}{\mathrm{deg}}^{2}$$, and $2\,{\rm
events/yr}$ event below $10\text{\,}{\mathrm{deg}}^{2}$. As a test, to isolate
the effect of the choice of population, we have performed further runs using
${\cal R}_{0,\rm
BNS}=$365\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$, ${\rm
SNR}\geq 8$ and the TaylorF2 waveform,383838However keeping the cut frequency
of the waveform to $2f_{\rm ISCO}$, rather than using $4f_{\rm ISCO}$ as in
Ronchini et al. (2022). so that the only difference is now in the model for
the redshift evolution of the population. Under these conditions, in our model
the population consists of $3.5\text{\times}{10}^{5}$ BNS mergers per year, to
be compared with $9.0\text{\times}{10}^{5}$ in Ronchini et al. (2022), which
is higher by a factor $\simeq 2.6$. Out of these mergers, we find that ET
alone would detect $8.6\text{\times}{10}^{4}$ events/yr, to be compared with
$1.4\text{\times}{10}^{5}$ events/yr in Ronchini et al. (2022), which is now
larger only by a factor $\simeq 1.7$. This is due to the fact that the excess
in the population used in Ronchini et al. (2022), compared to ours, is
especially significant at large $z$ where, however, BNSs are more difficult to
detect. Then, selecting the events with
$\Delta\Omega_{90\%}<$1\text{\,}{\mathrm{deg}}^{2}$$,
$10\text{\,}{\mathrm{deg}}^{2}$ and $100\text{\,}{\mathrm{deg}}^{2}$ we find,
respectively, $(1,\,16,\,298)$ events/yr, to be compared with $(2,\,10,\,370)$
in Ronchini et al. (2022). Given the remaining difference in the redshift
distribution of the population, we can then conclude that the results are
fully consistent, once we make the same assumptions on local rate, and we use
the SNR threshold and the same waveform. This underlines that, within current
uncertainties of the rate and on the population model, and different choices
for technical aspects such as waveform, SNR, etc, the results for BNSs can
still vary sensibly. A comparison with the results reported in Borhanian &
Sathyaprakash (2022) for their ECS configuration with our results for ET+2CE
also show that our results broadly consistent, once taken into account the
different rates, SNR threshold and network configuration.393939In the same
sense, our results are also broadly consistent with previous results presented
in Chan et al. (2018); Maggiore et al. (2020); Grimm & Harms (2020); Nitz &
Dal Canton (2021). They also agree with those presented in Pieroni et al.
(2022) (which are limited to BBHs and set a threshold ${\rm SNR}\geq 20$);
indeed a number of cross–checks between our codes were performed prior to
publication of Pieroni et al. (2022).

Technical differences. Among the various technical aspects of the code, one
worth mentioning here is the difference in the treatment of ill–conditioned
Fisher matrices. In the present work these are discarded, on the basis of a
criterion on the inversion error discussed in Sect. 2.2 and App. C. A similar
strategy is adopted in GWBENCH, where the criterion is rather based on the
condition number, i.e. the ratio of the largest to the smallest eigenvalue. In
contrast, in GWFISH, ill–conditioned Fisher matrices are conditioned and
inverted (i.e., after passing to a diagonal basis, the lines and columns
corresponding to the eigenvector with almost vanishing eigenvalue are deleted,
and the remaining reduced–order matrix is inverted) (Harms et al., 2022). This
leads to treating a larger number of events with the Fisher matrix technique;
however, the vanishing, or nearly vanishing, of an eigenvalue means that the
linear approximation on which the Fisher matrix technique is based becomes
invalid, and one must worry about the effect of non–linear contributions; the
basis in parameter space that diagonalizes the likelihood at quadratic order,
in general, will no longer diagonalize it when cubic and higher–order terms
are included; so the parameter that, at quadratic order, has vanishing
eigenvalue and is discarded, beyond quadratic order could have non–trivial
correlations with all other parameters, possibly leading to larger
marginalized errors on them, compared to the answer obtained from the
quadratic approximation, i.e. from the Fisher matrix. For this reason, we have
chosen to discard events with ill–conditioned Fisher matrices. In App. C we
show a comparison of different techniques. Note that, however, GWFAST offers
different choices, including the conditioning of the singular values, and the
final choice is left to the user (see Iacovelli et al. (2022) for a detailed
description).

As we see from Tab. 3, for BBHs the number of events that we have discarded
because the Fisher matrix was ill–conditioned is about $4\%$ of the total
number of detected events at ET, and just $0.3\%$ for ET+2CE. Therefore, at
the population level, the effect is quite marginal. For BNSs we see from Tab.
4 that the number of events discarded is slightly higher but still consistent,
$4.5\%$ at ET and $1.3\%$ at ET+2CE. The effect of different treatments of the
problem could then in principle be relevant. On the other hand, discarding
ill–conditioned matrices is not a perfect solution either, since it can
potentially introduce some bias; for instance, events with $\cos\iota$ close
to one are preferentially discarded, since in this limit the derivative with
respect to $\iota$ vanishes. Basically, we are facing here an intrinsic
limitations of the Fisher matrix technique, and it is reassuring that codes
employing different approaches reach consistent conclusions. At present, the
uncertainty in the results is completely dominated by the astrophysical
uncertainties.

##  Appendix C Comparison among different ways to quantify uncertainty from
the Fisher matrix

Parameter | Prior range
---|---
$\mathcal{M}_{c},\ d_{L}$ |  $(0,\,+\infty)$
$\eta$ |  $(0,\,0.25]$
$\phi,\ \Phi_{c}$ |  $[0,\,2\pi]$
$\theta,\ \iota,\ \psi,$ |  $[0,\,\pi]$
$t_{c}$ |  $[0,\,1]$
$\chi_{i,z}$ $(i=\\{1,2\\})$ |  $[-1,\,1]$
$\Lambda_{i}$ $(i=\\{1,2\\})$ |  $[0,\,+\infty)$
Table 6: Summary of the prior ranges used when sampling from the posterior
probability.

In this appendix we provide a comparison among different methods of extracting
the forecasted uncertainty on the waveform parameters from the Fisher matrix,
which motivates our choice for this paper as discussed in Sect. 2.2. As a test
case we use the sample of BNSs, which is generally the most prone to showing
large condition numbers. In particular, discarding matrices with condition
number higher than the inverse machine precision ($10^{-15}$ in our case)
would lead to having only a few events left, which does not allow a
statistical study of the population unless the even more stringent and less
realistic assumption of discarding some parameters (e.g. neglecting the
presence of spins and tidal deformabilities) is made. The options that we
consider are (see Iacovelli et al. (2022) for a detailed description of the
methods and their implementation in GWFAST):

  1. a)

we invert the Fisher matrix by a Cholesky decomposition, compute the inversion
error as defined in Eq. (12), and exclude events with inversion error less
than a chosen threshold $\epsilon_{\rm max}$. In particular, we compare the
choices $\epsilon_{\rm max}=${10}^{-3}$$ and $\epsilon_{\rm
max}=$5\text{\times}{10}^{-2}$$;

  2. b)

we invert all Fisher matrices by using the singular–value decomposition (SVD),
and exclude from the inversion the singular values below a threshold of
$10^{-10}$, to ensure that the numerical error does not propagate to
well–measured parameters (Harms et al., 2022). Note that the inversion error
is not defined in this case since the singular values of the matrix are
regularized before inversion;

  3. c)

We impose physically–motivated priors on top of the likelihood in Eq. (11),
draw samples from the corresponding posterior, and quantify the resulting
forecasted uncertainty from the samples. This procedure is detailed below and
does not require any inversion of the Fisher matrix, hence it is not sensitive
to numerical instabilities and encodes the correct boundaries for the waveform
parameters. On the other hand, this technique is computationally too expensive
to be used for all catalogs; we only resort to it in this case for a
validation of the forecasts.

Let us describe the method used for the third option. The posterior
probability for source parameters ${\bf\it\theta}$ given data $s$ is defined
as
$p({\bf\it\theta}\,|\,s)\propto\pi({{\bf\it\theta}})\,\mathcal{L}(s\,|\,{\bf\it\theta})$
where the likelihood $\mathcal{L}(s\,|\,{\bf\it\theta})$ is defined in Eq.
(11) and $\pi({{\bf\it\theta}})$ is the chosen prior. We use flat priors with
the boundaries summarized in Tab. 6 and also fix
$\langle{\delta\theta}^{i}\rangle=0$, i.e. we assume that the posterior mean
coincides with the true value, that we denote by ${\bf\it\mu}$. This is the
same assumption made when quantifying the uncertainty directly from the Fisher
matrix. Written in the form of Eq. (11), the evaluation of the likelihood does
not require the explicit inversion of the Fisher Matrix, and one can extract
samples from this distribution. For each event in the simulated catalogue of
detections, we draw samples from the corresponding posterior by the following
procedure:

  1. 1)

for variables whose $1\sigma$ uncertainty as forecasted by the Fisher analysis
is larger than the prior range, we add a Gaussian prior with standard
deviation equal to the prior range; this has only the purpose of improving the
efficiency and numerical stability, and such prior is undone during the
sampling procedure;

  2. 2)

We sample from the likelihood Eq. (11) as follows. First, we generate a vector
${\bf\it z}=(z_{1},\ ...,\ z_{D})$ of independent, standard normal variables
(with $D$ being the dimensionality of the Fisher matrix). Then, we find the
Cholesky decomposition ${\bf\it C}$ of the Fisher matrix ${\bf\it\Gamma}$,
such that ${\bf\it C}{\bf\it C}^{T}={\bf\it\Gamma}$. Note that this
decomposition might not be found if the Fisher matrix is ill–conditioned,
hence the need of the Gaussian prior discussed in the previous step, which
ensures regularization. Then, the desired samples are given by
${\bf\it\mu}+{\bf\it X}$ where ${\bf\it X}$ is the solution of the system
${\bf\it C}^{T}{\bf\it X}={\bf\it z}$, and ${\bf\it\mu}$ the vector of the
true values;

  3. 3)

We use rejection sampling to discard samples outside the prior range, as well
as to undo the Gaussian prior previously imposed to regularize the Fisher
matrix.

The uncertainty is finally quantified by computing the standard deviation of
the samples. The only caveat concerns the use of Eq. (28) to determine the sky
localization, since its derivation assumes a multivariate Gaussian
distribution, while the use of priors can spoil this assumption. Hence we
compute the $90\%$ sky localization by projecting the samples on the sphere
with a HEALPix pixelization (Zonca et al.,
2019),404040https://healpix.sourceforge.io. locating the pixels within the
$90\%$ localization region, and computing the corresponding sky area. This
procedure is exact, but can require a very extensive and expensive sampling
for events whose localization is not well determined. A second caveat is that
the number of pixels used must be adapted on an event–by–event basis. In order
to reduce the computational cost, we only extract $2000$ samples/event, which
can make this procedure underestimate the sky localization for poorly
localized sources. However, this is enough for an accurate computation for
well–localized events, in which case this provides an important numerical
confirmation of the validity of Eq. (28). As an alternative check, we compute
the sky localization by binning the samples in $(\theta,\phi)$, finding the
bins within the $90\%$ localization region, and computing the corresponding
area on the sphere. This procedure does not require a large number of samples,
but it is sensitive to the number of bins used, which is a quantity that has
to be fixed _a priori_ , and whose mis–specification can lead to wrong
estimates. In summary, there is no fast and $100\%$ accurate way to compute
the localization region, but comparing the techniques listed above we will
show that the different techniques agree very well for events with
sufficiently small ($\lesssim 10^{3}\ \rm deg^{2}$) localization region.

Fig. 25 and Fig. 26 summarize our results. We show the cumulative distribution
of the forecasted $1\sigma$ errors on the detector–frame chirp mass, symmetric
mass ratio, luminosity distance, $90\%$ sky localization and spins along the
$z$ axis, which are the most interesting for the discussion. Fig. 25 shows the
results for ET alone, while Fig. 26 for ET and two CE. From these figures we
can draw a number of interesting conclusions.

First, we see that for the left tails of the distributions, corresponding to
well–resolved events, the result obtained with the SVD inversion and
regularization (green curve) and the one obtained setting a threshold on the
inversion error of $\epsilon_{\rm max}=$5\text{\times}{10}^{-2}$$ (orange
curve) agree almost perfectly for all variables except the spin, in which case
the regularization of the singular values leads to severely underestimating
the error. We conclude that degeneracies among the spin parameters are those
that mostly affect the condition number and lead to the presence of very small
singular values, which are consequently regularized. Hence, the regularization
procedure should not be adopted if one is interested in forecasting errors on
spin parameters. Secondly we note that, adopting a threshold on the inversion
error which is too stringent, such as $\epsilon_{\rm max}=${10}^{-3}$$ (blue
curves), can lead to missing the tail of the best resolved events for some
parameters, in particular the chirp mass. On the contrary, the choice
$\epsilon_{\rm max}=$5\text{\times}{10}^{-2}$$ ensures that all the lower
tails of the distributions match perfectly both the distribution obtained with
SVD and the one obtained by sampling (red curves). This motivates the choice
of this threshold in this paper.

Finally, we comment on the errors obtained with the sampling procedure (red
curves). As expected, the left tail of the distributions, corresponding to
very–well measured parameters, match the results obtained by the inversion of
the Fisher matrix. In general, the use of priors affects differently different
parameters and different configurations. For the mass parameters, we note a
general improvement, which can be explained by the fact that most of the
detected events lie very close to the threshold value $\eta=0.25$ for the
symmetric mass ratio. This leads the finite prior range $\eta\leq 0.25$ to
play a significant role in reducing the forecasted uncertainty; in turn, it
also affects the constraint on the chirp mass, due to the correlation of this
parameter with $\eta$. Similarly, the relative error on the luminosity
distance is reduced at ET alone by the use of a prior because this reduces the
allowed range for the inclination angle $\iota$, whose degeneracy with $d_{L}$
contributes significantly to the error on the latter. For ET in combination
with two CE, this effects is mitigated because $\iota$ is generally better
constrained without the need of a prior. As for the spins, these are generally
poorly constrained and subject to strong degeneracies among the two
components; this explains the improvement given by the prior. We conclude with
the sky localization. The red curve corresponds to the results obtained by the
binning procedure described above, while we additionally show in violet the
results obtained by the HEALPix projection. As anticipated, both methods match
well the forecast obtained by the Fisher matrix alone for
$\Delta\Omega_{90\%}\lesssim 10^{3}\ \rm deg^{2}$, providing a confirmation of
Eq. (28). The use of priors has a larger effect on the forecasts for ET alone,
because for this configuration a precise localization is more difficult. In
this sense, the Fisher forecast provides a conservative estimation of the
localization capabilities. The method based on HEALPix (violet lines) shows
some disagreement for localization regions larger than $\sim 10^{3}\ \rm
deg^{2}$, due to the limited number of samples available, which do not
sufficiently cover all the pixels. Drawing more samples would lead to better
convergence also of the upper tail. This results confirm that the use of of
Eq. (28) is reliable for all events with sufficiently small localization
regions.

Figure 25: Comparison of cumulative distributions of the errors on the
detector–frame chirp mass, symmetric mass ratio, luminosity distance, $90\%$
sky localization and spin components along the $z$ axis, for BNS observed at
ET, with different choices to quantify the uncertainty. The default choice for
this paper, $\epsilon<0.05$, corresponds to the orange curve.  Figure 26: As
in Fig. 25, for a network of ET and two CE.

##  Appendix D Comparison with realistic neutron star equations of state

The equation of state (EoS) of dense matter above the nuclear density is still
unknown, thus we cannot yet reliably relate the mass of a neutron star to its
adimensional tidal deformability. In the main text we made the agnostic choice
of sampling the tidal deformabilities of the simulated BNS systems from a
uniform distribution with broad limits, while here we explore the impact of
assuming a common equation of state for all BNS systems. In particular, we
consider the so–called APR4 (Akmal et al., 1998) and ALF2 (Alford et al.,
2005) EoS, as representative of soft and hard EoS, respectively. The former
thus predicts more compact stars with respect to the latter, leading to
smaller values of the adimensional tidal deformability parameter for the same
mass, which seems to be favoured by the observation of the BNS event GW170817
(Abbott et al., 2017a). The results, for the ET + 2CE network only, are shown
in Fig. 27, where we compare sources having the same parameters apart from the
tidal deformability, which in the case of the green curve is sampled uniformly
in the interval $[0,\,2000]$, while for the others is computed from the mass
of the two components assuming one of the chosen EoS.414141The computation of
$\Lambda(M)$ for the two EoS has been performed using the
Tolman–Oppenheimer–Volkoff equations solver implemented in LALSuite (LIGO
Scientific Collaboration, 2018). Notice that the uniform distribution used in
the main text for the BNS masses, ranging up to
$2.5\text{\,}\mathrm{M_{\odot}}$, is not suitable for performing this test,
given that it exceeds the maximum masses of the chosen EoS.424242For the APR4
EoS the maximum mass is $2.21\text{\,}\mathrm{M_{\odot}}$, while for ALF2 it
is $2.09\text{\,}\mathrm{M_{\odot}}$. To perform this analysis we thus adopted
two independent Gaussian distribution for the masses of the components, ${\cal
N}(1.33,0.09)$ in units of $\mathrm{M_{\odot}}$. From the plot it is apparent
that different choices for the tidal deformability distribution have
negligible impact on the distributions of the parameter errors apart from the
tidal deformability itself. In particular, we find the stiff EoS to produce
systems whose relative error on the adimensional tidal deformability parameter
is smaller, as compared to the soft. This is in line with the fact that a
stiffer EoS predicts higher values for the adimensional tidal deformability
parameter, so the effect on the waveform is larger, which makes the estimation
easier. This also provides a consistency check of our code. We thus conclude
that our results for BNS systems are robust with respect to different choices
for the adimensional tidal deformability distribution.

Figure 27: Histograms of the parameter errors for BNSs with the ET+2CE network
configuration, using the waveform model IMRPhenomD_NRTidalv2, and sampling the
tidal deformability uniformly in the interval $[0,\,2000]$ (red line), or
assuming for all BNSs a common equation of state, chosen to be APR4 (violet)
or ALF2 (orange).

##  Appendix E Comparison among different waveform models

We report here the comparison of the results for the SNR and parameter
estimation errors for the three classes of sources, with waveform models
different from the ones adopted in the main text. Again comparison is only
made using the ET + 2CE detector network.

For BNS and NSBH systems, we compare the results obtained using the full
inspiral–merger–ringdown models IMRPhenomD_NRTidalv2 and IMRPhenomNSBH,
respectively, to the inspiral–only restricted PN model TaylorF2 with terms up
to 3.5 order (Buonanno et al., 2009b; Ajith, 2011; Mishra et al., 2016), also
including in it tidal effects (which enter at order 5 and 6 in the PN
expansion) (Wade et al., 2014), and a cut at twice the Innermost Stable
Circular Orbit frequency, $f_{\rm ISCO}$. As it is apparent from Fig. 28 and
29, the change in the distribution of SNRs and extrinsic parameters errors,
i.e. the sky location, luminosity distance, inclination angle and polarisation
angle, is slight, as expected, since these parameters are not linked to the
waveform model itself in absence of higher modes, which are not expected to
give strong contributions for these classes of sources. The situation is
different for intrinsic parameters, i.e. the chirp mass, symmetric mass ratio,
spins and tidal deformability, whose impact on the GW signal is obviously
better captured by the full inspiral–merger–ringdown models, resulting in
significant improvements in the estimation.

Figure 28: Histograms of the parameter errors for BNSs with the ET+2CE network
configuration, using the full inspiral–merger–ringdown waveform model
IMRPhenomD_NRTidalv2 and the 3.5 restricted PN model TaylorF2 with tidal
effects. Figure 29: Histograms of the parameter errors for NSBH with the
ET+2CE network configuration, using the full inspiral–merger–ringdown waveform
model IMRPhenomNSBH and the 3.5 restricted PN model TaylorF2 with tidal
effects.

## References

  * Abbott et al. (2017a) Abbott, B., et al. 2017a, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Abbott et al. (2019a) —. 2019a, Phys. Rev. X, 9, 011001, doi: 10.1103/PhysRevX.9.011001
  * Abbott et al. (2019b) —. 2019b, Phys. Rev. X, 9, 031040, doi: 10.1103/PhysRevX.9.031040
  * Abbott et al. (2016) Abbott, B. P., et al. 2016, Phys. Rev. Lett., 116, 061102, doi: 10.1103/PhysRevLett.116.061102
  * Abbott et al. (2017b) —. 2017b, Astrophys. J. Lett., 848, L13, doi: 10.3847/2041-8213/aa920c
  * Abbott et al. (2019c) —. 2019c, Astrophys. J. Lett., 882, L24, doi: 10.3847/2041-8213/ab3800
  * Abbott et al. (2020a) —. 2020a, The Astrophysical Journal Letters, 892, L3, doi: 10.3847/2041-8213/ab75f5
  * Abbott et al. (2020b) —. 2020b, Living Rev. Rel., 23, 3, doi: 10.1007/s41114-020-00026-9
  * Abbott et al. (2020c) Abbott, R., et al. 2020c, Phys. Rev. Lett., 125, 101102, doi: 10.1103/PhysRevLett.125.101102
  * Abbott et al. (2020d) —. 2020d, Astrophys. J. Lett., 900, L13, doi: 10.3847/2041-8213/aba493
  * Abbott et al. (2021a) —. 2021a, Phys. Rev. X, 11, 021053, doi: 10.1103/PhysRevX.11.021053
  * Abbott et al. (2021b) —. 2021b.  https://arxiv.org/abs/2111.03606
  * Abbott et al. (2021c) —. 2021c, The population of merging compact binaries inferred using gravitational waves through GWTC-3.  https://arxiv.org/abs/2111.03634
  * Abbott et al. (2021d) —. 2021d.  https://arxiv.org/abs/2112.06861
  * Abbott et al. (2021e) —. 2021e.  https://arxiv.org/abs/2111.03604
  * Abbott et al. (2021f) —. 2021f, Astrophys. J. Lett., 913, L7, doi: 10.3847/2041-8213/abe949
  * Abbott et al. (2021g) —. 2021g, The Astrophysical Journal Letters, 915, L5, doi: 10.3847/2041-8213/ac082e
  * Ajith (2011) Ajith, P. 2011, Phys. Rev. D, 84, 084037, doi: 10.1103/PhysRevD.84.084037
  * Ajith et al. (2006) Ajith, P., Hewitson, M., & Heng, I. S. 2006, Class. Quant. Grav., 23, S741, doi: 10.1088/0264-9381/23/19/S13
  * Akmal et al. (1998) Akmal, A., Pandharipande, V. R., & Ravenhall, D. G. 1998, Phys. Rev. C, 58, 1804, doi: 10.1103/PhysRevC.58.1804
  * Alford et al. (2005) Alford, M., Braby, M., Paris, M., & Reddy, S. 2005, The Astrophysical Journal, 629, 969, doi: 10.1086/430902
  * Antonelli et al. (2021) Antonelli, A., Burke, O., & Gair, J. R. 2021, Mon. Not. Roy. Astron. Soc., 507, 5069, doi: 10.1093/mnras/stab2358
  * Barack & Cutler (2004) Barack, L., & Cutler, C. 2004, Phys. Rev. D, 69, 082005, doi: 10.1103/PhysRevD.69.082005
  * Bavera et al. (2022) Bavera, S. S., Fishbach, M., Zevin, M., Zapartas, E., & Fragos, T. 2022.  https://arxiv.org/abs/2204.02619
  * Belczynski et al. (2022) Belczynski, K., Doctor, Z., Zevin, M., et al. 2022.  https://arxiv.org/abs/2204.11730
  * Belgacem et al. (2019a) Belgacem, E., Dirian, Y., Foffa, S., et al. 2019a, JCAP, 1908, 015, doi: 10.1088/1475-7516/2019/08/015
  * Belgacem et al. (2018a) Belgacem, E., Dirian, Y., Foffa, S., & Maggiore, M. 2018a, Phys. Rev. D, 97, 104066, doi: 10.1103/PhysRevD.97.104066
  * Belgacem et al. (2018b) —. 2018b, Phys. Rev. D, 98, 023510, doi: 10.1103/PhysRevD.98.023510
  * Belgacem et al. (2019b) Belgacem, E., et al. 2019b, JCAP, 1907, 024, doi: 10.1088/1475-7516/2019/07/024
  * Berti et al. (2005) Berti, E., Buonanno, A., & Will, C. M. 2005, Phys. Rev. D, 71, 084025, doi: 10.1103/PhysRevD.71.084025
  * Berti et al. (2016) Berti, E., Sesana, A., Barausse, E., Cardoso, V., & Belczynski, K. 2016, Phys. Rev. Lett., 117, 101102, doi: 10.1103/PhysRevLett.117.101102
  * Biscoveanu et al. (2022) Biscoveanu, S., Callister, T. A., Haster, C.-J., et al. 2022, Astrophys. J. Lett., 932, L19, doi: 10.3847/2041-8213/ac71a8
  * Borhanian (2021) Borhanian, S. 2021, Class. Quant. Grav., 38, 175014, doi: 10.1088/1361-6382/ac1618
  * Borhanian et al. (2020) Borhanian, S., Dhani, A., Gupta, A., Arun, K., & Sathyaprakash, B. 2020, Astrophys. J. Lett., 905, L28, doi: 10.3847/2041-8213/abcaf5
  * Borhanian & Sathyaprakash (2022) Borhanian, S., & Sathyaprakash, B. S. 2022, Listening to the Universe with Next Generation Ground-Based Gravitational-Wave Detectors.  https://arxiv.org/abs/2202.11048
  * Bradbury et al. (2018) Bradbury, J., Frostig, R., Hawkins, P., et al. 2018, JAX: composable transformations of Python+NumPy programs, 0.2.5.  http://github.com/google/jax
  * Broekgaarden & Berger (2021) Broekgaarden, F. S., & Berger, E. 2021, Astrophys. J. Lett., 920, L13, doi: 10.3847/2041-8213/ac2832
  * Broekgaarden et al. (2021) Broekgaarden, F. S., Berger, E., Neijssel, C. J., et al. 2021, Mon. Not. Roy. Astron. Soc., 508, 5028, doi: 10.1093/mnras/stab2716
  * Buonanno et al. (2009a) Buonanno, A., Iyer, B., Ochsner, E., Pan, Y., & Sathyaprakash, B. S. 2009a, Phys. Rev. D, 80, 084043, doi: 10.1103/PhysRevD.80.084043
  * Buonanno et al. (2009b) Buonanno, A., Iyer, B. R., Ochsner, E., Pan, Y., & Sathyaprakash, B. S. 2009b, Phys. Rev. D, 80, 084043, doi: 10.1103/PhysRevD.80.084043
  * Cardoso et al. (2016) Cardoso, V., Hopper, S., Macedo, C. F. B., Palenzuela, C., & Pani, P. 2016, Phys. Rev. D, 94, 084031, doi: 10.1103/PhysRevD.94.084031
  * Chan et al. (2018) Chan, M. L., Messenger, C., Heng, I. S., & Hendry, M. 2018, Phys. Rev. D, 97, 123014, doi: 10.1103/PhysRevD.97.123014
  * Chase et al. (2022) Chase, E. A., O’Connor, B., Fryer, C. L., et al. 2022, Astrophys. J., 927, 163, doi: 10.3847/1538-4357/ac3d25
  * Chassande-Mottin et al. (2019) Chassande-Mottin, E., Leyde, K., Mastrogiovanni, S., & Steer, D. A. 2019, Phys. Rev. D, 100, 083514, doi: 10.1103/PhysRevD.100.083514
  * Chatterji et al. (2006) Chatterji, S., Lazzarini, A., Stein, L., et al. 2006, Phys. Rev. D, 74, 082005, doi: 10.1103/PhysRevD.74.082005
  * Chruślińska (2022) Chruślińska, M. 2022.  https://arxiv.org/abs/2206.10622
  * Cornish & Larson (2003) Cornish, N. J., & Larson, S. L. 2003, Phys. Rev. D, 67, 103001, doi: 10.1103/PhysRevD.67.103001
  * Cowperthwaite et al. (2019) Cowperthwaite, P. S., Villar, V. A., Scolnic, D. M., & Berger, E. 2019, Astrophys. J., 874, 88, doi: 10.3847/1538-4357/ab07b6
  * Cutler (1998) Cutler, C. 1998, Phys. Rev. D, 57, 7089, doi: 10.1103/PhysRevD.57.7089
  * Cutler & Flanagan (1994) Cutler, C., & Flanagan, E. E. 1994, Phys. Rev. D, 49, 2658, doi: 10.1103/PhysRevD.49.2658
  * De Luca et al. (2020) De Luca, V., Franciolini, G., Pani, P., & Riotto, A. 2020, JCAP, 06, 044, doi: 10.1088/1475-7516/2020/06/044
  * De Luca et al. (2021a) —. 2021a, JCAP, 05, 003, doi: 10.1088/1475-7516/2021/05/003
  * De Luca et al. (2021b) —. 2021b, JCAP, 11, 039, doi: 10.1088/1475-7516/2021/11/039
  * Dietrich et al. (2019) Dietrich, T., Samajdar, A., Khan, S., et al. 2019, Phys. Rev. D, 100, 044003, doi: 10.1103/PhysRevD.100.044003
  * Dominik et al. (2013) Dominik, M., Belczynski, K., Fryer, C., et al. 2013, Astrophys. J., 779, 72, doi: 10.1088/0004-637X/779/1/72
  * Evans et al. (2021) Evans, M., et al. 2021.  https://arxiv.org/abs/2109.09882
  * Ezquiaga (2021) Ezquiaga, J. M. 2021, Phys. Lett. B, 822, 136665, doi: 10.1016/j.physletb.2021.136665
  * Ezquiaga & Holz (2022) Ezquiaga, J. M., & Holz, D. E. 2022.  https://arxiv.org/abs/2202.08240
  * Farr et al. (2019) Farr, W. M., Fishbach, M., Ye, J., & Holz, D. 2019, Astrophys. J. Lett., 883, L42, doi: 10.3847/2041-8213/ab4284
  * Farrow et al. (2019) Farrow, N., Zhu, X.-J., & Thrane, E. 2019, Astrophys. J., 876, 18, doi: 10.3847/1538-4357/ab12e3
  * Finke et al. (2021) Finke, A., Foffa, S., Iacovelli, F., Maggiore, M., & Mancarella, M. 2021, JCAP, 08, 026, doi: 10.1088/1475-7516/2021/08/026
  * Finke et al. (2022) —. 2022, Phys. Dark Univ., 36, 100994, doi: 10.1016/j.dark.2022.100994
  * Finn (1992) Finn, L. S. 1992, Phys. Rev. D, 46, 5236, doi: 10.1103/PhysRevD.46.5236
  * Fishbach et al. (2021) Fishbach, M., Doctor, Z., Callister, T., et al. 2021.  https://arxiv.org/abs/2101.07699
  * Franciolini (2021) Franciolini, G. 2021, PhD thesis, Geneva U., Dept. Theor. Phys., doi: 10.13097/archive-ouverte/unige:156136
  * Freise et al. (2009) Freise, A., Chelkowski, S., Hild, S., et al. 2009, Classical and Quantum Gravity, 26, 085012, doi: 10.1088/0264-9381/26/8/085012
  * García-Quirós et al. (2020) García-Quirós, C., Colleoni, M., Husa, S., et al. 2020, Phys. Rev. D, 102, 064002, doi: 10.1103/PhysRevD.102.064002
  * Giacobbo & Mapelli (2018) Giacobbo, N., & Mapelli, M. 2018, Monthly Notices of the Royal Astronomical Society, 480, 2011, doi: 10.1093/mnras/sty1999
  * Goldberg et al. (1967) Goldberg, J. N., Macfarlane, A. J., Newman, E. T., Rohrlich, F., & Sudarshan, E. C. G. 1967, Journal of Mathematical Physics, 8, 2155, doi: 10.1063/1.1705135
  * Goncharov et al. (2022) Goncharov, B., Nitz, A. H., & Harms, J. 2022, Utilizing the null stream of Einstein Telescope, arXiv.  https://arxiv.org/abs/2204.08533
  * Graham et al. (2020) Graham, M., et al. 2020, Phys. Rev. Lett., 124, 251102, doi: 10.1103/PhysRevLett.124.251102
  * Grimm & Harms (2020) Grimm, S., & Harms, J. 2020, Phys. Rev. D, 102, 022007, doi: 10.1103/PhysRevD.102.022007
  * Gupta et al. (2022) Gupta, P. K., Puecher, A., Pang, P. T. H., et al. 2022.  https://arxiv.org/abs/2205.01182
  * Gürsel & Tinto (1989) Gürsel, Y., & Tinto, M. 1989, Phys. Rev. D, 40, 3884, doi: 10.1103/PhysRevD.40.3884
  * Hall & Evans (2019) Hall, E. D., & Evans, M. 2019, Class. Quant. Grav., 36, 225002, doi: 10.1088/1361-6382/ab41d6
  * Harms et al. (2022) Harms, J., Dupletsa, U., Banerjee, B., et al. 2022.  https://arxiv.org/abs/2205.02499
  * Harry & Fairhurst (2011) Harry, I. W., & Fairhurst, S. 2011, Phys. Rev. D, 83, 084002, doi: 10.1103/PhysRevD.83.084002
  * Hernandez Vivanco et al. (2019) Hernandez Vivanco, F., Smith, R., Thrane, E., et al. 2019, Phys. Rev. D, 100, 103009, doi: 10.1103/PhysRevD.100.103009
  * Hild et al. (2011) Hild, S., et al. 2011, Class. Quant. Grav., 28, 094013, doi: 10.1088/0264-9381/28/9/094013
  * Himemoto et al. (2021) Himemoto, Y., Nishizawa, A., & Taruya, A. 2021, Phys. Rev. D, 104, 044010, doi: 10.1103/PhysRevD.104.044010
  * Hinderer et al. (2010) Hinderer, T., Lackey, B. D., Lang, R. N., & Read, J. S. 2010, Phys. Rev. D, 81, 123016, doi: 10.1103/PhysRevD.81.123016
  * Husa et al. (2016) Husa, S., Khan, S., Hannam, M., et al. 2016, Phys. Rev. D, 93, 044006, doi: 10.1103/PhysRevD.93.044006
  * Iacovelli (2022) Iacovelli, F. 2022, WF4Py: Gravitational waves waveform models in pure Python language, v1.0.0, Zenodo, doi: 10.5281/zenodo.7060240
  * Iacovelli & Mancarella (2022) Iacovelli, F., & Mancarella, M. 2022, gwfast: a Fisher information matrix Python package for GW detector networks, v1.0.1, Zenodo, doi: 10.5281/zenodo.7060236
  * Iacovelli et al. (2022) Iacovelli, F., Mancarella, M., Foffa, S., & Maggiore, M. 2022, Astrophys. J. Suppl., To appear.  https://arxiv.org/abs/2207.06910
  * Ivezić et al. (2019) Ivezić, v., et al. 2019, Astrophys. J., 873, 111, doi: 10.3847/1538-4357/ab042c
  * Jaranowski et al. (1998) Jaranowski, P., Królak, A., & Schutz, B. F. 1998, Phys. Rev. D, 58, 063001, doi: 10.1103/PhysRevD.58.063001
  * Jin-Ping et al. (2021) Jin-Ping, Z., Shichao, W., Yuan-Pei, Y., et al. 2021, Astrophys. J., 917, 24, doi: 10.3847/1538-4357/abfe5e
  * Kalaghatgi et al. (2020) Kalaghatgi, C., Hannam, M., & Raymond, V. 2020, Phys. Rev. D, 101, 103004, doi: 10.1103/PhysRevD.101.103004
  * Kalogera et al. (2021) Kalogera, V., et al. 2021.  https://arxiv.org/abs/2111.06990
  * Karathanasis et al. (2022) Karathanasis, C., Mukherjee, S., & Mastrogiovanni, S. 2022.  https://arxiv.org/abs/2204.13495
  * Keelin (2016) Keelin, T. W. 2016, Decision Analysis, 13, 243, doi: 10.1287/deca.2016.0338
  * Khan et al. (2016a) Khan, S., Husa, S., Hannam, M., et al. 2016a, Phys. Rev. D, 93, 044007, doi: 10.1103/PhysRevD.93.044007
  * Khan et al. (2016b) —. 2016b, Phys. Rev. D, 93, 044007, doi: 10.1103/PhysRevD.93.044007
  * Kinugawa et al. (2014) Kinugawa, T., Inayoshi, K., Hotokezaka, K., Nakauchi, D., & Nakamura, T. 2014, Mon. Not. Roy. Astron. Soc., 442, 2963, doi: 10.1093/mnras/stu1022
  * Lattimer & Schramm (1976) Lattimer, J. M., & Schramm, D. N. 1976, ApJ, 210, 549, doi: 10.1086/154860
  * Leyde et al. (2022) Leyde, K., Mastrogiovanni, S., Steer, D. A., Chassande-Mottin, E., & Karathanasis, C. 2022.  https://arxiv.org/abs/2202.00025
  * Li et al. (2022) Li, Y., Heng, I. S., Chan, M. L., Messenger, C., & Fan, X. 2022, Phys. Rev. D, 105, 043010, doi: 10.1103/PhysRevD.105.043010
  * LIGO Scientific Collaboration (2018) LIGO Scientific Collaboration. 2018, LIGO Algorithm Library - LALSuite, free software (GPL), doi: 10.7935/GT1W-FZ16
  * 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
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, Ann. Rev. Astron. Astrophys., 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Madau & Fragos (2017) Madau, P., & Fragos, T. 2017, Astrophys. J., 840, 39, doi: 10.3847/1538-4357/aa6af9
  * Maggiore (2007) Maggiore, M. 2007, Gravitational Waves. Vol. 1: Theory and Experiments, Oxford Master Series in Physics (Oxford University Press)
  * Maggiore et al. (2020) Maggiore, M., et al. 2020, JCAP, 2003, 050, doi: 10.1088/1475-7516/2020/03/050
  * Mancarella et al. (2022) Mancarella, M., Genoud-Prachex, E., & Maggiore, M. 2022, Phys. Rev. D, 105, 064030, doi: 10.1103/PhysRevD.105.064030
  * Mandel & Broekgaarden (2022) Mandel, I., & Broekgaarden, F. S. 2022, Living Rev. Rel., 25, 1, doi: 10.1007/s41114-021-00034-3
  * Mandel & Fragos (2020) Mandel, I., & Fragos, T. 2020, Astrophys. J. Lett., 895, L28, doi: 10.3847/2041-8213/ab8e41
  * Mandel & Smith (2021) Mandel, I., & Smith, R. J. E. 2021, Astrophys. J. Lett., 922, L14, doi: 10.3847/2041-8213/ac35dd
  * Margossian (2018) Margossian, C. C. 2018, CoRR, abs/1811.05031
  * Mastrogiovanni et al. (2021a) Mastrogiovanni, S., Haegel, L., Karathanasis, C., Magana-Hernandez, I., & Steer, D. 2021a, JCAP, 02, 043, doi: 10.1088/1475-7516/2021/02/043
  * Mastrogiovanni et al. (2021b) Mastrogiovanni, S., Leyde, K., Karathanasis, C., et al. 2021b, Phys. Rev. D, 104, 062009, doi: 10.1103/PhysRevD.104.062009
  * McKernan et al. (2019) McKernan, B., Ford, K. E. S., Bartos, I., et al. 2019, Astrophys. J. Lett., 884, L50, doi: 10.3847/2041-8213/ab4886
  * Meacher et al. (2016) Meacher, D., Cannon, K., Hanna, C., Regimbau, T., & Sathyaprakash, B. S. 2016, Phys. Rev. D, 93, 024018, doi: 10.1103/PhysRevD.93.024018
  * Mink & King (2017) Mink, S. E. d., & King, A. 2017, Astrophys. J. Lett., 839, L7, doi: 10.3847/2041-8213/aa67f3
  * Mishra et al. (2016) Mishra, C. K., Kela, A., Arun, K. G., & Faye, G. 2016, Phys. Rev. D, 93, 084054, doi: 10.1103/PhysRevD.93.084054
  * Ng et al. (2022) Ng, K. K. Y., Franciolini, G., Berti, E., et al. 2022.  https://arxiv.org/abs/2204.11864
  * Ng et al. (2021a) Ng, K. K. Y., Vitale, S., Farr, W. M., & Rodriguez, C. L. 2021a, Astrophys. J. Lett., 913, L5, doi: 10.3847/2041-8213/abf8be
  * Ng et al. (2021b) Ng, K. K. Y., Chen, S., Goncharov, B., et al. 2021b.  https://arxiv.org/abs/2108.07276
  * Nissanke et al. (2010) Nissanke, S., Holz, D. E., Hughes, S., Dalal, N., & Sievers, J. L. 2010, Astrophys. J., 725, 496, doi: 10.1088/0004-637X/725/1/496
  * Nitz & Dal Canton (2021) Nitz, A. H., & Dal Canton, T. 2021, Astrophys. J. Lett., 917, L27, doi: 10.3847/2041-8213/ac1a75
  * Palmese et al. (2021a) Palmese, A., Bom, C. R., Mucesh, S., & Hartley, W. G. 2021a.  https://arxiv.org/abs/2111.06445
  * Palmese et al. (2021b) Palmese, A., Fishbach, M., Burke, C. J., Annis, J. T., & Liu, X. 2021b, Astrophys. J. Lett., 914, L34, doi: 10.3847/2041-8213/ac0883
  * Pannarale et al. (2015) Pannarale, F., Berti, E., Kyutoku, K., Lackey, B. D., & Shibata, M. 2015, Phys. Rev. D, 92, 084050, doi: 10.1103/PhysRevD.92.084050
  * Perna et al. (2016) Perna, R., Lazzati, D., & Giacomazzo, B. 2016, Astrophys. J. Lett., 821, L18, doi: 10.3847/2041-8205/821/1/L18
  * Pieroni et al. (2022) Pieroni, M., Ricciardone, A., & Barausse, E. 2022.  https://arxiv.org/abs/2203.12586
  * Pizzati et al. (2022) Pizzati, E., Sachdev, S., Gupta, A., & Sathyaprakash, B. 2022, Phys. Rev. D, 105, 104016, doi: 10.1103/PhysRevD.105.104016
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, The Astronomical Journal, 156, 123
  * Prince et al. (2002) Prince, T. A., Tinto, M., Larson, S. L., & Armstrong, J. W. 2002, Phys. Rev. D, 66, 122002, doi: 10.1103/PhysRevD.66.122002
  * Punturo et al. (2010) Punturo, M., et al. 2010, Class. Quant. Grav., 27, 194002, doi: 10.1088/0264-9381/27/19/194002
  * Qin et al. (2018) Qin, Y., Fragos, T., Meynet, G., et al. 2018, Astron. Astrophys., 616, A28, doi: 10.1051/0004-6361/201832839
  * Raidal et al. (2019) Raidal, M., Spethmann, C., Vaskonen, V., & Veermäe, H. 2019, JCAP, 02, 018, doi: 10.1088/1475-7516/2019/02/018
  * Rakhmanov (2006) Rakhmanov, M. 2006, Class. Quant. Grav., 23, S673, doi: 10.1088/0264-9381/23/19/S05
  * Regimbau et al. (2017) Regimbau, T., Evans, M., Christensen, N., et al. 2017, Phys. Rev. Lett., 118, 151105, doi: 10.1103/PhysRevLett.118.151105
  * Regimbau & Hughes (2009) Regimbau, T., & Hughes, S. A. 2009, Phys. Rev. D, 79, 062002, doi: 10.1103/PhysRevD.79.062002
  * 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., et al. 2019, Bull. Am. Astron. Soc., 51, 035.  https://arxiv.org/abs/1907.04833
  * Rodriguez et al. (2013) Rodriguez, C. L., Farr, B., Farr, W. M., & Mandel, I. 2013, Phys. Rev. D, 88, 084013, doi: 10.1103/PhysRevD.88.084013
  * Ronchini et al. (2022) Ronchini, S., Branchesi, M., Oganesyan, G., et al. 2022.  https://arxiv.org/abs/2204.01746
  * Rossi et al. (2020) Rossi, A., et al. 2020, Mon. Not. Roy. Astron. Soc., 493, 3379, doi: 10.1093/mnras/staa479
  * Rozner & Perets (2022) Rozner, M., & Perets, H. B. 2022, Astrophys. J., 931, 149, doi: 10.3847/1538-4357/ac6d55
  * Sagiv et al. (2014) Sagiv, I., et al. 2014, Astron. J., 147, 79, doi: 10.1088/0004-6256/147/4/79
  * Samajdar et al. (2021) Samajdar, A., Janquart, J., Van Den Broeck, C., & Dietrich, T. 2021, Phys. Rev. D, 104, 044003, doi: 10.1103/PhysRevD.104.044003
  * Santoliquido et al. (2021a) Santoliquido, F., Mapelli, M., Giacobbo, N., Bouffanais, Y., & Artale, M. C. 2021a, Mon. Not. Roy. Astron. Soc., 502, 4877, doi: 10.1093/mnras/stab280
  * Santoliquido et al. (2021b) —. 2021b, Monthly Notices of the Royal Astronomical Society, 502, 4877, doi: 10.1093/mnras/stab280
  * Sathyaprakash et al. (2019) Sathyaprakash, B. S., et al. 2019.  https://arxiv.org/abs/1903.09277
  * Schutz (2011) Schutz, B. F. 2011, Class. Quant. Grav., 28, 125023, doi: 10.1088/0264-9381/28/12/125023
  * Schutz & Sathyaprakash (2020) Schutz, B. F., & Sathyaprakash, B. S. 2020.  https://arxiv.org/abs/2009.10212
  * Singer et al. (2016) Singer, L. P., et al. 2016, Astrophys. J. Lett., 829, L15, doi: 10.3847/2041-8205/829/1/L15
  * Smith et al. (2021) Smith, R., et al. 2021, Phys. Rev. Lett., 127, 081102, doi: 10.1103/PhysRevLett.127.081102
  * Spergel et al. (2015) Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv e-prints, arXiv:1503.03757.  https://arxiv.org/abs/1503.03757
  * Spyromilio et al. (2008) Spyromilio, J., Comerón, F., D’Odorico, S., Kissler-Patig, M., & Gilmozzi, R. 2008, The Messenger, 133, 2
  * Taylor et al. (2012) Taylor, S. R., Gair, J. R., & Mandel, I. 2012, Phys. Rev. D, 85, 023535, doi: 10.1103/PhysRevD.85.023535
  * Usman et al. (2019) Usman, S. A., Mills, J. C., & Fairhurst, S. 2019, Astrophys. J., 877, 82, doi: 10.3847/1538-4357/ab0b3e
  * Vallisneri (2008) Vallisneri, M. 2008, Phys. Rev. D, 77, 042001, doi: 10.1103/PhysRevD.77.042001
  * van Son et al. (2022) van Son, L. A. C., de Mink, S. E., Callister, T., et al. 2022, Astrophys. J., 931, 17, doi: 10.3847/1538-4357/ac64a3
  * Vitale & Chen (2018) Vitale, S., & Chen, H.-Y. 2018, Phys. Rev. Lett., 121, 021303, doi: 10.1103/PhysRevLett.121.021303
  * Wade et al. (2014) Wade, L., Creighton, J. D. E., Ochsner, E., et al. 2014, Phys. Rev. D, 89, 103012, doi: 10.1103/PhysRevD.89.103012
  * Weaver et al. (2022) Weaver, J. R., et al. 2022, Astrophys. J. Supp., 258, 11, doi: 10.3847/1538-4365/ac3078
  * Wen & Chen (2010) Wen, L., & Chen, Y. 2010, Phys. Rev. D, 81, 082001, doi: 10.1103/PhysRevD.81.082001
  * Wen & Schutz (2005) Wen, L., & Schutz, B. F. 2005, Class. Quant. Grav., 22, S1321, doi: 10.1088/0264-9381/22/18/S46
  * Wong & Li (2022) Wong, I. C. F., & Li, T. G. F. 2022, Phys. Rev. D, 105, 084002, doi: 10.1103/PhysRevD.105.084002
  * Wong et al. (2021) Wong, I. C. F., Pang, P. T. H., Lo, R. K. L., Li, T. G. F., & Van Den Broeck, C. 2021.  https://arxiv.org/abs/2105.09485
  * Zhao & Wen (2018) Zhao, W., & Wen, L. 2018, Phys. Rev. D, 97, 064031, doi: 10.1103/PhysRevD.97.064031
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source Software, 4, 1298, doi: 10.21105/joss.01298