# Systematic Biases in Estimating the Properties of Black Holes
Due to Inaccurate Gravitational-Wave Models

Arnab Dhani<EMAIL_ADDRESS>Max Planck Institute for Gravitational
Physics (Albert Einstein Institute), Am Mühlenberg 1, Potsdam 14476, Germany
Sebastian H. Völkel  Max Planck Institute for Gravitational Physics (Albert
Einstein Institute), Am Mühlenberg 1, Potsdam 14476, Germany  Alessandra
Buonanno  Max Planck Institute for Gravitational Physics (Albert Einstein
Institute), Am Mühlenberg 1, Potsdam 14476, Germany  Department of Physics,
University of Maryland, College Park, MD 20742, USA  Hector Estelles  Max
Planck Institute for Gravitational Physics (Albert Einstein Institute), Am
Mühlenberg 1, Potsdam 14476, Germany  Jonathan Gair  Max Planck Institute for
Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, Potsdam
14476, Germany  Harald P. Pfeiffer  Max Planck Institute for Gravitational
Physics (Albert Einstein Institute), Am Mühlenberg 1, Potsdam 14476, Germany
Lorenzo Pompili  Max Planck Institute for Gravitational Physics (Albert
Einstein Institute), Am Mühlenberg 1, Potsdam 14476, Germany  Alexandre
Toubiana  Max Planck Institute for Gravitational Physics (Albert Einstein
Institute), Am Mühlenberg 1, Potsdam 14476, Germany

###### Abstract

Gravitational-wave (GW) observations of binary black-hole (BBH) coalescences
are expected to address outstanding questions in astrophysics, cosmology, and
fundamental physics. Inference of BBH parameters relies on waveform models,
and realizing the full discovery potential of upcoming LIGO-Virgo-KAGRA
observing runs and new ground-based facilities (such as the Einstein Telescope
and Cosmic Explorer) hinges on the accuracy of these waveform models. Using
linear-signal approximation methods and Bayesian analysis, we start to assess
our readiness for what lies ahead using two state-of-the-art quasi-circular,
spin-precessing models: SEOBNRv5PHM and IMRPhenomXPHM. We find that systematic
biases increase with the spin of the BH, with parameter biases being
approximately 6 to 8 times likelier, if the primary-spin magnitude exceeds 0.5
compared to when it is less than 0.5. Additionally, we ascertain that current
waveforms can accurately recover the distribution of masses in the LVK
astrophysical population, but not spins. Upon exploring the broader parameter
space of BHs, we find that systematic biases increase with detector-frame
total mass, binary asymmetry, and spin-precession, with a majority of such
binaries incurring parameter biases, extending up to redshifts $\sim 3$ in
future detectors. Furthermore, we examine three “golden” events characterized
by mass ratios of approximately 6 to 10, significant spin magnitudes
($0.6\mbox{--}0.9$), and high precession, evaluating how systematic biases may
affect their scientific outcomes. Our findings reveal that current waveforms
fail to enable the unbiased measurement of the Hubble-Lemaître parameter and
sky localization from loud signals, even for current detectors. Moreover,
highly asymmetric systems within the lower BH mass-gap exhibit biased
measurements of the secondary-companion mass, which impacts the physics of
both neutron stars and formation channels. Similarly, we deduce that the
primary mass of massive binaries ($>60M_{\odot}$) will also be biased,
affecting supernova physics. Future progress in analytical calculations and
numerical-relativity simulations, crucial for calibrating the models, must
target regions of the parameter space with significant biases to develop more
accurate models. Only then can precision GW astronomy fulfill the promise it
holds.

##  I Introduction

Almost a decade ago, the first observation of a gravitational wave (GW) from
the coalescence of two black holes marked an important milestone in the
history of GW astronomy [1]. Since then, the LIGO-Virgo-KAGRA (LVK)
Collaboration [2, 3, 4] has detected more than 90 compact binary mergers [5,
6, 7], and independent research groups [8, 9, 10, 11, 12, 13] have discovered
additional events. Thus, GWs have become a novel tool to explore the Universe.
The observed signals have been used to measure the mass and spin distributions
of BHs and neutron stars, their formation channels, and the co-evolution of
their properties with that of the Universe [14, 15]. BNS mergers have improved
the bounds on the nuclear equation of state and the maximum allowed mass of a
NS [16, 17, 18]. The mass distributions have been employed to constrain the
observed lower and predicted upper-mass gaps and other features in the mass
spectrum. In conjunction with the electromagnetic (EM) counterparts observed
for GW170817, or together with available galaxy catalogs, they have also been
used to constrain the Hubble-Lemaître parameter ($H_{0}$) [19, 20]. GW
measurements have also probed general relativity (GR) as the fundamental
theory of gravity [21, 22, 23].

Improvements in the sensitivity of current GW detectors and proposed next-
generation (XG) ground-based observatories like the Einstein Telescope (ET)
and Cosmic Explorer (CE) [24, 25, 26], will significantly increase the
observational volume, and with it the number of GW sources. For instance, a
network of XG detectors will observe _every_ stellar-origin BBH merger and
most BNS mergers across the observable Universe [27]. A number of studies have
explored in detail the extent to which various science objectives can be
accomplished [28, 29, 30, 31, 32, 27]. With $\mathcal{O}(100)$ observations by
the LVK Collaboration, and $\mathcal{O}(10^{5})$ promised detections with XG
detectors, meaningful inferences on the properties of the astrophysical
distribution of BBHs will constrain more and more the underlying distribution
of main sequence stars and their evolution. EM observations in our galaxy
indicate that stellar-origin BHs have masses above $5\,\rm M_{\odot}$.
However, these observations may be biased by properties that are unique to our
galaxy. Similarly, the pair-instability supernova (PISN) process is expected
to suppress BH formation in the mass range $\sim 50\mbox{--}120\,\rm
M_{\odot}$ [33, 34]. Confident detections of BBHs in these mass ranges would
pose challenges to stellar-evolution models, as well as, constrain the
$\prescript{12}{}{C}(\alpha,\gamma)\prescript{16}{}{O}$ reaction rate that
drives the PISN process [35]. Gravitational-wave astronomy can also determine
the cosmological evolution of the Universe. In particular, by combining many
GW signals, it will contribute to resolving the $H_{0}$ tension, and provide
new constraints on structure formation. On the other hand, loud individual
events carry a lot of information too. Individual “golden” BBHs can also
resolve the Hubble-Lemaître tension [36]. Finally, precision tests of GR can
be derived from high SNR observations [37, 38, 39, 40, 41, 42, 43]. All of
these scientific objectives are vulnerable to false positives arising from
waveform inaccuracies.

The source properties are estimated from the GWs via Bayesian inference using
waveform models predicted by GR. Since there is no complete, closed-form
analytic solution for the gravitational waveform of a compact-binary
coalescence, various approximate and numerical methods have been developed to
describe the GW signal faithfully. The effective-one-body (EOB) waveforms [44,
45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58] can combine and resum
several perturbative results, such as post-Newtonian (PN), post-Minkowskian
(PM) and gravitational self-force information for the conservative and
dissipative dynamics, with physically motivated ansatze for the merger, and BH
perturbation theory for the ringdown. They are made highly accurate through
calibration to numerical relativity (NR) simulations [59, 60, 61]. Fast and
accurate inspiral-merger-ringdown phenomenological (IMRPhenom) models [62, 63,
64, 65, 66, 67] are built fitting EOB, PN and NR waveforms. NR simulations
give the most accurate representation of a GW signal although they are still
limited by numerical truncation error [68, 69, 70], imperfect outer boundary
conditions [71, 72, 73] and issues with GW extraction and extrapolation [74,
75]. Moreover, NR simulations are not available in the entire parameter space,
and are limited in length due to their high computational cost. NR surrogate
models (NRSur) [76, 77, 78, 79] are constructed by directly interpolating NR
waveforms, where available.

Figure 1: The GW strains for a BBH system with parameters given in Table 1
(_Binary 1_) at the LIGO-Livingston detector of the _O5_ network. The black
curve is the injected signal SEOBNRv5PHM, the green curve is the template
IMRPhenomXPHM, evaluated for the injection parameters, and time shifted and
global-phase rotated to maximize their overlap with the signal, while the
brown curve is the template IMRPhenomXPHM evaluated at the maximum likelihood
values obtained using a Bayesian analysis. The reference for the time axis,
$t=0$, is taken to be the peak of the GW multipole $h_{22}$ of the signal.

Thanks to advancement in GW modeling since the discovery of GW150914 [80],
waveform models have been sufficiently accurate to analyze most signals in the
LVK GW Transient Catalogs (GWTC) [80, 81]. In Ref. [82], the authors used the
absolute value of the difference between waveform models to quantify the
accuracy of a given pair of models, finding that a few high signal-to-noise
ratio (SNR) events in GWTC-3 and GWTC-2.1. fail their criterion. They also
find that parameter estimation of such events show greater inconsistencies. A
reanalysis of the GWTC-3 catalog by Ref. [83] finds that the NRSur7dq4 model
recovers noticeably different parameters compared to LVK analyses using
IMRPhenomXPHM and SEOBNRv4PHM waveform models for $\sim 20\%$ of the events
where NRSur7dq4 model can be used 111Note that some differences are likely to
be attributed to sampler issues rather than waveform systematics.. A
hypermodel approach to identify waveform systematics has also been carried out
on the 13 heaviest GW events from the GWTC-3 catalog. In this approach,
waveform models are treated as parameters and directly sampled over, yielding
a direct probability for each waveform model. The authors do not find any
waveform model to be preferred except for 3 events which are marred by data
quality issues [84]. Recently, there have also been efforts to marginalize
over waveform modeling uncertainties [85, 86]. Other studies have found that
even relatively low SNR events could be affected by systematic biases if they
lie in a region of the parameter space where calibration with NR is sparse.
This would include binaries that are asymmetric, eccentric, have large spin
magnitudes, and/or have precessing orbits [51, 82].

With increasing detector sensitivity and number of detections, the median SNR
of the observed population of binaries, as well as the likelihood of detecting
a binary from a region of the parameter space where waveform inaccuracies are
greater, will increase. While statistical uncertainties decrease with
increasing SNR, systematic biases are independent of the signal power. Several
studies have explored the validity of waveform models for the parameter
estimation of quasi-circular binary black hole (BBH) mergers in upgraded and
XG detectors, mainly focusing on the biases for individual events [87, 81, 82,
223], with Ref. [81] also showing the inferred distribution of the primary
mass to be biased. While the negligence of subdominant modes can significantly
bias the parameter estimation of individual events [88, 89], a recent study
indicated that such biases do not affect the inference of the LVK-like
astrophysical distribution of BBHs [90]. Other studies have focused on
waveform systematics in the presence of eccentricity [91, 92], matter effects
and spin-precession [93, 94, 95, 96, 97]. Recent studies have also explored
the effect of truncation errors in NR simulations employing finite
differencing methods and concluded that current simulations are not accurate
enough for highly asymmetric binaries and binaries whose orbits are inclined
with respect to the line of sight [98, 99]. However, state-of-the-art waveform
models, such as SEOBNRv5PHM and IMRPhenomXPHM, are calibrated to the
Simulating-eXtreme-Spacetimes (SXS) Collaboration waveforms, which employ
spectral methods and the effect of truncation errors on these waveforms have
not been explored systematically. An indistinguishability criterion [100] has
also been used as an easy-to-compute metric to determine accuracy requirements
of waveforms [101, 81]. However, this measure has been found to be very
conservative. Reference [102] proposes a correction to it to improve the
reliability of the measure.

We illustrate the effect of waveform mismodelling in Fig. 1, using a BBH with
parameters given in Table 1 222We refer the reader to Sec. II.1, Sec. II.2,
and Sec. III.2 for discussions on GW parameters, waveform models, and
maximization of overlaps between waveforms.. We show the multipolar, spin-
precessing GW strains in the LIGO-Livingston detector from the SEOBNRv5PHM
waveform model [51] as signal (black curve) and the IMRPhenomXPHM waveform
model [65] as template (brown and green curves). For the green curve, we fix
the polarization angle and time at coalescence, by maximizing its overlap
against the SEOBNRv5PHM signal. We employ the LIGO-Virgo detectors assuming
the sensitivity of the upcoming fifth observing (_O5_) run. If the green GW
strain faithfully represented the signal (black), they would perfectly match
throughout the coalescence. However, this is not the case, the amplitude
modulations are different during the long inspiral, and in particular during
the late inspiral, merger and ringdown. Furthermore, while the signal and the
template phases match during the early inspiral, there is significant
dephasing near the late inspiral, merger and ringdown. In Fig. 1, we also show
the IMRPhenomXPHM template (brown curve) evaluated at the maximum-likelihood
parameters (obtained through a Bayesian analysis). It has a much better match
to the SEOBNRv5PHM signal even during the late inspiral, merger and ringdown.
This best match is obtained at the expense of introducing a bias in the
parameters, notably the total mass, mass ratio, and the spin-precession
parameter are biased by $\sim 3\%$, $\sim 6\%$ and $\sim 13\%$, respectively.
The brown curve also has an associated brown band representing the measurement
errors at 90% credible interval in the GW parameters, but it is barely visible
to the naked eye, illustrating that this uncertainty, which represents the
estimated statistical uncertainty from instrumental noise is much smaller than
the waveform difference between the signal and the template evaluated at the
best-fit parameters. As previously stated, this inconsistency manifests itself
as biased parameter estimation, which could affect the various science
objectives.

In this work, we start to quantify the systematic biases that can be expected
in future observing runs with current facilities and XG detectors using the
SEOBNRv5PHM and IMRPhenomXPHM waveform models, which are employed for
parameter-estimation studies of BBHs by the LVK Collaboration. Both models are
valid for quasi-circular binaries and incorporate subdominant spherical
harmonics and spin-precession effects. While it would be ideal to quantify the
biases of each of these models against the true GR signal, it is infeasible to
do it everywhere in the parameter space since NR waveforms are not available.
We leave to a future study the use of NR waveforms as synthetic signals where
available. Throughout this paper, we instead generate GW signals using
SEOBNRv5PHM, considering these to represent the true signal, and analyze them
using IMRPhenomXPHM. However, there is a drawback to this approach. If both
the waveform models deviate in a similar way from the true GR signal, the
present analysis would predict small biases even when the true bias is large.
This is especially true since the two waveform models are not completely
independent. IMRPhenomXPHM uses the SEOBNR waveforms (although from a previous
version, i.e., SEOBNRv4) for calibration in parts of the parameter space where
there is a dearth of NR simulations. This is precisely the regions where
systematic biases are expected to be more common. In this sense, our analysis
is a conservative assessment of the prevalence of systematic biases.

