# Forecasting the detection capabilities of third–generation
gravitational–wave detectors using GWFAST

Francesco Iacovelli Département de Physique Théorique and Gravitational Wave
Science Center,
Université de Genève, 24 quai Ernest Ansermet, 1211 Genève 4, Switzerland
Michele Mancarella Département de Physique Théorique and Gravitational Wave
Science Center,
Université de Genève, 24 quai Ernest Ansermet, 1211 Genève 4, Switzerland
Stefano Foffa Département de Physique Théorique and Gravitational Wave Science
Center,
Université de Genève, 24 quai Ernest Ansermet, 1211 Genève 4, Switzerland
Michele Maggiore Département de Physique Théorique and Gravitational Wave
Science Center,
Université de Genève, 24 quai Ernest Ansermet, 1211 Genève 4, Switzerland

###### Abstract

We introduce GWFAST, a novel Fisher–matrix code for gravitational–wave
studies, tuned toward third–generation gravitational–wave detectors such as
Einstein Telescope (ET) and Cosmic Explorer (CE). We use it to perform a
comprehensive study of the capabilities of ET alone, and of a network made by
ET and two CE detectors, as well as to provide forecasts for the forthcoming
O4 run of the LVK collaboration. We consider binary neutron stars, binary
black holes and neutron star–black hole binaries, and compute basic metrics
such as the distribution of signal–to–noise ratio (SNR), the accuracy in the
reconstruction of various parameters (including distance, sky localization,
masses, spins and, for neutron stars, tidal deformabilities), and the redshift
distribution of the detections for different thresholds in SNR and different
levels of accuracy in localization and distance measurement. We examine the
expected distribution and properties of ‘golden events’, with especially large
values of the SNR. We also pay special attention to the dependence of the
results on astrophysical uncertainties and on various technical details (such
as choice of waveforms, or the threshold in SNR), and we compare with other
Fisher codes in the literature. In the companion paper Iacovelli et al. (2022)
we discuss the technical aspects of the code. Together with this paper, we
publicly release the code GWFAST 111https://github.com/CosmoStatGW/gwfast, and
the library WF4Py 222https://github.com/CosmoStatGW/WF4Py implementing
state–of–the–art gravitational–wave waveforms in pure Python.

††software: JAX (Bradbury et al., 2018), astropy (Price-Whelan et al., 2018),
LALSimulation (LIGO Scientific Collaboration, 2018)

tablenum SIXtablenum

ET-0141A-22

##  1 Introduction

The discoveries made in the last few years by the LIGO and Virgo gravitational
wave (GW) detectors have opened a new window on the Universe. After the first
historic detection of the GWs from a binary black hole (BBH) coalescence in
September 2015 (Abbott et al., 2016), and the first detection in 2017 of a
binary neutron star (BNS) coalescence and the identification and follow–up of
its electromagnetic counterpart (Abbott et al., 2017a, b), which opened the
field of multi–messenger astronomy, GW observations have become a routine,
with detections made on a weekly basis. After three observing runs (O1–O3),
the current catalog of GW detections contains about 90 events, mostly BBHs,
but include also two BNSs and two neutron star–black hole (NSBH) binaries
(Abbott et al., 2021a, b). These discoveries are already starting to have a
remarkable impact on astrophysics, cosmology and fundamental physics
including, e.g., studies of population properties of compact binary systems
(Abbott et al., 2021c), tests of General Relativity (Abbott et al., 2021d),
the determination of the speed of GWs to one part in $10^{15}$ (Abbott et al.,
2017b), first constraints on the expansion history of the Universe using GWs
(Abbott et al., 2021e), or constraints on the neutron star equation of state
(Hernandez Vivanco et al., 2019).

Currently LIGO and Virgo, now joined by KAGRA and forming the LVK
collaboration, are further improving their sensitivity and preparing for their
fourth observational run, O4, scheduled for 2023, and eventually for the O5
run, currently scheduled for 2026–2028. However, the community has long been
preparing the jump toward ‘third–generation’ (3G) detectors, Einstein
Telescope (ET) in Europe (Punturo et al., 2010; Hild et al., 2011) and Cosmic
Explorer (CE) in the US (Reitze et al., 2019; Evans et al., 2021). Thanks to
an increase by one order of magnitude in sensitivity and a significant
enlargement of the bandwidth toward both low and high frequencies, these
detectors will have an extraordinary potential for discoveries. In particular,
for coalescing binaries, they will allow us to make a huge jump in the
distances reached, in the number of detections, in the accuracy of signal
reconstructions, and in the range of masses that can be explored, providing
data that have the potential of triggering revolutions in fundamental physics,
in cosmology and in astrophysics. Building on many previous works, the science
case for 3G detectors is being systematically investigated (see Maggiore et
al. (2020); Kalogera et al. (2021) and references therein). The activity
around 3G detectors recently received a significant boost from the inclusion
of ET in the ESFRI Roadmap, the roadmap of large European scientific
infrastructures,333See https://www.esfri.eu/latest-esfri-news/new-ris-
roadmap-2021. resulting in an acceleration of activities around ET and, very
recently, the ET Collaboration has formally been created at the XII ET
Symposium.444See https://indico.ego-gw.it/event/411/. On the US side, activity
around Cosmic Explorer has started more recently, and has led to the Cosmic
Explorer Horizon Study document (Evans et al., 2021).

A crucial aspect, in order to assess the scientific potential of 3G detectors,
is the development of codes that allow us to forecast the performance of these
detectors, in terms of some rather general “metrics”, such as the distribution
of signal–to–noise ratio (SNR) for coalescing binaries, and the accuracy in
the reconstruction of their parameters; in particular, distance, sky
localization, masses, spins and, for neutron stars, quantities (such as the
tidal deformability) that are sensitive to their inner structure. It is also
very important to understand the reach in redshift of these detectors and how
the distribution in SNR and parameter reconstruction depends on redshift. The
reference tool for this kind of investigations is an approach based on the
Fisher matrix. Despite some of its well–known limitations, that we will review
in Sect. 2, this is at the moment the only computationally feasible way of
performing parameter estimation on large populations, such as the
$10^{4}-10^{5}$ BBHs and BNSs per year that, as we will discuss below, 3G
detectors are expected to detect. This tool also allows us to compare the
performance of different configurations of a given detector, and of different
detector networks, which is necessary in order to take informed decisions on
individual detector designs and on optimal network configurations.

For this reason, a number of parameter estimation codes tuned to 3G detectors
have been developed recently, in particular GWBENCH (Borhanian, 2021) and
GWFISH (Harms et al., 2022) [see also Chan et al. (2018); Grimm & Harms
(2020); Nitz & Dal Canton (2021); Li et al. (2022); Pieroni et al. (2022)]. In
this paper we present a novel parameter estimation code, GWFAST, also meant
mainly for application to 3G detectors. The existence of several different
parameter estimation codes is a welcome, and in fact necessary, feature; these
codes can contribute to taking decisions on detectors which are meant to
dominate the scientific landscape for decades (and require the huge financial
and human resources typical of Big Science), so they must be extremely
reliable and have undergone cross–checks between different groups, that
developed different codes independently. In this spirit, we have undergone a
process of cross–checking between GWBENCH, GWFISH and GWFAST, which are,
arguably, the most complete and advanced codes currently available, finding
broad consistency.555These checks are being performed in the context of the
activities of the Observational Science Board (OSB) of ET, which is in charge
of developing the Science Cases and the technical tools relevant for ET. See
https://www.et-gw.eu/index.php/observational-science-board for a repository of
papers relevant for ET, produced in the context of the OSB activities. Each of
these codes has different technical implementations. We discuss the technical
aspects of our code, together with reliability and performance tests, in the
companion paper Iacovelli et al. (2022), while here we briefly summarize them
in Sect. 3.3. One relevant aspect is that the problem of computing Fisher
matrices for a large catalog of independent events is clearly parallel.
Parallelization techniques can be used, but one limitation of current
Python–based implementations is that they still usually compute the SNRs and
Fisher matrices serially on each node/CPU. In our implementation, we are able
to ‘vectorize’ the evaluation of the Fisher matrices even on a single CPU (_on
top of_ the use of parallel computing), resulting in a gain in computational
speed, which is at the origin of the name GWFAST. This is also due to the
implementation of the waveforms in Python, which motivated the development and
release of the library WF4Py. Furthermore, this allows us to make use of
automatic differentiation to compute derivatives, which is a technique
alternative to finite differences leading to faster and more robust
evaluations. Note that GWFAST also supports an interface with the LIGO
Algorithm Library, LAL (LIGO Scientific Collaboration, 2018), which allows the
use of all waveforms available in this library.

In this paper we discuss the results obtained for BBHs, BNSs and NSBHs with
our code, and we compare them, in particular, with those reported for BBHs and
BNSs in Borhanian & Sathyaprakash (2022) (using GWBENCH), and for BNSs in
Ronchini et al. (2022) (which uses GWFISH). The comparison also needs to take
care of different assumptions on astrophysical populations, that reflect our
current uncertainties, as well as of different technical choices (waveforms,
thresholds on SNR, detector network configurations, etc.). Our results will
then contribute to giving an overall picture of the capabilities of 3G
detectors. Compared in particular to Borhanian & Sathyaprakash (2022), our
work is more focused on ET. To avoid a proliferation of plots, and of lines in
each plot, in this work we will only report the results for two 3G
configurations: ET alone, and a network of ET together with two CE detectors
(ET+2CE), and we will compare them with the expectations for the forthcoming
LVK–O4 run. 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. Our code,
however, can be used to study a large variety of networks.

The paper is organized as follows. In Sect. 2 we recall the Fisher matrix
formalism, as well as its limitations. In Sect. 3 we discuss the modelization
of the GW signal. In particular, since BNSs can stay in the bandwidth of 3G
detectors for as long as hours or a day, it is important to include the effect
of the Earth’s rotation on the signal, which can be exploited to improve the
sky localization of the source. In Sect. 4 we discuss the assumptions that we
make on the astrophysical populations of BBHs, BNSs and NSBHs, and we present
the detectors that we will study, computing their horizon (and the distance at
which $50\%$ of the events are detected) for different type of sources. Sect.
5 contains the bulk of our results for parameter estimation of BBHs, BNSs and
NSBHs. Further results are included in the Appendices, including, in App. B, a
comparison with the results presented in Borhanian & Sathyaprakash (2022) and
in Ronchini et al. (2022).

##  2 Formalism for parameter estimation

###  2.1 Fisher Information Matrix

In this section we recall the basic formalism and interpretation of the use of
the Fisher matrix formalism in GW parameter estimation. We refer to Cutler &
Flanagan (1994); Vallisneri (2008); Rodriguez et al. (2013) for comprehensive
treatments. We assume that the time–domain signal in a GW detector can be
written as the superposition of an expected signal $h_{0}$ and stationary,
Gaussian noise $n$ with zero mean:

$s(t)=h_{0}(t)+n(t)\;.$ (1)

The statistical properties of the noise are encoded in the one–sided power
spectral density, defined by

$\langle\tilde{n}^{*}(f)\,\tilde{n}(f^{\prime})\rangle=\frac{1}{2}\delta(f-f^{\prime})S_{n}(f)\;,$
(2)

where the tilde denotes a temporal Fourier transform. This determines an inner
product for any two time–domain signals $a(t)$, $b(t)$:

$\left(a\,|\,b\right)=4\Re{\int_{0}^{\infty}\differential{f}\,\frac{\tilde{a}^{*}(f)\,\tilde{b}(f)}{S_{n}(f)}}\;.$
(3)

Under the assumption of stationarity and Gaussianity, from Eq. (2) we have
that the variance of a Fourier mode with frequency $f$ is $S_{n}(f)/2$, so the
probability that the noise has a given distribution $n_{0}$ can be written in
terms of Eq. (3) as $p(n_{0})\propto\exp{-\left(n_{0}\,|\,n_{0}\right)/2}.$
Together with Eq. (1), this results in the following likelihood for a data
realisation $s$ conditioned on the waveform parameters ${\bf\it\theta}$:

$\mathcal{L}(s\;|\;{\bf\it\theta})\propto\exp{-\left(s-h({\bf\it\theta})\,|\,s-h({\bf\it\theta})\right)/2}\;.$
(4)

Note that ${\bf\it\theta}$ refers to the parameters of the template waveform
model, which may be different from the ones of the actual signal $h_{0}$, that
we denote by ${\bf\it\theta}_{0}$. Using Eq. (3) we can also express the
signal–to–noise ratio (SNR) of the true signal as

$\text{SNR}=\left(h_{0}\,|\,h_{0}\right)^{\nicefrac{{1}}{{2}}}\;.$ (5)

The Fisher Information Matrix (FIM) for the likelihood in Eq. (4) is defined
as:

$\Gamma_{ij}\equiv-\evaluated*{\langle\partial_{i}\partial_{j}\log\mathcal{L}(s\,|\,{\bf\it\theta})\rangle_{n}}_{{\bf\it\theta}={\bf\it\theta}_{0}}=\left(h_{i}\,|\,h_{j}\right)\;,$
(6)

where $h_{i}\equiv\partial_{i}h$, and the notation $\langle\,\dots\rangle_{n}$
denotes an average over noise realizations with fixed parameters. The last
equality is a consequence of the property
$\langle\left(a\,|\,n\right)\left(n\,|\,b\right)\rangle_{n}=\left(a\,|\,b\right)$
that follows from Eq. (2). Intuitively, near a maximum of the likelihood, the
latter is approximated by a multivariate Gaussian with covariance
$\Gamma_{ij}^{-1}$. To understand the correct interpretation of this statement
it is useful to consider an expansion of the template signal around the
waveform with true parameters’ values ${\bf\it\theta}_{0}$ as

$h({\bf\it\theta})=h_{0}+h_{i}\,\delta\theta^{i}+\dots\qquad\;,$ (7)

where $\delta\theta^{i}\equiv\theta^{i}-\theta_{0}^{i}$. This approximation,
where only first derivatives of the signal are included, is known as the
_linearized signal approximation_ (LSA). The LSA likelihood is obtained by
inserting Eq. (7) in Eq. (4):

$\mathcal{L}(s\,|\,{\bf\it\theta})\propto\exp[-\frac{1}{2}\left(n\,|\,n\right)+\delta\theta^{i}\left(n\,|\,h_{i}\right)-\frac{1}{2}\delta\theta^{i}\delta\theta^{j}\left(h_{i}\,|\,h_{j}\right)\Big{]}\;.$
(8)

It can be shown that the LSA is equivalent to the limit of large SNR
(Vallisneri, 2008).666Formally, the LSA likelihood is equivalent to the
leading order term in a series expansion in $1/\text{SNR}$. Two
interpretations of the meaning of the FIM appearing in Eq. (8) are possible:

  * •

From a _frequentist_ point of view, one can compute the maximum likelihood
(ML) estimator $\hat{{\bf\it\delta\theta}}$ and its covariance over noise
realizations at fixed parameters, which yields

$\displaystyle\hat{\delta\theta}^{i}$
$\displaystyle={\left(h_{i}\,|\,h_{j}\right)}^{-1}\left(n\,|\,h_{j}\right)\,,$
(9)
$\displaystyle\langle\hat{\delta\theta}^{i}\hat{\delta\theta}^{j}\rangle_{n}$
$\displaystyle=\left(h_{i}\,|\,h_{j}\right)^{-1}\,.$

This shows that _the inverse of the FIM is equal to the covariance of the
frequentist maximum–likelihood estimator in the LSA/large SNR limit and for
Gaussian noise_.777It can also be shown that the ML estimator is unbiased,
i.e. the expectation value of ${\bf\it\theta}$ is ${\bf\it\theta}_{0}$.

  * •