To quantify the systematics of the aforementioned waveform models in a wide
range of applications, we utilize Bayesian analysis as well as the linear-
signal approximation (LSA). The former is the most reliable tool to obtain the
posterior distribution for a GW signal, but computationally expensive. The
latter allows for computational efficiency, but approximate the predictions
for the posterior properties, including systematic biases, and should become a
good approximation only at large SNR. We use the LSA to study biases for BBH
populations, and a wider parameter space, which is not feasible with
conventional Bayesian methods. We consider three detector networks comprising
of the current LIGO-Virgo network at design sensitivity (_O5_), a planned
network where the current LIGO detectors are upgraded to improved sensitivity
(_A#_), and a XG network comprising of two CE and an ET. The BBH populations
we consider follow the LVK-like distributions where the binary masses are
distributed as determined by LVK while the spins are assumed to be
isotropically oriented and distributed uniformly in magnitude. We do this to
allow for a wider range of spins. The binaries extend up to a redshift of 3
following the Madau-Dickinson star formation rate (SFR) [103]. Following this,
we embark on a parameter-exploration study where we consider large redshifted
total masses of $200\rm\,M_{\odot}$, asymmetric systems with inverse mass
ratios going up to 30, and highly spin-precessing systems. We study these as-
yet unobserved regions of the parameter space in anticipation of future
observations. We also consider three distinct prototypes of BBH mergers, which
hold great potential for various science objectives, but are non-trivial to
model due to precession or large mass ratio. The details of these three
“golden” binary systems can be found in Table 1.

The paper is organized as follows. In Sec. II, we introduce the main
characteristics of the GW signal and its parameters, the waveform models that
we use, and the detector networks in which signals are simulated. In Sec. III,
we describe our methodology comprising of Bayesian analysis and LSA. We point
out the importance of having consistent parameter definitions across waveform
models, and its impact on the systematic bias, where we show a comparison of a
Bayesian analysis with the estimates from LSA. We also discuss the limitations
of the LSA for parameter estimation (notably the Fisher information matrix)
and biases. The study of systematic biases in the LVK-like BBH population and
a hierarchical Bayesian inference on parameter distributions, reweighted to
the LVK population, is reported in Sec. IV. A much broader study across the
binary parameter space with particular focus on massive, highly asymmetric and
spin-precessing binaries is reported in Sec. V. A ramification on the
different science applications for GWs can be found in Sec. VI where we study
selected GW events or “golden” binaries. The discussion and conclusion can be
found in Sec. VII. In Appendix A, we illustrate the effect of nonuniform-
parameter definitions across waveform models on the estimates of the
systematic bias through a toy model. In Appendix B, we discuss the effect of
different harmonics of the EOB model starting at different frequencies. In
Appendix C, we discuss the effect of the starting frequency of the analysis on
parameter estimation and systematic biases. In Appendix D, we provide a
complimentary plot to Fig. 6 by reporting the dependence of ratio of
systematic bias to statistical error as a function of the SNR. In Appendix E,
we show the effect of the SNR threshold on the distribution of the population
parameters. In Appendix F, we report the bias horizon for the $\chi_{1}$
parameter of the exploratory binaries of Sec. V.

##  II Gravitational-wave parameters, models and detectors

###  II.1 Gravitational-wave parameters

We are interested in estimating the properties of quasi-circular, spin-
precessing BBHs observed with current and future ground-based detector
networks. The GW strain emitted by such binaries is characterized by 15
parameters. The parameters intrinsic to the source are the component masses,
$m_{i}$ 333We adopt the convention $m_{1}\geq m_{2}$. and the dimensionless
spin vectors, $\bm{\chi}_{i}=\bm{S}_{i}/m_{i}^{2}$ $(i=1,2)$. The position of
the binary is described by its luminosity distance, $D_{L}$, and the
coordinates on the plane of the sky, $(\alpha,\delta)$. The orientation of the
binary is described by the polar angle, $\iota$, and the azimuthal angle,
$\varphi$, to the observer in the source frame [104] at the reference
frequency, $f_{\rm ref}$, which we set to $f_{\rm ref}=20\rm\,Hz$ throughout
this paper. Finally, the relative contribution of the two gravitational
polarizations, $h_{+}(t)$ and $h_{\times}(t)$, is described by the
polarization angle, $\psi$, while the reference for the time is given by the
coalescence time, $t_{c}$. With these definitions, the GW strain can be
expressed as

$\begin{split}h(t)=&h_{+}(t;m_{z,i},\bm{\chi}_{1,2},D_{L},\iota,\varphi,t_{c})F_{+}(\alpha,\delta,\psi)\\\
&+h_{\times}(t;m_{z,i},\bm{\chi}_{1,2},D_{L},\iota,\varphi,t_{c})F_{\times}(\alpha,\delta,\psi)\,,\end{split}$
(1)

where $F_{+,\times}(\alpha,\delta,\psi)$ are the antenna pattern functions
[105, 106]. The detector- and source-frame masses, $m_{i}^{z}$ and $m_{i}$,
respectively, are related by $m_{z,i}=m_{i}(1+z)$ with $z$ being the redshift
of the source. A superscript on any mass parameter indicates that it is
detector frame while its absence indicates it is source frame. The parameters
$D_{L}$ and $z$ are related for a given cosmological model, which we take to
be the one from Planck18 [107]. The two GW polarizations can be decomposed in
the basis of $-2$ spin-weighted spherical harmonics,
$\prescript{}{-2}{Y}_{lm}$, as

$h_{+}(t)-ih_{\times}(t)=\sum_{l=2}^{\infty}\sum_{m=-l}^{+l}\prescript{}{-2}{Y}_{lm}(\iota,\varphi)h_{lm}(t)$
(2)

where $h_{lm}(t)$ are the GW multipoles and $\varphi=\pi/2-\phi_{\rm ref}$.

It is often helpful to express the GW signal in terms of parameters that are
combinations of the component masses and spins, either because they appear in
such combinations in PN expressions or because they are conserved up to
certain PN orders. In particular, the chirp mass, $\mathcal{M}_{c}$, and the
symmetric mass ratio, $\nu$, are defined by
$\mathcal{M}_{c}=(m_{1}m_{2})^{3/5}/M^{1/5}$ and $\nu=(m_{1}m_{2})/M^{2}$,
respectively, where $M=m_{1}+m_{2}$ is the total mass. The effective spin,
$\chi_{\rm eff}$, and spin-precession, $\chi_{\rm p}$, parameters are given by

$\displaystyle\chi_{\rm eff}$
$\displaystyle=\frac{m_{1}\chi_{1z}+m_{2}\chi_{2z}}{m_{1}+m_{2}}\,,$ (3a)
$\displaystyle\chi_{\rm p}$
$\displaystyle=\frac{1}{B_{1}m_{1}^{2}}\max{(B_{1}m_{1}^{2}\chi_{1,\perp},B_{2}m_{2}^{2}\chi_{2,\perp})}\,,$
(3b)

where $B_{1,2}=2+3\,m_{2,1}/m_{1,2}$, and $\chi_{i\perp}$ and $\chi_{iz}$ are
the magnitudes of the projection of $\bm{\chi}_{i}$ on to the orbital plane
and perpendicular to it, respectively. While alternative definitions of
$\chi_{\rm p}$ have been proposed [108, 109], the GW–Bayesian-analysis package
Bilby that we use to analyze simulated signals uses Eq. (3b).

When transforming to a spherical coordinate system with the $z$-axis
perpendicular to the instantaneous orbital angular momentum, the tilts of the
two spin vectors with respect to the $z$-axis, $\theta_{i}$, are given by

$\cos\theta_{i}=\frac{\chi_{iz}}{\chi_{i}}\,,$ (4)

where $\chi_{i}\equiv|\bm{\chi}_{i}|$ are the magnitudes of the dimensionless
spin vectors. The relative angle between them in the orbital plane is
parameterized by $\phi_{12}=\phi_{1}-\phi_{2}$ where $\phi_{1,2}$ are the
azimuthal angles of the two spin vectors in the spherical coordinates.
Finally, the direction of the total angular momentum, $\bm{J}$, in the plane
perpendicular to the orbital angular momentum, $\bm{L}$, at some reference
time is given by the parameter $\phi_{\rm JL}$. Since
$\bm{J}=\bm{L}+\bm{S}_{1}+\bm{S}_{2}$, $\phi_{\rm JL}$ also defines the
direction of the total spin vector in the orbital plane. The total angular
momentum also defines the angle, $\theta_{\rm JN}$, which gives the
orientation of the total angular momentum vector relative to the line-of-
sight, $\bm{N}$, of the observer. The angle $\theta_{\rm JN}$ can be expressed
in terms of the inclination angle, $\iota$, at the reference frequency,
$f_{\rm ref}$, e.g., through Eq. (C9) of Pratten  _et al._ [65]. In summary,
the waveform depends on the following 15 parameters:

$\bm{\vartheta}=\\{\mathcal{M}_{c},\nu,\chi_{i},\cos\theta_{i},\phi_{\text{JL}},\phi_{12},\alpha,\delta,D_{\text{L}},t_{\text{c}},\theta_{\text{JN}},\psi,\phi_{\rm
ref}\\}.$ (5)

In the following, we will use a bold-face variable, like $\bm{\vartheta}$, to
describe a set of parameters and regular-face variable, like $\vartheta$, to
describe a particular parameter in the set.

###  II.2 Waveform models

We consider two state-of-the-art, quasi-circular spin-precessing waveform
models incorporating subdominant spherical harmonics — SEOBNRv5PHM and
IMRPhenomXPHM. The GW modes $(l,m)\neq(2,2)$ are important both for detection
[110, 111, 112], where their non-inclusion leads to a loss of signal power for
asymmetric binaries and inclined orbits, and parameter estimation, where these
modes can break degeneracies between various parameters and improve the
measurement accuracy [88, 89, 113, 114].

The SEOBNRv5PHM waveforms contain the spherical harmonics
$(l,|m|)=(2,2),(2,1),(3,3),(3,2)$, and $(4,4),(4,3),(5,5)$ in the coprecessing
frame. However, in this paper we do not include the $(l,m)=(5,5)$ mode. The
IMRPhenomXPHM waveforms include the $(l,|m|)=(2,2),(2,1),(3,3),(3,2),(4,4)$
modes in the coprecessing frame 444We note that, since our work started, there
have been a few important updates on phenomenological models [67, 115], which
included NR calibration to the precessing sector, a more faithful ringdown
model, and improvements to the spin-precessing equations. However, we do not
expect that our results would change substantially, if we used those new
waveform models.. Both waveform models can be used for a wide range of mass
ratios, as well as BH-spin magnitudes up to the maximal values. However, only
the aligned-spin sectors of both waveform models were calibrated to NR
simulations, and their accuracy has been assessed only in regions of parameter
space where NR is available.

In this work, we consider a signal generated using the SEOBNRv5PHM model to be
the true GW signal and analyze it using IMRPhenomXPHM as the template model.
For the Bayesian analyses of this paper, this is done because IMRPhenomXPHM is
quicker to evaluate due to it being a frequency-domain model while SEOBNRv5PHM
is a time-domain model and it is slower. Furthermore, the computational
efficiency of IMRPhenomXPHM can be improved by utilizing the multibanding
approach [116] while no such analogous methods exist for time-domain models.
For the Fisher–information-matrix analysis discussed later, we find
instabilities in the numerical derivatives of the SEOBNRv5PHM waveform with
respect to the GW parameters, for some regions of the parameter space, and
hence restrict ourselves to computing derivatives of the IMRPhenomXPHM model.
We expect to address this issue in the future.

###  II.3 Detector networks

Figure 2: The amplitude–spectral-density curves of the various detectors used
in this paper (see footnote 5). The curves labeled by A+ and V+ denote the
design sensitivity of the LIGO and Virgo detectors, respectively, which form
part of the network of the fifth observing run (O5), while A# refers to the
LIGO detectors at upgraded sensitivity. The next-generation observatories are
the ET and CE with the baseline network for the latter consisting of a 20 km
and a 40 km detector.

The current detectors are expected to achieve design sensitivity in the next
few years, during the fifth observing (O5) run, and continue operating till
the end of the decade [117]. It is anticipated that the detectors would
undergo major upgrades thereafter and operate until next-generation detectors
come online or even in tandem with them. Since plans for future detector
networks have not yet been finalized, a number of studies have explored the
capabilities of different combinations of detector configurations to
understand what the optimal design is for various science goals [27, 29, 30,
118, 28]. In this work, three GW detector networks, consisting of the current
detectors at design and upgraded sensitivity, and proposed future detectors,
are considered to emulate a highly probable observing scenario for the coming
decades 555 The A+, V+, ET, and CE sensitivity curves in this work are those
used in Ref. [27] while the A# sensitivity curve is taken from
https://dcc.ligo.org/LIGO-T2300041/public .. These are enumerated below:

  * •

_O5_ network: This is comprised of the advanced LIGO detectors located at
Hanford and Livingston, and the advanced Virgo detector operating at design
sensitivities, A+ and V+, respectively.

  * •

_A#_ network: In this configuration, the LIGO detectors operate at upgraded A#
sensitivity while the Virgo detector continues to operate at design
sensitivity, V+.

  * •

_XG_ network: This comprises of three proposed XG observatories consisting of
the baseline 40 km and 20 km CE in the United States, and an ET in Europe.

The power spectral density (PSD) of the individual detectors are shown in Fig.
2.

##  III Statistical Methods

This section describes the data-analysis methods used in this paper. In Sec.
III.1 we start by describing the Bayesian framework for analyzing GW signals
and lay out the different choices of priors and frequency bands used for the
different networks. Thereafter, in Sec. III.2, we introduce the LSA for the
likelihood, and recount the Fisher information matrix (FIM) [105] method for
estimating measurement errors. Following that, in Sec. III.3, we elucidate the
computation of biases (or systematic errors) under the LSA [100, 119], and
emphasize the importance of minimizing the mismatch between (i.e., aligning)
the signal and template for reliable estimates of the bias. As an example, we
compare the posterior distributions for a chosen binary system, as obtained
from a full Bayesian analysis with the estimates from the LSA. Specifically,
we point out the differences if the bias is computed without aligning the two
waveforms. Finally, in Sec. III.4 we discuss the hierarchical Bayesian method,
which we employ to understand the impact of biases on the inference of the
properties of the BBH population.

###  III.1 Bayesian analysis

The posterior probability distribution on the parameters of the waveform
model, $\bm{\vartheta}$, given the observational data $d$, and the hypothesis
(model description) $\mathcal{H}$, is obtained using Bayes’ theorem,

$p(\bm{\vartheta}|d,\mathcal{H})=\frac{p(d|\bm{\vartheta},\mathcal{H})\,p(\bm{\vartheta}|\mathcal{H})}{p(d|\mathcal{H})}\,,$
(6)

where $p(\bm{\vartheta}|\mathcal{H})$ is the prior probability distribution,
$p(d|\bm{\vartheta},\mathcal{H})$ is the likelihood function, and
$p(d|\mathcal{H})$ is the evidence of the hypothesis $\mathcal{H}$. If one is
interested solely in parameter estimation, and not in model selection, the
latter serves as a normalization constant, and can be discarded.

For a detector with stationary, Gaussian noise, the likelihood function for
the data given the parameters $\bm{\vartheta}$ is defined as

$\ln
p(d|\bm{\vartheta})=\left<d-h(\bm{\vartheta})|d-h(\bm{\vartheta})\right>\,,$
(7)

where we define the noise-weighted inner product as

$\left<h_{1}|h_{2}\right>=4\Re\left[\int_{f_{\text{low}}}^{f_{\text{high}}}\frac{h_{1}(f)\times
h_{2}^{*}(f)}{S_{h}(f)}df\right]$ (8)

with $S_{h}(f)$ being the noise PSD, and $f_{\text{low}}$ and
$f_{\text{high}}$ are the minimum and maximum frequency in the detectors’
bandwidth. This inner product also defines the optimal, matched-filtering SNR
in a detector, $\rho_{n}$, by

$\rho_{n}^{2}=\left<h(\bm{\vartheta})|h(\bm{\vartheta})\right>.$ (9)

The total SNR is $\rho^{2}=\sum_{n=1}^{N}\rho_{n}^{2}$, where $N$ is the
number of detectors in the network. We note that all our injections are
noiseless which corresponds to averaging over multiple noise realizations.

While the current detectors’ sensitivity is limited to a minimum frequency of
20 Hz, at design sensitivity and with further upgrades, they are expected to
reach a low-frequency sensitivity of 10 Hz. Meanwhile, XG observatories are
aiming to further this improvement to 5 Hz. Therefore, the minimum frequency
for the _O5_ and _A#_ networks are assumed to be $f_{\text{low}}=10\,\rm Hz$,
while for _XG_ detectors, it is taken as $f_{\text{low}}=5\,\rm Hz$. On the
other hand, the maximum frequency is kept the same for all three networks at
$f_{\text{high}}=1024\,\rm Hz$. This does not limit the analysis whatsoever
since all the BBH systems considered in Sec. VI merge at much lower
frequencies.

As we mentioned earlier, the signal is generated using the SEOBNRv5PHM model
with the same starting frequency as the analysis — $f_{\text{low}}=10\,\rm Hz$
for _O5_ and _A#_ ; $f_{\text{low}}=5\,\rm Hz$ for _XG_. Since SEOBNRv5PHM is
a time-domain waveform model, this $f_{\text{low}}$ refers to the starting
frequency of the $(l,m)=(2,2)$ mode. Subdominant harmonics with
$m^{\prime}\neq 2$ start at higher frequencies given by
$f_{lo}^{m^{\prime}}={m^{\prime}}f_{lo}/2$. For instance, in _O5_ and _A#_
networks, the $m^{\prime}=3$ modes start at 15 Hz while the $m^{\prime}=4$
modes start at 20 Hz. In Appendix B, we show that this choice does not affect
our results. This is because of the minimal additional information contained
in the missing frequencies compared to the rest of the signal.

To simulate and analyze the GW signals in Sec. VI, we use the publicly
available Bilby package [120, 121], which incorporates the nested sampler
dynesty [122], interfaced through the Bilby-pipe wrapper. Initially, a
14-dimensional GW parameter space is sampled using the dynesty sampler with a
distance-marginalized likelihood. The full posterior probabilities are then
reconstructed using semi-analytic methods [123, 124].

All the detectors used in this study have an L-shaped interferometer
configuration except the ET, which is proposed to have a triangular
configuration. However, the Bilby-pipe wrapper is limited to L-shaped
interferometer configurations. Consequently, the ET telescope is assumed to be
L-shaped in Sec. VI. Our conclusions remain unaffected as the interferometer’s
shape has no significant impact on the science cases discussed here [30].

We make standard choices for the priors for all the parameters [125]. The
priors for the component masses are taken to be uniform, and the spins are
assumed to be isotropic in direction and uniform in magnitude. For the
distance we choose the prior $\propto d_{L}^{2}$, corresponding to a uniform
in comoving volume distribution at low redshift. We assume that the binary’s
position in the sky and the inclination of its orbit in the coprecessing frame
are random. Therefore, we assign uniform priors on $\alpha$, $\cos\delta$, and
$\cos\theta_{\rm JN}$ across their domains. The other extrinsic parameters,
namely, the polarization angle, coalescence time, and coalescence phase are
also taken to be uniform in their respective ranges.

###  III.2 Linear-signal approximation for measurement errors, systematic
biases and alignment

#### Measurement errors

The evaluation of the posterior probability distribution, as described in the
previous section, is computationally expensive. This makes the estimation of
the measurement accuracies and systematic biases for large number of sources
computationally prohibitive using the Bayesian method. An inexpensive
approximate method is the LSA, which we now briefly introduce.

To estimate the parameter-estimation errors, the waveform model is expanded to
linear order in the parameters around the maximum likelihood (best-fit)
values, $\bm{\vartheta}_{\text{bf}}$. This results in a Gaussian likelihood
distribution whose covariance, $C_{ij}$, is given by the inverse of the FIM,
$C_{ij}=\Gamma_{ij}^{-1}$, which takes the form [105, 119],

$\Gamma_{ij}\equiv\left<\frac{\partial
h}{\partial\bm{\vartheta}^{i}}\middle|\frac{\partial
h}{\partial\bm{\vartheta}^{j}}\right>\bigg{|}_{\bm{\vartheta}=\bm{\vartheta}_{\text{bf}}}.$
(10)

The marginalized one-dimensional errors are then given by the diagonal
elements, $\Delta\bm{\vartheta}^{i}=\sqrt{C_{ii}}$ 666$C_{ii}$ is the i-th
element of $C_{ij}$. The approximation holds for large SNR.

We use the publicly available package GWBENCH [126] to calculate the
measurement errors. GWBENCH is an easy-to-use FIM analysis tool for ground-
based detectors that implements finite difference derivatives to estimate the
approximate measurement errors. Other recent FIM analysis codes for compact-
binary coalescences are GWFAST [32] and GWFish [127]. The waveform models are
internally referenced from the LALSuite [128] libraries. While the
IMRPhenomXPHM model is directly present in LALSuite, the SEOBNRv5PHM model is
interfaced through the pySEOBNR package [129] within LALSuite. For
IMRPhenomXPHM, the default model in LALSuite implements a multibanding
approach [116] for faster waveform computation. However, we turn this off in
our FIM analysis because we found that the output of the last frequency bin
has some randomness associated to it. This is harmless in a Monte Carlo
sampling of the likelihood since the amplitude in that frequency bin is
subdominant and does not contribute to the integral of Eq. (7). However, a FIM
analysis involves taking waveform derivatives with respect to binary
parameters and the randomness manifests as a delta-function which dominates
the integral in Eq. (10).

#### Systematic biases

A further assumption in the FIM formalism is that there are no mismodeling
errors, that is, the signal is accurately represented by the model waveform
and errors are only due to a measurement process using detectors with finite
sensitivity. In reality, we do not know the true GW waveform and use various
approximate models to faithfully represent the true signal. As such, there is
a source of error arising from a difference between the signal and the
waveform model used to represent the signal (template). As a result, the
parameters that maximize the likelihood are biased from the true parameters of
the GW signal by $\delta\bm{\vartheta}^{i}$. This mismodeling error,
henceforth called bias $\delta\bm{\vartheta}^{i}$, is given by [100,
119]777This expression first appeared in Ref. [100], but it is often referred
to as the Cutler-Vallisneri formula after a later paper [119], which was the
first to explore its implications. ,

$\delta\bm{\vartheta}^{i}=C^{ij}\,\left<\partial_{j}h|\delta
h\right>\big{|}_{\bm{\vartheta}=\bm{\vartheta}_{\text{bf}}},$ (11)

at the leading order, where $\delta h=h_{\text{s}}-h$ with $h_{\text{s}}$
being the true signal. In practice, the true signal is not known, so this
formula can only be used if the true signal is replaced by some fiducial
reference model, here taken to be SEOBNRv5PHM.

#### Waveform alignment

We now discuss a few subtleties in the use and applicability of Eq. (11) for
the estimation of biases. Note that the bias is directly proportional to the
waveform difference, $\delta h$. In part due to different conventions for some
extrinsic parameters, $\delta h$ can be artificially large when evaluated at
the same value of all parameters, but can be significantly reduced by changing
the values of certain extrinsic parameters, such as the global phase and time
shift, while keeping the intrinsic parameters fixed. In Appendix A we describe
a toy model that illustrates how a simple time shift can cause biases in
physical parameters to become large. Since Eq. (11) is derived under the LSA,
large waveform differences stretch the formula beyond its domain of validity
resulting in unreliable estimates. However, large uncertainties in the
extrinsic parameters are typically not problematic for scientific applications
of GW observations, so if, by changing only a subset of the extrinsic
parameters, we can bring the waveform difference back into the range of
validity of the LSA, this should be done to improve the accuracy of the
inferred results.

Waveform-accuracy studies in the literature that use $\delta h$ as a metric to
quantify waveform differences, and estimate expected biases, have typically
followed this approach and minimized $\delta h$ over the extrinsic parameters
[130, 82] (alignment). However, to the best of our knowledge, many studies
employing Eq. (11) to estimate the bias either neglect this aspect and naively
use the difference between waveform models to estimate the systematic bias, or
at least do not discuss it. The incorrect use leads to unreasonably large
estimated biases, particularly for the luminosity distance. Therefore, we
describe here how we implement the alignment in the bias formula.

Using Eq. (8), we define the unfaithfulness or mismatch between two waveforms
$h_{1}$ and $h_{2}$ as

$\mathcal{M}=\min_{\bm{\lambda}}\left\\{1-\frac{\left<h_{1}|h_{2}\right>}{\sqrt{\left<h_{1}|h_{1}\right>\left<h_{2}|h_{2}\right>}}\right\\}$
(12)

where the minimization is done on a subset of the binary’s parameters that we
denote $\bm{\lambda}$. For nonprecessing waveform models employing only the
dominant quadrupolar mode, $\bm{\lambda}=\\{\psi,t_{c}\\}$. In this case,
$\psi$ is degenerate with $\phi_{\rm ref}$, so we need to consider only one of
them. On the other hand, since we are considering spin-precessing waveform
models, $\bm{\lambda}=\\{\psi,t_{c},\phi_{\rm ref},\phi_{\rm JL}\\}$ where
$\phi_{\rm JL}$ is a rotation of the in-plane spin angles. For spin-precessing
waveform models, some studies have chosen to minimize the mismatch over the
reference frequency instead of in-plane spin rotations [131, 132]. However, in
this study, we choose to optimize the mismatch by rotating the in-plane spin
components [65, 51], thus keeping the reference frequency fixed at $f_{\rm
ref}=20\rm\,Hz$.

Starting with the set of parameters $\bm{\vartheta}$, we find the parameters
$\overline{\bm{\lambda}}$ that minimizes $\mathcal{M}$ in Eq. (12) for the
detector network being considered. The minimization over the polarization
angle $\psi$ is done analytically, while the coalescence time $t_{c}$ is
optimized by convolving the two waveforms utilizing the convolution theorem
[106, 133, 134]. The reference phase $\phi_{\rm ref}$ and in-plane spin
rotations $\phi_{\rm JL}$ are optimized numerically by using standard
optimization algorithms. Having found the parameters that minimize Eq. (12),
$\overline{\bm{\lambda}}$, we have a new set of parameters
$\overline{\bm{\vartheta}}_{\text{bf}}$, where the parameters
$\bm{\lambda}=\\{\psi,t_{c},\phi_{\rm ref},\phi_{\rm JL}\\}$ have been
replaced by the values obtained through Eq. (12). This is the set of
parameters that we use to compute the FIM, as well as, the $\delta h$ in Eq.
(11). Therefore, the alignment procedure modifies Eq. (11) to

$\delta\bm{\vartheta}^{i}=C^{ij}(\overline{\bm{\vartheta}}_{\text{bf}})\,\left<\partial_{j}h(\overline{\bm{\vartheta}}_{\text{bf}})|h_{\text{s}}(\bm{\vartheta})-h(\overline{\bm{\vartheta}}_{\text{bf}})\right>\,.$
(13)

If the parameters $\bm{\lambda}$ are uncorrelated to the other binary
parameters, the bias formula Eq. (11) should give $\delta\bm{\lambda}=0$.
However, in general, that is not the case and, therefore,
$\delta\bm{\lambda}\neq 0$. Thus the total bias for $\bm{\lambda}$ is
$\Delta\bm{\lambda}=\delta\bm{\lambda}+\left(\bm{\lambda}_{s}-\overline{\bm{\lambda}}\right)$
where $\bm{\lambda}_{s}-\overline{\bm{\lambda}}$ is the difference between the
parameters $\bm{\lambda}$ of the fiducial signal and those obtained after the
optimization procedure.

Note that we chose to modify the template in Eq. (11) following the
optimization procedure Eq. (12). Under the LSA, we are free to modify the
signal evaluating the template at the fiducial parameters. However, we notice
a slightly better agreement of the bias with full Bayesian results when
modifying the template. This is because the Bayesian analyses are performed
using the fiducial parameters as the values of the synthetic-injected signal,
and we find the systematic bias to be more sensitive to small changes in the
injected values compared to the measurement errors.

Lastly, we note that $\mathcal{M}$ could already be close to the minimum for
certain pairs of waveform models at a given set of parameters out of the box.
In such cases, the optimization procedure will have minimal effect on the
total bias and one could simply use the bias formula as it is. However, the
total bias $\Delta\bm{\lambda}$ would be the same regardless of whether one
chooses to do the initial optimization or not even though the output of the
bias formula will not be. For the same reason, the net bias does not depend
sensitively on the precision of the optimization routine as the bias formula
compensates for it. Therefore, it is more prudent to compare the net biases
rather than the optimized values following the initial minimization. We have
verified our minimization procedure using a brute-force 4D minimization
algorithm and find that while $\overline{\bm{\lambda}}$ is slightly different
between the two minimization routines, $\Delta\bm{\lambda}$ remains the same.

In the following, we calculate the systematic biases with and without the
optimization procedure outlined above. We compare these estimates with the
full Bayesian-analysis results for a subset of the parameters, and find
agreement with the Bayesian analysis when using the optimization procedure Eq.
(12).

###  III.3 Comparing Bayesian and linear-signal approximation analyses

Figure 3: Comparison of the posterior distributions for the chirp mass,
symmetric mass ratio, luminosity distance, and primary spin magnitude for
_Binary 1_ with parameters given in Table 1 and in the _O5_ detector network.
The distributions obtained from a Bayesian parameter estimation using Bilby
are shown in green. The estimates from the LSA with and without the
minimization procedure (12) (alignment) are shown in orange and brown,
respectively. The black cross-hairs show the true injected value. The
parameter estimation is performed by injecting a SEOBNRv5PHM signal and
recovering it with the IMRPhenomXPHM waveforms. The Bayesian posteriors are
accurately represented by the LSA when the alignment is enforced.

As a representative case to compare results between the Bayesian method, and
the LSA with and without the optimization procedure of Eq. (12), we consider a
BBH with parameters given in Table 1, and denoted as _Binary 1_ in the _O5_
detector network. In Fig. 3, we show the posterior distributions for selected
parameters, namely, $\mathcal{M}_{c}$, $\nu$, $D_{L}$, and $\chi_{1}$, using
the Bayesian analysis and the LSA estimate for the errors and biases computed
from Eq. (10) and Eq. (11), respectively. The full Bayesian posterior
estimates are shown in green. The LSA posterior distributions are
multidimensional normal distributions centered at the biased value with the
covariance matrix given by Eq. (10). The curves in orange show the
distributions when the biases are estimated by minimizing the mismatch between
the waveforms (see Eq. (13)), while the ones in brown are the estimates when
the minimization is not performed (as it is typically done in the literature;
Eq. (11)). Since the covariance matrix is approximately the same in a
neighborhood, the posterior widths are similar. However, the predicted bias
differs substantially between the two procedures. The effect on the estimation
for the distance bias is especially noticeable, with the traditional method
predicting a $\approx 50\%$ bias when it is unbiased in actuality. This would
be of particular importance for cosmology studies where the traditional
estimates done in the literature would be overly pessimistic. The contours
show the 90% credible intervals of the parameters.

We now briefly discuss the validity of the LSA. Even though we find excellent
agreement between the LSA and the Bayesian analysis for this fiducial case, it
is important to keep in mind that the estimates are approximate. Particularly,
both the FIM and the bias formula (Eq. (13)) are derived under the assumption
that a waveform model can be expanded linearly in its parameters. While the
FIM approximation improves with increasing SNR, with higher order
contributions scaling as $\mathcal{O}(1/\acs{SNR})$, the bias is independent
of the SNR both in the linear approximation and the full likelihood. For the
LSA, this can be easily gauged from the bias equation, Eq. (13), which is
independent of the distance and/or simple scaling of the PSD. For the Bayesian
analysis, one can conclude from Eq. (7) that a simple scaling of the PSD will
not affect the stationary points. Therefore, the point in the parameter space
where the likelihood peaks remains constant. This means that the error in the
bias computation is also constant. In addition, note that we are interested in
the bias in units of the statistical errors. Hence, while the measurement
becomes better with improving sensitivity (or larger SNR), the error in the
systematic bias estimated using LSA becomes more important. A priori it is
difficult to know the range of the sweet spot where both approximations hold.
However, the event shown in Fig. 3 has an $\acs{SNR}\sim 75$ and we also
observe similar agreement in the _A#_ network where the event has an
$\acs{SNR}\sim 220$ (see Table 1) prompting us to make the reasonable
assertion that the LSA is most trustworthy for such ranges of the SNR. We were
not able to directly compare the Bayesian results in the XG network with the
LSA estimates because as we explain in Sec. III.1, the former assumed an
L-shaped interferometer for ET while the latter was done using a triangular ET
configuration. We would also like to stress that one would expect the LSA to
hold when the mismatch between two waveform models is not too large. For the
case illustrated above, we find the mismatch, $\mathcal{M}\sim 3\%$. However,
the binaries that are considered in Sec. IV and Sec. V can have much larger
mismatches and a more detailed analysis is required to quantify the validity
of the LSA as a function of the mismatch which is beyond the scope of this
study.

###  III.4 Hierarchical Bayesian analysis

We now discuss the method, which we employ in Sec. IV to understand the impact
of the biases on the inference of the properties of the BBH population. Given
a set of $N_{\rm obs}$ observed data $\\{d_{i}\\}$, we can estimate the
underlying distribution of parameters that generated it through a hierarchical
Bayesian analysis. We denote by $\vartheta$ ($\subset\bm{\vartheta}$) the set
of parameters, whose distribution we wish to infer. Assuming a form for the
number density of observed events, $\frac{{\rm d}N}{{\rm
d}\vartheta}(\Lambda)$, that depends on _hyperparameters_ $\Lambda$, the
posterior on the latter is given by [135, 136]

$p(\Lambda|\\{d_{i}\\})\propto\pi(\Lambda)e^{-N(\Lambda)}\prod_{i=1}^{N_{\rm
obs}}\int\frac{{\rm d}N}{{\rm
d}\vartheta}(\Lambda)\frac{p(\vartheta|d_{i})}{\pi_{\text{PE}}(\Theta)}{\rm
d}\vartheta,$ (14)

where $p(\vartheta|d_{i})$ is the single-event posterior,
$\pi_{\text{PE}}(\vartheta)$ is the prior used for parameter estimation,
$\pi(\Lambda)$ is the prior on the hyperparameters, and $N(\Lambda)$ is the
total number of events, defined as

$N(\Lambda)=\int\frac{{\rm d}N}{{\rm d}\vartheta}(\Lambda){\rm d}\vartheta.$
(15)

In the analysis of real data, the above equation must be modified to include
selection effects. We interpret ${\rm d}N/{\rm d}\vartheta$ as the rate
density of the full population, and modify the argument of the exponential to
$p_{\rm det}(\Lambda)N(\Lambda)$, where $p_{\rm det}(\Lambda)$ is the
probability of detection of a source, averaged over the population model. In
the analysis performed here, we instead approximate selection as a hard cut on
the intrinsic SNR of the source. In this model, selection is now defined on
the source parameters, not the data, and the above equation can be used
directly, but ${\rm d}N/{\rm d}\vartheta$ must now be interpreted as the rate
density in this observed portion of the population. This approach, which is
common in the literature, ignores the fuzziness at the detection horizon that
arises from instrumental noise, but will give quantitatively reliable and
unbiased results, provided the data is simulated from the same model. In Eq.
(14), we use the proportionality symbol instead of the equality one because we
have omitted numerical factors that depend on the observed data $\\{d_{i}\\}$,
but not on $\Lambda$, i.e., the individual event evidences and the overall
model evidence. These factors are required to perform model selection, but are
unimportant when the goal is to obtain the posterior distribution on
$\Lambda$.

We perform a hierarchical Bayesian analysis for each source parameter
separately, i.e., $\chi_{1}$, $q$, $\cos\theta_{1}$, and $\mathcal{M}_{c}$, so
that ${{\rm d}N}/{{\rm d}\vartheta}$ is a one-dimensional function. This
yields optimistic measurements for the number densities as compared to the
full inference, but allows us to have a quick assessment of the impact of
systematic biases on population inference. Adopting the approach of Toubiana
_et al._ [137], we describe the number density of observed events, ${{\rm
d}N}/{{\rm d}\vartheta}(\Lambda)$, as a piece-wise linear function. The
extremities of the $\vartheta$ range over which we perform the inference are
fixed, and determined by the minimum and maximum samples present in the data.
Thus, our hyperparameters are: the values of the number densities at the
extremities, the number of knots, their positions and the value of the number
density at the knots. The number density at any point is then obtained by
linear interpolation. We stress that the number of knots is a free parameter
of the model, and is inferred by using a reversible-jump Markov chain Monte-
Carlo algorithm [138]. In this way, the complexity of the model is determined
by the data itself.

For a given detector network, we perform population inference on a mock
catalog with systematic biases and on one without, generated as follows.

  1. 1.

We draw the parameters $\vartheta_{0}$ from the population model described in
Sec. IV and select those with SNR above a given threshold.

  2. 2.

We compute the measurement error and the systematic bias for all observable
events using the LSA, as described in Sec. III.2.

  3. 3.

For the catalog with systematic biases, we shift the true parameters by
$\delta\vartheta_{bf}$ to obtain the biased parameters, $\vartheta_{bf}$.

  4. 4.

For each event $\vartheta_{i}$, we attribute a measurement error $\sigma_{i}$
drawn randomly among the set of computed measurement errors, allowing for
replacement.

  5. 5.

We draw a noisy measurement $\vartheta_{n,i}$ of each event from a Gaussian
centered at $\vartheta_{i}$ ($\vartheta_{bf,i}$ for the biased catalog), with
standard deviation given by the error drawn in step 4.

Under the LSA, the posterior distribution on $\vartheta$ is a truncated
Gaussian:

$p(\vartheta|\vartheta_{n,i})=\frac{2\exp\left[-\frac{1}{2}\frac{(\vartheta-\vartheta_{n,i})^{2}}{\sigma_{i}^{2}}\right]}{\sqrt{2\pi}\sigma_{i}\left[\operatorname{erf}\left(\frac{\vartheta_{\rm
max}-\vartheta_{n,i}}{2\sigma_{i}}\right)+\operatorname{erf}\left(\frac{\vartheta_{n,i}-\vartheta_{\rm
min}}{2\sigma_{i}}\right)\right]},$ (16)

where $\vartheta_{\rm min}$ and $\vartheta_{\rm max}$ are the boundaries of
the prior domain on $\vartheta$. The purpose of the randomization of the
errors (step 4) is to remove the dependency on $\vartheta$ from the standard
deviation entering the posterior distribution. If we were to use the
corresponding value predicted by the FIM for each event, we would have to
account for the complicated dependency of $\sigma$ on $\vartheta$, and the
posterior on $\vartheta$ would no longer be a Gaussian, requiring to go to
beyond quadratic order in the LSA. Moreover, $\sigma$ would also depend on the
remaining parameters in $\bm{\vartheta}$, and, by performing the inference on
a single parameter, we would not be accounting for this dependence correctly,
making our analysis not self-consistent. However, we observe that, for
$\vartheta=q$ or $\chi_{1}$ or $\cos\theta_{1}$, the amount by which the
estimated uncertainty in the parameter varies over the range of our priors is
small, so we expect our procedure to yield realistic results for those
parameters. Step 5 is crucial to make sure our mock catalog is self-consistent
from the statistical point of view. Working in the so-called zero-noise
approximation is valid for performing parameter estimation on single mock
events, because it is a fair realization of the noise in the detector. On the
other hand, having zero noise for all events is no longer a fair realization,
and would be valid only if all the events were perfectly measured. Note that,
in steps 3 and 5, we allow the biased parameters and the noisy ones to be
outside of the prior range. The rationale is that those steps are meant to
mimick the behavior of the likelihood function in the presence of systematic
biases and noise, which, as a function, does not contain information on the
physically allowed range of a given parameter. The posterior in turn is
truncated to the prior range, as done explicitly in Eq. (16).

In the hierarchical Bayesian analysis, we take the parameter estimation prior
$\pi_{\text{PE}}(\vartheta)$ to be flat in $\vartheta$. Thus, each of the
integrands in Eq. (14) is the product of a Gaussian with a piece-wise linear
function, and we can perform the integration analytically. This allows us to
evade problems related to having an insufficient number of samples when
performing Monte-Carlo integration [139] and speeds up the analysis.

##  IV Systematic biases in the BBH population

In this section, we study the effect of systematic biases on a BBH population.
We use the GWTC-3 results [5, 14] only for the distribution of masses. We
explore the impact of systematic biases on this LVK-like population,
considering the three detector networks introduced in Sec. II.3. We also do a
hierarchical Bayesian inference of the underlying population where we re-
weight our population distribution to the current LVK distribution of
astrophysical BBHs.

###  IV.1 LVK-like population

We simulate $10^{5}$ binaries in each of the detector networks described in
Sec. II.3 up to a redshift $z=3$ using the SEOBNRv5PHM waveform model.
Following Borhanian and Sathyaprakash [27], this is around the expected number
of BBH mergers per year. We choose a network SNR threshold of 12 for detection
and use this to identify the subset of simulated binaries that are in the
population observed by each network.

The redshift distribution for the population is drawn from a probability
distribution given by

$p(z)\propto\frac{dV_{c}}{dz}\frac{1}{1+z}\psi(z),$ (17)

where $dV_{c}/dz$ is a comoving volume element per unit redshift and $\psi(z)$
is the SFR which is taken to be [103]

$\psi(z)=0.015\frac{(1+z)^{2.7}}{1+[(1+z)/2.9]^{5.6}}\rm
M_{\odot}yr^{-1}Mpc^{-3}.$ (18)

The distribution of the primary source-frame mass is assumed to follow the
Power Law + Peak model of Abbott  _et al._ [14] with the parameters fixed to
their maximum likelihood values. For completeness and ease of reference, we
elucidate the model here. The Power Law + Peak model is given by

$\begin{split}p(&m_{1}|\lambda_{\rm peak},\alpha_{m},m_{\rm min},m_{\rm
max},\mu_{m},\sigma_{m},\delta_{m})\propto\\\ &[(1-\lambda_{\rm peak}){\rm
PL}(m_{1}|\alpha_{m},m_{\rm max})+\lambda_{\rm
peak}\mathcal{N}(m_{1}|\mu_{m},\sigma_{m})]\\\ &\times S(m_{1}|m_{\rm
min},\delta_{m}),\end{split}$ (19)

where $\lambda_{\rm peak}$ gives the weight of the peak component, ${\rm
PL}(m_{1}|\alpha_{m},m_{\rm min},m_{\rm max})$ is a normalised power-law
distribution with spectral index $\alpha_{m}$ and truncated to the range
$[m_{\rm min},m_{\rm max}]$, $\mathcal{N}(m_{1}|\mu_{m},\sigma_{m})$ is a
normalised Gaussian distribution and finally, $S(m_{1}|m_{\rm
min},\delta_{m})$ is a smoothing function defined by

$S(m|m_{\rm min},\delta_{m})=\begin{cases}&0;\ {\rm if}\;m<m_{\rm min},\\\
&[f(m-m_{\rm min},\delta_{m})+1]^{-1};\\\ &\ {\rm if}\;m_{\rm min}\leq m\leq
m_{\rm min}+\delta_{m},\\\ &1;\ {\rm if}\;m>m_{\rm
min}+\delta_{m},\end{cases}$ (20)

with

$f(m,\delta_{m})=\exp\left(\frac{\delta_{m}}{m}+\frac{\delta_{m}}{m-\delta_{m}}\right).$
(21)

The maximum likelihood values for the fit to GWTC-3 [14] were $\lambda_{\rm
peak}=0.02$, $\alpha_{m}=-3.5$, $m_{\rm min}=4.8M_{\odot}$, $m_{\rm
max}=83M_{\odot}$, $\mu_{m}=34M_{\odot}$, $\sigma=1.9M_{\odot}$ and
$\delta_{m}=5.4M_{\odot}$. The mass ratio distribution is modeled using a
power-law with a smoothing function and takes the form

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

where the maximum likelihood value for the spectral index $\beta_{q}=0.76$.

Figure 4: The distribution of component masses for a population of $10^{5}$
BBHs following the astrophysical distribution as determined by the LVK
Collaboration.

The analyses performed on GWTC-3 [14] suggest a broad distribution for the
spin magnitude, peaking around $0.2$ and falling off to 0 for large spins.
However, we are particularly interested in estimating systematic biases for
large-spin systems, so we draw the magnitude uniformly between 0 and 1. The
analyses on systematic effects are performed on this population, in particular
we compute the measurement uncertainties and systematic biases within the LSA
for these parameters, but when performing the hierarchical Bayesian analysis,
we generate the mock catalog by performing importance sampling on this flat
population to obtain a spin distribution in agreement with the results of the
DEFAULT model of Abbott  _et al._ [14]. Finally, we assume the spins’
orientation to be distributed isotropically, which is in qualitative agreement
with the results on GWTC-3.

The location and orientation of the binary in the plane of the sky is assumed
to be randomly distributed. Therefore, the declination angle, $\delta$, right
ascension, $\alpha$, and inclination angle, $\iota$, follow the distributions
$\cos\delta\in\mathcal{U}[-1,1]$, $\alpha\in\mathcal{U}[0,2\pi]$, and
$\cos\theta_{\rm JN}\in\mathcal{U}[-1,1]$, respectively. The polarization
angle, $\psi$, and the coalescence phase, $\phi_{c}$, are also drawn from a
uniform distribution, $\psi,\phi_{c}\in\mathcal{U}[0,2\pi]$.

We report the SNR distribution of the $10^{5}$ binaries simulated in the three
detector networks in Fig. 5. We shade the region below the SNR threshold of 12
in grey. The tail of the three distributions exhibit the $\propto\rho^{-4}$
dependence of the rate of mergers per unit redshift in accordance with the
uniform in comoving-volume distribution of sources in the nearby universe. On
the other hand, the peak of the distribution correlates with the peak of the
SFR, while the initial slope depicts the first generation of stars following
the dark ages. While we are limited by the sensitivity of the current detector
networks and their upgrades in our ability to observe GWs from the mergers of
the first stellar-origin BBHs, the _XG_ network will enable us to study BBH
mergers immediately following reionization.

Figure 5: The distribution of the SNRs of the $10^{5}$ BBHs simulated in the
three detector networks computed using the IMRPhenomXPHM model. The SNR
distribution is similar when using the SEOBNRv5PHM model. We indicate in
shaded gray the region with network-SNR threshold below 12. Figure 6: The
ratio between the systematic bias and statistical errors,
$\delta\vartheta/\Delta\vartheta$, for source-frame chirp mass, symmetric mass
ratio, primary spin magnitude, and luminosity distances as a function of the
mismatch between SEOBNRv5PHM and IMRPhenomXPHM waveform models for a
population of BBH mergers as observed by the LVK. A network SNR threshold of
12 was imposed on the $10^{5}$ binaries in the population resulting in $\sim
1800$, $\sim 8100$, $\sim 99000$ detected events in the _O5_ (top), _A#_
(middle), and _XG_ (bottom) networks, respectively. The SNRs of the binaries
are depicted using the different colorscales. The outlier events in the top
and middle panels having the largest mismatches are the heaviest and most
distant events with redshifted total mass $>400M_{\odot}$ and redshift $>2$.

###  IV.2 Systematic bias

In the following, we discuss, at first, the systematic biases for the
individual events in the LVK-like population. Then, we carry out a
hierarchical Bayesian inference of the population distributions by reweighting
the parameter distributions to the LVK distribution. Finally, we determine the
type of binaries more likely to be biased, which motivates our analysis of
Sec. V.

Figure 7: Inferred distribution of $\chi_{1}$, $\cos\theta_{1}$, and $q$ for
different detector networks, in the case where the parameters are measured
perfectly (black), where there is measurement error but no bias (brown), and
where there is measurement error and bias (green). The solid lines show the
mean value in a given bin and the shaded areas the $90\%$ confidence
intervals. The “no measurement error and no bias” case does not correspond to
any realistic scenario and is there to show the underlying distribution of
observed events, accounting for Poisson errors. For $\chi_{1}$, when including
biases, we observe a shift towards smaller values together with a broadening
of the distribution, particularly for _A#_ and _XG_. For $\cos\theta_{1}$, we
observe an excess of events at both ends and for $q$ at $\sim 1$. The
distributions in the _XG_ case for those two parameters are less biased
because the bias for individual events is smaller in _XG_.

#### Variation with mismatch

Armed with the biases in the GW parameters and their measurement errors, we
report in Fig. 6 their ratio, $|\delta\vartheta/\Delta\vartheta|$, as a
function of the mismatch, $\mathcal{M}$, between SEOBNRv5PHM and IMRPhenomXPHM
for the chirp mass, symmetric mass ratio, primary spin magnitude, and
luminosity distance. We recall that the biases are computed using Eq. (13).
The SNR for the binaries of the population are portrayed using a colorscale
with lighter colors representing smaller SNR and vice-versa. A value of
$|\delta\vartheta/\Delta\vartheta|>1$ indicates that systematic biases are
larger than the typical size of statistical errors. A common feature for all
the parameters is a direct correlation between
$|\delta\vartheta/\Delta\vartheta|$ and the mismatch. This is intuitive
because a larger mismatch implies a greater difference between the two
waveform models and, therefore, larger biases assuming the measurement errors
do not vary significantly with the mismatch which we find to be broadly true
for the population. On the other hand, the colorscale shows that the loudness
of a signal is not a guarantee for a dominant systematic effect with quieter
signals exhibiting significant systematic biases particularly when the
mismatch is greater. This is especially true for the _O5_ and _A $\\#$_
networks. We provide a complementary plot of
$|\delta\vartheta/\Delta\vartheta|$ as a function of the SNR in Fig. 24 in the
appendix for the interested reader.

A few events appear as outliers in the figure with large mismatches but
extremely small $|\delta\vartheta/\Delta\vartheta|$ for the _O5_ and _A $\\#$_
networks. These are the heaviest and most distant events with redshifted total
masses $>400M_{\odot}$ and redshifts $>2$. As such, these signals are
extremely short consisting of only the merger-ringdown. The FIM in these cases
is close to singular because the signal does not contain much information
resulting in extremely large errors. This suppresses the ratio for every
parameter except the luminosity distance in which case the ratio is directly
proportional to the mismatch and hence large. The LSA approximation is not
reliable for such cases and we should resort to full Bayesian analyses.
Nevertheless, we include these binaries in the figure to show their existence
in the population observable by these two networks.

Figure 6 also shows that the number of binaries with biased parameters as a
fraction of the detected population increases with improving detector
sensitivity. This can be simply understood as due to the independence of the
overall scale of the PSD in estimating the parameter bias (see Eq. (11)) while
the covariance is inversely proportional to it. While an improving detector
sensitivity leads to an increase in the total number of binaries that are
significantly biased, the increase in the biased fraction has to do with the
finite number of stellar origin BBH mergers in the universe. Note that because
of different PSD shapes, interferometer designs, and minimum frequencies, the
three rows of Fig. 6 are not simply shifted versions of one another.
Nevertheless, even for the _XG_ network, only a minority of events are biased
with the biased fraction ranging from $10-25\%$ depending on the parameter.
For detectors of the current generation, the biased fraction is even smaller
with only $\sim 2\%$ and $\sim 2.5\%$ of binaries significantly biased for the
_O5_ and _A $\\#$_ networks respectively. This suggests that biases in
parameter estimation will only be of importance for extraordinary individual
events rather than for inferring general characteristics of the population. We
check this more carefully in the next section.

It is also important to realize that a larger value of
$|\delta\vartheta/\Delta\vartheta|$ does not necessarily mean larger value of
the absolute bias. For instance, as we will see in Sec. VI, the _XG_ network
quite often has smaller absolute biases due to improved low frequency
sensitivity where waveform models agree to a greater extent. However, the
improved sensitivity reduces the measurement error more than the decrease in
the systematic bias resulting in a larger value of
$|\delta\vartheta/\Delta\vartheta|$. We will see the effect of this in the
next section.

#### Inferred distributions

Figure 7 shows the inferred number density of events for $\chi_{1}$,
$\cos\theta_{1}$, and $q$ for the three detector networks. Solid lines
indicate the mean of the number density and colored bands the $90\%$
confidence intervals. The results shown in black correspond to an unphysical
scenario where all observed events are perfectly measured. In this case, the
only source of uncertainty is the Poisson error due to the finite number of
events. The black curves serve as guidance to indicate the underlying
distribution on which we perform the inference. We show in brown the case with
measurement error and no bias, following the procedure outlined in Sec. III.4.
The brown and black shaded areas overlap, with the red encompassing the blue
most of the time due to the inclusion of measurement errors, indicating that
our procedure yields unbiased results. Thus, differences between the cases
without (in brown) and with biases (in green) are due to waveform systematics.

We recall that, when performing the hierarchical Bayesian analysis, we
resample the flat $\chi_{1}$ distribution into the distribution inferred by
the LVK Collaboration [14] by means of importance sampling. For O5, the biased
and non-biased distributions are mostly compatible. However, for _A#_ and _XG_
, we observe that, when including bias, the inferred $\chi_{1}$ distribution
is broadened, with the peak being shifted to lower values. The broadening is a
consequence of the occurrence of large biases, while the shift happens because
the systematic bias typically increases with $\chi_{1}$ (see Fig. 8 and the
associated discussion). Events with large spin are more shifted, with a small
preference for shifts towards lower $\chi_{1}$, than events with small spins.
From the astrophysical point of view, the shift of the peak is rather
negligible, but the tail of high spin events would lead to an overestimation
of the number of BHs with high spins, by up to two orders of magnitude,
potentially challenging formation scenarios.

Figure 8: The percentage of binaries in each $\chi_{1}$ bin with
$|\delta\vartheta/\Delta\vartheta|>1$ for different parameters. The top,
middle, bottom panels show the percentages for the _O5_ , _A#_ , and _XG_
networks, respectively.

For the tilt angle, we observe an excess at the ends due to events with more
precession ($\cos\theta_{1}\sim 0$) being more biased than those with aligned
spins. We note that the _XG_ population is less biased than the _A#_ one. This
is a consequence of the improvement at detector at low frequencies, which
increases the proportion of inspiral signal that is observed, where waveform
models agree best, and yielding a less biased estimate of precession effects.
We observe a similar behavior for the $q$ distribution. Asymmetric events are
more biased than nearly-equal-mass ones, shifting the overall distribution to
$q\sim 1$. As for $\cos\theta_{1}$, the _XG_ population is less biased than
the _A#_ one. At first glance, this might seem in contradiction with Fig. 6,
which shows that the ratio between the systematic bias and statistical error
on $\nu$ tends to be larger in the _XG_ case. However, as further discussed in
Sec. VI, when comparing the full posteriors, in many cases the result in the
_XG_ case is closer to the true value than in the _A#_ case. The ratio between
the systematic bias and statistical error is larger for _XG_ because the
measurement error decreases more than the bias (in relative terms), but both
errors decrease and the fact that the absolute bias is smaller ends up
reducing the bias at the level of the population inference. We have also
performed hierarchical Bayesian analysis on $\mathcal{M}_{c}$ and found no
bias at the level of the population, as expected given that this parameter is
typically little biased.

Finally, let us stress again that those results were obtained using the LSA
for the measurement error and the systematic bias. As explained in Sec. III.4,
we do allow the biased estimate of the parameters to be outside of the
physical range, with the idea that this would mimic the likelihood behavior:
it is reasonable that the likelihood of an event $\chi_{1}\sim 0$ seems to
peak at $\chi_{1}=-0.1$, and when performing parameter estimation we would
observe a truncated distribution due to the physical prior. However, in some
cases, our formula predicts biases that are orders of magnitude outside of the
physical range (e.g $\chi_{1}\sim-10$), most likely indicating that the LSA
should not be trusted. Indeed, the LSA relies on the quadratic approximation
to the likelihood, which should hold only in a region around the peak of the
likelihood, with a better agreement at high SNRs. Thus, the reason for our
estimates of the bias with _O5_ (_A#_) being so much larger than with _XG_
might also be due to the invalidity of the LSA in some cases. However,
observing more of the inspiral certainly contributes, as discussed in more
detail in Sec. VI. Overall, we expect the results shown here for _XG_ to be
the most reliable.

#### Which binaries are likely biased?

Having explored the effect of waveform systematics on the full detected LVK-
like population and studied the inferred population properties via
hierarchical Bayesian inference, we turn our attention to the subset of
binaries with significant parameter biases. In the following, we identify the
properties of binaries that have a greater susceptibility to systematic
biases. To accomplish this, we explore the dependence of the systematic bias
as a direct function of the binary parameters. In Fig. 8, we show the
percentage of binaries in each $\chi_{1}$ bin with
$|\delta\vartheta/\Delta\vartheta|>1$ for various parameters such as the chirp
mass, symmetric mass ratio, primary spin magnitude, and luminosity distance.
It is immediately clear that the number of binaries with biased parameters
increases with increasing $\chi_{1}$. Notice that for current detectors and
its upgrades, only a tiny fraction of binaries ($\lesssim 1\%$) have biased
parameters when $\chi_{1}<0.4$. Even for highly spinning binaries, the biased
percent is less than 10%. In contrast, we observe that a relatively large
fraction ($10-25\%$) of the binaries have biased parameters in the _XG_
network even when $\chi<0.4$.

Before ending this section, we remark that the results for the systematic
biases of the LVK-like population have been obtained by comparing two state-
of-the-art quasi-circular, spin-precessing, multipolar waveform models.
However, when assessing the accuracy of the waveform models, more robust and
definitive results can be achieved when comparing models to NR waveforms. We
plan to carry out such a study in the near future, although it will be limited
by the number of NR waveforms and their length.

##  V Systematic biases across binary parameter space

In Sec. IV, the general properties of systematic bias across an LVK-like BBH
population were explored. Here, we consider an agnostic BBH population in
order to explore a wider region of the binary parameter space and identify the
regions with greater susceptibility to systematic biases.

Figure 9: Comparison of the LVK-like population and distribution of
exploratory binaries.

We sample uniformly in the total redshifted mass,
$M^{z}\in\mathcal{U}(10,200)[M_{\odot}]$, and inverse mass ratio,
$1/q\in\mathcal{U}(1,30)$. However, we impose a constraint on the mass of the
lighter object, $m_{2}\geq 5M_{\odot}$, and only select those binaries that
satisfy this constraint. This results in a non-uniform distribution in the two
masses. Regardless, in this section, our interest is not in any particular
distribution of parameters but rather with the coverage of the parameter
space. Exploring the region of large inverse mass-ratio is interesting since
the waveform models considered here are expected to differ more in this part
of the parameter space due to differences in the models’ calibration. For
SEOBNRv5PHM, SXS NR simulations for aligned-spin systems at $1/q\geq 15$ and a
non-spinning simulation at $1/q\geq 30$ were included in the calibration of
the model, and second-order gravitational self-force information was
incorporated, improving the reliability of the model at large $1/q$ [52]. An
example of the different behaviour of the models for large $1/q$ is shown in
Figure 21 of Pompili  _et al._ [50], which illustrates the differences in
parameter recovery for a large inverse mass-ratio NR simulation is shown
between the aligned-spin versions of SEOBNRv5PHM and IMRPhenomXPHM, with
SEOBNRv5PHM being more reliable in the recovery of the parameters.

Due to the lack of calibration to precessing NR waveforms in both waveform
models, we expect waveform models to have greater differences for large
$\chi_{\rm p}$ values as evident from the mismatch plot of Fig. 10 in Ramos-
Buades  _et al._ [51]. Therefore, it is of interest to understand the behavior
of systematic biases in these parts of parameter space. Thus, we create a
sample of binaries having a uniform distribution in $\chi_{\rm p}$ (the LVK-
like population, in turn, disfavors large values of $\chi_{\rm p}$). To do
this, we generate a large set of samples for the spin magnitudes and tilt
angles from the precessing prior and retain a subset of these samples such
that the resulting distribution in $\chi_{\rm p}$ is uniform. This selection
procedure has a negligible effect on the distribution of spin magnitude and
tilt of the secondary companion while giving greater weight to large and in-
plane spin for the primary companion. We draw 50,000 binaries using this
procedure to cover the binary parameter space. A comparison in the
$1/q-\chi_{\rm eff}$ and $1/q-\chi_{\rm p}$ planes between the LVK-like and
the agnostic population is shown in Fig. 9.

The distributions of all other parameters are the same as for the LVK-like
population, except for the distance, which is kept fixed at $D_{L}=235\,\rm
Mpc$ while computing the measurement errors and biases. However, errors have a
simple scaling with the distance for given redshifted masses, while the biases
remain unaffected, and we use it to obtain results at other distances.

#### Bias horizon

Figure 10: The scatter shows the distribution of the 50,000 binaries in the
parameter space represented in Fig. 9, while the color scale shows the
distance to which the $D_{L}$ parameter is biased ($\delta D_{L}/\Delta
D_{L}\geq 1$) for the _O5_ (left) and _A#_ (right) networks. Systematic biases
become less important if a binary is at a larger distance since measurement
precision decreases with distance. Therefore, a large bias horizon signifies
that a given parameter ($D_{L}$ in this case) is measured well enough for
systematic biases to be important even at such large distances. Notice that
the binaries are biased up to a greater distance in the _A#_ network compared
to the _O5_ network due its greater sensitivity, and the resultant improvement
in measurement precision. Figure 11: Same as Fig. 10 for the _XG_ network.
BBHs observed with the _XG_ network are biased up to a greater distance than
those observed with either the _O5_ or _A#_ networks, due its greater
sensitivity, and the resultant improvement in measurement precision, with a
majority of the binaries having a bias horizon $\geq 25\,\rm Gpc$ ($z\approx
3$) beyond which stellar-origin BBHs are not expected to exist.

We compute the biases and the measurement errors on the parameter set
$\bm{\vartheta}$ for the 50,000 binaries considered in this section under the
LSA as was done in Sec. IV. The ratio $\delta\vartheta/\Delta\vartheta$ of
systematic errors to statistical errors is a function of the distance to a
binary through the dependence of the statistical error $\Delta\vartheta\sim
D_{L}$. We exploit this to calculate the distance at which the ratio
$\delta\vartheta/\Delta\vartheta=1$ for any given parameter. Since the
measurement errors for a binary with given redshifted parameters increase with
its distance, systematic biases will become less important the farther the
binary is located. Hence, the distance at which
$\delta\vartheta/\Delta\vartheta=1$ is a measure of the _bias horizon_ , i.e.,
the maximum distance up to which systematic biases dominate statistical
errors. Given that the biases and errors for each binary parameter are
different, the bias horizons for different GW parameters are also different.

We show the bias horizon for the $D_{L}$ parameter in Fig. 10 for the _O5_
(left) and _A#_ (right) networks, while that for the _XG_ network is shown in
Fig. 11. The distribution of the 50,000 binaries in the 4D space given by
$\\{M^{z},1/q,\chi_{\rm eff},\chi_{\rm p}\\}$ is illustrated by projecting
them into 2D subspaces. The bias horizon for each binary is shown by the color
bar. The bias horizon increases with increasing detector sensitivity implying
that systematic biases will be prevalent up to greater distances. In
particular, the majority of the binaries in the _XG_ network have $D_{L}$ bias
horizons exceeding $25\,\rm Gpc$ ($z\approx 3$), the distance around which the
first stars formed. This is qualitatively different from the conclusions of
Sec. IV.2, in particular Fig. 6, where all the binaries have $z\leq 3$ but
$\sim 75\%$ of them are not systematic-error dominated. This is due to the
different distributions of the parameters in this section compared to the LVK-
like population. Particularly, the majority of the binaries considered in this
section are highly asymmetric with large redshifted masses. Moreover, a
uniform distribution in $\chi_{\rm p}$ leads to a large fraction of highly
spinning binaries. We gather from Fig. 8 that systematic biases are more
prevalant for such binaries. Note that the boundary in the $\chi_{\rm
eff}-\chi_{\rm p}$ space is physical and is a result of the maximum value of
the spin magnitude being 1. When comparing the biases for a given binary for
different detector networks, it is important to keep in mind that a larger
value of the ratio $\delta\vartheta/\Delta\vartheta$ for a better detector
network need not be necessarily due to a bigger $\delta\vartheta$, but rather
a much more precise measurement (see Sec. IV.2 for further discussion).

We also observe that the importance of systematic biases depend on the
parameter space inhabited by the binaries. For instance, from the $\chi_{\rm
eff}-\chi_{\rm p}$ space, it is clear that binaries with small spins have
smaller $D_{L}$ bias horizon compared to binaries with large $\chi_{\rm eff}$
and/or $\chi_{\rm p}$, with the binaries lying on the parameter space boundary
having the largest $D_{L}$ bias horizon. Similarly, one can also observe that
for positive $\chi_{\rm eff}$ the bias horizon at large $\chi_{\rm p}$ is
higher than for $\chi_{\rm eff}<0$. This can be intuited from Fig. 13 of
Pompili  _et al._ [50] which shows that the mismatch is larger for positive
$\chi_{\rm eff}$ compared to negative $\chi_{\rm eff}$. Highly precessing
binaries with aligned spins have greater $D_{L}$ bias horizon compared to
similarly highly precessing but with anti-aligned spins. Along the same lines,
we also observe that binaries with large total masses and inverse mass ratios
have larger $D_{L}$ bias horizons. Note that we intentionally do not talk
about properties relating to the distribution of binaries in the parameter
space since they do not follow any physically motivated parameter
distributions. We report the bias horizon for $\chi_{1}$ in Figs. 26 and 27 of
the Appendix. They broadly show the same dependence across the parameter space
although, the quantitative values of the bias horizon are different for each
binary parameter.

##  VI Impact of systematics on the science of individual events

Until now, we have discussed the effect of systematic biases for the LVK-like
population in Sec. IV, and explored systematic biases across parameter space
within the LSA in Sec. V. In recent years a number of studies have emphasized
the science objectives that can be accomplished using GWs in the near future
[140, 30, 29]. In this section, we consider a few of those science
applications. We then handpick three binaries (see Table 1) with very relevant
science potential and discuss the effects of systematic biases on various
science objectives. The majority of the results in this section are obtained
using a full Bayesian analysis except where pointed out.

###  VI.1 Science objectives

In the following, we introduce the science applications that will be
considered in this section.

  1. 1.

_Cosmology_ : There exists a tension, at the level of $4.4\sigma$, between the
value of the Hubble-Lemaître parameter, $H_{0}$, measured at high redshift
from the cosmic microwave background [107], and measured using the local
distance ladder comprising Cepheid variables and type Ia supernova [141]. The
emergence of GW astronomy provides an avenue for an independent measurement of
the Hubble-Lemaître parameter, which could importantly contribute to resolving
this tension. Indeed, GWs have already provided multiple independent
measurements, albeit not yet at an accuracy to resolve the tension [19, 20,
142, 143, 144].

A measurement of $H_{0}$ requires both the luminosity distance and redshift of
a source to be estimated. GW observations provide the former, but additional
data or assumptions are required to provide the latter. For GW observations
accompanied with an EM counterpart, primarily binaries containing NSs, the
redshift information is provided by spectroscopic/photometric observations of
the host galaxy [145, 146]. In binary neutron star (BNS) or neutron-
star–black-hole (NSBH) observations, NS tides can be used to provide an
independent redshift measurement [147, 148, 149] and, thereby, infer $H_{0}$
[150, 151, 152, 153]. Alternatively, features in the mass distribution of
compact binaries can be exploited [154, 155, 156, 157]. Finally, a statistical
measurement of the redshift using galaxy catalogues [158, 159, 36, 160, 161,
162] or galaxy cross-correlation techniques [163, 164, 165] can be employed.

Here, the focus is on the last technique, which can be applied to BBH systems,
but relies on an accurate measurement of the distance and sky position. An
indirect effect on cosmological inference due to inaccurate determination of
the mass distribution will be discussed in the following subsections.

Events that have the smallest volumetric uncertainties are the most
informative systems for statistical measurements using galaxy catalogues [20,
160]. This can be intuitively understood as follows: if there is a single
galaxy in the localization volume of a GW event, and one assumes that the
event originated in a galaxy, then there is unit probability that the
identified galaxy hosted the event and the redshift is known as well as the
galaxy redshift. On the other hand, if the volumetric localization of the GW
event is poor and there are a large number of galaxies that are potential
hosts, the redshift distribution would essentially be uniform, getting
contributions from each possible host. All three events studied here are prime
candidates for such a method due to their asymmetrical component masses that
make distance estimates more precise than analogous comparable mass mergers.
This is because of a greater contribution from subdominant harmonics that
break the distance-inclination degeneracy. Spin precession also helps in
breaking this degeneracy since it mixes different modes in the inertial frame.

  2. 2.

_Lower mass gap_ : The nature of compact objects with masses between
$2\rm\,M_{\odot}-3\rm\,M_{\odot}$ can have wide-ranging consequences in
fundamental physics — from the physics of nuclear matter to primordial BH
formation mechanisms — and astrophysics — from hierarchical formation
probabilities to the proportion of rapidly spinning NSs. The LVK Collaboration
has observed two events where one of the components of the binary
unambiguously lies in this mass range — GW190814 [166] and GW200210_092254
[5]. Tidal effects on GW waveforms in highly asymmetric mergers hosting
candidate NSs, such as the _Binary 1_ system, are minimal. Therefore, indirect
constraints on the nature of the secondary and the binary’s formation history
are derived indirectly from their mass and spin measurements.

The existence of ultra-heavy NSs has consequences for the nuclear equation of
state (EoS) at a few times the nuclear saturation density [167] with the
possibility of nontrivial structures in the speed-of-sound relation and
related phase transition phenomena [168, 169] or rapidly rotating NSs
stabilized against collapse by its rotation [170, 171, 172, 173]. On the other
hand, BHs in this mass range will inform the primordial BH formation scenarios
[174, 175, 176] and hierarchical mergers in dense environments [177, 178, 179,
180, 181]. It is also proposed that the secondary gains mass due to accretion
either prior to a supernova explosion [182] or following it [183]. In all of
these scenarios, accurate measurements of the mass and spin are essential to
further the storyline.

  3. 3.

_PISN mass-gap_ : In the mass range of about $50\mbox{--}120\rm\,M_{\odot}$,
there is expected to be a dearth of stellar origin BHs. This is because main
sequence stars with masses heavier than $\sim 120\rm\,M_{\odot}$ have core
temperatures that facilitate electron-positron pair production which lead to a
decrease in radiation pressure in the core of the star, causing explosive
oxygen burning, and a resultant disruption of the entire star. This process,
known as a PISN process, does not leave behind a remnant, thereby producing a
dearth of BHs above $\sim 60\rm\,M_{\odot}$. However, if the mass of the main
sequence star is greater than $\sim 250\rm\,M_{\odot}$, all the heavy elements
undergo photo-disintegration, first to alpha particles and then further. This
reduces the radiation pressure causing the star to implode forming a BH with
mass greater than $\sim 120\rm\,M_{\odot}$.

The determination of the boundaries of this mass-gap can inform the physics of
PISN, such as the $\prescript{12}{}{C}(\alpha,\gamma)\prescript{16}{}{O}$
reaction rate [35, 184], and the role of stellar rotation [185]. Other
astrophysical processes such as mass reversal, mass growth due to accretion,
and hierarchical mergers can result in the formation of BHs that populate this
mass gap.

  4. 4.

_Spin Morphology_ : The spin distribution of BHs in binaries provides crucial
information on their formation channels. Upcoming observing runs of LVK
detectors and XG observatories will measure the spins of compact binaries to
ever greater precision. This will help to constrain the spin distribution of
astrophysical BBH populations and their formation channels. For instance,
binaries formed via isolated evolution tend to have their spins aligned with
the orbital angular momentum, while those formed dynamically are likely to
have an isotropic spin distribution. Similarly, hierarchical formation is
expected to produce larger spins compared to stellar collapse. However, given
that the spin measurements are expected to be precise, it is crucial for them
to be accurate as well to allow unbiased inference of the source properties of
the underlying population. We have already discussed the indirect effects of
the spin measurement on the inference of the nature of the secondary component
of the _Binary 1_ system in Sec. VI.3.2. In the following, we discuss the
systematic biases on the spin for another binary system.

The origin of massive BBHs, particularly those filling the upper mass-gap, can
be traced using their effective spin, $\chi_{\rm eff}$, and spin-precession,
$\chi_{\rm p}$, parameters. A hierarchical formation mechanism leads to large
component spins since the remnant of the previous merger is expected to be
spinning.

While the $\chi_{\rm eff}$ and $\chi_{\rm p}$ parameters can broadly inform
and differentiate between an isolated and dynamical formation channel, an
accurate measurement of the tilts of the two spin vectors with respect to the
orbital angular momentum, $\theta_{1}$ and $\theta_{2}$, and their relative
orientations in the orbital plane, $\phi_{12}$, provide detailed knowledge of
the formation mechanisms and spin distributions of the BBH merger population.
Precessing binaries exist in different spin morphologies [186, 187, 188, 189]
due to spin-orbit resonances [190]. Therefore, precessing binaries can exist
in subpopulations characterized by their spin morphology depending on their
tilt angles at formation [191, 192, 193, 194, 195, 196]. Recent efforts have
gone in to better understanding the spin-precession dynamics [197] and the
ability of current detectors to measure spin-precession effects on the
waveform [198]. There have also been efforts to probe whether the BBHs
detected by the LVK Collaboration are associated with a particular spin
morphology [199] as well as studies on the capability of current detectors at
improved sensitivities to distinguish different spin morphologies [200, 201].

  5. 5.

_Remnant quantites_ : The properties of the remnant BH following a BBH merger
can be determined from its binary parameters. However, the nonlinear merger
makes a fully analytical calculation intractable. As such, estimates of the
final mass ($M_{f}$) and spin ($\chi_{f}$) of the remnant include information
from NR. Several studies have proposed fits for the remnant quantities for
non-precessing binaries [202, 203]. In the precessing case, the final mass is
found to agree very well using a non-precessing formula but the same does not
hold true for the final spin where the in-plane spin components are important
[204, 205]. Some simple arguments to include in-plane contributions have also
been proposed in the literature [206, 207, 208]. Such arguments have been used
to augment the non-precessing final spin estimates [202, 203]. A surrogate
model for the final spin using NR has also been proposed for moderate mass
ratios [205]. In this paper, the reported final mass estimates are the average
of the non-precessing fits in Refs. [202, 203]. The final spin is computed by
averaging the estimates of Refs.[202, 203, 208] with in-plane spin
augmentations for the non-precessing fits of Refs.[202, 203].

Such estimates of the remnant quantities are important not only for modeling
the ringdown in inspiral-merger-ringdown waveforms, but also for probing the
dynamics of the merger and the nature of the remnant. The ringdown of a Kerr
BH is described by quasi-normal modes whose complex frequencies are determined
from the properties of the final stationary BH. Therefore, the ringdown signal
can be used to estimate the mass and spin of the final remnant independently.
The remnant quantities can be affected by deviations from GR which can modify
the radiated energy and angular momentum as well as an exotic remnant which
will have a modified spectrum, resulting in inconsistent estimates from the
ringdown signal and NR inspired fits. Such consistency tests have been
performed for GW150914 [21, 209] and GW190521 [210, 211]. With improving
detector sensitivities, such tests form an important part of null hypothesis
tests of GR.

  6. 6.

_Maximal BH spins_: Obtaining accurate spin measurements from astrophysical
BHs is a non-trivial endeavor for which different electromagnetic approaches
exist. In the context of stellar-mass BHs, it is possible to apply the
continuum fitting method or reflection spectroscopy to X-ray observations
[212, 213] . However, modeling the astrophysical environment is non-trivial
and can introduce systematic effects. GW measurements provide a novel way to
measure BH spins, and due to their vacuum environment, might be easier and
more robust to model. While the spin distribution of BBH systems can provide
valuable information about their formation channels, an individual event with
very high spins, especially if close to extremal Kerr, would be of great
scientific interest. For instance, due to the cosmic censorship hypothesis, no
naked singularities should exist, which implies that BHs should not spin above
$a>1$. Moreover, quasi-normal modes of a rapidly rotating BH become long-
lived, which could lead to turbulence phenomena [214].

###  VI.2 Handpicked binary black holes

We pick three asymmetric and precessing systems out of which one has low total
mass and two have large total masses. These systems have some precedence in
the current LVK-TCs, though they are not the most common events. We pick them
due to their science potential and to explore the parts of the parameter space
that are expected to be considerably affected by systematic biases. The
parameters of these systems, their SNRs in the different detector networks,
and their science objectives are listed in Table 1. We now discuss the
properties of these systems and their analogs in the LVK-TCs.

  * •

_Binary 1_ : highly asymmetric, spin-precessing, low total-mass binary. This
system is modeled after GW190814 and has masses, distance, and inclination
compatible with the measured values [166]. While GW190814 has no measurable
spin precession, the system we consider is highly precessing. This is a
reasonable choice because a dynamical capture or hierarchical merger are
widely accepted as possible formation channels for GW190814 and these channels
can produce highly precessing binaries. The merger rate of such systems is
estimated to be $7^{+16}_{-6}\,\rm Gpc^{-3}\,yr^{-1}$. Hence, even a
pessimistic merger rate of $1\,\rm Gpc^{-3}\,yr^{-1}$ equates to more than 400
such mergers within a redshift of 3 each year. With planned and future
observatories, a handful of these mergers will be observed with large SNRs.

  * •

_Binary 2_ : rather asymmetric, nonprecessing, high total-mass binary. Among
the LVK observations, the candidate events GW190403_051519 [6] and
GW200208_222617 [5] have properties similar to this system. Both the
candidates are asymmetric binaries with large effective spins indicating that
the spins are aligned to the orbital angular momentum. The median of the
posterior distribution of the primary spin magnitude for GW190403_051519 is
$\chi_{1}=0.89$ while the probability that the primary spin of GW200208_222617
is $\chi_{1}>0.8$ is 51%. The primary components for both candidates also have
a considerable probability of being in the upper mass-gap. While these events
were marginal detections with moderate astrophysical significance, their
science potential is immense and confident detections of similar systems in
the future will be invaluable. From a modeling perspective, it is also an
interesting region of the parameter space to probe, because different state-
of-the-art waveform models have significant mismatches for binaries with large
aligned spins.

  * •

_Binary 3_ : rather asymmetric, spin-precessing, high total-mass binary. This
system has the same parameters as _Binary 2_ except for the spins which are
now misaligned with respect to the orbital angular momentum. The origin and
associated formation channels for such systems through stellar evolution is
highly uncertain [215, 216, 217, 218, 219]. Nonetheless, the current
observations provide sufficient evidence for the existence of such systems
and, therefore, it is essential to understand whether current waveform models
have the necessary accuracy to study such mergers.

|  _Binary 1_ : highly asymmetric, spin-precessing, low total mass |  _Binary 2_ : rather asymmetric, nonprecessing, high total mass |  _Binary 3_ : rather asymmetric, spin-precessing, high total mass
---|---|---|---
Parameters | $m_{1}[M_{\odot}]$ | 23.2 | 61.8 | 61.8
$m_{2}[M_{\odot}]$ | 2.6 | 9.5 | 9.5
$\chi_{1}$ | 0.7 | 0.9 | 0.9
$\chi_{2}$ | 0.4 | 0.8 | 0.3
| $\theta_{1}[^{\circ}]$ | 40 | 0 | 140
| $\theta_{2}[^{\circ}]$ | 40 | 0 | 120
| $\chi_{\rm eff}$ | 0.51 | 0.89 | -0.43
| $\chi_{\rm p}$ | 0.45 | 0 | 0.77
Network SNR | _O5_ | 75.3 | 222 | 119
_A#_ | 137 | 405 | 219
_XG_ | 1040 | 3150 | 2490
Science cases | Cosmology | ✓ | ✓ | ✓
Lower mass gap | ✓ | ✗ | ✗
PISN mass gap | ✗ | ✓ | ✓
Spin morphology | ✓ | ✗ | ✓
Remnant quantities | ✓ | ✓ | ✓
Maximally spinning secondary | ✗ | ✓ | ✓
Table 1: We list the properties of the three systems that we study using
Bayesian analysis in Sec: VI. The top row lists the intrinsic parameters of
the systems, the middle row enumerates the SNR of the systems in the three
detector networks used, and the last row illustrates the science applications
associated with each system.

Figure 12: We show the bias horizon for the most relevant parameters of the
three systems (from top to bottom panels) for the three different networks.
The colored lines indicate up to which distance the systematic error of a
parameter is outside the 90 % credible interval of the posterior and thus
describes biased parameters. The circles at the end of each line indicate the
bias horizon, which is the distance where it is at 90 %. The number at the end
of each line reports the SNR at that distance. Note that there is a cutoff at
25 Gpc (lines without circle), and that we only explicitly show SNRs up to 20.

For all three events we perform a full Bayesian analysis placing them at a
luminosity distance of 235 Mpc. While BBH mergers at this distance are
certainly possible and have been observed, they are few in number and the
majority of events will originate from larger distances. On the other hand,
the SNR for events that are sufficiently far away will be small and the
systematic biases will be inconsequential compared to the measurement errors.
It is therefore desirable to know the maximum distance up to which a given
science objective will be affected due to biases in parameter estimation.
However, it is infeasible to repeat the full parameter estimation calculation
for multiple distances. For that reason, we model the posterior distribution
using its median and covariance. This is equivalent to the LSA approximation
in the large SNR limit. In this limit, the bias is independent of the strength
of the signal while the covariance increases with the square of the distance.
We compute the bias horizon for the three sources and these are shown in Fig.
12. A cursory look at the SNRs required for a parameter to be biased reveals
that biases could be present for much smaller SNRs than what is simulated
here. Below, we will discuss the implications of these results in the context
of various science applications. The posteriors for events with low SNR are
not well approximated by a Normal distribution. The simple scaling argument
employed here fails for such cases. Therefore, for parameters for which the
SNRs the bias horizon is $<20$, we do not quote the exact value but rather
denote it with an inequality sign. Similarly, we do not show the exact
projected distances if they are larger than $25\,\rm Gpc$ ($z\sim 3$). Note
that there is an implicit assumption in this scaling argument that the
detector-frame masses are constant. For large distances, the source-frame
masses would be materially smaller. Therefore, the binaries considered may no
longer be appropriate for science objectives related to the source-frame mass.
If, in turn, these binaries are placed at a larger distance, while keeping the
source-frame quantities constant, they will have a larger detector-frame total
mass and, consequently, fewer GW cycles in the detectable band. Since the
differences between various waveform models are greater closer to the merger,
one expects that the bias will be larger for such systems compared to the
scaled binaries. For instance, it can be seen from Figs. 10 and 11 that the
bias horizon is generally greater for larger total detector-frame masses. In
this sense, the reported scaled numbers can be considered conservative
estimates.

Figure 13: The distance posteriors for the _Binary 1_ , _Binary 2_ , and
_Binary 3_ systems (from top to bottom) in the three detector networks
considered in this study, The shaded region shows the 90% credible interval.
The contours in the inset show the 90% credible region for the corresponding
sky position measurement. The true values are denoted by black lines.

###  VI.3 Impact on science objectives

In the rest of the section, we investigate the effect of systematic bias on
various science objectives. Even though we isolate different science
objectives to discuss them individually, this is an artificial separation
since the science is interconnected and so are the biases. Therefore, some of
the same results may be discussed in different sections in slightly different
ways. For instance, a biased measurement of the maximum NS mass has not only
consequences for nuclear physics, but also for cosmology. Similarly, different
representations of the same parameter space may be helpful in highlighting
different aspects of the science. For instance, we refer to the magnitude and
tilt of the spin vector when discussing lower mass-gap events, effective spin
and spin-precession parameters when discussing astrophysical formation
channels of heavy BBHs, and the relative orientations of the spin vectors when
talking about spin morphologies. They are all different slices of the same
spin space. However, the different parameterizations are useful for telling
and highlighting different aspects of the full story.

####  VI.3.1 Cosmology

The luminosity distance and sky localization posteriors for the three
simulated systems in the different detector networks are shown in Fig. 13. An
illustration of an expected distribution of galaxies in a volume uncertainty
region is depicted as blue ellipses. The volume uncertainty includes galaxies
with redshift between $z_{\rm min}=H_{0,{\rm min}}/D_{L,{\rm max}}$ and
$z_{\rm max}=H_{0,{\rm max}}/D_{L,{\rm min}}$, where $H_{0,{\rm
min}}=35\rm\,km\,s^{-1}Mpc^{-1}$ and $H_{0,{\rm
max}}=140\rm\,km\,s^{-1}Mpc^{-1}$, and $D_{L,{\rm min}}$ and $D_{L,{\rm max}}$
are the edges of the 90% credible interval of the $D_{L}$ posterior. The sky
patch is the size of the inset panel that contains the sky location posterior.
The areal density of luminous L* galaxies in the local Universe is taken to be
$0.07\,\rm deg^{-2}$ at a distance of 100 Mpc [124].

While the localisation volume of the _Binary 1_ system (upper panel) is
accurately measured for current detectors at design sensitivity and their next
upgrade, both the distance and the sky position are biased for the _XG_
network. With a sky position measurement precision of $<0.1\rm\,deg^{2}$ at
90% credibility, atmost a single galaxy is expected to be in a volumetric cone
up to the distance to the event [36]. An inaccurate GW measurement will
completely miss the host galaxy in this specific example or result in the
identification of the wrong host in general. From Fig. 12, it is clear that
such a system will give a biased distance estimate for sources at distances up
to $\sim 35$ times farther, having an SNR $\gtrsim 30$ in the XG network.

For the _Binary 2_ system (middle panel), while _O5_ and _A#_ networks
accurately recover the sky position again, the distance is biased for all
three networks. A bias in the distance measurement towards larger values will
cause the $H_{0}$ measurement to be systematically biased towards smaller
values. As in the case of the previous system, the precision in the sky
position measurement for the _XG_ network implies the existence of a single
galaxy in the volume uncertainty region on average. Figure 12 tells us that
the distance bias is important in _O5_ and _A#_ networks for SNRs $\sim 70$
which implies distances of $\sim 800$ Mpc and $\sim 1300$ Mpc, respectively,
while the distance at which $|\delta\vartheta/\Delta\vartheta|=1$ for the _XG_
network is $\sim 5700$ Mpc. Even at a distance which is $\sim 7$ ($\sim 4$)
times the maximum distances in _O5_ (_A#_), the event in _XG_ will have an SNR
of 130.

The measurement of both the distance and sky position is inaccurate for the
_Binary 3_ system (lower panel) for all three detector networks under
consideration. However, the distance bias, particularly for the _XG_ network,
is smaller than the previous systems. On the other hand, the sky position is
biased even for _O5_ and _A#_ networks. Due to the smaller distance biases, we
see from Fig. 12 that the bias horizon for the _XG_ network is also smaller
resulting in a very loud signal of SNR $\sim 400$. Meanwhile, since the SNRs
in the three detector networks for this system are smaller than the _Binary 2_
system, an SNR of $\sim 35$ in the _O5_ network is sufficient for a material
distance bias.

Given that a two percent measurement of the Hubble constant will resolve the
Hubble tension, all the three events in the _XG_ network can single-handedly
resolve the tension while the _Binary 2_ and the _Binary 3_ event can do so
even in an _A#_ network, at the simulated distance of $D_{L}=235\rm\,Mpc$.

####  VI.3.2 Lower mass gap

Figure 14: Corner plot showing the posterior distributions for the mass and
dimensionless spin magnitude of the secondary companion of the _Binary 1_
system in the three detector networks in which the system is simulated. The
shaded region in the 1D posteriors and the contours in the 2D space denote the
90% credible region. The black lines show the true injected value. A smaller
minimum frequency for the _XG_ network ($f_{\rm low}=5\rm\,Hz$) compared to
the _O5_ and _A#_ networks ($f_{\rm low}=10\rm\,Hz$) results in a smaller bias
for the mass.

In Fig. 14, the recovered distributions for the mass and dimensionless spin
magnitude of the secondary companion of the _Binary 1_ system in the three
detector networks are shown. Since this is a highly asymmetric merger, the
spin of the secondary is poorly measured as is clear from the $\chi_{2}$
posteriors in the _O5_ and _A#_ networks. Nevertheless, the measurement
precision improves significantly in the XG network, both due to its broad
sensitivity improvement and a smaller minimum frequency. However, the measured
value is materially inaccurate predicting a much lower value than the
injection. This would affect inference on hierarchical formation channels and
accretion-induced mass growth of a NS since both predict a large component
spin. Instead, primordial BH formation scenarios that predict small spins
would be favoured [176]. The measurement bias in the component mass would also
have impacts on constraints on the speed-of-sound relation in NSs [168, 169]
and inference on the effect of rotation on NS radii [173].

More directly, the upper edge in the mass distribution of astrophysical NSs
can be used to determine the redshift and, thereby, estimate $H_{0}$ [154,
155, 156, 157]. Hence, an inaccurate determination of the edge of the mass
distribution could also bias cosmological parameter estimation. This effect
can be quantified using some simple calculations. Assume that the NS mass
distribution is uniform and this system lies at the edge of that distribution.
It can be shown that the error in the determination of the upper edge of the
mass distribution is given by $\Delta m_{\rm max}={\rm
max}\left[\sqrt{\sigma_{m}\,R_{o}/(N-1)},R_{o}/(N-1)\right]$, where
$R_{o}=(m_{{\rm max},o}-m_{{\rm min},o})$, $\sigma_{m}$ is the typical mass
measurement uncertainty for systems near the upper edge, $N$ is the number of
observations and $m_{{\rm max},o}$ and $m_{{\rm min},o}$ are the maximum and
minimum observed NS masses, respectively [220, 221]. The posterior peak of the
secondary companion in _A#_ network is $\sim 0.08M_{\odot}$ away from the true
value. It is easy to calculate from the above equation that with $N\sim 20$
observations and $m_{\rm min}=1M_{\odot}$, the measurement error in
determining the upper edge of the mass distribution is $\Delta m_{\rm max}\sim
0.08M_{\odot}$.

The merger rate of a _Binary 1_ binary within a distance of 235 Mpc is 1 every
3 years. Several studies have estimated that upcoming GW observing runs are
expected to detect tens of BNS mergers per year [27, 32]. Therefore, a
systematic bias in the mass and spin measurement can bias the inference of NS
properties and cosmological parameters in the near future.

####  VI.3.3 PISN mass-gap

Figure 15: The posterior distributions of the primary mass for the _Binary 2_
(solid lines) and _Binary 3_ (dashed lines) systems in the three detector
networks. The true injected value is shown by the vertical black line. The
filled regions depict the 90% credible interval.

We report the measurements of the primary mass for the _Binary 2_ and _Binary
3_ systems in Fig. 15 in the three detector networks. While the biases in the
aligned spin binary are less than in the precessing case, the measured values
for both are significantly different from the injected value.

Furthermore, knowledge of the BH mass spectrum, particularly the edge of the
mass distribution, can be used for cosmological inference as with the NS mass
distribution. Following similar calculations and considering the simplified
case where the BH mass distribution follows a power-law, $p(m)\propto
m^{\alpha_{m}}$, with power-law index $\alpha_{m}=-3.5$ and sharp cutoffs at
$m_{\rm min}=10M_{\odot}$ and $m_{\rm max}=65M_{\odot}$, it can be shown that
the error in the determination of the upper cutoff is given by,

$\Delta m_{\rm max}=\frac{(m_{{\rm max},o}^{\alpha_{m}+1}N-m_{{\rm
min},o}^{\alpha_{m}+1})^{\frac{1}{\alpha_{m}+1}}}{(N-1)^{\frac{1}{\alpha_{m}+1}}}-m_{{\rm
max},o},$ (23)

until $\Delta m_{\rm max}\sim\sigma_{m}$, the individual event mass-
measurement uncertainty. This uncertainty for $N=10^{3}$ events is then
$\Delta m_{\rm max}\approx 3\,M_{\odot}$. From Fig. 15, it can be observed
that the systematic error becomes dominant, even for the _O5_ network which is
expected to observe $\mathcal{O}(10^{3})$ events every year.

####  VI.3.4 Spin morphology

Figure 16: The posterior probability distribution of the $\chi_{\rm eff}$
(left) and $\chi_{\rm p}$ (right) parameters for the _Binary 3_ system in the
three detector networks. The true value of the injection is shown by the black
vertical line. The filled regions show the 90% credible interval.

In Fig. 16, the posterior distributions of $\chi_{\rm eff}$ and $\chi_{\rm p}$
are depicted for the _Binary 3_ system in the three detector networks. It is
immediately clear that the posterior distributions for both spin parameters
and in all three detector networks are far from the true value. However, it is
also noticeable that the absolute value of the bias is larger for the current
networks and their upgrades compared to the future _XG_ network. Spin
precession measurements are enabled by low-frequency sensitivity since the
modulations in the GW amplitude due to spin precession occur on these
timescales. Moreover, different waveform models agree to a greater degree at
lower frequencies because this regime is closely informed by PN calculations
in the various models. Therefore, a smaller minimum frequency and a better
low-frequency sensitivity enable a more precise and accurate measurement of
the $\chi_{\rm eff}$ and $\chi_{\rm p}$ parameters in the _XG_ network. Even
so, the parameters are significantly biased.

The posterior distributions of the parameters characterising the spin
morphology for the _Binary 1_ and _Binary 3_ system are reported in Figs. 17
and 18. Let us analyze the _Binary 1_ system first. The tilt angles determine
which morphology the binary falls into. It is observed that both the tilt
angles are recovered inaccurately. While the median bias for $\theta_{2}$ is
greater than that for $\theta_{1}$, with the median value of $\theta_{2}$
being $1.5-2$ times the injected value, the poor measurement accuracy due to
the highly asymmetric and low total mass nature of the binary results in the
median value of $\theta_{1}$ being farther away from the true value, when
expressed as a multiple of _sigma_ (the measurement error). The parameter
$\phi_{12}$, which characterizes the spin morphology, is also heavily biased
with the best-fit median values $2-3$ times the injected value depending on
the detector network.

Now, considering the _Binary 3_ system, it is observed that the measurement
accuracy increases due to the greater total mass and a resultant higher SNR.
Similar to the _Binary 1_ system, the systematic bias in the $\theta_{1}$ is
smaller in absolute terms compared to $\theta_{2}$, whose median inferred
value is $0.4-0.6$ the injected value. Interestingly, $\phi_{12}$ is less
biased for the _O5_ and _A $\\#$_ networks. The estimates in the _XG_ network
are not only more biased but also completely different from the _O5_ and _A
$\\#$_ networks revealing the sensitivity of the measurement to the minimum
frequency. Specifically, the measured value of $\phi_{12}$ in the _XG_ network
is consistent with $0^{\circ}$ while the injected value is around
$200^{\circ}$.

Figure 17: The posterior distributions of the tilt angles, $\theta_{1}$ and
$\theta_{2}$, and the relative in-plane spin angle, $\phi_{12}$, for the
_Binary 1_ (left) and _Binary 3_ (right) systems. The black cross-hairs show
the true injected value. For the _Binary 1_ system, on the left panel, the
recovered distributions in all three networks are shown while for the _Binary
3_ system, on the right panel, those for the _O5_ and _A $\\#$_ networks are
reported. The probability distributions in the _XG_ network are very different
and are, therefore, presented in Fig. 18.

Figure 18: The posterior distributions of the tilt angles, $\theta_{1}$ and
$\theta_{2}$, and the relative in-plane spin angle, $\phi_{12}$, for the
_Binary 3_ system in the _XG_ network. The black cross-hairs show the true
injected value. The result for this case is shown separately from Fig. 17
because it is significantly different from the other networks and
incorporating them in the same figure distorts its appearance.

####  VI.3.5 Remnant quantities

Figure 19: The final mass (left) and final spin (right) posteriors for the
_Binary 1_ , _Binary 2_ , and _Binary 3_ systems (from top to bottom) in the
three detector networks considered in this study. The true values are denoted
by black lines.

We report the inferred remnant quantities for the _Binary 1_ , _Binary 2_ ,
and _Binary 3_ systems (top to bottom) in the three detector networks in Fig.
19. Unsurprisingly, there is significant bias in most of the cases. Since we
reported the biases in various mass and spin quantities in the preceding
sections, the biases in the remnant quantities are expected. It can also be
observed from the figure that the remnant properties of high-mass systems are
better measured compared to the low-mass system. Consequently, the biases for
_Binary 2_ and _Binary 3_ as a multiple of the measurement error is larger.

####  VI.3.6 Maximal BH spin

Figure 20: The probability distributions of the secondary dimensionless spin
magnitude $\chi_{2}$ in the _Binary 2_ (left) and _Binary 3_ (right) systems
in the _O5_ and _A $\\#$_ networks. The injected signal is generated using
SEOBNRv5PHM while the template model is IMRPhenomXPHM. While the estimates are
biased for both the systems and in both the detector networks, the more
interesting feature is that the recovered $\chi_{2}$ distributions imply a
maximally spinning secondary. We do not show the posterior in the _XG_ network
because they do not rail against the boundary, albeit still being biased.

In Fig. 20, the posterior distributions on the secondary dimensionless spin
magnitude $\chi_{2}$ are reported for the _Binary 2_ (left panel) and _Binary
3_ (right panel) systems. It is observed that $\chi_{2}$ is biased for both
systems. While the effect of $\chi_{2}$ on the GW waveform is subdominant due
to the high mass asymmetry and smaller magnitude compared to $\chi_{1}$, which
results in a poorer measurement of this parameter compared to $\chi_{1}$, the
significantly biased recovery hints at systematic differences in the modeling
of the effects of the parameter on the GW waveform. Perhaps, more
interestingly, the secondary is measured to be maximally spinning. This is of
huge importance because such a measurement in a real event would be
revolutionary.

As shown in Fig. 10, the measurement of $\chi_{2}$ will be biased even for low
SNR, albeit the prior will start becoming important at such SNRs.

##  VII Discussion and Conclusion

In this work we have studied systematic biases arising in state-of-the-art BBH
waveform models used by the LVK Collaboration, in particular the quasi-
circular, spin-precessing, multipolar IMRPhenomXPHM and SEOBNRv5PHM models.
With increasing detector sensitivities of the current LVK network, its future
upgrades, and XG detectors ahead, quantifying waveform systematics becomes
central in achieving the promising science goals of GW astronomy. Unbiased
parameter estimation is crucial for a range of applications, from individual
events to the entire BBH population. Moreover, it is also vital for precision
tests of GR as the prevailing theory of gravity.

Throughout this work, we have assumed that the true signal can be represented
by SEOBNRv5PHM, and modeled it with IMRPhenomXPHM. Although both models do not
represent exact solutions, using them in this way for injection-recovery
studies allows one to explore a wide range of the BBH parameter space, for
which accurate NR simulations are not yet available.

To quantify the bias in parameter estimation, we utilized statistical tools
from LSA and full Bayesian analysis (see Sec. III). While LSA is approximate
and relies on large-SNR events, it can forecast results for a large number of
events, which is what we are interested in. Here, the two main methods are the
FIM to approximate measurement errors and the bias formula to predict the
systematic errors. The full Bayesian analysis is computationally expensive but
more reliable, which makes it suitable for understanding selected events in
detail. With these methods, we studied biases and performed hierarchical
inference on the BBH population in Sec. IV. We explored vast parts of the BBH
parameter space in Sec. V and investigated in detail the impact of systematics
on the science cases of individual events in Sec. VI. In the following, we
summarize our main results.

Our first result is mainly on the “use and abuse” of the widely used bias
formula based on the LSA. By comparing with full Bayesian results, we
explicitly demonstrated that the direct application of the bias formula
without a “waveform alignment” procedure, as described in Eq. (12), can yield
unreliable results. To our knowledge, this has not been discussed in the
literature yet, and it is not explained in the main references for the bias
formula [100, 119]. We find good agreement with full Bayesian results only
after performing alignment and proper bookkeeping of the adjusted parameters.
When not correcting for it, one would predict biases in the BBH population
that could overestimate the actual bias by two orders of magnitude. This is
particularly important when comparing waveforms of different families. Even
after such caretaking, it is important to keep in mind that the LSA is an
approximate estimate of the biases and statistical errors. We include
discussions on the validity of the estimates and possible shortcomings of the
method.

Regarding the LVK-like BBH population studied in Sec. IV, we find that the
vast majority of events measured with the _O5_ and _A $\\#$_ networks are not
biased, around $2.5\%$ of the events, whereas, for _XG_ detectors, up to
$25\%$ of the events are expected to show biases. This is most clearly
reported in Fig. 6, which relates the mismatch of each event with the ratio of
systematic to statistical error. Although the bias relative to the statistical
error is larger for _XG_ detectors, the absolute bias per-se is smaller thanks
to the detectors’ improvement at low frequencies, which allows to better
observe the signal in the inspiral regime, where waveform models agree best.
Thus, when combining all events through a hierarchical Bayesian analysis on
the observed population, we find that the impact of systematic biases is more
pronounced for _A $\\#$_ than for _XG_ detectors. As illustrated in Fig. 7,
the most biased parameter is the magnitude of the spin of the primary, with
the recovered distribution peaking at smaller values, while exhibiting a tail
at large spin magnitudes that is absent in the injected population. The latter
feature, in particular, could erroneously guide astrophysical formation
scenarios into explaining the existence of a sizeable number of BHs with large
spins. We stress that the bias estimates obtained with the LSA might be overly
pessimistic, since it relies on the quadratic approximation to the likelihood,
which holds only in the high SNR regime. Thus, the impact on the population