From a _Bayesian_ perspective one can compute the mean and variance of the
posterior probability $p({\bf\it\theta}\,|\,s)\propto
p({\bf\it\theta})\mathcal{L}({\bf\it\theta}\,|\,s)$. Using a flat prior
$p({\bf\it\theta})\propto\text{const.}$, we obtain

$\displaystyle\langle{\delta\theta}^{i}\rangle$
$\displaystyle=\frac{\int\differential{\bf\it\theta}\,{\delta\theta}^{i}\,p({\bf\it\theta}\,|\,s)}{\int\differential{\bf\it\theta}\,p({\bf\it\theta}\,|\,s)}={\left(h_{i}\,|\,h_{j}\right)}^{-1}\left(n\,|\,h_{j}\right)\,,$
(10) $\displaystyle\langle{\delta\theta}^{i}{\delta\theta}^{j}\rangle$
$\displaystyle=\langle\big{(}{\delta\theta}^{i}-\langle{\delta\theta}^{i}\rangle\big{)}\big{(}{\delta\theta}^{j}-\langle{\delta\theta}^{j}\rangle\big{)}\rangle=\left(h_{i}\,|\,h_{j}\right)^{-1}\,.$

In the above equations, $\langle\,\dots\rangle$ denotes an average over
parameter realizations at fixed data. Hence, _the inverse of the FIM is also
equal to the covariance of the Bayesian posterior probability distribution of
the true parameters ${\bf\it\theta}_{0}$ for a given experiment with data $s$,
assuming: a flat prior, the LSA/large SNR limit and Gaussian noise._888In this
case, note that an additional uncertainty might be present, i.e. the fact that
the posterior might not peak at the true parameters ${\bf\it\theta}_{0}$ due
to contribution of the noise. In this sense, one should talk about uncertainty
rather than error (Vallisneri, 2008; Rodriguez et al., 2013). However, under
the LSA, this bias has expectation value zero under repeated realizations of
the noise with covariance $\Gamma_{ij}^{\mkern-4.0mu-1}$. Note also that this
coincidence between frequentist ML estimator and posterior mean is in general
not true beyond the LSA (Vallisneri, 2008).

The Bayesian point of view is the most useful to compare to results of an
actual parameter estimation as well as to incorporate the effects of priors on
the source parameters. In this case, one can explicitly re–write the LSA
likelihood (omitting factors that do not depend on ${\bf\it\delta\theta}$) as

$\mathcal{L}(s\,|\,{\bf\it\theta})\propto\exp\left[-\frac{1}{2}\left(h_{i}\,|\,h_{j}\right)\big{(}\delta\theta^{i}-\langle{\delta\theta}^{i}\rangle\big{)}\big{(}\delta\theta^{j}-\langle{\delta\theta}^{j}\rangle\big{)}\right]\,.$
(11)

On the other hand, in absence of real data, as is the case when forecasting
parameter estimation capabilities of future experiments, it makes sense also
to consider the frequentist approach which makes statements about ensembles of
possible data realizations (Cutler & Flanagan, 1994).

Finally, another important interpretation of the FIM concerns the Cramer–Rao
bound: the inverse FIM gives the lower bound for the covariance of any
unbiased estimator of the true parameters ${\bf\it\theta}_{0}$ under different
noise realization (hence, it is a frequentist error estimate). However, its
interpretation is subtle (Vallisneri, 2008) and this bound does not translate
in a bound on the variance of the Bayesian posterior (Rodriguez et al., 2013).

###  2.2 Singularities

The most relevant issues when computing the FIM are the accurate computation
of the derivatives and the presence of singularities, which can result in
instability of the numerical inversion. The computational aspects as
implemented in GWFAST are described in Sect. 3.3, while here we discuss the
interpretation of singular or ill–conditioned matrices.

If the FIM is exactly singular, the likelihood has one or more directions in
parameter space along which it is exactly constant, corresponding to null
eigenvalues. If the LSA approximation was exact, one could conclude that one
or more combination of parameters are impossible to constrain and discard
them. In particular, one can discard the combination of parameters
corresponding to the null singular values of a singular–value decomposition.
However, this is true only in the LSA, as higher–order contributions may cure
the singularity, and using this approach in the FIM does not necessarily lead
to realistic forecasts of the parameter estimation performance, since only
truly degenerate combinations should be discarded.

A more common situation to face in practice is the case where the matrix is
not exactly singular but _ill conditioned_ , i.e. the inverse of the condition
number (defined as the ratio of the largest to the smallest eigenvalues) is
comparable or smaller than machine precision. This can lead to large
amplifications of inversion errors, resulting in a covariance that may be
inaccurate even at the 100% level. Condition numbers might be improved by
numerical techniques (see Sect. 3.3 for the implementation in GWFAST), but the
physical interpretation is that the directions in parameter space
corresponding to small eigenvalues need large variations of the parameters to
produce relevant changes in the waveform (at least comparable to the noise
level). Hence, the LSA might no longer be sufficient to accurately describe
the likelihood over all the range of interest (Vallisneri, 2008), even if the
SNR at the true parameters’ values is high enough. We will discuss this point
in Sect. 3.3. One can resort to some regularization of the singularities, for
example truncating the singular values of the singular–value decomposition to
the minimum allowed numerical precision; however this has not a direct link to
actual generalizations of the FIM that might cure the singularity, such as the
inclusion of higher order terms (or of priors, see Sect. 2.3). In general,
when using the FIM to forecast parameter estimation, one must be aware that
the bad conditioning of the matrix is pointing to some possible breakdown of
the LSA on the likelihood surface irrespective of the specific regularization
technique used.

In this work we prefer not to resort to regularization of the singular values
nor to discard combinations of parameters corresponding to small eigenvalues.
The reason is that this might lead to underestimating the errors on the
waveform parameters that are most affected by degeneracies, unless careful
checks and comparisons with other methods are performed.

One possible alternative approach corresponds to discarding all the matrices
with condition number larger than the inverse machine precision. This is the
most conservative option, which might on the other hand lead to discarding
many events for which a sensible inversion of the Fisher matrix could still be
obtained. A somewhat intermediate possibility is to discard all the events for
which the inversion error of the FIM is larger than a given threshold. In this
work [as, e.g., in Berti et al. (2005), see their App. B], we adopt the latter
strategy and quantify the inversion error as

$\epsilon={||\Gamma\cdot\Gamma^{-1}-\mathbb{1}||}_{\rm
max}=\max_{ij}|(\Gamma\cdot\Gamma^{-1}-\mathbb{1})_{ij}|\;.$ (12)

The chosen threshold for this work is $\epsilon_{\rm
max}=$5\text{\times}{10}^{-2}$$.999Note that the v1 arXiv version of this
paper we used $\epsilon_{\rm max}=$1\text{\times}{10}^{-3}$$. This is at the
origin of some minor differences in the results. In App. C we detail the
choice of this value and compare the results with other methods to obtain
forecasts on the uncertainty, including a direct sampling from the likelihood
in Eq. (11) which does not rely on the inversion of the Fisher matrix. In App.
B we compare to other choices in the literature and discuss the impact of
different approaches on the population of sources studied in this paper.

In a Bayesian context, the use of priors can cure singularities and is the
most realistic way to proceed, as priors are actually used in parameter
estimation. We discuss the role of priors in the next section, while their
effect is also studied in App. Cfv.

###  2.3 Role of priors

In a Bayesian parameter estimation problem, prior distributions are expected
to have a role whenever their information content is more restrictive than the
one in the likelihood, typically either because the prior changes
significantly in the region of non–vanishing likelihood or because hard
boundaries are present. If the LSA is valid, the presence of priors can play
an important role in curing possible singularities of the FIM. In practice,
the FIM formalism allows to treat in a simple way only the case of a Gaussian
prior
$p(\mathbf{\theta})\propto\exp{-P_{ij}\delta\theta^{i}\delta\theta^{j}}$, in
which case adding the prior amounts to substitute
$\left(h_{i}\,|\,h_{j}\right)\mapsto\left(h_{i}\,|\,h_{j}\right)+P_{ij}$ in
Eq. (10)–(11), i.e. to add the prior matrix to the FIM.

In a GW parameter estimation problem, we have to deal with the fact that some
parameters might require priors that are far from Gaussian, in particular to
incorporate information on the physical range of some of them, such as the
symmetric mass ratio $\eta$ (which is constrained to be in the range
$[0;\,0.25]$), the luminosity distance $d_{L}$ (constrained to be positive),
and the angles. In Markov chain Monte Carlo (MCMC) analyses, priors on the
angles might not be used if not for increasing speed, since the full
likelihood carries information on the periodicity of these variables. In
contrast, the Fisher formalism does not have any information about the
periodicity of the likelihood with respect to some parameters, hence the use
of a prior can be a more realistic choice. For angular variables, a crude but
simple approximation can be the use of a Gaussian prior of width $2\pi$. The
situation is more complicated for other parameters, in particular $\eta$ and
$d_{L}$. When using a FIM analysis to forecast overall trends in future
experiments rather than concentrating on single events, one can check that the
majority of the events have a predicted $1\sigma$ contour that does not exceed
the physical boundary. For more realistic estimates, an exact prior can be
included by explicitly drawing samples from the likelihood in Eq. (11), using
rejection sampling to account for the prior, and estimating the posterior
covariance from the remaining samples. We show the effect of adopting this
procedure, and compare it to different inversion methods, in App. C.

###  2.4 Limits of applicability

---
Figure 1: Scatter plot of the condition numbers of the Fisher matrices as a
function of the system inclination angle, computed for BBH events detected at
a network consisting of ET+2CE, using the IMRPhenomHM waveform model. In the
left panel we report on the color bar the inversion errors associated with the
events and in the central panel the SNRs and in the right panel the relative
errors on the luminosity distance. The horizontal dashed line shows the
inverse machine precision limit, of ${10}^{15}$.

We have seen that the applicability of the FIM relies on the LSA/high SNR
limit. From the discussion in Sect. 2.2, however, it is clear that such an
approximation should hold within all the likelihood region of interest – for
example, the 1$\sigma$ or 2$\sigma$ contours. For this to be the case, not
only the SNR at the true values should be high, but of particular relevance is
the problem of conditioning: a high condition number might signal a breakdown
of the LSA in the region of interest (Finn, 1992; Vallisneri, 2008).

One instructive example concerns the inclination angle of the source with
respect to the observer, $\iota$. The amplitudes of the two GW polarisations
depend differently on this parameter and their joint measurement would allow
disentangling $\iota$ from the distance to the source. However, it is for
nearly face–on binaries ($\iota\sim 0$) that the signal is louder, but in this
case the two polarisations are nearly equal (Nissanke et al., 2010; Schutz,
2011; Usman et al., 2019), which leads to a strong degeneracy of the
inclination angle with the distance to the source. Thus, an event observed
face–on will be louder than the same system observed in an inclined
configuration, but its FIM will be ill–conditioned. In the limit of exactly
face–on systems, $\iota=0$, and ignoring the presence of higher modes in the
GW signal, the FIM will even be exactly degenerate, since the first derivative
of the signal with respect to $\iota$ vanishes exactly in this limit. We refer
to Fig. 1 for an illustration based on the population studied in this work.
This shows that, for events with nearly face–on configuration, we can be in
the situation where the SNR is large but the FIM is ill–conditioned (Cutler &
Flanagan, 1994).101010In these cases, distance estimates will be particularly
affected, since the luminosity distance has a strong correlation with the
inclination, and the FIM can yield inaccurate predictions (Cutler & Flanagan,
1994). It is however possible to extend the approximation of the full
posterior beyond linear order for the marginal posterior in the subspace
$(\iota,d_{L})$ (Cutler & Flanagan, 1994; Chassande-Mottin et al., 2019).

In practise, assessing the validity of the LSA on all the likelihood surface
without explicitly calculating higher–order corrections or resorting to
explicit Monte Carlo analyses is a subtle problem. Vallisneri (2008) proposed
a “maximum mismatch criterion” to determine the validity of the LSA, which
consists in sampling the $1\sigma$ likelihood surface, compute the difference
between the waveforms on this surface and at the true value, $\Delta
h=h({\bf\it\theta}_{1\sigma})-h_{0}$ and computing the ratio $r$ of the LSA
likelihood to the full likelihood, which can be shown to be

$|\log(r)|=\frac{1}{2}\,\left(\theta^{i}h_{i}-\Delta
h\,|\theta^{j}h_{j}-\Delta h\,h_{j}\right)\;.$ (13)

Then, a signal can be considered linear if some threshold is satisfied by this
mismatch, for example by requiring that $|\log(r)|<0.1$ for $90\%$ of the
points on the $1\sigma$ likelihood surface (Vallisneri, 2008; Rodriguez et
al., 2013). This is an internal consistency criterion rather than a proof of
validity. However, explicit comparison with MCMC showed that, while this is
indeed a sufficient condition for the LSA, many systems showing good agreement
between the FIM and a full MCMC analysis failed this test (Rodriguez et al.,
2013). In summary, there seems not to be a conclusive, flexible and
computationally cheap test to assess the validity of the FIM over all the
surface of interest. In this work we apply the threshold on the inversion
error of the FIM described in Sect. 2.2 and App. C. However, one must be aware
of such limitations of the FIM approach. We believe these issues to be a
further motivation for the development and comparison of different
implementations of the FIM technique, in order to assess the robustness of the
predictions to different choices.

##  3 Modeling the gravitational–wave signal

The response of a detector to a GW signal emitted by a coalescing binary
system is given by a linear combination of the two polarisations
$(h_{+},h_{\times})$, obtained by their projection on the detector arms by
suitable “antenna pattern functions” $(F_{+},F_{\times})$ that depend on the
source position and polarisation angle, as well as the location, orientation
and shape of the detector, which we denote collectively by ${\bf\it\lambda}$.
In full generality, we denote the parameters of the waveform by
${\bf\it\theta}=\\{{\cal
M}_{c},\eta,d_{L},\theta,\phi,\iota,\psi,t_{c},\Phi_{c},\chi_{1,x},\chi_{2,x},\chi_{1,y},\chi_{2,y},\chi_{1,z},\chi_{2,z},\Lambda_{1},\Lambda_{2}\\}$
[see e.g. Maggiore (2007)], where ${\cal M}_{c}$ denotes the detector–frame
chirp mass, $\eta$ the symmetric mass ratio, $d_{L}$ the luminosity distance
to the source, $\theta$ and $\phi$ are the sky position coordinates, defined
as $\theta=\pi/2-\delta$ and $\phi=\alpha$ (with $\alpha$ and $\delta$ right
ascension and declination, respectively), $\iota$ the inclination angle of the
binary with respect to the line of sight, $\psi$ the polarisation angle,
$t_{c}$ the time of coalescence, $\Phi_{c}$ the phase at coalescence,
$\chi_{i,c}$ the dimensionless spin of the object $i=\\{1,2\\}$ along the axis
$c=\\{x,y,z\\}$ and $\Lambda_{i}$ the dimensionless tidal deformability of the
object $i$ (which is present only for systems containing a NS). Instead of
$\Lambda_{1},\Lambda_{2}$, we will actually use (Wade et al., 2014)

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

which have the advantage that $\tilde{\Lambda}$ is the combination that enters
the inspiral waveform at 5 PN, while $\delta\tilde{\Lambda}$ first enters at 6
PN.

In the time domain, the signal of the quadrupole mode is given by

$h(t,{\bf\it\theta},{\bf\it\lambda})=A_{+}(t,{\bf\it\theta})F_{+}(t,{\bf\it\theta},{\bf\it\lambda})\cos{\Phi(t,{\bf\it\theta})}+A_{\times}(t,{\bf\it\theta})F_{\times}(t,{\bf\it\theta},{\bf\it\lambda})\sin{\Phi(t,{\bf\it\theta})}\,.$
(15)

The amplitudes $A_{+,\times}$ and the phase $\Phi$ are obtained from a
waveform model. We describe the models used in Sect. 3.1. In particular, in
this work and in GWFAST we work in the frequency domain. In this case, when
computing the full signal in Eq. (15), an important point is that the
low–frequency sensitivity of 3G detectors, and in particular of ET, makes it
possible to observe the inspiral phase of low–mass events, such as BNSs, for
several hours to possibly one day. In those cases the pattern functions evolve
in time during the detection due to the change in the relative position of the
source and the detector. Since we work in frequency rather than time domain,
it is important to correctly account for this time evolution when Fourier
transforming the signal. We describe this in detail in Sect. 3.2.

###  3.1 Waveform models

We adopt Fourier domain full inspiral–merger–ringdown models, tuned on
Numerical Relativity (NR) simulations. Our code can be adapted to a large
variety of waveforms. As reference waveforms, for BBH, BNS and NSBH systems,
we will use, respectively:

IMRPhenomHM

(London et al., 2018; Kalaghatgi et al., 2020) this is a recent model, tuned
for BBH systems (in quasi–circular orbits) with non–precessing spins, which
takes into account the quadrupole of the signal and the sub–dominant modes
$(l,m)=(2,1),\ (3,2),\ (3,3),\ (4,3),\ {\rm and}\ (4,4)$. The contribution of
these higher modes is of fundamental importance for parameter estimation,
since they can break the degeneracy between the luminosity distance and
inclination angle. In fact, each mode depends on a different combination of
sines and cosines of the inclination angle $\iota$, through the spin–weighted
spherical harmonics ${}_{-2}Y^{lm}$ (Goldberg et al., 1967), while the
fundamental mode (2,2), only depends on the cosine, leading to a degeneracy
with $d_{L}$;

IMRPhenomD_NRTidalv2

(Dietrich et al., 2019) this model is an extension of IMRPhenomD (Husa et al.,
2016; Khan et al., 2016a), which also accounts for tidal effects in BNS
systems. In fact, the two neutron stars in a coalescing binary, differently
from black holes, will deform when getting closer to each other, and this
leaves clear signatures on the waveform, which have to be accurately modelled
and tuned to specific NR simulations;

IMRPhenomNSBH

(Pannarale et al., 2015; Dietrich et al., 2019), this model can describe the
signals coming from the merger of a NS and a BH, which have very distinctive
features. In particular, the mass ratios of these systems can be much higher
than that of BNSs and BBHs, so this model is tuned up to $q\sim 100$, and it
also accounts for tidal effects, since the NS will deform getting closer to
the BH. Moreover, this model is built and tuned to account for the fact that,
when it is close enough to the BH, the NS can either plunge into it or get
disrupted forming an accretion torus111111A torus can be formed also if the NS
plunges into the BH, resulting in a mildly disruptive merger, which is also a
case taken into account in the tuning of IMRPhenomNSBH. (Lattimer & Schramm,
1976), two scenarios that result in very different features on the waveforms.

We will then perform a comparison between the results obtained using different
waveform models. Note that, in addition to the above waveform models, GWFAST
includes Python implementations of the models TaylorF2_RestrictedPN and
IMRPhenomD, as well as a wrapper to all waveforms contained in the LIGO
Algorithm Library, LAL (note, however, that when using the latter, which are
implemented in C, it is no longer possible to exploit Python vectorization).

###  3.2 Detector response and effect of Earth’s rotation

In this work and in GWFAST we use the expression of the ‘pattern functions’
which takes into account the rotation of the Earth given in Jaranowski et al.
(1998). We collectively denote the parameters characterizing the detector by
${\bf\it\lambda}=\\{\lambda,\ \varphi,\ \gamma,\ \zeta\\}$, where $\zeta$ is
the angle between the detector’s arms (e.g.
$90\text{\,}\mathrm{\SIUnitSymbolDegree}$ for an L–shaped detector), $\lambda$
and $\varphi$ are the detector’s latitude and longitude, respectively, and
$\gamma$ is the angle formed by the arms bisector and East. We then have for
the pattern functions:

$\displaystyle F_{+}(t;\theta,\phi,\psi,{\bf\it\lambda})$
$\displaystyle=\sin{\zeta}\,\big{[}a(t;\theta,\phi,{\bf\it\lambda})\cos{2\psi}+b(t;\theta,\phi,{\bf\it\lambda})\sin{2\psi}\big{]}\,,$
(16) $\displaystyle F_{\times}(t;\theta,\phi,\psi,{\bf\it\lambda})$
$\displaystyle=\sin{\zeta}\,\big{[}b(t;\theta,\phi,{\bf\it\lambda})\cos{2\psi}-a(t;\theta,\phi,{\bf\it\lambda})\sin{2\psi}\big{]}\,,$

with (omitting from now on the explicit dependence on the parameters)

$\displaystyle a(t)$
$\displaystyle=\dfrac{1}{16}\sin{2\gamma}(3-\cos{2\lambda})(3+\cos{2\theta})\cos{[2(\phi-\varphi-2\pi
f_{\oplus}t)]}-\dfrac{1}{4}\cos{2\gamma}\sin{\lambda}(3+\cos{2\theta})\sin{[2(\phi-\varphi-2\pi
f_{\oplus}t)]}$ (17)
$\displaystyle\quad\,+\dfrac{1}{4}\sin{2\gamma}\sin{2\lambda}\sin{2\theta}\cos{(\phi-\varphi-2\pi
f_{\oplus}t)}-\dfrac{1}{2}\cos{2\gamma}\cos{\lambda}\sin{2\theta}\sin{(\phi-\varphi-2\pi
f_{\oplus}t)}$
$\displaystyle\quad\,+\dfrac{3}{4}\sin{2\gamma}\cos^{2}{\lambda}\sin^{2}{\theta}\,,$
$\displaystyle b(t)$
$\displaystyle=\cos{2\gamma}\sin{\lambda}\cos{\theta}\cos{[2(\phi-\varphi-2\pi
f_{\oplus}t)]}+\dfrac{1}{4}\sin{2\gamma}(3-\cos{2\lambda})\cos{\theta}\sin{[2(\phi-\varphi-2\pi
f_{\oplus}t)]}$
$\displaystyle\quad\,+\cos{2\gamma}\cos{\lambda}\sin{\theta}\cos{(\phi-\varphi-2\pi
f_{\oplus}t)}+\dfrac{1}{2}\sin{2\gamma}\sin{2\lambda}\sin{\theta}\sin{(\phi-\varphi-2\pi
f_{\oplus}t)}\,,$

where $f_{\oplus}\simeq$1\text{\,}{\mathrm{day}}^{-1}$$ is the Earth’s
rotational frequency.

The effect of the Earth’s rotation on the signal consists of an amplitude
modulation, due to the variation of the pattern functions with time as in Eq.
(16), a phase modulation arising for the same reason, due to the fact that the
pattern functions relative to the two polarisations evolve differently in
time,121212The phase modulation is implicit when writing the two components of
the signal separately, as in Eq. (15), but becomes apparent by re–writing it
as
$h(t)=A(t)\cos\\{\Phi(t)-\arctan{[{A_{\times}F_{\times}(t)}/{(A_{+}F_{+}(t))}}]\\}$
with
$A(t)=[(A_{\times}F_{\times}(t))^{2}+(A_{+}F_{+}(t))^{2}]^{\nicefrac{{1}}{{2}}}$.
and a Doppler contribution to the phase due to the relative motion between the
source and the detector (Cutler, 1998; Cornish & Larson, 2003). In time
domain, the Doppler contribution can be conveniently expressed as a
time–dependent shift of the time variable [see, e.g. Sect. 7.6.2 of Maggiore
(2007) and Wen & Chen (2010)]. We work with signals in frequency domain. To
compute the Fourier transform one can adopt the Stationary Phase
Approximation, that applies if the change in the amplitude during a cycle is
much slower than the corresponding change in the phase. This is the case for
each of the two terms in the sum in Eq. (15),131313The Stationary Phase
Approximation is usually adopted when the pattern functions do not depend on
time, in which case the condition is satisfied, as shown e.g. in Maggiore
(2007). This remains true even including the time dependence in Eq. (16),
given that the frequency of Earth’s rotation is much smaller than the
frequency of the gravitational wave when it enters the detectors band. so,
neglecting for the moment the time shift due to the Doppler effect, we get

$\displaystyle\tilde{h}_{+,\times}(f)$ $\displaystyle=\int
A_{+,\times}(t)F_{+,\times}(t)\cos{\Phi(t)}\,e^{2\pi ift}\differential{t}$
(18)
$\displaystyle\simeq\dfrac{1}{2}A_{+,\times}(t^{\ast})F_{+,\times}(t^{\ast})\left(\dfrac{2\pi}{\ddot{\Phi}(t^{\ast})}\right)^{\nicefrac{{1}}{{2}}}\exp{i\left[2\pi
ft^{\ast}\,-\,\Phi(t^{\ast})-\dfrac{\pi}{4}\right]}\,,$

where the stationary point $t^{\ast}(f)$ is determined by the condition $2\pi
f=\dot{\Phi}(t^{\ast})$. To lowest order in the post–Newtonian (PN) expansion,
this gives

$t^{\ast}(f)=t_{c}-\dfrac{5}{256}\left(\dfrac{G{\cal
M}_{c}}{c^{3}}\right)^{-\nicefrac{{5}}{{3}}}(\pi
f)^{-\nicefrac{{8}}{{3}}}\times\left[1+\order{(\pi fGM_{\rm
tot}/c^{3})^{\nicefrac{{2}}{{3}}}}\right]\,,$ (19)

where $t_{c}$ is the time of coalescence and $M_{\rm tot}$ the total mass of
the system. Note, however, that in our analysis we compute the time to
coalescence at 3.5 PN order, as in Buonanno et al. (2009b), Eq. (3.8b). The
amplitude modulation is reflected in the fact that the Fourier transformed
signal contains the expressions of $F_{+}$ and $F_{\times}$ evaluated at the
stationary point $t^{\ast}(f)$, according to Eq. (18) (see also Zhao & Wen
(2018)). An illustration of this effect is given in Fig. 2, for a BNS signal
analogous to GW170817.

We next consider the Doppler effect due to the rotation of the Earth around
its axis. Its contribution to the calculation of $t^{\ast}(f)$ and to the
amplitude of the Fourier transform on the right–hand side of Eq. (18) is
totally negligible, since the frequency scale $1\,{\rm day}^{-1}$ is very
small compared to the frequencies relevant for ground–based detectors, so we
only need to consider its effect in the phase. This is encoded in the time
dependence of the time delay corresponding to the travel time of the signal
from the origin of the reference frame, i.e. the center of the Earth, to the
detector. This time delay in time domain results in a phase shift in Fourier
domain, known as location phase factor, given by

$\phi_{L}=2\pi f\Delta t_{L}=-2\pi f\dfrac{R_{\oplus}}{c}\hat{{\bf\it
s}}(\theta,\phi)\dotproduct\hat{{\bf\it d}}(\lambda,\varphi,t^{\ast}(f))\,,$
(20)

where $R_{\oplus}\simeq$6371\text{\,}\mathrm{km}$$ is the Earth’s radius,
$\hat{{\bf\it s}}(\theta,\phi)$ and $\hat{{\bf\it d}}(\lambda,\varphi)$ are
the unit vectors pointing to the source and the detector, respectively, with
the second one being evaluated at the stationary point $t^{\ast}(f)$ when
including Earth’s rotation, i.e.

$\displaystyle\hat{{\bf\it s}}(\theta,\phi)$
$\displaystyle=\Big{(}\sin{\theta}\cos{\phi},\ \sin{\theta}\sin{\phi},\
\cos{\theta}\Big{)}\,,$ (21) $\displaystyle\hat{{\bf\it
d}}(\lambda,\varphi,t^{\ast}(f))$
$\displaystyle=\Big{(}\cos{\lambda}\cos{(\varphi+2\pi f_{\oplus}t^{\ast})},\
\cos{\lambda}\sin{(\varphi+2\pi f_{\oplus}t^{\ast})},\
\sin{\lambda}\Big{)}\,.$

Figure 2: Comparison between the amplitude of a BNS signal analogous to
GW170817, as seen by a single instrument of the ET triangular detector (whose
sensitivity curve is shown in green) accounting for the effect of Earth’s
rotation (violet curve) and without accounting for it (orange curve). The
dashed vertical lines indicate the amount of time left before coalescence and
the amplitude is computed using the waveform model IMRPhenomD_NRTidalv2. For
comparison, we also report in blue a representative LIGO Livingston
sensitivity curve during the second observing run, available at
https://dcc.ligo.org/LIGO-G1801952/public.

The effect of this location phase is of particular relevance in the case of a
network of detectors, since the difference in the arrival time of the signal
between different detectors gives a fundamental information to localize the
source. The Doppler effect due to the rotation of the Earth around the Sun,
which is relevant for instance for LISA, is negligible for the signals
expected at ground–based detectors.

In summary, the signal in the frequency domain takes the form

$\begin{split}\tilde{h}(f)=\Big{[}F_{+}\big{(}t^{\ast}(f);{\bf\it\theta},{\bf\it\lambda}\big{)}\,A^{(F)}_{+}(f,{\bf\it\theta})&+F_{\times}\big{(}t^{\ast}(f);{\bf\it\theta},{\bf\it\lambda}\big{)}\,A^{(F)}_{\times}(f,{\bf\it\theta})\Big{]}\times\\\
&\exp{i\left[2\pi
ft_{c}+\phi_{L}\left(t^{\ast}(f)\right)-\Phi_{c}+\Phi^{(F)}(f;{\bf\it\theta})-\frac{\pi}{4}\right]}\,,\end{split}$
(22)

where the quantities $A^{(F)}_{+,\times}(f,{\bf\it\theta})$ and
$\Phi^{(F)}(f;{\bf\it\theta})$ are the output of the waveform model in
frequency domain [which already includes the factor
$(2\pi/\ddot{\Phi}(t^{\ast}))^{\nicefrac{{1}}{{2}}}$ in Eq. (18)], and
$\Phi_{c}$ the value of the phase at coalescence, which we explicitly
separated from $\Phi^{(F)}(f;{\bf\it\theta})$ for clarity.

When including the contribution of higher modes, the signal can be expressed
as a superposition of multipoles of spin–weighted spherical harmonics
(Goldberg et al., 1967) of weight $-2$, whose Fourier transform can be
computed independently. To compute the frequency–domain signal one can thus
proceed as before, obtaining

$\begin{split}\tilde{h}(f,{\bf\it\theta},{\bf\it\lambda})=\Big{[}F_{+}(t^{\ast}(f),{\bf\it\theta},{\bf\it\lambda})\,\tilde{h}_{+}(f,{\bf\it\theta})&+F_{\times}(t^{\ast}(f),{\bf\it\theta},{\bf\it\lambda})\,\tilde{h}_{\times}(f,{\bf\it\theta})\Big{]}\times\\\
&\exp{i\left[2\pi
ft_{c}+\phi_{L}\left(t^{\ast}(f)\right)-\Phi_{c}-\frac{\pi}{4}\right]}\,,\end{split}$
(23)

where the quantities $\tilde{h}_{+,\times}$ are now the output of the waveform
model, computed in frequency domain [see e.g. García-Quirós et al. (2020)] as

$\displaystyle\tilde{h}_{+}(f,{\bf\it\theta})$ $\displaystyle=\sum_{\ell\geq
2}\sum_{m=0}^{\ell}\dfrac{1}{2}\left[{}_{-2}Y_{\ell
m}(\iota)+(-1)^{\ell}_{-2}Y_{\ell-m}^{\ast}(\iota)\right]A^{(F)}_{\ell
m}(f,{\bf\it\theta})\exp{i\Phi^{(F)}_{\ell m}(f,{\bf\it\theta})},$ (24)
$\displaystyle\tilde{h}_{\times}(f,{\bf\it\theta})$
$\displaystyle=\sum_{\ell\geq
2}\sum_{m=0}^{\ell}\dfrac{i}{2}\left[{}_{-2}Y_{\ell
m}(\iota)-(-1)^{\ell}_{-2}Y_{\ell-m}^{\ast}(\iota)\right]A^{(F)}_{\ell
m}(f,{\bf\it\theta})\exp{i\Phi^{(F)}_{\ell m}(f,{\bf\it\theta})}\,,$

with $A^{(F)}_{\ell m}$ and $\Phi^{(F)}_{\ell m}$ being the amplitude and
phase of the $(\ell,m)$ mode in frequency domain.

Another noteworthy aspect is that, when using a triangular–shaped detector,
like ET (see Sect. 4.2 for a brief description of its design), it is not
needed to compute the signal in each interferometer, since the sum of all the
signals will vanish by geometrical reasons (Freise et al., 2009). In fact, the
difference among the signals observed in each instrument will arise only from
the pattern functions, that depend on the orientation of the detector as

$F_{n}\propto
C_{1}\sin\left(2\gamma+n\dfrac{2\pi}{3}\right)+C_{2}\cos\left(2\gamma+n\dfrac{2\pi}{3}\right);\qquad
n\in\\{1,2,3\\},$ (25)

where $\gamma$ denotes the orientation of the first interferometer (the angle
between east and the bisector of the first arm in the expressions we use) and
the terms $C_{1}$ and $C_{2}$ collect all the other terms appearing in the
pattern functions. We here assumed the three interferometers to be equal and
co–located, which is a sensible approximation for ground–based detectors. It
can then be trivially shown that $F_{1}+F_{2}+F_{3}=0$, thus the signal in one
of the instruments can simply be obtained from the signals in the other two,
reducing the amount of needed calculations by one third.141414 This is true in
general for a detector forming a closed loop of generic shape, again assuming
the various interferometers to be equal and co–located: the generalisation of
Eq. (25) to a detector with $N_{\rm ifo}$ interferometers is given by
$F_{n}\propto C_{1}\sin\left(2\gamma+{2\pi n}/{N_{\rm
ifo}}\right)+C_{2}\cos\left(2\gamma+{2\pi n}/{N_{\rm ifo}}\right)$ with
$n\in\\{1,2,\dots,N_{\rm ifo}\\},$ and computing the total signal one obtains
$\sum_{n=1}^{N_{\rm ifo}}F_{n}=(C_{1}\sin 2\gamma+C_{2}\cos
2\gamma)\sum_{n=1}^{N_{\rm ifo}}\cos\left(n\dfrac{2\pi}{N_{\rm
ifo}}\right)+(C_{1}\cos 2\gamma-C_{2}\sin 2\gamma)\sum_{n=1}^{N_{\rm
ifo}}\sin\left(n\dfrac{2\pi}{N_{\rm ifo}}\right)=0\,,$ since the two sums
trivially vanish. This is the basis of the so–called null–stream, that is the
data stream obtained by summing the signals in the instruments of a
closed–loop detector, and offers outstanding capabilities for analysing data
(Gürsel & Tinto, 1989; Wen & Schutz, 2005; Ajith et al., 2006; Chatterji et
al., 2006; Rakhmanov, 2006; Harry & Fairhurst, 2011; Regimbau et al., 2012;
Schutz & Sathyaprakash, 2020; Wong et al., 2021; Wong & Li, 2022). In
particular, as it has recently be shown in Goncharov et al. (2022) , it can be
used to eliminate glitches contaminating the detector data stream, and to get
unbiased estimates of the power spectral density (PSD). Notice also that the
null–stream in ET corresponds to the so–called T channel of the LISA space
interferometer (Prince et al., 2002).

###  3.3 Technical aspects

We here briefly outline some technical choices made to implement the above
features in GWFAST, which are described in more detail, and tested, in the
companion paper Iacovelli et al. (2022). GWFAST is a pure Python code that
fully exploits the vectorization capabilities of this language, and is able to
rapidly get signal–to–noise ratios and Fisher matrices for large catalogs of
events. In particular, we entirely re–wrote in Python, in fully vectorized
form, the waveform models IMRPhenomD, IMRPhenomD_NRTidalv2, IMRPhenomHM and
IMRPhenomNSBH, available in C in the LIGO Algorithm Library, LAL (LIGO
Scientific Collaboration, 2018). Together with this paper and GWFAST, we also
release the open–source library WF4Py containing the Python implementation of
the waveforms. The agreement with LAL is excellent, at the level of $\sim
10^{-14}$ in the inspiral part and $\sim 10^{-5}$ in the worst case in the
merger–ringdown. The difference in the latter case is entirely due to
interpolation routines needed to compute the waveform in the ringdown phase.
GWFAST anyway also implements an interface with LAL, which makes possible to
easily use all the waveforms available in that library.

When using pure Python waveforms, derivatives in GWFAST are computed using a
mixture of analytical differentiation and automatic differentiation. The
latter is a technique alternative to finite–difference, that allows us to
compute derivatives up to machine precision without issues of convergence due
to the step–size choice in finite difference, based on a decomposition of the
function in elementary functions, see Margossian (2018) for a review. In
particular, GWFAST makes use of the implementation in the JAX package
(Bradbury et al., 2018). This further has the advantage of allowing
vectorization of the computation of derivatives, thus fully exploiting the
capabilities of Python. For this to be efficient, a vectorized implementation
of the waveforms is needed, which motivated the development of WF4Py. Finally,
differentiation with respect to the parameters $d_{L},\ \theta,\ \phi,\
\iota,\ \psi,\ t_{c}$ and $\Phi_{c}$, which do not depend on the waveform
model, is implemented analytically, to further gain in speed and accuracy. As
additional checks of reliability, we verified that, for these parameters, the
analytical results and the result obtained by JAX agree at machine precision
($10^{-15}$ in our case). We also implemented the FIM formalism in an
independent code in Wolfram Mathematica in the case of the
TaylorF2_RestrictedPN waveform model, for which the calculation of all
derivatives can be more easily performed analytically. The agreement on the
diagonal elements is never worse than $10^{-4}$, including differences in the
integration routines in the two languages. We refer to Iacovelli et al. (2022)
for more details of these tests.

When resorting instead to the LAL waveforms, the computation of the
derivatives is performed using finite difference techniques, as implemented in
the numdifftools library,151515https://pypi.org/project/numdifftools/. with an
adaptive step–size, while the differentiation with respect to the parameters
$d_{L},\ \theta,\ \phi,\ \iota,\ \psi,\ t_{c}$ and $\Phi_{c}$, is still
performed analytically.

Coming to the inversion of the FIM, each row and column is normalized to the
square root of the diagonal entries, so that the resulting matrix has ones on
the diagonal and the remaining elements in the interval $[-1,\,1]$ (Harms et
al., 2022). This transformation is applied again after the inversion of the
resulting matrix to obtain the inverse of the original FIM. The inversion
itself is done by means of the Cholesky decomposition, which amounts to
express a (hermitian, positive–definite) matrix as a product of a lower
triangular matrix and its conjugate transpose. The inversion of a triangular
matrix is an easier task than that of a full matrix, which improves the
inversion.161616There is a small sub–sample of matrices which may be not
positive–definite due to the presence of very small eigenvalues that can
assume small negative values due to numerical fluctuations, in which case the
Cholesky decomposition cannot be found. For those matrices, GWFAST resorts by
default to a singular–value decomposition for the inversion. In any case, the
inversion error for these events is always larger than the threshold adopted,
so they are discarded.  Other methods supported by GWFAST are discussed in
Iacovelli et al. (2022). The inversion makes use of the Python library mpmath
for precision arithmetic.

##  4 Applications to current and future ground–based detectors

In this section, we use GWFAST to study the detection and parameter estimation
capabilities of current and future observatories for the three kind of GW
sources that have been detected so far, namely, BBHs, BNSs, and NSBHs. The
goal is to give realistic forecasts based on updated population models. For
the short term, we focus on the forthcoming O4 run of the LIGO–Virgo–KAGRA
(LVK) collaboration. For the long term, we consider a single Einstein
Telescope (ET) observatory, and a network “ET+2CE”, made by ET, located in
Europe, and two Cosmic Explorers (CE), located in the US.

###  4.1 Populations

We here describe our baseline assumptions for the populations of the three
kinds of compact binary systems considered in the analysis, (astrophysical)
BBHs, BNSs and NSBHs, also summarised in Tab. 1. The functional form and
numerical values of the parameters for all distributions are reported in App.
A.

BBH

We adopt a source–frame mass and spin distribution calibrated on the latest
LVK results (Abbott et al., 2021c), using the Power Law + Peak profile for the
former, and the Default model for the latter, and assume that the mass and
spin distributions do not evolve with redshift. For the local rate, the value
inferred from the GWTC–3 catalog is ${\cal R}_{0,{\rm
BBH}}=17^{+10}_{-6.7}\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$.171717This is the median
value with $90\%$ c.l. error, inferred using the Power Law + Peak
distribution. It is not explicitly given in Abbott et al. (2021c) (a plot of
${\cal R}_{{\rm BBH}}(z)$ is anyhow shown in their Fig. 13), but is available
in the associated data release at
https://zenodo.org/record/5655785#.YnUnPS8QN70, inside the file
PowerLawPeakObsOneTwoThree.json. We will then adopt ${\cal R}_{0,{\rm
BBH}}=$17\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ as our
reference value. Given the still limited redshift range of the detected
events, the rate distribution in redshift has large uncertainties, and only
its power-law behavior at low redshift has been constrained. Since 3G
detectors will cover a much broader redshift range, for which a power–law
behavior is not realistic, we choose to adopt a Madau–Dickinson profile (Madau
& Dickinson, 2014) in which the low–end slope is fixed to the LVK value and
the other parameters assume typical values used in the literature (Madau &
Fragos, 2017), see App. A for details. With these choices, we find that the
number of BBHs coalescing in one year, out to $z=20$, is $N_{\rm
BBH}\simeq$7.5\text{\times}{10}^{4}$$.

Several caveats are here in order. Beside the uncertainty in the local rate
(see Mandel & Broekgaarden (2022) for a comprehensive review of the
uncertainties from both observation and theory), there is an uncertainty on
the BBH mass function, for which the LVK data already provide some
information, but which can still vary significantly; even more important are
the uncertainties on the redshift evolution of the merger rate (Dominik et
al., 2013; Santoliquido et al., 2021a; Rozner & Perets, 2022; Chruślińska,
2022) and of the mass and spin distributions. For the redshifts of interest at
3G detectors, these distributions cannot be significantly constrained by
current data, while theoretical modelizations still have large uncertainties.
Our choice of neglecting any redshift dependence in the mass and spin
distributions is the simplest one, and is consistent with Abbott et al.
(2021c), that find no evidence for redshift dependence in the range of
redshifts currently explored by 2G detectors. For the broader range of
redshifts that will be accessible to 3G detectors, this is not expected to
continue to hold, see e.g. Fishbach et al. (2021); van Son et al. (2022);
Belczynski et al. (2022) for the redshift dependence of the mass distribution,
and Qin et al. (2018); Biscoveanu et al. (2022); Bavera et al. (2022) for the
redshift dependence of the spin distribution. Our choice of
redshift–independent mass and spin distributions should therefore be
considered only as dictated by simplicity, and will likely have to be modified
as the observational and theoretical understanding improve. It should also be
stressed that, already for BBHs of astrophysical origin, different formation
channels have merger rates with different redshift dependence; in particular,
BBHs whose progenitors were population III stars have a redshift dependence of
the merger rate which is sensibly different from the one that we have assumed,
and can extend to redshift $z\sim 20$ and beyond, see Kinugawa et al. (2014);
Ng et al. (2021a) and references therein. Furthermore, BBHs with a primordial
origin would have a completely different merger rate, that increases
monotonically with redshift as $R_{\rm
PBH}\propto[t(z)]^{-\nicefrac{{34}}{{37}}}$ up to $z=\order{${10}^{3}$}$
(Raidal et al., 2019; De Luca et al., 2020, 2021a), see also Franciolini
(2021) for recent review.

BNS

The knowledge of the population of merging BNS binaries is still limited
compared to BBH systems, given the very small number of GW detections of this
class of sources. Following Abbott et al. (2021f, c), we assume that the
(source–frame) masses of neutron stars in merging binaries have a flat
distribution in the interval $[1,2.5]\,$\mathrm{M_{\odot}}$$. The other
distribution commonly used in the literature is a Gaussian for each of the two
NSs, such as $\mathcal{N}(\mu=1.33,\sigma=0.09)$ (with masses in units of
$\mathrm{M_{\odot}}$), which comes from galactic electromagnetic (EM)
observations (Farrow et al., 2019). The local rate for BNS mergers inferred
from GW observations is quite sensitive to the assumptions made for the mass
function distribution, and it is therefore important to choose a rate
consistent with the assumed mass distribution. The value inferred from the
GWTC–3 catalog, assuming a flat mass distribution, is ${\cal R}_{0,\rm
BNS}=105.5^{+190.2}_{-83.0}\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ (Abbott et al., 2021c). In
the following, we will then use as reference value ${\cal R}_{0,\rm
BNS}=$105.5\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$. However,
the uncertainty on this number is still quite large and, depending on the
astrophysical modelization used, values of ${\cal R}_{0,\rm BNS}$ in the range
$(10-1700)\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ are consistent with current
observations (Abbott et al., 2021c). We will emphasize this uncertainty
whenever we quote numbers that depend on the local BNS rate. With our choices,
the number of BNSs coalescing in one year, out to $z=20$, is $N_{\rm
BNS}\simeq$1\text{\times}{10}^{5}$$.

Given the expected small values for the spins of these objects, we sample
their components aligned with the orbital angular momentum independently and
uniformly in the range $[-0.05,\,0.05]$. Another parameter characterising NSs
is the adimensional tidal deformability, $\Lambda$, which strongly depends on
the equation of state (EoS) of dense matter above the nuclear density, and is
still largely unknown (Hinderer et al., 2010). We thus make the agnostic
choice of sampling $\Lambda$ for each component uniformly in the range
$[0,\,2000]$, which is compatible with current observations (Abbott et al.,
2019a, 2020a). In App. D we compare with the results obtained by using some
specific NS equations of state. For the redshift distribution of these systems
we assume a Madau–Dickinson profile convolved with a time delay distribution
$P(t_{d})\propto 1/t_{d}$ with a minimum time delay of
$20\text{\,}\mathrm{Myr}$, which is again a standard choice in the literature
[see e.g. Regimbau et al. (2012); Belgacem et al. (2019a)] and, again, no
redshift evolution for the mass distribution.

NSBH

The population of this class of sources is the most uncertain, with only two
GW observations by the time of writing (Abbott et al., 2021g), and no
electromagnetic observation. For the NS we adopt a Gaussian distribution for
the masses, $\mathcal{N}(1.33,0.09)$ (with masses in units of
$\mathrm{M_{\odot}}$), which seems a good approximation to realistic
astrophysical scenarios [see e.g. Fig. 14 of (Broekgaarden et al., 2021)], an
aligned spin uniformly distributed in the interval $[-0.05,\,0.05]$, and an
adimesional tidal deformability parameter uniformly sampled in $[0,\,2000]$.
Both from observations and astrophysical simulations, there seems to be a
preference for low–mass BHs appearing in this systems, also with small spins.
We thus adopt for the BH mass distribution the fitting function provided in
Eq. (6) of Jin-Ping et al. (2021), which is tuned on the astrophysical
simulations of Giacobbo & Mapelli (2018), and shows a peak around $m_{\rm
BH}\sim$8\text{\,}$Msun{\mathrm{M_{\odot}}}$$$, while for the BH aligned spin
component we assume a Gaussian distribution $\mathcal{N}(0,0.15)$ as suggested
in Jin-Ping et al. (2021), which is consistent with the values found in the
analysis of GW200105 and GW200115. The rate distribution is the most
uncertain, as can be seen for example from the various plots and discussion in
Santoliquido et al. (2021b); Broekgaarden et al. (2021). We thus adopt the
same distribution used in the BNS case, i.e. a Madau–Dickinson profile
convolved with a time delay distribution $P(t_{d})\propto 1/t_{d}$ with a
minimum time delay of $20\text{\,}\mathrm{Myr}$. With these choices, the
number of merger per year out to $z=20$ is $N_{\rm
NSBH}\simeq$4.5\text{\times}{10}^{4}$$. We will then study the dependence of
our results on these choices for the BH mass distribution and the event rate,
comparing with some of the models presented in Broekgaarden et al. (2021).

The assumptions for the remaining parameters are common to all the three kinds
of sources: the non–aligned spin components are set to 0, the sky position
parameters $\theta$ and $\phi$ are sampled uniformly over the whole sphere,
the cosine of the inclination angle is sampled uniformly in the interval
$[-1,\,1]$, while the polarisation angle in the interval $[0,\,\pi]$, and
coalescence phase in $[0,\,2\pi]$, and the time of coalescence is sampled
uniformly in time over a 10 year period.181818What really matters in the
analysis is anyway the Greenwich Mean Sidereal Time, GMST, associated to the
GPS time, which, expressed in days, can only range between 0 and 1. The
luminosity distances are computed from the redshifts assuming a flat
$\Lambda$CDM cosmology with Planck18 parameters (Planck Collaboration et al.,
2020), using the astropy.cosmology package (Price-Whelan et al., 2018).

Parameter | BBH | BNS | NSBH
---|---|---|---
$m_{1}$ | Power Law + Peak (Abbott et al., 2021c) | uniform in $[1,\,2.5]$ $\mathrm{M_{\odot}}$ |  (Jin-Ping et al., 2021) Eq. (6)
$m_{2}$ |  $\mathcal{N}(1.33,0.09)$ in $\mathrm{M_{\odot}}$
$z$ | Madau–Dickinson (Madau & Dickinson, 2014) | Madau–Dickinson + $P(t_{d})\propto 1/t_{d}$, $t_{d,{\rm min}}=$20\text{\,}\mathrm{Myr}$$
$d_{L}$ | computed from $z$ assuming Planck18 flat $\Lambda$CDM (Planck Collaboration et al., 2020)
$\chi_{1,z}$ | Default (Abbott et al., 2021c) | uniform in $[-0.05,\,0.05]$ |  $\mathcal{N}(0,0.15)$
$\chi_{2,z}$ | uniform in $[-0.05,\,0.05]$
$\chi_{x},\chi_{y}$ | 0
$\Lambda_{1}$ | 0 | uniform in $[0,\,2000]$ | 0
$\Lambda_{2}$ | uniform in $[0,\,2000]$
$\theta$ |  $\cos(\theta)$ uniform in $[-1,\,1]$
$\phi$ | uniform in $[0,\,2\pi]$
$\iota$ |  $\cos(\iota)$ uniform in $[-1,\,1]$
$\psi$ | uniform in $[0,\,\pi]$
$t_{c}$ | uniform in 10 yr
$\Phi_{c}$ | uniform in $[0,\,2\pi]$
Table 1: Summary of the distributions assumed for all the parameters of the
three classes of compact binary systems considered, refer to the the text for
more details. For BBHs and BNSs the index 1 refer to the heaviest component
and the index 2 to the lightest. For NSBH systems the index 1 refers to the BH
and the index 2 to the NS.

###  4.2 Detector networks

Detector | arms length | latitude $\lambda$ | longitude $\varphi$ | orientation $\gamma$ | arms aperture $\zeta$ | shape | duty cycle
---|---|---|---|---|---|---|---
CE 1 | $40\text{\,}\mathrm{km}$ | $46.5\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $-119.4\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $171\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $90\text{\,}\mathrm{\SIUnitSymbolDegree}$ | L | 85%
CE 2 | $20\text{\,}\mathrm{km}$ | $30.6\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $-90.8\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $242.7\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $90\text{\,}\mathrm{\SIUnitSymbolDegree}$ | L | 85%
ET | $10\text{\,}\mathrm{km}$ | $40.5\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $9.4\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $0\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $60\text{\,}\mathrm{\SIUnitSymbolDegree}$ | Triangle | 85%
LIGO H1 | $4\text{\,}\mathrm{km}$ | $46.5\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $-119.4\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $171\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $90\text{\,}\mathrm{\SIUnitSymbolDegree}$ | L | 70%
LIGO L1 | $4\text{\,}\mathrm{km}$ | $30.6\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $-90.8\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $242.7\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $90\text{\,}\mathrm{\SIUnitSymbolDegree}$ | L | 70%
Virgo | $3\text{\,}\mathrm{km}$ | $43.6\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $10.5\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $115.6\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $90\text{\,}\mathrm{\SIUnitSymbolDegree}$ | L | 70%
KAGRA | $3\text{\,}\mathrm{km}$ | $36.4\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $137.3\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $15.4\text{\,}\mathrm{\SIUnitSymbolDegree}$ | $90\text{\,}\mathrm{\SIUnitSymbolDegree}$ | L | 70%
Table 2: Summary of the positions, orientations, angle between arms, shapes
and duty cycles of the detectors used in our analysis. The orientation denotes
the angle between the bisector of the arms (the first arm in the case of a
triangle) and East.

Our analysis is carried out for three different networks of detectors, which
are also summarised in Tab. 2:

LVK O4

This network consists of the four L–shaped ground–based GW detectors which are
operative at the time of writing, namely the two LIGO detectors, in the U.S.,
in the sites of Hanford and Livingston, Virgo, in Italy, near Cascina, and
KAGRA, in Japan, near Hida. The expected PSDs for all the detectors can be
downloaded from https://dcc.ligo.org/LIGO-T2000012/public. To forecast the
capabilities of the O4 observational run we use the Advanced LIGO sensitivity
with a BNS range of $190\text{\,}\mathrm{Mpc}$, the Advanced Virgo sensitivity
with a BNS range of $120\text{\,}\mathrm{Mpc}$, and the KAGRA sensitivity with
a BNS range of $10\text{\,}\mathrm{Mpc}$, which corresponds to current
expectations for the best sensitivities that could be reached during
O4.191919See https://observing.docs.ligo.org/plan/. For Virgo, current
expectations are rather of a maximum range of $115\text{\,}\mathrm{Mpc}$, but
we use the publicly available PSD, that still corresponds to
$120\text{\,}\mathrm{Mpc}$. It should also be stressed that initial O4
sensitivities will be much lower, with expected ranges of
$160\text{\,}\mathrm{Mpc}$ for LIGO, $80\text{\,}\mathrm{Mpc}$ for Virgo and
1–$3\text{\,}\mathrm{Mpc}$ for KAGRA. To be more realistic, we produce our
results assuming an uncorrelated 70% duty cycle for each detector, as
suggested in Abbott et al. (2020b).

ET

Einstein Telescope (ET) is a proposed third generation GW detector, to be
built in Europe. A candidate site is in the municipality of Lula, in Sardinia,
Italy, and we use this location for definiteness. Very similar results would
be obtained choosing the candidate site in the Meuse–Rhine Euroregion, across
the borders of the Netherlands, Belgium and Germany (basically, the only
difference is in the effect of the Earth’s rotation on the localization
capability of BNSs, that slightly improves increasing the latitude).
Differently from current detectors, ET will be an underground detector and has
a triangular design, consisting of three nested interferometer with
$10\text{\,}\mathrm{km}$ long arms forming angles of
$60\text{\,}\mathrm{\SIUnitSymbolDegree}$. Each interferometer actually has a
xylophone design, meaning that each arm will consist of two separate
instruments, one optimised for high frequencies and one for low frequencies.
For our purposes they can be treated effectively as three equal co–located
detectors rotated by $120\text{\,}\mathrm{\SIUnitSymbolDegree}$ with respect
to each other, neglecting the distance between each arm.202020See
http://www.et-gw.eu/index.php/relevant-et-documents for an collection of
documents on the ET design study and Science Case. The official ET PSD
(corresponding to what was formerly called the ET–D design) can be downloaded
at https://apps.et-gw.eu/tds/?content=3&r=14065. We also assume an
uncorrelated duty cycle of 85% for each arm of the detector, as in Ronchini et
al. (2022).

ET + 2CE

Beside ET alone, we will study a network of three 3G ground–based detectors,
made by ET in Europe, and by two Cosmic Explorer detectors in the U.S. The CE
detectors are planned to be L–shaped interferometers, on the surface (rather
than underground, as ET). Two main configurations are being investigated, one
in which both detectors have $40\text{\,}\mathrm{km}$ arms, and one that
consists of a detector with $40\text{\,}\mathrm{km}$ arms and another with
$20\text{\,}\mathrm{km}$ arms (Evans et al., 2021). The latter, beside its
baseline design (that, for the given length, maximizes the range to compact
binaries), can also be occasionally tuned to have a better sensitivity to the
post–merger phase of BNSs. The network with a $40\text{\,}\mathrm{km}$
detector and a tunable $20\text{\,}\mathrm{km}$ detector appears to maximize
the science output and is currently the reference CE configuration (Evans et
al., 2021). We will then study the case in which one CE detector has
$40\text{\,}\mathrm{km}$ arms and the other has $20\text{\,}\mathrm{km}$ arms
(which we will set in its baseline configuration, that maximizes the range to
compact binaries), and we will refer to the network made by ET and these two
CE detectors as “ET+2CE”.212121The most recent PSDs of CE can be found at
https://dcc.cosmicexplorer.org/CE-T2000017/public. Given the current
uncertainty about their construction sites, for definiteness we assume the
$40\text{\,}\mathrm{km}$ interferometer to be located and oriented as the LIGO
Hanford detector while the $20\text{\,}\mathrm{km}$ instrument as LIGO
Livingston; these are not expected to be the actual locations but, at the
level of the present analysis, this will not be very important, as long as
their relative distance is comparable to the one assumed here. Another option
under consideration, in which one of the two CE detectors could rather be
placed in Australia, would of course lead to better angular localization. In
our analysis we will not consider the two CE detectors alone (except in the
plots showing the range, in Fig. 4 and 5 below), but always a network
consisting of them and ET. Again, we assume an uncorrelated 85% duty cycle for
the two CE detectors and each arm of ET.

The noise spectral densities of the various detectors considered are shown in
Fig. 3. We take into account that ET is made of three nested interferometers,
with an opening angle of $60\text{\,}\mathrm{\SIUnitSymbolDegree}$. In the
${\rm SNR}^{2}$, this gives a factor $3\times(\sqrt{3}/2)^{2}=9/4$, i.e. a
factor $3/2$ in the SNR. Therefore, to compare with the sensitivity of a
single L–shaped interferometer, we multiply the ET sensitivity by a factor
$2/3$. In Fig. 4 we show, for LVK during the O4 run, for ET, and for a network
of two CE (without ET), the corresponding horizon distance to BBHs (i.e., the
maximum distance to which a BBH optimally oriented and with optimal sky
location can be detected), requiring ${\rm SNR}\geq 8$ for detection, where
SNR is the signal–to–noise ratio of the detector network. The lower edge of
the shaded bands gives the distance to which about $50\%$ of the BBHs can be
detected. This has been obtained requiring that the SNR, averaged over sky
position and inclination, is above our threshold value, which is a proxy for a
more accurate computation obtained performing a sampling of a population of
events for each mass bin. The left panel shows the case of equal mass and
non–spinning binaries, and is analogous to a well–known plot presented in Hall
& Evans (2019).222222Actually, the plot in Hall & Evans (2019) refers to two
CEs both of $40\text{\,}\mathrm{km}$, while we consider the
$40\text{\,}\mathrm{km}$+$20\text{\,}\mathrm{km}$ configuration. Note also
that we use the IMRPhenomHM waveform, while the plot in Hall & Evans (2019)
was obtained using IMRPhenomD (we thank Evan Hall for providing this
information). However, in the equal–mass case, higher modes give a negligible
contribution. In the central panel we show the result for non–spinning
binaries with a mass ratio $q=m_{1}/m_{2}=5$, so in this case the fact that we
include higher modes matters for the $50\%$ detection range (but not for the
horizon, which is obtained for $\iota=0$, in which case the higher modes
vanish). Note that, for a given total mass, the horizon decreases by
increasing the mass ratio, simply because, for fixed total mass, the chirp
mass ${\cal M}_{c}$ decreases as the mass ratio moves away from $q=1$, and the
amplitude in the inspiral phase is proportional to ${\cal
M}_{c}^{\nicefrac{{5}}{{3}}}$. In the right panel we also turn on the spins,
taking aligned spins for the two BHs, choosing for definiteness
$\chi_{1,z}=\chi_{2,z}=0.8$, while still keeping $q=5$. We see that the
parallel spins have the effect of raising the horizon distance again, because
of their repulsive effect in the inspiral phase, which delays the merger
(conversely, anti-parallel spins accelerate the merger and lower the horizon
distance).

Figure 3: Noise Amplitude Spectral Densities, ASDs, of the various detectors
considered. To have a uniform comparison between different geometries, the
single–detector ET sensitivity curve has been rescaled by a factor 2/3 to take
into account that, having a triangular design, ET actually consists of 3
nested instruments with an opening angle of
$60\text{\,}\mathrm{\SIUnitSymbolDegree}$. Figure 4: Detection horizons for
different monochromatic populations of BBH sources for the LVK detector
network during O4, for ET alone, and for a network of two CE detectors
(without ET), taken one with 40 km arms and one with 20 km arms. The solid
lines show the maximum redshift out to which a binary with optimal sky
location and inclination can be detected, while the shaded bands represent the
maximum $z$ out to which a binary could be observed after averaging the SNR
over the sky, which is a proxy for the redshift out to which $50\%$ of the
population with given masses and spins could be detected. We set a detection
threshold ${\rm SNR}=8$. In the left panel we show the result for equal mass
non–spinning binaries, in the central panel for non–spinning binaries with a
mass ratio $q=5$, and in the right panel for binaries with a mass ratio $q=5$
and aligned spin components $\chi_{1,z}=\chi_{2,z}=0.8$. We use the
IMRPhenomHM waveform. Figure 5: As in Fig. 4 for NSBH sources. The NS mass is
fixed to $m_{\rm NS}=$1.33\text{\,}\mathrm{M_{\odot}}$$ while we vary the BH
mass in the astrophysically motivated range $[3,\,30]\,$\mathrm{M_{\odot}}$$
(Broekgaarden et al., 2021). In the left panel we show the results for
non–spinning binaries neglecting tidal effects; in the central panel for
binaries having a BH with aligned spin component $\chi_{{\rm BH},z}=0.8$, and
in the right panel we consider both the effect of a spinning BH with
$\chi_{{\rm BH},z}=0.8$ and a tidal deformability parameter for the NS
$\Lambda=360$, which is consistent with the APR4 (Akmal et al., 1998) soft
equation of state. We use the waveform model IMRPhenomNSBH.

For BNS, taking them of equal mass and spinless, which is a very good
approximation for actual systems, and neglecting the small effect of tidal
deformability, the horizon and $50\%$ detection range can also be read from
the left panel of Fig. 4, using a value of $M_{\rm tot}$, such as $M_{\rm
tot}\simeq 2.7$, appropriate to typical BNSs.

For NSBHs, instead, the mass ratio is very different from one, and the BH spin
is not necessarily small. In Fig. 5 we show the results for NSBH using a
horizontal scale for the total mass appropriate to current expectations for
these systems, obtained fixing the NS mass to $m_{\rm
NS}=$1.33\text{\,}\mathrm{M_{\odot}}$$ and varying the BH mass in the
astrophysically motivated range $[3,\,30]\,$\mathrm{M_{\odot}}$$ (Broekgaarden
et al., 2021), and using the IMRPhenomNSBH waveform. In the left panel we show
the case of a spinless BH, and we set to zero the parameter $\Lambda$ that
describes the tidal deformability of the NS. In the central panel we turn on
the BH spin (setting, for definiteness, $\chi_{{\rm BH},z}=0.8$) and, in the
right panel, we also turn on the NS tidal deformability. We see that, for
NSBHs, over the relevant range of masses, the horizon curves are much more
flat than for BBHs. This, first of all, is simply due to the much smaller
range that we have taken for $M_{\rm tot}$, based on current expectations for
NSBH systems. Furthermore, the chirp mass is related to the total mass $M_{\rm
tot}$ and to the mass ratio $q=m_{1}/m_{2}$ by ${\cal
M}_{c}=[q/(1+q)^{2}]^{\nicefrac{{3}}{{5}}}M_{\rm tot}$. Therefore, for $q\gg
1$, as is the case for NSBHs (recall that, with our conventions, $m_{1}\geq
m_{2}$, so $q\geq 1$), we have ${\cal M}_{c}\simeq
q^{-\nicefrac{{3}}{{5}}}M_{\rm tot}$ and a given change in $M_{\rm tot}$
results in a smaller change in ${\cal M}_{c}$, which is the mass scale that
characterizes the amplitude. We also see that, for CE, increasing $M_{\rm
tot}$, beyond some value in the considered range the horizon starts to
decrease; this is due to the fact that, increasing $M_{\rm tot}$, the merger
takes place at lower frequencies, so the signal is moved toward the region
where the sensitivity of CE degrades faster, compared to ET.

##  5 Results

We now show the results of our analysis for the three populations of sources,
as seen by the different detector networks and with the assumptions discussed
above. We first compute, for each event, the network SNR, defined by

${\rm SNR}^{2}=\sum\nolimits_{i}{\rm SNR}^{2}_{i}\,,$ (26)

where the sum runs over the detectors, and the matched filter signal–to–noise
ratios ${\rm SNR}_{i}$ of the individual detectors are given in Eq. (5) or,
more explicitly,

${\rm SNR}^{2}_{i}=4\int_{f_{\rm min}}^{f_{\rm
cut}}\dfrac{|\tilde{h}_{(i)}(f)|^{2}}{S_{n,i}(f)}\differential{f}\,.$ (27)

Here $\tilde{h}_{(i)}(f)$ denotes the GW signal as detected by the $i^{\rm
th}$ detector in the network (thus including the contribution of the pattern
function of the detector and the location phase factor), $S_{n,i}(f)$ is the
one–sided noise PSD of the $i^{\rm th}$ detector, $f_{\rm min}$ denotes the
minimum frequency of the adopted frequency grid, which we set at
$2\text{\,}\mathrm{Hz}$ for ET, $5\text{\,}\mathrm{Hz}$ for CE and
$10\text{\,}\mathrm{Hz}$ for 2G detectors, and $f_{\rm cut}$ is the maximum
frequency of the grid, which depends on the source characteristics and the
adopted waveform model.232323For TaylorF2 this is set to twice the the binary
Innermost Stable Circular Orbit frequency, $f_{\rm
ISCO}=1/(2\pi\,6\sqrt{6}\,GM_{\rm tot}/c^{3})$, while for the chosen full
inspiral–merger–ringdown waveforms, as in LALSimulation, we set the cut
frequency at $(GM_{\rm tot}/c^{3})\,f_{\rm cut}=0.2$. After the SNRs have been
determined, we perform a Fisher analysis restricting to the events having a
network ${\rm SNR}\geq 12$ (in the case of ET, the network SNR is obtained
combining the contributions of the three arms). We will also compare the
results for the number of detections and horizon distances with the results
obtained with a network ${\rm SNR}\geq 8$, while not performing a Fisher
matrix analysis in this case, since it becomes less reliable for such low
values of the SNR.

We further discard the signals with an inversion error of the Fisher matrix
bigger than $5\times 10^{-2}$ (see Sect. 2.2). The Fisher matrices are
computed according to Eq. (6). For a network of detectors, the total Fisher
matrix is just the sum of the Fisher matrices computed for each detector. For
BBHs the set of parameters is ${\bf\it\theta}=\\{{\cal M}_{c},\ \eta,\ d_{L},\
\theta,\ \phi,\ \iota,\ \psi,\ t_{c},\ \Phi_{c},\ \chi_{1},\ \chi_{2}\\}$,
where we denoted the aligned spin components, $\chi_{1,z}$ and $\chi_{2,z}$,
simply as $\chi_{1}$ and $\chi_{2}$, to simplify the notation. For BNS and
NSBH systems, ${\bf\it\theta}$ also includes the tidal deformability
parameters $\tilde{\Lambda}$ and $\delta\tilde{\Lambda}$, defined in Eq. (14),
where, in the case of NSBH, the parameter $\Lambda_{1}$ corresponding to the
BH must be set to zero (recall, from Tab. 1, that for NSBH we use the
convention that the index $i=1$ always refers to the BH). After the inversion
of the Fisher matrix, we compute the sky localization area for the events
according to the definition (Barack & Cutler, 2004; Wen & Chen, 2010)

$\Delta\Omega_{{\rm
X}\%}=-2\pi|\sin\theta|\mathchoice{{\hbox{$\displaystyle\sqrt{\left(\Gamma^{-1}\right)_{\theta\theta}\,\left(\Gamma^{-1}\right)_{\phi\phi}-\left(\Gamma^{-1}\right)_{\theta\phi}^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=9.30444pt,depth=-7.44359pt}}}{{\hbox{$\textstyle\sqrt{\left(\Gamma^{-1}\right)_{\theta\theta}\,\left(\Gamma^{-1}\right)_{\phi\phi}-\left(\Gamma^{-1}\right)_{\theta\phi}^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=9.30444pt,depth=-7.44359pt}}}{{\hbox{$\scriptstyle\sqrt{\left(\Gamma^{-1}\right)_{\theta\theta}\,\left(\Gamma^{-1}\right)_{\phi\phi}-\left(\Gamma^{-1}\right)_{\theta\phi}^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=6.53888pt,depth=-5.23112pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\left(\Gamma^{-1}\right)_{\theta\theta}\,\left(\Gamma^{-1}\right)_{\phi\phi}-\left(\Gamma^{-1}\right)_{\theta\phi}^{2}\,}$}\lower
0.4pt\hbox{\vrule height=5.03888pt,depth=-4.03113pt}}}\
\ln{\left(1-\dfrac{{\rm X}}{100}\right)}\,,$ (28)

where X denotes the confidence level. The result of this expression is in
units of steradian, so a further multiplication by
$($180\text{\,}\mathrm{\SIUnitSymbolDegree}$/\pi)^{2}$ is needed to get the
estimation in the usual ${\mathrm{deg}}^{2}$ units. We will give our results
in terms of $\Delta\Omega_{90\%}$.

###  5.1 Binary black holes

We first focus on BBH systems, which are on average the loudest signals that
can be observed in the frequency band of current and future ground–based
detectors. As recalled in Sect. 4.1, the most recent estimate for the local
rate from the GWTC–3 catalog is ${\cal R}_{0,{\rm
BBH}}=17^{+10}_{-6.7}\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ (Abbott et al., 2021c), so
we use the value ${\cal R}_{0,{\rm BBH}}=17\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ as our reference value for
normalizing the redshift distribution of the merger rate discussed in App. A.
In particular, this means that our merger rate as a function of redshift,
${\cal R}_{{\rm BBH}}(z)$, will be consistent with that shown in Fig. 13 of
Abbott et al. (2021c). For these systems we simulate a population of $N_{\rm
BBH}=$7.5\text{\times}{10}^{4}$$ sources out to $z=20$ that, as mentioned in
Sect. 4.1, using this value for ${\cal R}_{0,{\rm BBH}}$ and our choice for
the redshift distribution of the merger rate, corresponds to the number of
systems coalescing in one year.

Figure 6: Histograms of the SNR and parameter errors for BBHs for ET alone,
using the waveform model IMRPhenomHM, which includes sub–dominant modes of the
signal, and IMRPhenomD, which only contains the quadrupole. For the SNR, in
each bin we show the fraction of events, normalized to the total number of
sources in our sample. For the parameters of the waveform, we limit to the
$\tilde{N}_{\rm det}$ events that pass the cut ${\rm SNR}\geq 12$, and that
furthermore admit a reliable inversion of the Fisher matrix, and we show the
fraction of events, normalized to $\tilde{N}_{\rm det}$. Figure 7: As in Fig.
6, for the ET+2CE network. We here further show the results obtained without
imposing the 85% uncorrelated duty cycle for the detectors.

At the methodological level, it is important first of all to understand the
effect of the waveform model used and, in particular, for heavy systems such
as BBHs, the role of higher modes. We then begin by comparing the results
obtained using IMRPhenomHM, which is our baseline waveform for BBHs and, as
discussed in Sect. 3.1, includes several higher modes, with the results that
we get analysing the same catalog of sources with the waveform model
IMRPhenomD (Husa et al., 2016; Khan et al., 2016a). This is a full
inspiral–merger–ringdown model, tuned to binaries with non–precessing spins
but, differently from IMRPhenomHM, it only contains the dominant quadrupole
mode of the signal. The results are shown in Fig. 6 for ET alone and in Fig. 7
for ET+2CE. In the upper left panels of these figures we show, for the two
waveforms, the distribution in SNR of the whole sample of events, displayed as
a fraction of events with respect to the total number of events in the sample,
i.e. the $7.5\text{\times}{10}^{4}$ sources that we have simulated. We then
restrict parameter estimation to the events that pass the cut ${\rm SNR}\geq
12$, which is our criterion for detection and, for these, we perform the
Fisher matrix analysis. As discussed in Sect. 2.2, for some of these events
the Fisher matrix is ill–conditioned and its inversion can lead to large
amplification of numerical errors, resulting in a covariance matrix that can
be inaccurate even at the $100\%$ level, and our strategy is to simply discard
these events, accepting this as a limitation of the Fisher matrix formalism.
We denote by $N_{\rm det}$ the number of events detected in our sample, i.e.
those that pass the cut ${\rm SNR}\geq 12$, and by $\tilde{N}_{\rm det}$ the
number of detected event that, furthermore, have a Fisher matrix that can be
inverted reliably, according to our criterion discussed in Sect. 2.2. In the
plots showing the distribution of the errors on the parameters we only include
these $\tilde{N}_{\rm det}$ events. The corresponding panels in Fig. 6 and 7,
as well as all similar plots in the following, show the distribution of these
events, as a fraction normalized to $\tilde{N}_{\rm det}$.

From Fig. 6 we see that, for ET alone, while the distribution of SNR is
unaffected by the inclusion of higher modes, parameter estimation is
significantly improved. This is especially remarkable for the angular
resolution and for the error on the luminosity distance, where the inclusion
of higher–order modes allows us to break the distance–inclination degeneracy,
but the estimation of all parameters is significantly improved, thanks to the
better description offered by the waveform model. For all parameters shown,
the inclusion of higher modes has the effect of increasing the fraction of
events for which accurate parameter reconstruction is possible, and also of
cutting the long tails corresponding to events with large errors, as
especially evident in $\Delta\iota$ and $\Delta\eta$. For ET+2CE we see from
Fig. 7 that a rather similar pattern emerges, except for the angular
localization which, in this case, is completely dominated by the
triangulation, rather than by the accuracy of the waveform. Another
interesting information that emerges from Fig. 7 is the role of the duty
cycle. The dashed lines show the results obtained assuming a $100\%$ duty
cycle, while the solid line use a more realistic independent duty cycle of
$85\%$ for each detector. We see that the duty cycle has a large effect on the
tails of the distribution of $\Delta d_{L}/d_{L}$, and of
$\Delta\Omega_{90\%}$. This is natural, since these are the quantities that
are more directly sensitive to the network SNR, which decrease when one or
more detectors in the network is down. This highlights the importance of
having detectors with a high duty cycle, or of having a network of more
detectors.

In the following, we therefore focus only on the results obtained using
IMRPhenomHM. While a histogram of detection fractions contains the information
in its most raw form, the corresponding cumulative distributions give
information that is more condensed and often easier to interpret (and less
sensitive to the specific random realization used, when small numbers are
involved). In Fig. 8 we then show, for LVK–O4, ET alone, and ET+2CE, the
corresponding cumulative distributions function (CDF). For the SNR, we show
the cumulative distribution of the fraction of events, where the fraction is
obtained normalizing the events in a given bin to the total number of events
that we have generated; for the other panels, involving parameter estimation,
we show the cumulative distribution of the fraction of events obtained
normalizing the events in a given bin to $\tilde{N}_{\rm det}$, i.e. to the
total number of detected events which, furthermore, admit a reliable inversion
of the Fisher matrix. For the SNR we actually plot $1-{\rm CDF}$, so the
vertical axis gives the fraction of events with signal–to–noise ratio larger
than a given value, while for the parameters we show the CDF, so the fraction
of events with error smaller than a given value.

Figure 8: Cumulative distributions of the detection efficiency, for the SNRs
and for the error on the parameters, for BBH signals with the three considered
network configuration, using the waveform model IMRPhenomHM.

We find that ET alone can detect $67\%$ of the BBH population, while a network
of ET and 2CE could detect $93\%$ of the full population of sources. With our
reference value for the local rate, ${\cal R}_{0,{\rm
BBH}}=$17\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$, a
threshold in signal–to–noise ratio ${\rm SNR}\geq 12$, and our choices for the
BBH mass function and merger rate distribution in redshift, discussed in Sect.
4.1, this means, for ET alone, that the number of events detected in one year
is $N_{\rm det}\simeq$5.1\text{\times}{10}^{4}$$ while, for ET+2CE, the number
of BBH events detected per year is $N_{\rm
det}\simeq$7.0\text{\times}{10}^{4}$$. Using, more generally, the currently
$90\%$ c.l. allowed range ${\cal R}_{0,{\rm
BBH}}=[10.3,27]\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ for the local rate, the
number of detections per year for ET is in the range
$[$3.1\text{\times}{10}^{4}$,$8.0\text{\times}{10}^{4}$]$ and, for ET+2CE, is
in the range $[$4.2\text{\times}{10}^{4}$,$1.1\text{\times}{10}^{5}$]$. For
LVK in the O4 observational run, we find instead 86 detections per year (which
raise to 141 assuming a duty cycle of $100\%$). Our result for the O4 run is
perfectly consistent with the forecast in Abbott et al. (2020b). Given the
logarithmic scale, the lines in Fig. 8 (in particular, on this scale, the blue
lines that refers to LVK–O4) drop vertically to zero when the accuracy
required on a parameter becomes so small that no more system, in our sample of
detections, satisfies it. The precise value of the parameters when this
happens, as well as the discrete steps apparent in these curves, depend of
course on the specific random realization of our sample.

When performing parameter estimation, we restrict as usual to the
$\tilde{N}_{\rm det}$ detected events that pass our criterion on the inversion
of the Fisher matrix; for ET we find that, with $N_{\rm
det}\simeq$5.1\text{\times}{10}^{4}$$ detection, we have $\tilde{N}_{\rm
det}\simeq$4.8\text{\times}{10}^{4}$$, corresponding to $96\%$ of the detected
events; for ET+2CE, with $N_{\rm det}\simeq$7.0\text{\times}{10}^{4}$$,
$99.2\%$ of the detected events pass the criterion on the inversion, so we
have $\tilde{N}_{\rm det}\simeq$6.9\text{\times}{10}^{4}$$; for LVK–O4,
$N_{\rm det}=\tilde{N}_{\rm det}=86$. We expect that most of the events that
we have discarded because their Fisher matrix is ill–conditioned correspond to
cases in which an analysis based on the full likelihood, rather than on the
Fisher matrix approximation, would anyhow return large errors on the
parameters, so most of the detected events that we have discarded should only
contribute to the tails of the distributions, corresponding to large parameter
errors.

The cumulative detection fraction is a useful metric to evaluate the potential
of GW detectors, in particular to appreciate how a sample of detections is
representative of the whole population, and also has the advantage of being
independent of the value chosen for the local rate, which makes easier the
comparison between different papers in the literature. However, the full
potential of 3G detectors is better appreciated by showing, instead, the
cumulative distribution of the absolute number of detected events, since the
overall normalizations, in particular between 2G and 3G detectors, are very
different. Indeed, in some of the panels of Fig. 8, such as that for $\Delta
d_{L}/d_{L}$, the blue curve for LVK–O4 partly stays above the green curve for
ET alone. This, of course, does not mean that LVK–O4 has a better sensitivity
to $d_{L}$ than ET; rather, it is simply due to the fact that ET sees many
more events, much further away, and some of these events have worse
resolution, so have the effect of decreasing the fraction of detected events
which are accurately measured. In Fig. 9 we show the same plots as in Fig. 8,
but now in terms of the cumulative number of events, rather than the
cumulative fraction of events; for the panel on parameter estimation this is
obtained multiplying the cumulative detection fraction curves for LVK–O4, ET
and ET+2CE, by the respective values of $\tilde{N}_{\rm det}$, while, for the
SNR, we multiply by number of simulated events, $N_{\rm BBH}$. We can then
better appreciate how significantly even ET alone improves on LVK–O4. In
particular, now, for $\Delta d_{L}/d_{L}$, the green curve for ET alone is
well above the blue curve for LVK–O4: for instance, while LVK–O4 is expected
to detect only $\order{10}$ events per year with an accuracy on $d_{L}$ better
than $10\%$, ET alone will detect $\order{10^{3}}$ events per year with
$\Delta d_{L}/d_{L}<10\%$, of which $\order{10^{2}}$ per year will have
$\Delta d_{L}/d_{L}<4\%$. It is also important, in particular for applications
to multi–messenger astronomy and to cosmology, that a single ET improves over
LVK–O4 even on sky localization accuracy. This was not obvious a priori since
a single detector, compared to a network of four widely separated detectors,
cannot exploit triangulation. Nevertheless, the increase in sensitivity of ET,
and therefore its capability to reconstruct all parameters of the signal,
provides a significant improvement, compared to LVK–O4, even in the number of
detected events with a given angular resolution, as we see from the panel for
$\Delta\Omega_{90\%}$ in Fig. 9: for instance, we see that ET alone can reach
an accuracy on $\Delta\Omega_{90\%}$ better than
$100\text{\,}{\mathrm{deg}}^{2}$ on about 2000 events per year (compared to
about 45 events/yr for LVK–O4); better than $10\text{\,}{\mathrm{deg}}^{2}$ on
about 160 events per year (compared to about 15 events/yr for LVK–O4); and can
even reach $\Delta\Omega_{90\%}\leq$1\text{\,}{\mathrm{deg}}^{2}$$ on a few
events per year. An ET+2CE network, combining the sensitivity of 3G detectors
with the long baselines for triangulation, provides further remarkable
improvement on source localization, with about 3400 BBH/yr localized to better
than $1\text{\,}{\mathrm{deg}}^{2}$, and the very best events localized to
less than ${10}^{-2}\text{\,}{\mathrm{deg}}^{2}$.

Figure 9: Cumulative distributions of the number of detections per year, for
the SNRs and for the error on the parameters, for BBH signals with the three
considered network configuration, using the waveform model IMRPhenomHM, and
the value ${\cal R}_{0,{\rm
BBH}}=$17\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ for the
local rate. Except for the panel showing the distribution of the SNR, only the
events whose Fisher matrix could be reliably inverted are included.

Another crucial information, for understanding the scientific potential of 3G
detectors, is provided by the redshift distribution of the detected events,
for different values of the threshold in SNR, or for different selection cuts
on the accuracy that can be obtained on some especially important parameters.
In the left panel of Fig. 10 we show the redshift distributions of the
detected BBHs at ET alone, selected on the basis of the SNR, and we compare it
to the redshift distribution of the total population. Beside showing the
distribution of detected events, with detection defined from the criterion
${\rm SNR}\geq 12$, we also show the distribution in redshift of the events
that satisfy a less stringent cut on the SNR, that we take for definiteness
${\rm SNR}\geq 8$, as well as the events that pass some very high thresholds
on the SNR, while in the central and right panels we show the redshift
distribution of the observed events (defined from the condition ${\rm SNR}\geq
12$), which furthermore pass stringent cuts on the error on the luminosity
distance (central panel) or sky localization (right panel). In this figure,
the upper panels give the fraction of events per redshift bin, and the lower
panels give the corresponding cumulative distributions. A lower threshold on
the SNR can be useful in multi–messenger observations. In general, BBHs are
not expected to have an electromagnetic counterpart; however, should a
counterpart be present for any of those events, the temporal association
between the gravitational and electromagnetic signals could be used to lower
the SNR, while maintaining a high statistically significance.242424BBH systems
could trigger detectable EM emissions if their merger happen in dense
environments (Perna et al., 2016; Mink & King, 2017), but these are expected
to be more difficult to detect, see Palmese et al. (2021b) and references
therein for a recent discussion. For instance, a mechanism that has been
investigated, for producing an electromagnetic counterpart to a BBH, is that
two BHs could merge in the accretion disk of an AGN. The final BH then
receives a kick and moves at high speed through the accretion disk, emitting a
flare (McKernan et al., 2019). In Graham et al. (2020) has been proposed that
the flare ZTF19abanrhr, detected by the Zwicky Transient Facility, could
indeed be associated to the event GW190521, described in Abbott et al. (2020c,
d). See also Mastrogiovanni et al. (2021a); Finke et al. (2021) for the
possibility of testing modified GW propagation with the association of
GW190521 with ZTF19abanrhr, which could also remove possible objections to the
association, related to the difference in their luminosity distance. Another
reason for using a lower SNR, specific to ET, is that, as we have discussed in
Sect. 3.1, the three–arms configuration of ET allows us to construct a null
stream, where the GW signal cancels. This can provide a powerful veto against
glitches and other non–Gaussian noise (Goncharov et al., 2022), that allows us
to dig more confidently into events with a lower SNR. However, one must keep
in mind that the Fisher matrix analysis becomes less reliable for such low
values of the SNR. On the opposite side, events with very large SNR, such as
${\rm SNR}\geq 100$, or events with especially good accuracy on parameter
reconstruction, are obviously interesting for precision physics, such as tests
of GR or measurements of $H_{0}$. In Sect. 5.4 we will discuss these ‘golden
binaries’ further. For the moment, we just observe from Fig. 10 that, with the
detection criterion ${\rm SNR}\geq 12$, ET alone can detect basically $100\%$
of the population up to $z\sim 1$, and about $67\%$ of the BBH population out
to $z\sim 20$; furthermore, about $1\%$ of the population up to $z\sim 10$ is
detected with ${\rm SNR}\geq 100$. Fig. 11 shows the same results for the
ET+2CE network. In this case, basically $100\%$ of the BBH population out to
$z\simeq 2$ is detected, and $93\%$ up to $z\sim 20$.252525Individual BBHs can
be detected to much higher redshifts. In particular, observing a BBH at
$z\gtrsim 30$ , with sufficient accuracy on the redshift reconstruction, would
allow us to establish with good confidence that the event is due to a BBH with
a primordial, rather than astrophysical, origin (De Luca et al., 2021b, a; Ng
et al., 2021b, 2022). The central and right panels show the redshift
distribution of events that pass the detection cut ${\rm SNR}\geq 12$, and
furthermore satisfy more stringent requirement on the accuracy that can be
obtained on $d_{L}$ or on sky localization. We will further comment on them in
Sect. 5.4.

|  |
---|---|---
Figure 10: Redshift distributions of the BBHs observed at ET alone in 1 yr, 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 BBH population in the astrophysical model that we have assumed. We set ${\cal R}_{0,{\rm BBH}}=$17\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). 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 BBH events in our sample, $N_{\rm BBH}=7.5\times 10^{4}$.  |  |
---|---|---
Figure 11: As in Fig. 10, for ET+2CE.

In Fig. 12 we present scatter plots showing the correlations of the errors on
quantities that are relevant in cosmological and astrophysical applications of
GWs. In particular, in the left panels we show the scatter of detections in
the plane $(\Delta\Omega_{90\%},\,\Delta d_{L}/d_{L})$ (with a color scale
giving information on the redshift). We see that, for an ET + 2CE network,
about 36% of the events will have a localization region smaller than
$10\text{\,}{\mathrm{deg}}^{2}$, while also having an accuracy on the
luminosity distance $\Delta d_{L}/d_{L}\lesssim 10^{-1}$, even at $z\gtrsim
1$. This would have an important impact on late–time cosmology, in particular
for the possibility of determining the Hubble constant and constraining the
phenomenon of “modified GW propagation” (Belgacem et al., 2018a, b, 2019b) by
correlating GWs and galaxy catalogs [see Finke et al. (2021); Abbott et al.
(2021e); Palmese et al. (2021a) for the most recent applications]. A small
fraction of the events, of order ${10}^{-3}$ (which still, for ET+2CE, and
with our assumptions on the population and the local rate, corresponds to
about 90 events per year), could even have a localization region
$\Delta\Omega_{90\%}\lesssim${10}^{-1}\text{\,}{\mathrm{deg}}^{2}$$ and
$\Delta d_{L}/d_{L}\lesssim${10}^{-2}\text{\,}$$, which could allow a unique
host galaxy identification even in absence of a direct EM counterpart
(Borhanian et al., 2020), either from a galaxy catalog or from a dedicated
follow–up.

---
Figure 12: Scatter plots of the observed BBH population at the various
networks. In the first row we show the results for the LVK network in O4, in
the middle 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, with the color code representing the redshift, the
central panel the distribution of errors on the luminosity distance and
detector–frame chirp mass, again with the color code representing the
redshift, and the right panel the distribution of errors on the source–frame
chirp mass and symmetric spin $\chi_{s}=(\chi_{1}+\chi_{2})/2$, with the color
code representing the total source–frame mass. The numbers reported on the
dashed lines refer to the fraction of observed events lying inside the
corresponding region.

Very significant results are obtained also with ET alone, despite the
limitations in the localization with a single detector, with a fraction of
events of order $10^{-3}$ (which, in the case of ET alone, and with our
reference value for the local rate, corresponds to about 90 events per year)
having at the same time a sky localization better than
$10\text{\,}{\mathrm{deg}}^{2}$ and an error on the luminosity distance
$\Delta d_{L}/d_{L}\lesssim${10}^{-1}$$. As we mentioned above, ET alone will
see about $10^{3}$ BBH events with $\Delta d_{L}/d_{L}\lesssim${10}^{-1}$$ and
about $160$ BBH events with
$\Delta\Omega_{90\%}\leq$10\text{\,}{\mathrm{deg}}^{2}$$. We see that, among
these, about $90$ satisfy both criteria, $\Delta
d_{L}/d_{L}\lesssim${10}^{-1}$$ and
$\Delta\Omega_{90\%}\leq$10\text{\,}{\mathrm{deg}}^{2}$$, and will therefore
be particularly useful for cosmological studies.

For the current network of detectors in the O4 run we find that it should
observe $\sim 5$ BBHs per year with localization
$\Delta\Omega_{90\%}\lesssim$10\text{\,}{\mathrm{deg}}^{2}$$ and $\Delta
d_{L}/d_{L}\lesssim${10}^{-1}$$, including some sources at $z\gtrsim 0.3-0.4$.
This would allow substantial improvements in the constraints on the Hubble
constant within $\sim 2$ years of data taking, using the correlation with
galaxy catalogs. By comparison, in the GWTC–2 and GWTC–3 catalogs, which
include the detections up to the O3 run, the best localized BBH event,
GW190814, has a sky localization of $19\text{\,}{\mathrm{deg}}^{2}$, and only
two more BBH events, GW190412 ($21\text{\,}{\mathrm{deg}}^{2}$) and
GW200208$\\_$130117 ($30\text{\,}{\mathrm{deg}}^{2}$) have
$\Delta\Omega_{90\%}\leq$30\text{\,}{\mathrm{deg}}^{2}$$ (Abbott et al.,
2021a, b).

In the central panels of Fig. 12 we show our forecasts for the joint errors on
the detector–frame chirp mass and the luminosity distance, again with a color
scale giving information on the redshift. An accurate measurement of
detector–frame masses, together with information on the distribution of
source–frame masses, can be used to break the mass–redshift degeneracy in the
GW waveform and constrain the underlying cosmological model together with the
compact binary population (Farr et al., 2019; Ezquiaga, 2021; Mastrogiovanni
et al., 2021b; Abbott et al., 2021e; Mancarella et al., 2022; Leyde et al.,
2022; Karathanasis et al., 2022).

Finally, in the right column of Fig. 12, we show the expected errors
attainable on the source–frame chirp mass and symmetric spin parameter,
$\chi_{s}=(\chi_{1}+\chi_{2})/2$, with a color code that provides information
on the total source–frame mass. These three parameters are of fundamental
importance to characterise the astrophysical BBH population (Abbott et al.,
2021c) and also to disentangle it from a possible primordial component (De
Luca et al., 2020). Again, for both combinations of parameters, these results
are very promising, showing the potential of the next generation of GW
detectors.

For ease of readability, some selected results presented in this section are
summarised in Tab. 3.

BBH
---
Network | Detected | Analysed | $\rm SNR\geq 100$ | $\Delta d_{L}/d_{L}\leq 10\%$ |  $\Delta\Omega_{90\%}\leq$10\text{\,}{\mathrm{deg}}^{2}$$
LVK–O4 | 86 | 86 | 0 | 8 | 15
ET | 50607 | 48456 | 885 | 1120 | 161
ET+2CE | 69799 | 69610 | 5384 | 45331 | 30889
Table 3: A selection of results from the analysis of the
$7.5\text{\times}{10}^{4}$ BBHs (corresponding to the full population in about
$1\text{\,}\mathrm{yr}$ with our choices for the parameters) at the considered
networks. Here, and in similar tables below, the column labeled “Analysed”
reports the number of events that passed our detection threshold (${\rm
SNR}\geq 12$) and for which, furthermore, we could get a reliable inversion of
the Fisher matrix, according to our criterion explained in Sect. 2.2.

###  5.2 Binary neutron stars

We next present our results for the population of BNS systems. The current
best estimates for the local merger rate, ${\cal R}_{0,\rm BNS}$, is obtained
from the GWTC–3 catalog (Abbott et al., 2021c), but still suffers from a large
uncertainty. Using a rather specific set of assumptions (in particular, a mass
distribution for NS in binaries flat between $1$\mathrm{M_{\odot}}$$ and
$2.5$\mathrm{M_{\odot}}$$), Abbott et al. (2021c) finds ${\cal R}_{0,\rm
BNS}=105.5^{+190.2}_{-83.0}\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$, consistent with the value
${\cal R}_{0,\rm BNS}=320^{+490}_{-240}\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$ inferred from the GWTC–2
catalog using the same mass distribution (Abbott et al., 2021f), but still
with a median value lower by a factor $\sim 3$. However, these numbers are
quite model dependent. Three different models have been investigated in
(Abbott et al., 2021c) and, taking the lowest $5\%$ and highest $95\%$
credible interval out of all three models, ${\cal R}_{0,\rm BNS}$ is inferred
to be in the range $(10-1700)\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$. In the following, we will
use as reference value ${\cal R}_{0,\rm
BNS}=$105.5\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$, which is
the median value obtained for the flat mass distribution that we are using,
but it is important to keep in mind the uncertainty on this number, which
could be up to a factor $\sim 10$ lower, or up to a factor $\sim 16$ higher.
We then simulate a population of $N_{\rm BNS}=${10}^{5}$$ BNS out to $z=20$,
which, for ${\cal R}_{0,\rm
BNS}=$105.5\text{\,}{\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$, and our
assumed redshift dependence for the merger rate, corresponds to the number of
coalescences in about one year.

Figure 13: Cumulative distributions of the event fraction for BNS signals as
function of SNR (for LVK–O4, ET and ET+2CE) and as a function of parameter
errors (for ET and ET+2CE), using the waveform model IMRPhenomD_NRTidalv2.
Figure 14: As in Fig. 13, in terms of the total number of events rather than
detection fraction, using ${\cal R}_{0,\rm BNS}=105.5\,{\rm Gpc}^{-3}{\rm
yr}^{-1}$. Except for the panel showing the distribution of the SNR, only the
events whose Fisher matrix could be reliably inverted are included.

For this value of the local rate, for the current LVK network of detectors
during the O4 run and accounting for the duty cycle, we find that only about 4
events per year should detected, consistently with the forecast in Abbott et
al. (2020b). Even assuming the current detectors to be operational 100% of the
time during O4, the number of detections just raises to 5. Given the smallness
of this sample, for LVK–O4 we do not plot the results for the parameter
estimation of the corresponding events, since they would strongly depend on
the particular random realization considered. Therefore, except for the
cumulative distribution of the SNR, we only show the results for ET and for
ET+2CE.

|  |
---|---|---
Figure 15: Redshift distributions of the BNSs 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 BNS population in the astrophysical model that we have assumed. We set ${\cal R}_{0,{\rm BNS}}=$105.5\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 uncertainly on ${\cal R}_{0,{\rm BNS}}$, which currently can be in the range $(10-1700)\leavevmode\nobreak\ ${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$, one should keep in mind that the absolute number can still change by a factor $\order{10}$ or more. 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 BNS}=${10}^{5}$$. |  |
---|---|---
Figure 16: As in Fig. 15, for ET+2CE.

For 3G detectors, with our threshold ${\rm SNR}\geq 12$ we find $N_{\rm
det}\simeq$9.6\text{\times}{10}^{3}$$ detections per year for ET alone
(corresponding to $9.6\%$ of the population) and $N_{\rm
det}\simeq$4.0\text{\times}{10}^{4}$$ detections per year for ET + 2CE
(corresponding to $40.1\%$ of the population). Varying ${\cal R}_{0,\rm BNS}$
in the currently allowed range $(10-1700)\,{\rm Gpc}^{-3}{\rm yr}^{-1}$, while
keeping fixed the redshift distribution of our population model, the number of
detections per year at LVK–O4 ranges between zero and 65, for ET varies in the
range $[$9.1\text{\times}{10}^{2}$,$1.6\text{\times}{10}^{5}$]$, and for
ET+2CE varies in the range
$[$3.8\text{\times}{10}^{3}$,$6.5\text{\times}{10}^{5}$]$. When performing
parameter estimation, we will restrict again to the $\tilde{N}_{\rm det}$
detected events that also pass our criterion for having a reliable inversion
of the Fisher matrix. We find $\tilde{N}_{\rm
det}\simeq$9.2\text{\times}{10}^{3}$$ for ET (corresponding to about $96\%$ of
the detected events), and $\tilde{N}_{\rm
det}\simeq$4.0\text{\times}{10}^{4}$$ (corresponding to about $99\%$ of the
detected events) for ET+2CE.

Our results for the SNR and for parameter estimation are shown in Fig. 13
(where we show the detection fractions) and in Fig. 14 (where, for parameter
estimation, we multiply the detection fractions at the different detector
networks by the respective values of $\tilde{N}_{\rm det}$, and for the SNR we
multiply by $N_{\rm BNS}$), while in Fig. 15 we show the redshift dependence
of various distributions at ET, similarly to Fig. 10 for the BBH case; in Fig.
16 we show the corresponding results for ET+2CE.

Compared to BBHs, the smaller mass of BNSs allows them to stay longer in the
detector band, but makes them less loud, thus more difficult to see,
especially at high redshifts. As we discussed in Sect. 3.2, thanks to the long
time spent in the detector bandwidth at ET, the BNS localization can be
improved by exploiting the rotation of the Earth. From the panel for
$\Delta\Omega_{90\%}$ in Fig. 13 and 14 we see that, for BNS, with our
fiducial values of the rate and our assumptions on the population, ET as a
single detector can reach an error of order $10\,${\mathrm{deg}}^{2}$$ for the
very best localized sources, and about $0.6\%$ of the detected BNS will have
$\Delta\Omega_{90\%}<${10}^{2}\text{\,}{\mathrm{deg}}^{2}$$.262626Observe that
the CDF for $\Delta\Omega_{90\%}$ in Fig. 13, in the case of ET, does not
saturate to one as $\Delta\Omega_{90\%}$ approaches $4\pi\,{\rm
sr}\simeq$4.1\text{\times}{10}^{4}\text{\,}{\mathrm{deg}}^{2}$$. This is due
to the fact that the Fisher matrix does not know about the physical constraint
that $\Delta\Omega$ must be smaller or equal than $4\pi$ and, for poorly
constrained events, can return a $90\%$ c.l. value larger than that. Indeed,
we see from the plot that, for BNSs at ET alone, only about $25\%$ of the
detected events (and that, furthermore, have passed our criterion on the
inversion of the Fisher matrix) satisfy $\Delta\Omega_{90\%}\leq 4\pi$. The
effect was already present, although less pronounced, for BBHs at ET alone, as
we can see from the panel for $\Delta\Omega_{90\%}$ in Fig. 8. Similar
considerations hold for $\iota$ and $\psi$, for which the Fisher matrix does
not know the constraints $\Delta\iota\leq\pi$ and $\Delta\psi\leq\pi$, and for
the dimensionless spin variables $\chi_{i}$, that range between $-1$ and $1$,
and so should satisfy $\Delta\chi_{i}\leq 2$. As a consequence, the
corresponding CDF, computed with the Fisher matrix, do not necessarily
saturate to one when these boundaries are reached. When comparing with other
results in the literature, it is important to observe that, given the current
large uncertainty on the BNS local rate and population properties, different
papers can use very different assumptions, which can easily lead to results
differing by about one order of magnitude. In particular, in App. B we compare
with the results found in Borhanian & Sathyaprakash (2022) and in Ronchini et
al. (2022), where we find that our results are consistent, once the
differences in the assumptions made (all legitimate, within the current
uncertainties), are taken into account.

It is interesting to observe that, despite the improvement brought by the
inclusion of the Earth’s rotation (which is only relevant for BNSs, given the
long duration of their signal) still, for ET as a single detector, the
localization accuracy of BNSs is not comparable to that of BBHs, that was
shown in the corresponding panel of Fig. 9. For instance, with a single ET
detector a fraction of about $3\text{\times}{10}^{-3}$ of the detections
(corresponding, for our fiducial value of the local rate and our population
model, to $\order{150}$ BBH/yr) can be localized to better than
$10\text{\,}{\mathrm{deg}}^{2}$, while, with our assumptions on the BNS local
rate and BNS population model, we only find 2 BNS that can be localized to
such an accuracy. This shows that the loudness of typical BBH signals is more
important, for sky localization, than the long time spent in the bandwidth by
BNS signals.

For a network ET+2CE, thanks to triangulation, we see again from Fig. 14 that
we can reach an accuracy below ${10}^{-1}\text{\,}{\mathrm{deg}}^{2}$ for the
best localized BNS systems, to be compared with the localization below
${10}^{-2}\text{\,}{\mathrm{deg}}^{2}$ that ET+2CE can reach for the best
localized BBHs, as we saw in Fig. 9. However, it is important to notice that,
for multi–messenger observations of BNSs and for the identification and
follow–up of the associated kilonova, already the angular resolution of a
single ET detector can be adequate. For instance, to understand the
nucleosynthesis spectra of the kilonova, for sources up to $z\simeq 0.3-0.4$
the best instrument is the Extremely Large Telescope
(ELT),272727https://elt.eso.org. planned to start observations in 2027, and
that could still be operational by the time of 3G GW detectors (Spyromilio et
al., 2008; Rossi et al., 2020). However, a direct pointing with ELT would
require arcsec localization, which is anyhow out of question for GW detectors.
The actual, strategy for sources at these moderate redshifts, is rather to use
telescopes with large field–of–view (FOV): these instruments, with a
localization of the GW events of the order of tens to hundreds of square
degrees, can indeed localize the source. Notice that such a GW resolution can
be given already by ET alone. In particular, the Vera Rubin Observatory’s LSST
(Ivezić et al., 2019) has a FOV of $9.6\text{\,}{\mathrm{deg}}^{2}$ and can
observe kilonovae up to $z\sim 0.1$ (with $5\%$ of the kilonovae observable up
to $z\simeq 0.4$) and several other instruments, such as ULTRASAT (Sagiv et
al., 2014), can reach $z\simeq 0.1$; the Nancy Roman Space Telescope (formerly
WFIRST) (Spergel et al., 2015) has the highest reach, being able to observe
$50\%$ of the kilonovae up to $z\sim 0.2-0.3$, with $5\%$ of the kilonovae
observable up to $z\simeq 1$, although its FOV,
$0.28\text{\,}{\mathrm{deg}}^{2}$, is not as large, so it requires sub–degree
localization, or an earlier localization by instruments with a larger FOV (see
Cowperthwaite et al. (2019); Chase et al. (2022); Ronchini et al. (2022) for
recent discussions). Once localized the kilonova with these large FOV
instruments, telescopes such as the ELT can perform a more detailed follow–up
of the source. For BNSs at larger redshifts, the identification of the
electromagnetic counterpart can only be made by ${\rm X}/\gamma$–ray
satellites with large FOV; a sky localization of order
$100\text{\,}{\mathrm{deg}}^{2}$, as can be provided by ET alone, is about
$1/10$ of the typical FOV of wide–field X–ray telescope, which in general is
larger than $1\text{\,}\mathrm{sr}$, so such a localization can already
provide sufficient information for the search and localization of the
$\gamma$– or X–ray counterpart; then, the ${\rm X}/\gamma$–ray satellites can
provide the arcmin localization needed to drive the ground–based follow–up
(Ronchini et al., 2022). For these multi–messenger studies, a lower threshold
on the SNR, such as ${\rm SNR}=8$, can be appropriate, since the statistical
significance will be enhanced by the temporal coincidence of the gravitational
and electromagnetic signals. From the upper–left panel of Fig. 15 we see that,
for ET alone and a threshold ${\rm SNR}=8$, the farthest detected BNSs in our
sample are at $z\sim 4.5$ (the farthest one happens to be at $z\simeq 4.77$
and has $M_{\rm tot}=$4.15\text{\,}\mathrm{M_{\odot}}$$ and nearly equal
masses). Higher values of the rate, or longer observation times, would lead to
the gradual appearance of more and more rare events. For ${\rm SNR}=8$, at ET
alone, an equal mass non–spinning BNS, with
$m_{1}=m_{2}=$1.35\text{\,}$Msun{\mathrm{M_{\odot}}}$$$, optimal inclination
and sky location could be detected up to $z_{\rm max}\simeq 3.5$, while a
heavier equal mass non–spinning BNS, say with
$m_{1}=m_{2}=$2.1\text{\,}$Msun{\mathrm{M_{\odot}}}$$$, optimal inclination
and sky location could be detected up to $z_{\rm max}\simeq 6.9$.282828Because
of the long duration of the BNS signal and the Earth movement, the sky
location changes in time. In practice, we find that the optimal sky location
is such that the BNS is near a maximum of the pattern function at the time of
merger. These estimates are in agreement with those in Ronchini et al. (2022).
For ET+2CE, the farthest detection in our sample, with ${\rm SNR}\geq 8$, is
at $z\simeq 8.8$ and has $M_{\rm tot}=$3.73\text{\,}\mathrm{M_{\odot}}$$ and a
mass ratio of $q\simeq 1.33$. An ideally oriented BNS with the same (rather
extreme) characteristics could in principle be detected up to $z\simeq 20$
(where, of course, NSs are not even expected to exist), while a more realistic
equal–mass non–spinning BNS, with
$m_{1}=m_{2}=$1.35\text{\,}$Msun{\mathrm{M_{\odot}}}$$$, could be detected up
to $z\simeq 11$.

---
Figure 17: Scatter plots of the observed BNS 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 redshift, and the right panel the
distribution of errors on the source–frame chirp mass and adimensional tidal
deformability 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.

When drawing conclusions from plots featuring the absolute number of BNS
detections, it is important to keep in mind that the estimate of ${\cal
R}_{0,{\rm BNS}}$ from current observations is still very uncertain, given the
limited sample of detected BNS coalescences, much more than for BBHs, for
which we already have a much larger sample of detections. As we have mentioned
above, current data, depending on the model used (e.g. on the assumed mass,
spin and redshift distributions) are consistent with values of ${\cal
R}_{0,{\rm BNS}}$ in the range $(10-1700)\leavevmode\nobreak\
${\mathrm{Gpc}}^{-3}\text{\,}{\mathrm{yr}}^{-1}$$. Therefore, as more and more
BNS detections will accumulate with future runs, and the BNS local rate will
become better constrained, the absolute scale on the vertical axis in Fig. 14
could still change, possibly by as much as a factor $\order{10}$. Furthermore,
a change in the rate is not just reabsorbed in a rescaling of the absolute
vertical scale of these plots. For a fixed observation time, such as the
standard reference choice of one year, the tails of the distributions, i.e.
rare and particularly interesting events, change. As an obvious example, if,
with a given choice of rate, in one year there is no event with, say, $z>4$,
this does not mean that, with a rate three times larger, there will be three
times zero events, which is still zero. Rarer and rarer events will appear as
we increase the rate or the observation time. Therefore, given the current
large uncertainty, when comparing different results for BNSs in the