# Anatomy of parameter-estimation biases in overlapping gravitational-wave
signals

Ziming Wang1,2, Dicong Liang2, Junjie Zhao3,4, Chang Liu1,2,5, Lijing Shao2,6
1 Department of Astronomy, School of Physics, Peking University, Beijing
100871, China  2 Kavli Institute for Astronomy and Astrophysics, Peking
University, Beijing 100871, China  3 Institute for Frontiers in Astronomy and
Astrophysics, Beijing Normal University, Beijing 102206, China  4 Department
of Astronomy, Beijing Normal University, Beijing 100875, China  5 Laboratoire
des 2 Infinis - Toulouse (L2IT-IN2P3), Université de Toulouse, CNRS, UPS,
F-31062 Toulouse Cedex 9, France  6 National Astronomical Observatories,
Chinese Academy of Sciences, Beijing 100012, China<EMAIL_ADDRESS>(LS)

###### Abstract

In future gravitational-wave (GW) detections, a large number of overlapping GW
signals will appear in the data stream of detectors. When extracting
information from one signal, the presence of other signals can cause large
parameter estimation biases. Using the Fisher matrix (FM), we develop a bias
analysis procedure to investigate how each parameter of other signals affects
the inference biases. Taking two-signal overlapping as an example, we show
detailedly and quantitatively that the biases essentially originate from the
overlapping of the frequency evolution. Furthermore, we find that the
behaviors of the correlation coefficients between the parameters of the two
signals are similar to the biases. Both of them can be used as
characterization of the influence between signals. We also corroborate the
bias results of the FM method with full Bayesian analysis. Our results can
provide guidance for the development of new PE algorithms on overlapping
signals, and the analysis methodology has the potential to generalize.

##  1 Introduction

Since the first gravitational-wave (GW) event was detected in 2015 [2],
nowadays GW observatories LIGO, Virgo, and KAGRA have detected nearly 90 GW
events, which greatly enhanced our understanding of fundamental physics,
astrophysics, and cosmology [1, 3, 4]. In the future, the next-generation (XG)
ground-based GW detectors, such as Cosmic Explorer (CE) [36, 37] and Einstein
Telescope (ET) [21, 41] will be constructed. Compared with the present GW
detectors, XG detectors improve the sensitivity by an order of magnitude, and
have a lower cut-off frequency [28], which leads to more detectable GW events
with longer signal duration. Therefore, in the data stream of XG detectors,
there will be a large number of overlapping signals (OSs) [35, 40, 34, 38,
22]. For space-based GW detectors, such as Laser Interferometer Space Antenna
(LISA) [5], Taiji [24] and TianQin [32, 20], a large number of OSs will also
be present [6].

Searching for GW signals from detector data and extracting information of GW
sources are the basis for GW science. The impact of OSs is mainly reflected in
two aspects. On the one hand, every OS comes from the superposition of several
independent signals, so the dimension of the parameter space is high.
Conducting parameter estimation (PE) on the whole OS requires considerable
computational time and resources. On the other hand, if one deals with every
single signal in the OS one by one, the existence of other signals will have
an impact on the PE of the target, usually reflecting as a systematic error,
or the PE bias more quantitatively. Janquart et al. [27] referred to these two
PE cases as the joint PE (JPE) and the single PE (SPE) respectively. If the
SPE bias of each signal is small, then we do not need to conduct complex JPE.
On the contrary, when the number of GW events is known, JPE is a better choice
when SPE biases are large. Compared to JPE, computing SPE is faster, and the
results of SPE can tell us whether we need JPE. Therefore, as the basis of
studying the whole PE technology of OS, it is necessary and useful to
investigate the behaviors of SPE first.

A number of publications [35, 40, 34, 38, 22, 27, 46] have studied how
overlapping will affect SPE. They consider different detectors, such as
Advanced LIGO (AdvLIGO), LIGO Voyager, ground-based XG detectors (CE and ET),
and space-based detector LISA [6] with different compact binary coalescences,
including binary black holes (BBHs) [38, 27], binary neutron stars (BNSs) [34,
46], and neutron star–black hole (NSBH) systems [40]. The GW parameters used
in these works are also different. In addition to the necessary
parameters—binary masses, luminosity distance, and merger time—Himemoto et al.
[22] considered the spin-orbit coupling and spin-spin coupling effects, and Hu
et al. [23] considered the parameterized post-Newtonian coefficients to test
general relativity and calculated the accumulating errors in the PE of OSs.

Qualitatively, when the merging times of different signals are close to each
other, the mutual influence between the signals is expected to be more
significant. In this case, the SPE results will have obvious biases, and we
call the corresponding OS the “biased OS”. In other words, the biased OSs are
real OSs in the sense of PE. On the contrary, if the OS can be decomposed into
individual signals and analyzed separately, we call it “unbiased OS”.
Quantitatively, Relton et al. [38] and Himemoto et al. [22] stated that for an
OS generated by two BBHs, it can be regarded as unbiased as long as the merger
time difference $\Delta t_{c}$ is larger than $0.1\,\rm{s}$. The critical
$\Delta t_{c}$ given by Pizzati et al. [34] is more conservative, which is
$\Delta t_{c}>0.5\,\rm{s}$. Samajdar et al. [40] argued that in most cases,
for the OS of two compact binaries the intrinsic parameters of each signal can
be extracted with negligible biases. These inconsistencies are mainly caused
by the different parameter choices. In addition, these works only pointed out
whether the simulated OSs are biased OSs or unbiased OSs. They summarized the
simulation results but did not systematically investigate the dependence of
bias on the parameters.

In this work, we analyze the SPE biases of OSs generated by two BBHs, which is
the basis to conduct PE for more complex OSs. Specifically, our work answers
the following questions:

  1. (i)

What is the origin of biases? In other words, why could there be biases when
conducting SPE for OSs?

  2. (ii)

How do the biases depend on the merger time of each component in the OS? Will
large biases arise when the merger times are close?

  3. (iii)

Given a target signal, what kind of signals are more likely to generate a
biased OS when overlapping with it?

As mentioned above, it is difficult to directly describe the behaviors of bias
in the whole parameter space. Therefore, we use the Fisher matrix (FM) method
to quickly forecast the SPE biases for numerous OSs. The FM method can be
regarded as an approximation of a full Bayesian analysis, and has been widely
used to forecast and estimate parameters in GW researches [29, 16, 8, 26, 51,
43, 31, 50]. In particular, the FM method has also been applied in OS analysis
[34, 22, 6].

As answers to the above questions, here we summarize our main conclusions:

  1. (i)

Biases originate from the overlapping of the frequency evolution of signals.
In the frequency domain, this is equivalent to a nearly constant phase
difference between signals in a certain frequency band. For unbiased OSs,
their component signals only overlap in the data stream in the time domain,
but can be decomposed into independent signals and analyzed separately in the
frequency domain.

  2. (ii)

A small merger time difference is only a necessary condition for large biases.
Although the bias is indeed larger when $\Delta t_{c}\sim 0$, there is a
significant asymmetry in the order of the mergers. When the binary with
lighter masses merges after the heavier one, the frequency evolution is easier
to overlap, leading to a large bias. In addition, when $|\Delta t_{c}|$
increases, biases decay to zero in an oscillating way rather than
monotonically. Moreover, the bias dependence on $\Delta t_{c}$ is different
for different parameters. The parameters can be divided into two categories
according to whether they enter in amplitude or phase in the waveform. Varying
$\Delta t_{c}$, when the bias for one category is close to extreme values, the
bias for the other category is almost zero.

  3. (iii)

Biased OSs are more likely to arise when the component masses, especially the
chirp masses, of individual signals are close. In this case, the frequency
evolution of the two signals is so similar that large biases will occur in a
wide range of $\Delta t_{c}$. In addition, the louder the other signal is, the
larger the SPE bias of the target signal is.

It is worth noting that all the masses in the above conclusions are defined in
the detector frame, since signals are overlapping in the detector data stream.
When considering a detector network, the biases will depend on the relative
locations between GW sources and detectors [38]. However, in each detector,
the biases’ behaviors are the same as we summarized above.

These conclusions are consistent with the intuitive understanding of the GW PE
process, but they are analyzed and justified in detail for the first time in
this work. In the attempt to obtain these conclusions, we establish an
explicit expression [see Eq. (3.2)] and its corresponding analyzing pipeline
for SPE bias. Benefiting from the rapidity of the FM method, this analyzing
pipeline allows one to discuss and analyze the dependence of bias on various
parameters easily.

In this paper, we also calculated the correlation coefficients between the
corresponding parameters of two signals using the FM method. In this work, we
detailedly calculate and compare the behaviors of correlation coefficients and
biases. Their similar explicit expressions [see Eq. (3.2) and Eq. (22)] lead
to their similar dependence on parameters. This is consistent with the
qualitative results by Himemoto et al. [22]. Large correlation coefficients
mean a strong indistinguishability and lead to a large bias in the SPE.
Finally, to verify the validity of FM, we conduct a full Bayesian analysis for
some representative OSs. For most OSs, the FM results are consistent with the
full Bayesian SPE results. In the large-bias cases, there are slight
differences between the FM method and the full Bayesian analysis. We justify
the use of the FM method to forecast biases, which will provide a reference
and guidance in the simulation of SPE and JPE for OSs.

The paper is organized as follows. Section 2 introduces the basic PE methods
for GWs, including the full Bayesian inference and the FM approximation. In
Section 3, we investigate the bias behaviors and analyze the origin of biases
in depth using the FM method. Section 4 calculates the correlation
coefficients of two signals and then compares them with the FM results. The
applicability of the FM method in analyzing biases is verified using full
Bayesian inference. Section 5 gives the summary and conclusion.

##  2 Parameter Estimation Methods

###  2.1 Basic concepts of Bayesian inference and FM approximation

The data $g(t)$ in a GW detector are the superposition of the GW strain
$h(t;\bm{\theta})$ and the detector noise $n(t)$. In PE, we want to extract
the GW parameters $\bm{\theta}$ from data $g(t)$. Using Bayes’ theorem, one
has the posterior distribution of $\bm{\theta}$,

$P(\bm{\theta}|g)\propto p(\bm{\theta})P(g|\bm{\theta})\,,$ (1)

where $p(\bm{\theta})$ is the prior, and $P(g|\bm{\theta})$ is called
likelihood, which describes the conditional probability of collecting the data
$g(t)$ given the parameters $\bm{\theta}$. Suppose that $n(t)$ is stationary
and follows a Gaussian distribution with a mean of zero, the likelihood can be
further expressed as [18]

$P(g|\bm{\theta})\propto e^{-\frac{1}{2}(g-h,\,g-h)}\,,$ (2)

where the inner product $(u\,,v)$ is defined for any two data streams, $u(t)$
and $v(t)$, as

$\displaystyle(u\,,v):=2\Re\int_{-\infty}^{\infty}\frac{u^{*}(f)v(f)}{S_{n}(|f|)}{\rm
d}f=4\Re\int_{0}^{\infty}\frac{u^{*}(f)v(f)}{S_{n}(f)}{\rm d}f\,,$ (3)

where $S_{n}(f)$ is the one-sided power spectral density of the noise; $u(f)$
and $v(f)$ are the Fourier transforms of $u(t)$ and $v(t)$, respectively.
Using $u(f)$ as example, the Fourier transform is defined by,

$u(f):=\int_{-\infty}^{\infty}u(t)e^{-2\pi ift}{\rm d}t\,.$ (4)

In addition, we define the whitened data stream $u_{\rm
w}(f):=2u(f)/\sqrt{S_{n}}$, which can be interpreted as the signal strength
with respect to the noise in the frequency domain. Once the GW data $g(t)$ is
obtained, one can construct the posterior distribution of the parameters
$\bm{\theta}$ by the above formulas. After that, we can use some sampling
techniques, such as the Markov-Chain Monte Carlo (MCMC) method [14, 15, 42]
and Nested Sampling [44, 45], to obtain the posterior distributions of
$\bm{\theta}$.

In principle, Bayesian inference gives the full posterior at the desired
precision, but it can also take considerable computational time and resources,
especially when the parameter space dimension is high. In contrast, the FM
method is one of the fastest methods to give the approximate posterior, and it
is widely used in the GW data analysis [48], where it reads

$F_{\alpha\beta}=\left(h_{,\alpha}\,,h_{,\beta}\right)\,,$ (5)

where $h_{,\alpha}\equiv\partial h/\partial\theta^{\alpha}$ is the derivative
of the waveform with respect to the $\alpha$-th parameter, computed at the
truth point $\tilde{\bm{\theta}}$ of the parameters. FM is a symmetric,
positive semi-definite matrix, whose inverse exists in most cases in GW
analysis [48, 49]. Denoting the inverse of FM as $\Sigma\equiv F^{-1}$, one
obtain

$\sigma_{\alpha}=\sqrt{\Sigma^{\alpha\alpha}}\,,\quad
c_{\alpha\beta}=\frac{\Sigma^{\alpha\beta}}{\sigma_{\alpha}\sigma_{\beta}}=\frac{\Sigma^{\alpha\beta}}{\sqrt{\Sigma^{\alpha\alpha}\Sigma^{\beta\beta}}}\,.$
(6)

Here, $\sigma_{\alpha}$ is the standard deviation of $\theta^{\alpha}$,
$c_{\alpha\beta}$ is the correlation coefficient between $\theta^{\alpha}$ and
$\theta^{\beta}$ [29]. The FM approximation, which uses Eqs. (5) and (6) to
obtain the errors and correlations of the parameters, is very fast to
calculate, as it requires the derivative of the waveform only at one point. FM
approximation is equivalent to the linearized-signal approximation [18, 48],
that is, the dependence of the waveform template on the parameters is linear,
which holds if and only if the GW signal is loud enough. In addition, when
using the FM approximation result (6), we have defaulted that the priors of
all parameters are uniform. In the frequentist view, the inverse matrix
$\Sigma$ is the covariance matrix of the maximum likelihood estimate (MLE)
values $\hat{\bm{\theta}}$ under linearized-signal approximation, and
$\sigma_{\alpha}$ can be interpreted as (the lower limit of) measurement
uncertainty due to the presence of noise [48].

###  2.2 PE methods of OSs

In the case of OSs, $h(t)$ should be the superposition of individual
independent GW signals, then

$g(t)=h(t;\tilde{\bm{\theta}})+n(t)=\sum_{j=1}^{m}h^{(j)}\big{(}t;\tilde{\bm{\theta}}^{(j)}\big{)}+n(t)\,,$
(7)

where $m$ is the number of GWs overlapping together,
$h^{(j)}\big{(}t;\tilde{\bm{\theta}}^{(j)}\big{)}$ is the strain of the $j$-th
signal, and $\tilde{\bm{\theta}}^{(j)}$ is the corresponding true parameters.
It should be noted that in practice the components of $g(t)$ are unknown; an
OS can contain any number of GWs and the template corresponding to each signal
can also be different [6] .

When conducting PE, one can use different models to fit the GW strain. If the
number of individual signals is known, and the parameters of all GWs are
investigated simultaneously, one should take a model
$h=h(t;\bm{\theta})=\sum_{j=1}^{m}h^{(j)}\big{(}t;\bm{\theta}^{(j)}\big{)}$,
with $\bm{\theta}=\\{\bm{\theta}^{(1)},\cdots,\bm{\theta}^{(m)}\\}$. This is
called JPE [27]. In the SPE case, only one signal, $h^{(1)}$, and its
corresponding parameters are considered. However, when conducting SPE, the
strain data are still determined by Eq. (7). The presence of
$\tilde{h}^{(2)},\cdots,\tilde{h}^{(m)}$ inevitably has an impact on the
estimation of $\bm{\theta}^{(1)}$, which leads to a bias in the SPE results.

For simplicity, we consider the likelihood function of SPE with the inclusion
of two GW signals overlapping,

$\displaystyle P\big{(}g(t)|\bm{\theta}^{(1)}\big{)}\propto\exp$
$\displaystyle\left\\{-\frac{1}{2}\Big{(}\tilde{h}^{(1)}+\tilde{h}^{(2)}+n-h^{(1)}\big{(}\bm{\theta}^{(1)}\big{)}\,,\right.$
$\displaystyle\left.\tilde{h}^{(1)}+\tilde{h}^{(2)}+n-h^{(1)}\big{(}\bm{\theta}^{(1)}\big{)}\Big{)}\right\\}\,.$
(8)

With Eq. (2.2), one can exhaustively discuss how one signal in the OS affects
the SPE of another signal. In principle, this requires exploring throughout
the parameter space and comparing the SPE result with the truth
$\tilde{\bm{\theta}}^{(1)}$ and $\tilde{\bm{\theta}}^{(2)}$. However, due to
the limited computational resources, one in practice can only select some
representative points for full Bayesian analysis, as was done in Refs. [40,
34, 38, 22]. It will be further discussed in Sec. 4.

Now, we use the FM method to derive the explicit expression of biases in SPE.
Using linearized-signal approximation and expanding the waveform at the truth
point, we have
$h^{(1)}\big{(}\bm{\theta}^{(1)}\big{)}\approx\tilde{h}^{(1)}+\tilde{h}^{(1)}_{,\alpha}\big{(}\theta^{(1)\alpha}-\tilde{\theta}^{(1)\alpha}\big{)}$,
where $\tilde{h}^{(1)}_{,\alpha}\equiv\partial
h^{(1)}\big{(}\bm{\theta}^{(1)}\big{)}/\partial\theta^{(1)\alpha}$. Note that
all quantities are computed at $\tilde{\bm{\theta}}^{(1)}$, and we have used
the Einstein summation convention. Substituting it into Eq. (2.2) and using
MLE, we obtain

$\hat{\theta}^{\alpha}-\tilde{\theta}^{\alpha}=\Sigma^{\alpha\beta}\big{(}\tilde{h}^{(1)}_{,\beta}\,,n\big{)}+\Sigma^{\alpha\beta}\big{(}\tilde{h}^{(1)}_{,\beta}\,,\tilde{h}^{(2)}\big{)}\,.$
(9)

Since the SPE is only for $\bm{\theta}^{(1)}$, the superscript “(1)” of
$\bm{\theta}$ is omitted for conciseness. To measure the deviation of the MLE
results from the true value in the statistical sense, we use the most common
approach, namely, the mean square error (MSE) [13], which is defined as the
expected value of the square of the error
$\hat{\theta}^{\alpha}-\tilde{\theta}^{\alpha}$

$\displaystyle\rm{MSE}_{\alpha}$
$\displaystyle:={\rm{E}}\left[\big{(}\hat{\theta}^{\alpha}-\tilde{\theta}^{\alpha}\big{)}^{2}\right]=\big{(}\Delta\theta^{\alpha}_{\rm
bias}\big{)}^{2}+\sigma_{\alpha}^{2}\,,$ (10)

where ${\rm{E}}\left[\cdot\right]$ means the ensemble average over the noise,
$\Delta\theta^{\alpha}_{\rm
bias}:={\rm{E}}[\hat{\theta}]-\tilde{\theta}^{\alpha}=\Sigma^{\alpha\beta}\big{(}\tilde{h}^{(1)}_{,\beta}\,,\tilde{h}^{(2)}\big{)}$
is called bias. Note this derivation has assumed stationary and Gaussian
noise, and used the relation ${\rm{E}}\left[(u\,,n)(n\,,v)\right]=(u\,,v)$
[18].

Equation (10) deserves a detailed discussion. Mathematically, the MSE can be
decomposed into the sum of the square of the bias, and the variance of the
estimator itself. The other signal, $\tilde{h}^{(2)}$, only contributes to the
bias part in the MSE via
$\big{(}\tilde{h}^{(1)}_{,\alpha}\,,\tilde{h}^{(2)}\big{)}$, while the
variance part has no relationship with $\tilde{h}^{(2)}$ and depends only on
the sensitivity of the waveform to the parameters (via
$\tilde{h}^{(1)}_{,\alpha}$) and the average strength of the noise $S_{n}$.

The MLE results are biased as long as
$\Sigma^{\alpha\beta}\big{(}\tilde{h}^{(1)}_{,\beta}\,,\tilde{h}^{(2)}\big{)}\neq
0$, however, this is not the condition for biased OSs. In practice, if the
bias $\Delta\theta^{\alpha}_{\rm bias}$ caused by $\tilde{h}^{(2)}$ is much
smaller than the statistical uncertainty $\Delta\theta^{\alpha}_{\rm stat}$
caused by $n$, we can not conclude the bias in the results at all, but instead
will attribute it to random errors caused by the noise or the sampling
algorithm. Therefore, it is more meaningful to define the (dimensionless)
reduced bias

$B_{\alpha}:=\frac{\Delta\theta^{\alpha}_{\rm
bias}}{\Delta\theta^{\alpha}_{\rm
stat}}=\frac{\Sigma^{\alpha\beta}\big{(}h^{(1)}_{,\beta}\,,h^{(2)}\big{)}}{\sqrt{\Sigma^{\alpha\alpha}}}\,.$
(11)

Note that in Eq. (11) only $\beta$ is summed. Also, we omit the tilde sign
since there is no essential difference whether the variable is $\bm{\theta}$
or $\tilde{\bm{\theta}}$. Once $h^{(1)}$ and $h^{(2)}$ are given, the reduced
bias $\bm{B}$ is determined by Eq. (11). Substituting $\bm{B}$ into the
expression of MSE, we obtain
${\rm{MSE}}_{\alpha}=\sigma^{2}\big{(}B_{\alpha}^{2}+1\big{)}$. Now, it is
easy to give a natural and simple criterion for biased OS: for a certain OS,
if there exists a component of $\bm{B}$ larger than $1$ or smaller than $-1$,
the OS is a biased OS; otherwise, the OS is an unbiased OS.

##  3 Bias in the Overlapping of Two Individual Signals

###  3.1 Parameter selection and detector configuration

Now we use the method developed in Sec. 2 to study the behaviors of $\bm{B}$.
Although the FM method quickly calculates the bias for given
${\bm{\theta}}^{(1)}$ and ${\bm{\theta}}^{(2)}$, exploring a high-dimensional
parameter space is still not an easy task. As mentioned earlier, the focus of
this work is not to simulate several OS PEs as realistically as possible, but
to discuss and explain the origin of SPE biases within a simple parameter
space.

When generating the OSs, we get the nonspinning stellar-mass BBH signals using
the IMRPhenomD waveform template for $h^{(1)}$ and $h^{(2)}$ [25, 30]. The
variable parameters are the chirp mass $\cal M$, the symmetry mass ratio
$\eta$, the luminosity distance $d_{L}$, and the merger time $t_{c}$. In
addition, we use the angle-averaged waveforms over the sky localization and
inclination, as well as the polarization angle of GWs. This corresponds to a
most conservative scenario, where two GW signals arrive at the detector from
the same direction with the same inclination and polarization. Taking into
account these angle parameters may reduce the biases [38], but this effect
originates more from geometry than the intrinsic properties of PE. As for the
merger phase $\phi_{c}$, some studies pointed out that the merging phase
difference can affect the appearance of biases significantly [34, 38], but no
further analysis has been done. Here we take the merger phase between two
signals $\Delta\phi_{c}=0$ temporarily. After establishing the analysis
process of biases, the $\Delta\phi_{c}\neq 0$ case can be easily incorporated
into our analysis, and how it affects the biases can also be clearly
explained. Finally, all parameters (especially the two mass parameters, ${\cal
M}$ and $\eta$) are defined in the detector frame, since the overlapping of
GWs occurs in the detector rather than at the GW source location.

As for the detector, we choose the designed sensitivity of
AdvLIGO111https://dcc.ligo.org/LIGO-T1800044/public. for the calculation of
the inner product. The reason for this selection is essentially the same as
[34]. On the one hand, although AdvLIGO is unlikely to detect OSs at present
[35, 38], this is due to the low total event rate in the detection horizon of
AdvLIGO. If two signals with close merging times are injected, a biased OS may
occur. Also, the effective duration of a signal in AdvLIGO is much shorter
than that in XG detectors, so the bias calculation is faster. On the other
hand, the shape of the AdvLIGO sensitivity curve is similar to those of XG
ground-based GW detectors (such as CE and ET), except that the sensitivity of
AdvLIGO is about one order of magnitude lower. This means that the behaviors
of an OS in AdvLIGO and XG detectors are similar, only differing roughly by a
constant factor. Therefore, we expect the properties of biases with AdvLIGO
sensitivity also apply to XG detectors. Moreover, as will be demonstrated in
the following sections, the behaviors of bias in AdvLIGO is consistent with
previous work, which considered OS in ET [22, 6]. This further increases the
confidence in our conclusions based on OS simulations in AdvLIGO.

As a short summary, when generating an OS, we have eight free parameters,
$\bm{\theta}=\\{{\cal M}^{(1)},\eta^{(1)},d_{L}^{(1)},t_{c}^{(1)},{\cal
M}^{(2)},\eta^{(2)},d_{L}^{(2)},t_{c}^{(2)}\\}$. The reduced bias
$\bm{B}=\\{B_{\cal M},B_{\eta},B_{d_{L}},B_{t_{c}}\\}$ is a $4$-dimensional
function defined on this $8$-dimensional space. For convenience, we use
indices $\alpha=\\{1,2,3,4\\}$ for $\\{{\cal M},\eta,d_{L},t_{c}\\}$
respectively. In the following, we will analyze the properties of $\bm{B}$ for
each parameter.

###  3.2 Luminosity distance, signal-to-noise ratio and merger time

In the IMRPhenomD waveform, the amplitude is inversely proportional to the
luminosity distance $d_{L}$, while the phase is independent of $d_{L}$. This
is relatively simple, so we analyze the dependence of $\bm{B}$ on
$d_{L}^{(1)}$ and $d_{L}^{(2)}$ first. For $d_{L}^{(2)}$, since
$h^{(2)}\propto 1/d_{L}^{(2)}$, according to Eq. (11) we have ${\bm{B}}\propto
1/d_{L}^{(2)}$. As for $d_{L}^{(1)}$, since for $\beta\neq 3$,

$h_{,\beta}^{(1)}\propto 1/d_{L}^{(1)}\,,\quad h_{,3}^{(1)}\propto
1/\big{(}d_{L}^{(1)}\big{)}^{2}\,,$ (12)

we have, for $\alpha,\,\beta\neq 3$,

$\Sigma^{\alpha\beta}\propto\big{(}d_{L}^{(1)}\big{)}^{2}\,,\quad\Sigma^{\alpha
3}\propto\big{(}d_{L}^{(1)}\big{)}^{3}\,,\quad\Sigma^{33}\propto\big{(}d_{L}^{(1)}\big{)}^{4}\,,$
(13)

from the definition of FM and its inverse [Eqs. (5) and (6)]. Substituting
them in Eq. (11), we find that
${\bm{B}}\propto\big{(}d_{L}^{(1)}\big{)}^{0}\big{(}d_{L}^{(2)}\big{)}^{-1}$.

The luminosity distance directly affects the amplitude of GWs, which in turn
affects the signal-to-noise ratio (SNR). For a GW signal $h$, its SNR is
$\rho=\sqrt{(h\,,h)}$. Since SNR is inversely proportional to $d_{L}$, we have
${\bm{B}}\propto\big{(}\rho^{(1)}\big{)}^{0}\big{(}\rho^{(2)}\big{)}^{1}$. At
first glance, this conclusion may seem strange, since the bias should be
smaller for louder $h^{(1)}$. Note that the PE uncertainty is also smaller for
louder signals, so the reduced bias $\bm{B}$ does not depend on the strength
of signal $h^{(1)}$. Similar results were mentioned by [6]. However, strictly
speaking, this conclusion holds only when $d_{L}$ varies or the sensitivity
has been consistently increased or decreased by an overall factor. Signals
with different internal parameters can have the same SNR, but lead to
different biases.

Intuitively, for the merger time, $\bm{B}$ should be related only to the
merger time difference $\Delta t_{c}=t_{c}^{(2)}-t_{c}^{(1)}$ of the two
signals. First, we decompose the frequency domain waveform into amplitude and
phase, $h(f)=A(f)e^{i\phi(f)}$. It is easy to know that only $\phi(f)$ depends
on $t_{c}$, and can be further written as $\phi(f)=\phi_{0}(f)-2\pi ft_{c}$,
where $\phi_{0}$ is the phase when $t_{c}=0$. On the other hand, $\Sigma$ is
only related to ${\cal M}^{(1)}$, $\eta^{(1)}$, and $d_{L}^{(1)}$, so $\bm{B}$
depends on $t_{c}^{(1)}$ and $t_{c}^{(2)}$ only through the inner product
$\big{(}h^{(1)}_{,\alpha}\,,h^{(2)}\big{)}$, which is written explicitly as

$\displaystyle\big{(}h^{(1)}_{,\alpha}\,,h^{(2)}\big{)}$
$\displaystyle=4\Re\int_{f_{\rm min}}^{{f_{\rm
max}}}\frac{A^{(1)}_{,\alpha}A^{(2)}-i\phi^{(1)}_{,\alpha}A^{(1)}A^{(2)}}{S_{n}(f)}e^{i\left(\phi_{0}^{(2)}-\phi_{0}^{(1)}-2\pi
f\big{(}t_{c}^{(2)}-t_{c}^{(1)}\big{)}\right)}{\rm d}f$
$\displaystyle=\int_{f_{\rm min}}^{{f_{\rm max}}}{\cal
A}^{c}_{\alpha}\cos\left(\Delta\phi_{0}-2\pi\Delta t_{c}f\right){\rm
d}f+\int_{f_{\rm min}}^{{f_{\rm max}}}{\cal
A}^{s}_{\alpha}\sin\left(\Delta\phi_{0}-2\pi\Delta t_{c}f\right){\rm d}f\,,$
(14)

where ${\cal A}^{c}_{\alpha}\equiv 4A^{(1)}_{,\alpha}A^{(2)}/S_{n}$, ${\cal
A}^{s}_{\alpha}\equiv 4\phi^{(1)}_{,\alpha}A^{(1)}A^{(2)}/S_{n}$, and
$\Delta\phi_{0}\equiv\phi_{0}^{(2)}-\phi_{0}^{(1)}$. The integration interval
$[f_{\rm min}\,,f_{\rm max}]$ is determined by the sensitivity curve of
AdvLIGO, independent of the OS parameters.

Mathematically, Eq. (3.2) is a modulated function integral, ${\cal
A}^{c}_{\alpha}$ and ${\cal A}^{s}_{\alpha}$ are the modulation amplitudes,
and the modulation phase $\Delta\phi_{0}$ is the phase difference between the
two signals at $t_{c}^{(1)}=t_{c}^{(2)}=0$ in the frequency domain. For the
IMRPhenomD template, $d_{L}$ appears only in the amplitude, so we have ${\cal
A}^{c}_{3}\neq 0$ and ${\cal A}^{s}_{3}=0$; $t_{c}$ appears only in the phase,
so ${\cal A}^{c}_{4}=0$ and ${\cal A}^{s}_{4}\neq 0$. The intrinsic parameters
${\cal M}$ and $\eta$ contribute to both amplitude and phase, but numerical
calculations show that ${\cal A}^{s}_{1/2}$ is about two orders of magnitude
larger than ${\cal A}^{c}_{1/2}$, namely ${\cal A}^{s}_{1/2}\gg{\cal
A}^{c}_{1/2}>0$ (see A). This is because that the current matched filtering
technique is more sensitive to the phase evolution of GWs. The phase
$\phi_{0}$ is only related to the internal parameters ${\cal M}$ and $\eta$,
and is generally a nonlinear function of $f$. As for $t_{c}^{(1)}$ and
$t_{c}^{(2)}$, they only contribute to the arguments of trigonometric
functions in the form of merger time difference $\Delta t_{c}$, which agrees
with our expectation. Moreover, the contribution to the phase is a linear term
of $f$.

According to the above discussion, we only need to explore a $5$-dimensional
parameter space, and the corresponding parameters are $\\{{\cal
M}^{(1)},\eta^{(1)},{\cal M}^{(2)},\eta^{(2)},t_{c}^{(2)}-t_{c}^{(1)}\\}$.
Therefore, unless specified, afterwards we fix $d_{L}^{(1)}=500\,{\rm{Mpc}}$,
$d_{L}^{(2)}=3000\,{\rm{Mpc}}$, and $t_{c}^{(1)}=0$.

###  3.3 Global dependence of biases on intrinsic parameters

The dependence of $\bm{B}$ on $d_{L}$ is only a simple rescale factor and does
not touch on the real origin of biased OSs. To further shrink the parameter
space, this subsection will take several sets of $\\{{\cal
M}^{(1)},\eta^{(1)},{\cal M}^{(2)},\eta^{(2)}\\}$ and study the generality of
the dependence of $\bm{B}$ on the internal parameters.

For $\\{{\cal M}^{(1)},\eta^{(1)}\\}$, we will take $5$ sets of parameters,
expressed in terms of binary masses, with the Solar mass ${\rm M}_{\odot}$ as
unit. The first set of parameters is intentionally taken as integers
$\big{(}m_{1}^{(1)},m_{2}^{(1)}\big{)}=(30,20)$. We call this configuration
the “Test Event.” For the other $4$ sets of parameters, we select the real GW
events GW150914, GW190602 175927, GW190814, and GW190924 021846\. GW150914 is
the first GW event. GW190602 175927, GW190814, and GW190924 021846 are the
events with the second largest chirp mass, the smallest mass ratio, and the
lightest chirp mass, respectively, in GWTC-1 and GWTC-2. Finally, for each set
of the above parameters, we calculate $\bm{B}$ for $\Delta t_{c}\in[t_{{\rm
min}}\,,t_{{\rm max}}]=[-0.1\,,0.1]\,{\rm s}$. Most of the biased OSs occur in
this range of $\Delta t_{c}$ [34, 38, 22], which is also verified by the
calculations below.

This subsection focuses on the dependence of $\bm{B}$ on internal parameters,
which means we need to choose some characteristic quantities to reflect the
distribution ${\bm{B}}$ in the range $\Delta t_{c}\in[-0.1\,,0.1]\,{\rm s}$.
In Fig. 1, we show the dependence of ${\bm{B}}$’s “max-max” and “max-mean”
values on $m_{1}^{(2)}$ and $m_{2}^{(2)}$ when
$\big{(}m_{1}^{(1)},m_{2}^{(1)}\big{)}=(30,20)$. To obtain the max-max (max-
mean) values, we first calculate the maximum (mean) of $|B_{\alpha}|$ over
$[t_{{\rm min}}\,,t_{{\rm max}}]$, then find the maximum value in the $4$
(dimensionless) components of $\bm{B}$. We take the absolute value of
$B_{\alpha}$ to identify biased OSs. In addition, to eliminate the effects of
signal strength on biases, we adjust $d_{L}^{(2)}$ so that $\rho^{(2)}=8$,
which is the detection threshold of the ground-based detectors. The results of
the other four sets of $m_{1}^{(1)}$ and $m_{2}^{(1)}$ are shown in B.

Figure 1: Dependence of the max-max value, max-mean value and average
asymmetry of ${\bm{B}}$ in the parameter space of $m_{1}^{(2)}$ and
$m_{2}^{(2)}$ for the Test Event. The location of
$\big{(}m_{1}^{(1)},m_{2}^{(1)}\big{)}$ is marked with green stars. The
magenta dashed lines indicate the values of
$\big{(}m_{1}^{(2)},m_{2}^{(2)}\big{)}$ for which ${\cal M}^{(2)}={\cal
M}^{(1)}$.

We find that for a wide range of $m_{1}^{(2)}$ and $m_{2}^{(2)}$, the max-max
value of $\bm{B}$ is greater than $1$. This means that for almost any two GW
signals, there exists a $\Delta t_{c}$ making the corresponding OS biased. For
the max-mean value, it is greater than $1$ only in the region around ${\cal
M}^{(2)}={\cal M}^{(1)}$. Therefore, in the average sense, there is a higher
probability of biased OS only when the chirp masses of the two signals are
close. It should be noted that since $\bm{B}$ is larger at $|\Delta t_{c}|\sim
0$ and decays to zero as $|\Delta t_{c}|$ increases, if the length of the
$\Delta t_{c}$ interval is shortened, e.g., to $|\Delta
t_{c}|\in[-0.05,0.05]\,{\rm s}$, the max-mean value will increase. However,
the global behaviors of $\bm{B}$ remain unchanged. Also, the max-max and max-
mean values involve the absolute size of the bias, and we fix the SNR
$\rho^{(2)}=8$.

In addition to the absolute magnitude of the biases, we are also interested in
their distribution with respect to $\Delta t_{c}$. One of the most important
questions is whether there is a significant difference in the biases when
$\Delta t_{c}>0$ and $\Delta t_{c}<0$, i.e., the dependence on the merger
order. First, we define ${\rm Asy}_{\alpha}$ as the difference between the
area enclosed by the horizontal axis $|B_{\alpha}|$ on $[-0.1,0]\,{\rm s}$ and
the area on $[0,0.1]\,{\rm s}$,

${\rm Asy}_{\alpha}\coloneqq\frac{\int_{0}^{t_{\rm max}}|B_{\alpha}|\,{\rm
d}\Delta t_{c}-\int^{0}_{t_{\rm min}}|B_{\alpha}|\,{\rm d}\Delta
t_{c}}{\int_{0}^{t_{\rm max}}|B_{\alpha}|\,{\rm d}\Delta
t_{c}+\int^{0}_{t_{\rm min}}|B_{\alpha}|\,{\rm d}\Delta t_{c}}\,.$

Then, we average ${\rm Asy}_{\alpha}$ over the parameters of
$\bm{\theta}^{(1)}$ to obtain the average asymmetry $\overline{\rm Asy}$.
$\overline{\rm Asy}$ is a function of internal parameters $\\{{\cal
M}^{(1)},\eta^{(1)},{\cal M}^{(2)},\eta^{(2)}\\}$. When $\overline{\rm
Asy}>0$, it means that the OS is more likely to be biased when $h^{(2)}$
merges after $h^{(1)}$. In addition, the closer $|\overline{\rm Asy}|$ is to
$1$, the stronger the dependence of the bias on the merger order. For example,
$\overline{\rm Asy}=0.8=(9-1)/(9+1)$ means that the average bias on cases
where $h^{(1)}$ merges first is $9$ times of that where $h^{(2)}$ merges
first. We calculate the average asymmetries of the above five sets of GW
events and put the results of the Test Event and the other four real GW events
in Fig. 1 and B respectively. There is a clear boundary between $\overline{\rm
Asy}>0$ and $\overline{\rm Asy}<0$. The boundary coincides with the ${\cal
M}^{(2)}={\cal M}^{(1)}$ curve. Roughly speaking, it is easier to generate a
biased OS when the signal with heavier masses merges first. In addition, we
found that the regions of small $|\overline{\rm Asy}|$ (the white regions in
the average asymmetry plot) are consistent with the regions of large bias (the
dark blue regions in the max-max, max-mean plots).

Table 1: Parameter configurations of the six cases of $h^{(2)}$, where $M=m_{1}+m_{2}$ is the total mass of binary. For $h^{(1)}$, we have $m_{1}=30.00\,{\rm M}_{\odot}$, $m_{2}=20.00\,{\rm M}_{\odot}$ and $d_{L}=500\,{\rm Mpc}$.  Name | $m_{1}\,({\rm M}_{\odot})$ | $m_{2}\,({\rm M}_{\odot})$ | $M\,({\rm M}_{\odot})$ | ${\cal M}\,({\rm M}_{\odot})$ | $\eta$ | $d_{L}\,({\rm Mpc})$ | $\rho$
---|---|---|---|---|---|---|---
Equal | 30.00 | 20.00 | 50.00 | 21.24 | 0.240 | 3000 | 5.5
Asymmetric | 27.00 | 18.00 | 45.00 | 19.11 | 0.240 | 3000 | 5.1
Asymmetric2 | 33.00 | 22.00 | 55.00 | 23.36 | 0.240 | 3000 | 5.9
Random | 40.00 | 34.00 | 74.00 | 32.08 | 0.248 | 3000 | 7.7
Symmetric and close | 33.00 | 18.20 | 51.20 | 21.15 | 0.229 | 3000 | 5.4
Symmetric and not close | 60.00 | 9.11 | 69.11 | 18.82 | 0.114 | 3000 | 3.8

###  3.4 Origin of biases

In Sec. 3.3, we vary
$\big{\\{}m_{1}^{(1)},m_{2}^{(1)},m_{1}^{(2)},m_{2}^{(2)}\big{\\}}$ and
summarize the generality of the behaviors of $\bm{B}$ with different internal
parameters. Now we present the mathematical characteristics and physical
interpretation of $B_{\alpha}\big{(}\Delta t_{c}\big{)}$ and draw the most
important conclusion in our work: the SPE biases originate from the
overlapping of frequency evolution processes of the two signals. To obtain
this, we fix the internal parameters and investigate the relationship between
$\bm{B}$ and $\Delta t_{c}$. For $h^{(1)}$, we take ${\bm{\theta}}^{(1)}$ as
the parameter in the Test Event. For $h^{(2)}$, we take six representative
sets of parameters, called Equal, Asymmetric, Asymmetric2, Random, Symmetric
and close, Symmetric and not close, according to the relationship between
${\bm{\theta}}^{(2)}$ and ${\bm{\theta}}^{(1)}$ and the behavior of
${\bm{B}}\big{(}\Delta t_{c}\big{)}$. The parameters are listed in Table 1. In
the main text, we only discuss the first two configurations to draw the main
conclusion. The other four additional configurations provide a more
comprehensive study in the parameter space, and we discuss them in the C.

Figure 2: The dependence of $\big{(}h^{(1)}_{,\alpha}\,,h^{(2)}\big{)}$ and
$|B_{\alpha}|=|\Delta\theta^{\alpha}_{\rm bias}|/\Delta\theta^{\alpha}_{\rm
stat}$ on $\Delta t_{c}$ for Equal and Asymmetric configurations. In the
shaded area, $|B_{\alpha}|<1$, which corresponds to unbiased OSs; the other
region corresponds to biased OSs. The magenta dashed lines mark some
representative $\Delta t_{c}$ for reference.

According to Eqs. (11) and (3.2), the reduced bias $B_{\alpha}$ is just a
linear combination of $\big{(}h^{(1)}_{,\alpha}\,,h^{(2)}\big{)}$, and the
combination coefficients only depend on ${\bm{\theta}}^{(1)}$. When $h^{(1)}$
is fixed, a large bias is roughly equivalent to a large
$\big{(}h^{(1)}_{,\alpha}\,,h^{(2)}\big{)}$, and the latter can be expressed
as the sum of the integrals of two modulated trigonometric functions, which is
easier to analyze. Thus, for each parameter configuration of $h^{(2)}$, we
vary $\Delta t_{c}$ within $[-0.1,0.1]\,{\rm s}$, calculate
$\big{(}h^{(1)}_{,\alpha}\,,h^{(2)}\big{)}$ and the corresponding
${B}_{\alpha}$, and display them in Fig. 2.

For the convenience of presentation, we define

$\displaystyle D_{\alpha}(f)\coloneqq$ $\displaystyle{\cal
A}^{c}_{\alpha}\cos\left(\Delta\phi_{0}-2\pi\Delta t_{c}f\right)+{\cal
A}^{s}_{\alpha}\sin\left(\Delta\phi_{0}-2\pi\Delta t_{c}f\right)\,,$
$\displaystyle R_{\alpha}(f)\coloneqq$
$\displaystyle\Sigma^{\alpha\beta}D_{\alpha}/\sqrt{\Sigma^{\alpha\alpha}}\,,$

which are the integrand functions of
$\big{(}h^{(1)}_{,\alpha},\,h^{(2)}\big{)}$ and ${B}_{\alpha}$, respectively.
According to the dependence of ${\cal A}^{c}$ and ${\cal A}^{s}$ on the
parameters, we have

$\displaystyle D_{1}$ $\displaystyle\approx{\cal
A}^{s}_{1}\sin\Delta\phi\,,\qquad D_{2}\approx{\cal
A}^{s}_{2}\sin\Delta\phi\,,$ (15) $\displaystyle D_{3}$
$\displaystyle\approx{\cal A}^{c}_{3}\cos\Delta\phi\,,\qquad D_{4}\approx{\cal
A}^{s}_{4}\sin\Delta\phi\,,$ (16)

where $\Delta\phi=\Delta\phi_{0}-2\pi\Delta t_{c}f$.

When $m_{1}^{(1)}$ and $m_{2}^{(1)}$ are fixed, the most special OS consists
of two GWs with the same internal parameters, i.e.,
$\big{(}m_{1}^{(1)},m_{2}^{(1)}\big{)}=\big{(}m_{1}^{(2)},m_{2}^{(2)}\big{)}$.
We call this configuration Equal. In this case, $\Delta\phi_{0}=0$, and Eq.
(3.2) degenerates to

$\displaystyle\big{(}h^{(1)}_{,\alpha}\,,h^{(2)}\big{)}=$
$\displaystyle\int_{f_{\rm min}}^{{f_{\rm max}}}{\cal
A}^{c}_{\alpha}\cos\left(2\pi\Delta t_{c}f\right){\rm d}f-\int_{f_{\rm
min}}^{{f_{\rm max}}}{\cal A}^{s}_{\alpha}\sin\left(2\pi\Delta
t_{c}f\right){\rm d}f\,.$

Thus, the modulation phase becomes a linear function of $f$, whose slope is
proportional to $\Delta t_{c}$.

In Fig. 2, $\big{(}h^{(1)}_{,3}\,,h^{(2)}\big{)}$ is an even function of
$\Delta t_{c}$, while the other three components behave like odd functions.
The $d_{L}$ only contributes to the amplitude part of the waveform, then
${\cal A}^{s}_{3}=0$. Therefore, $D_{3}(f;-\Delta t_{c})=D_{3}(f;\Delta
t_{c})$ is an even function of $\Delta t_{c}$, and remains an even function
after integration over $f$. Noting that $t_{c}$ only contributes to the phase,
we have ${\cal A}^{c}_{4}=0$, so $\big{(}h^{(1)}_{,4}\,,h^{(2)}\big{)}$ is an
odd function. For ${\cal M}$ and $\eta$, we have ${\cal A}^{s}\gg{\cal
A}^{c}$. So $D_{\alpha}$ is dominated by ${\cal A}^{s}$, and both
$\big{(}h^{(1)}_{,1}\,,h^{(2)}\big{)}$ and
$\big{(}h^{(1)}_{,2}\,,h^{(2)}\big{)}$ are almost odd functions. As for
$\left|B_{\alpha}\right|$, they are all basically even functions of $\Delta
t_{c}$. The slight asymmetry originates from the linear combination of
$\big{(}h^{(1)}_{,\alpha}\,,h^{(2)}\big{)}$. Physically, the intrinsic
parameters of the two signals are the same, and only the amplitudes are
different. Therefore, the relationship between biases and merger time
difference will inherit the symmetry.

Another thing is that the biases are large at $\Delta t_{c}\approx 0$, but
small when $|\Delta t_{c}|$ is far away from $0$. Obviously, the larger the
difference in the merger times is, the smaller the bias will be. To explain
this point quantitatively, we take some $\Delta t_{c}$ and draw the
corresponding $\Delta\phi(f)$, $D_{\alpha}$, and $R_{\alpha}$ in Fig. 3. When
$\Delta t_{c}=0$, $|B_{d_{L}}|>5$, while the other three parameters are almost
unbiased. In this case, $D_{\alpha}={\cal A}^{c}_{\alpha}$, and only ${\cal
A}^{c}_{3}$ significantly deviates from $0$. In the Equal configuration and
$\Delta t_{c}=0$, the two signals only have the difference in amplitude,
$h^{(1)}+h^{(2)}\propto h^{(1)}$. Therefore, the existence of $h^{(2)}$ is
equivalent to slightly increasing the amplitude of $h^{(1)}$, leading to a
large impact on the PE of $d_{L}$. The difference between $\Delta
t_{c}=2.1\,{\rm ms}$ and $\Delta t_{c}=-37\,{\rm ms}$ is only reflected in the
slope of $\Delta\phi$ with respect to $f$. For the case of small $\Delta
t_{c}$, $\Delta\phi$ changes slowly by $\pi$ in the range of $0\sim 200\,{\rm
Hz}$. For all parameters, $D_{\alpha}$ is significantly nonzero, which
eventually leads to large biases. However, in the $\Delta t_{c}=-37\,{\rm ms}$
case, the rapid increase of $\Delta\phi$ leads to violent oscillations of
$D_{\alpha}$, which can not generate large biases after integration. This
quantitatively explains why the biases are small for large $\Delta t_{c}$. At
last, it is worth noting that $\Delta\phi$ always contains the $2\pi\Delta
t_{c}f$ term that grows linearly with frequency, so this behavior can be
applied to all mass configurations.

Figure 3: The dependence of $\Delta\phi(f)$, $D_{\alpha}$, and $R_{\alpha}$ on
the frequency $f$ in the Equal configuration. The merger time differences are
also marked. Note that the vertical variable ranges, such as $\Delta\phi$, in
each subfigure are different. In the second row of each subfigure, the shaded
area represents $-\pi/2<\Delta\phi(f)<\pi/2$, while the hatched area means
$0<\Delta\phi(f)<\pi$. In the third row, the number after each parameter in
the legend is the maximal value of $|D_{\alpha}|$. In the forth row, the
number means the integration of $R_{\alpha}$ over $f$, i.e. $B_{\alpha}$.

In the Equal case, $\Delta\phi$ is completely determined by $\Delta t_{c}$,
but for general mass configurations, we have to consider the contribution of a
non-zero $\Delta\phi_{0}$. In the Asymmetric (Asymmetric2) configuration, we
decrease (increase) $m_{1}^{(2)}$ and $m_{2}^{(2)}$ simultaneously, while
keeping the mass ratio unchanged. Taking Asymmetric as an example, we
immediately find the asymmetry of $\bm{B}$ on $\Delta t_{c}$, i.e., the
significant dependence of the merging order. When $h^{(2)}$ merges after
$h^{(1)}$, the biases are generally large. Otherwise, the biases can be
ignored. Same as the Equal case, we also take some $\Delta t_{c}$ as examples
and show them in Fig. 4. Intuitively, we do expect the largest bias when two
signals are “closest.” However, from the $B_{\alpha}\big{(}\Delta
t_{c}\big{)}$ plot, we find that when $\Delta t_{c}=0$ the bias is not
maximum. In the following, we discuss how these larger biases arise, which
correspond to those OSs with $\Delta t_{c}>0$.

Figure 4: Same as Fig. 3, but for the Asymmetric configuration.

Mathematically, $D_{\alpha}$ is a modulated trigonometric function. The
modulation amplitude is ${\cal A}^{c}$ (${\cal A}^{s}$), and the phase
argument is $\Delta\phi=\Delta\phi_{0}-2\pi\Delta t_{c}f$. In order to
generate a significant non-zero integral, one needs a large modulation
amplitude and a stable $\Delta\phi(f)$. Since ${\cal A}^{c}$ (${\cal A}^{s}$)
is large in the low-frequency band and slowly decays to $0$ in the high-
frequency band, combined with the behavior of $\Delta\phi_{0}$, we can divide
the large bias cases into two categories according to the frequency band as
follows.

  1. (i)

Low-frequency band. In this case, $\Delta\phi_{0}$ changes drastically, but
adding a larger linear term can make $\Delta\phi_{0}$ stable in a small range.
Because ${\cal A}^{c}$ (${\cal A}^{s}$) is relatively large in the low-
frequency band, a large bias can also be integrated in a short frequency
range. In order to offset the drastic changes in $\Delta\phi_{0}$, the
corresponding $|\Delta t_{c}|$ is large when biased in this way.

  2. (ii)

High-frequency band. When $f$ is high, $\Delta\phi_{0}$ changes approximately
linearly, so it is easy to keep $\Delta\phi(f)$ to be a constant in a large
frequency range by adding a linear term. Although ${\cal A}^{c}$ (${\cal
A}^{s}$) is small at this time, a large bias can still be accumulated when
integrating in a wide frequency band where $\Delta\phi(f)$ is stable. Since
the change of $\Delta\phi_{0}$ in the high-frequency band is slow, the
corresponding $|\Delta t_{c}|$ is also small.

In Fig. 4, we give examples of the two cases. When $\Delta t_{c}=3.1\,{\rm
ms}$, due to the existence of the linear term $-2\pi\Delta t_{c}f$,
$\Delta\phi$ first increases and then decreases in the $50\sim 150\,{\rm Hz}$
band. The change of monotonicity implies that $\Delta\phi$ basically stays
around $\pi/2$ in the corresponding frequency band, leading to a large bias
after integration. As for $\Delta t_{c}=34.7\,{\rm ms}$, the monotonicity of
$\Delta\phi$ changes in the $30\sim 60\,{\rm Hz}$ band, producing a short
stable stage for $D_{\alpha}$ and $R_{\alpha}$. In this frequency band, both
${\cal A}^{c}$ and ${\cal A}^{s}$ also reach their maxima, so large biases can
occur.

To sum up, as long as there is a relatively stable stage in the phase argument
$\Delta\phi$, biased OSs are likely to arise. The location and even the
presence of the stable stage are determined by $\Delta t_{c}$. With this in
mind, we can explain why $B_{\alpha}$ oscillates and tends to $0$ with the
increase of $|\Delta t_{c}|$. Similar phenomenon has been noticed before, but
not explored in detail [22]. Taking $d_{L}$ as an example, the growth rate of
$\Delta\phi_{0}$ about $f$ gradually decreases, while the linear term
$-2\pi\Delta t_{c}f$ contributes a constant negative growth rate. Therefore,
increasing $\Delta t_{c}$ will cause the stable stage to shift left. On the
other hand, since $\Delta\phi_{0}$ is changing between $-\pi$ and $\pi$ (or
equivalently, between $0$ and $2\pi$), the change of the horizontal position
of the stable stage will be accompanied by the change of its vertical
position. This means that, the stable value of $\Delta\phi$, which is denoted
as $\phi_{\rm sta}$, will oscillate between $-\pi$ and $\pi$. For $d_{L}$, in
the stable stage, $D_{3}={\cal A}^{c}_{3}\cos\Delta\phi\approx{\cal
A}^{c}_{3}\cos\phi_{\rm sta}$. If $\phi_{\rm sta}\approx 0$ or $\pi$, we have
$D_{3}\approx 0$ and $\big{(}h^{(1)}_{,4}\,,h^{(2)}\big{)}\approx 0$ after
integration. As for $B_{3}$, it is a linear combination of
$\big{(}h^{(1)}_{,\alpha},h^{(2)}\big{)}$, and there still exists $\phi_{\rm
sta}$ making $B_{3}\approx 0$. Therefore, when $\Delta t_{c}$ increases from
$0$, the stable stage moves to the left, and $\phi_{\rm sta}$ oscillates
between $-\pi$ and $\pi$, so $B_{3}$ also oscillates. When $\Delta t_{c}$
further increases, the stable stage is very short, and the ${\cal A}^{c}$
(${\cal A}^{s}$) in this frequency band is small, then the cumulative bias is
small in the stable stage. Outside the stable stage, $\Delta\phi$ changes
rapidly with the frequency, and it is impossible to generate a large bias.
Finally, we have ${\bm{B}}\to 0$.

As for the other three parameters, their corresponding $D_{\alpha}$ is
dominated by ${\cal A}^{s}\sin\Delta\phi$. When $D_{3}=0$, the other
$D_{\alpha}$ corresponding to $\alpha=1,2,4$ will reach the extreme value.
Therefore, when the bias of $d_{L}$ is very small, the bias of the other three
parameters is close to the maximum and vice versa. This is consistent with the
results of numerical calculation, as shown in Fig. 2. Further, it is easy to
explain why the biases are generally small when $\Delta t_{c}<0$. Taking
$\Delta t_{c}=-13.1\,{\rm ms}$ as an example (see Fig. 4), both $-2\pi\Delta
t_{c}f$ and $\Delta\phi_{0}$ are monotonically increasing with respect to $f$,
and the phase argument changes more rapidly. So the oscillation of
$D_{\alpha}$ is more violent than that when $\Delta t_{c}>0$, and there is no
stable stage and large bias. This explains the strong dependence of
$B_{\alpha}$ on the merge order of the two signals.

Using the fact that biases are mainly generated at the stable stage, it is
easy to analyze the dependence of biases on $\phi_{c}$. The numerical results
of [34] indicate that when $|\Delta t_{c}|$ is small the bias is very
sensitive to the difference of merger phase
$\Delta\phi_{c}=\phi_{c}^{(2)}-\phi_{c}^{(1)}$.222Note that the merger phase
is originally defined in the time domain. One can even turn a BOS into an UOS
by adjusting $\Delta\phi_{c}$. Here we give an explanation. Adding the
parameter $\Delta\phi_{c}\neq 0$, the phase argument becomes
$\Delta\phi=\Delta\phi_{0}+\Delta\phi_{c}-2\pi\Delta t_{c}f$, and the stable
value becomes $\phi_{\rm sta}+\Delta\phi_{c}$. If $\phi_{\rm
sta}+\Delta\phi_{c}$ is close to $0$ or $\pi$, then the $D_{\alpha}$ of ${\cal
M}$, $\eta$, and $t_{c}$ are also very close to $0$. Even if there is a long
stable stage, it is difficult to generate a large bias. On the contrary, the
bias of $d_{L}$ will be close to the extremum. Similar analysis can be applied
to $\phi_{\rm sta}+\Delta\phi_{c}\approx\pm\pi/2$, where the bias of $d_{L}$
is small while biases of ${\cal M}$, $\eta$, and $t_{c}$ reach their extrema.

Figure 5: The whitened waveforms and frequency evolution diagrams of two
signals in the time domain for the Asymmetric configuration when $\Delta
t_{c}$ takes $3.1\,\,{\rm ms}$ and $34.7\,\,{\rm ms}$.

Through the above analysis, we have understood the origin of biases
mathematically. Physically, the biases come from the overlapping of the
frequency evolution processes of the two signals, that is, in the same time
period, the frequencies of the two signals are relatively close to each other.
As an example, in the Asymmetry configuration, we find that it is easier to
generate BOS when $h^{(2)}$ merges before $h^{(1)}$. The chirp frequency of
the lighter binary is higher, so the frequency evolution processes are easier
to “overlap” when it merges after the heavier one. Further, in Fig. 5, we show
the whitened waveform, $h_{\rm w}$, and the frequency evolution, $f(t)$, of
two signals in the time domain for some representative $\Delta t_{c}$. For the
case of $\Delta t_{c}=3.1\,\rm{ms}$ ($34.7\,\rm{ms}$), we note that the
frequency band $50\sim 150\,{\rm Hz}$ ($30\sim 60\,{\rm Hz}$), where the chirp
frequency of the two signals are relatively close, is also the frequency band
where the phase difference in frequency domain $\Delta\phi$ is stable. The
physical and mathematical explanations of the origin of biases are (and they
must be) consistent. Generally speaking, the frequency of GW changes rapidly
in the merger and ringdown stages, and only in the inspiral stage can there be
a long overlapping period of the frequency evolution processes. On the other
hand, in the inspiral phase, the frequency domain phase and the frequency
evolution have the relationship [12]

$\frac{{\rm d}\phi}{{\rm d}f}=-2\pi t(f)\,.$ (17)

Applying to $h^{(1)}$ and $h^{(2)}$, we have

$\frac{{\rm d}\Delta\phi}{{\rm d}f}\approx 0\Leftrightarrow\frac{{\rm
d}\phi^{(1)}}{{\rm d}f}\approx\frac{{\rm d}\phi^{(2)}}{{\rm
d}f}\Leftrightarrow t^{(1)}(f)\approx t^{(2)}(f)\,,$ (18)

which shows that the phase difference in the frequency domain keeps constant
in a certain frequency band and that the frequency evolution processes overlap
in the time domain are equivalent. They just explain the origin of bias from
different aspects. Both explanations are consistent with our intuitive
understanding of the matched filtering technology. If the frequency evolution
over time of the two signals is completely the same, it is impossible to
separate the OS into individual ones. The strong degeneracy will influence the
SPE, leading to large biases.

##  4 Other methods of analyzing biases

In this section, we verify the results and corresponding conclusions of biases
in Sec. 3, using the correlation coefficients and the full Bayesian inference.
We also point out the feasibility of identifying biased OSs and unbiased OSs
with FM in our analysis.

###  4.1 The correlation between two signals

In the previous section, we have used the SPE biases to characterize the
effect of signal overlapping. When large biases occur, it indicates that the
two signals have strong influences on the PEs of each other, i.e., the
correlation between the two signals is large. [34] briefly studied the
relationship between the correlation coefficient of the parameters of the two
signals, as well as $\Delta t_{c}$ under some specific configurations.
However, they did not discuss the relationship between the correlation
coefficients and the biases. [6] calculated some correlation coefficients and
biases, but the result is limited to a rough qualitative analysis due to the
high dimensionality of the parameter space. Here we calculate the correlation
coefficients between the parameters of the two signals by the FM method, and
focus on the relationship between the correlation coefficients and the biases.
Here we only consider the correlation coefficients between the corresponding
parameters of the two signals (for example, ${\cal M}^{(1)}$ and ${\cal
M}^{(2)}$), which represent the main correlation between the two signals. The
systematic analysis of the correlation coefficient between any two parameters
(e.g., between ${\cal M}^{(1)}$ and $\eta^{(2)}$) is left to future work.

When calculating the correlation coefficients of the parameters between the
two signals, the FM approximation should be applied to the JPE. The parameters
are $\\{\bm{\theta}^{(1)},\bm{\theta}^{(2)}\\}$, while the model takes
$h\big{(}\bm{\theta}^{(1)},\bm{\theta}^{(2)}\big{)}=h^{(1)}\big{(}\bm{\theta}^{(1)}\big{)}+h^{(2)}\big{(}\bm{\theta}^{(2)}\big{)}$.
According to Eq. (5), the FM is

$F=\begin{bmatrix}F^{(1)}&M\\\ M^{\intercal}&F^{(2)}\end{bmatrix}\,,$ (19)

where $F^{(1)}$ and $F^{(2)}$ are the FM calculated when the parameters are
taken as $\\{\bm{\theta}^{(1)}\\}$ and $\\{\bm{\theta}^{(2)}\\}$,
respectively. The inverses of them are noted as $\Sigma^{(1)}$ and
$\Sigma^{(2)}$. For the off-diagonal term, we have
$M_{\alpha\beta}=\big{(}\partial^{(1)}_{\alpha}h^{(1)}\,,\partial^{(2)}_{\beta}h^{(2)}\big{)}=\big{(}h^{(1)}_{,\alpha}\,,h^{(2)}_{,\beta}\big{)}$,
and $M^{\intercal}$ is the transposition of $M$. When $|\Delta t_{c}|$ is
large, the correlation between the two signals is small, so each element of
$M$ will be much smaller than those of $F^{(1)}$ and $F^{(2)}$. In this case,
the inverse of $F$ can be written as

$\Sigma=\begin{bmatrix}\Sigma^{(1)}&-\Sigma^{(1)}M\Sigma^{(2)}\\\
-\Sigma^{(2)}M^{\intercal}\Sigma^{(1)}&\Sigma^{(2)}\end{bmatrix}+{\cal
O}\big{(}M^{2}\big{)}\,.$ (20)

Ignoring the ${\cal O}\big{(}M^{2}\big{)}$ terms, the correlation coefficients
between the corresponding parameters of the two signals are

$c_{\theta^{(1)\alpha}\theta^{(2)\alpha}}=c_{\alpha,\alpha+4}=\frac{\Sigma^{(1)\alpha\beta}\big{(}\partial^{(1)}_{\beta}h^{(1)}\,,\partial^{(2)}_{\gamma}h^{(2)}\big{)}\Sigma^{(2)\gamma\alpha}}{\sqrt{\Sigma^{(1)\alpha\alpha}\Sigma^{(2)\alpha\alpha}}}\,.$
(21)

Note that $\alpha$ is not summed in the above equation, and both
$\\{\bm{\theta}^{(1)}\\}$ and $\\{\bm{\theta}^{(2)}\\}$ have $4$ independent
parameters. Correlation coefficients $c_{\alpha,\alpha+4}$ are linear
combinations of the inner products
$\big{(}h^{(1)}_{,\beta}\,,h^{(2)}_{,\gamma}\big{)}$, and the combination
coefficients are independent of $\Delta t_{c}$.

The expression of the correlation coefficients in Eq. (21) is very similar to
the expression of the biases in Eq. (11), except for the element in the inner
product. Therefore, we can use the analysis method in Sec. 3 to deal with
them. According to Fig. 2, we find $B_{\alpha}$ and
$\big{(}h^{(1)}_{,\alpha}\,,h^{(2)}\big{)}$ have similar behaviors with
respect to $\Delta t_{c}$, so are the positions of zeros and extrema.
Mathematically, this means that using $\Sigma^{(1)\alpha\beta}$ to linearly
combine $\big{(}h^{(1)}_{,\alpha}\,,h^{(2)}\big{)}$ is similar to using a
diagonal matrix. Applying this conclusion to Eq. (21), we only need to analyze
the dependence of $\big{(}h^{(1)}_{,\beta}\,,h^{(2)}_{,\gamma}\big{)}$ on
$\Delta t_{c}$.

Like in Eq. (3.2), we explicitly write the expression of
$\big{(}h^{(1)}_{,\beta}\,,h^{(2)}_{,\gamma}\big{)}$

$\displaystyle\big{(}h^{(1)}_{,\beta}\,,h^{(2)}_{,\gamma}\big{)}=$
$\displaystyle 4\int_{f_{\rm min}}^{f_{\rm
max}}\frac{A_{,\beta}^{(1)}A_{,\gamma}^{(2)}+A^{(1)}A^{(2)}\phi_{,\beta}^{(1)}\phi_{,\gamma}^{(2)}}{S_{n}}\cos\left(\Delta\phi_{0}-2\pi\Delta
t_{c}f\right){\rm d}f$ (22) $\displaystyle-4\int_{f_{\rm min}}^{f_{\rm
max}}\frac{A_{,\beta}^{(1)}A^{(2)}\phi_{,\gamma}^{(2)}+A^{(1)}A^{(2)}_{,\gamma}\phi_{,\beta}^{(1)}}{S_{n}}\sin\left(\Delta\phi_{0}-2\pi\Delta
t_{c}f\right){\rm d}f\,,$

where $\Delta\phi_{0}$ is the same as in Eq. (3.2). Since the linear
combinatorial effects of $\Sigma^{(1)}$ and $\Sigma^{(2)}$ are equivalent to
diagonal matrices, we only need to consider the contribution of the diagonal
elements. Take $\beta=\gamma=\alpha$, substitute them in Eq. (3.2) we have
$c_{\alpha,\alpha+4}\propto\big{(}h^{(1)}_{,\alpha}\,,h^{(2)}_{,\alpha}\big{)}$.
For the luminosity distance, it does not contribute to the phase, i.e.
$\phi_{,3}=0$. For the merger time, $A_{,4}=0$. As for ${\cal M}$ and $\eta$,
we have $A_{,\alpha}\ll A\phi_{,\alpha}$. Therefore, in the special case of
$\beta=\gamma=\alpha$, the integral of $\sin\Delta\phi$ in Eq. (22) can always
be insignificant, that is

$\big{(}h^{(1)}_{,\beta}\,,h^{(2)}_{,\gamma}\big{)}\approx 4\int_{f_{\rm
min}}^{f_{\rm
max}}\frac{A_{,\beta}^{(1)}A_{,\gamma}^{(2)}+A^{(1)}A^{(2)}\phi_{,\beta}^{(1)}\phi_{,\gamma}^{(2)}}{S_{n}}\cos\left(\Delta\phi_{0}-2\pi\Delta
t_{c}f\right){\rm d}f\,,$

which is the modulation integral of $\cos\Delta\phi$. Therefore, we can expect
the dependence of $c_{\alpha,\alpha+4}$ and $B_{3}$ (the bias of luminosity
distance) on $\Delta t_{c}$ to be almost the same.

Figure 6: The dependence of biases and correlation coefficients on $\Delta
t_{c}$ in the Symmetric and not close configuration. The behaviors of
$\big{(}h^{(1)}_{,\alpha}\,,h^{(2)}\big{)}$ and $1/\sqrt{1-c^{2}}$ are also
shown for comparison.

In Fig. 6, we show the dependence of
$\big{(}h^{(1)}_{,\alpha}\,,h^{(2)}\big{)}$, correlation coefficients
$c_{\alpha,\alpha+4}$, and $B_{\alpha}$ on the merger time difference $\Delta
t_{c}$. For the sake of brevity, we only show the results of two
configurations, Asymmetric and Symmetric and not close, which can clearly
reflect the oscillation of $B_{\alpha}$. When $|c|$ is small ($\lesssim 0.5$),
the behaviors of $c_{\alpha,\alpha+4}$ versus $\Delta t_{c}$ are highly
consistent with $B_{3}$, and the positions of extrema and zeros are also
similar. Correspondingly, for ${\cal M}$, $\eta$, and $t_{c}$, the zeros of
their biases basically coincide with the extrema of $c_{\alpha,\alpha+4}$.
When $|c|$ is large (such as the region around $\Delta t_{c}\sim 0$ in the
Asymmetric configuration), the oscillations of correlation coefficients and
biases are more severe, but their basic trends are still similar, namely, $c$
increases (decreases) when $\bm{B}$ increases (decreases).

In this work, we only discuss the dependence of the relative size of
correlation coefficients and biases on $\Delta t_{c}$. In fact, it is
incorrect to expect that there can not be any biased OS when $|c|$ is small.
For example, in the Symmetric and not close configuration, there is
$|c|\lesssim 0.4$, but the maximum of $B_{\alpha}$ is about $1.5$. For the
Asymmetric configuration, the maximum value of $|c|$ is close to $1$, but the
maximum of $B_{\alpha}$ is only slightly greater than $2$. After rescaling the
biases by the SNRs of the two configurations, their maxima of biases have only
small difference. On the other hand, when $|c|$ is large, the higher-order
term of $M$ in Eq. (20) cannot be ignored, and the relationship between
correlation coefficients and biases will become more complicated.

The correlation between any two parameters can also be analyzed using Eq.
(22). For example, if one wants to study the correlation between $d_{L}^{(1)}$
and $t_{c}^{(2)}$, taking $\beta=3$ and $\gamma=4$, only the
$A_{,\beta}^{(1)}A^{(2)}\phi_{,\gamma}^{(2)}$ term is nonzero. Thus, Eq. (22)
becomes the modulation integral of $\sin\Delta\phi$, so the correlation
coefficient will behave similarly to $B_{4}$. In addation, through simple
calculation we have
$c\propto\big{(}d_{L}^{(1)}\big{)}^{0}\big{(}d_{L}^{(2)}\big{)}^{0}$. This
means that the correlation coefficients do not depend on the strength of the
two signals, which is consistent with our intuitive understanding.

###  4.2 Full bayesian inference

Figure 7: The SPE results in the Equal configuration, rescaled by the
statistical uncertainty $\Delta\theta_{\rm stat}$ from FM. The maximum
likelihood points in the full Bayesian analysis and biases forecasted by the
FM are labeled “$+$” and “$\times$” respectively. For reference, the last
column shows the SPE results when the data stream contains only the louder
signal.

So far, we have used the FM method to discuss the biases in SPE. This requires
that the SNR of the target signal $h^{(1)}$ is large enough,333When
calculating the correlation coefficients with the FM method, the SNRs of both
$h^{(1)}$ and $h^{(2)}$ need to be large. and the biases are small enough to
make linearized-signal approximation valid [18, 48]. Although the results
obtained by FM in the previous sections are consistent with our physical
intuition, in principle, it is always necessary to verify the applicability of
the FM approximation. In this section, we select some representative points
for the $6$ group parameter configurations in Sec. 3.4, such as the extrema
and zeros of $B_{\alpha}$. We use Eq. (2.2) to perform full Bayesian PEs and
compare them with the results in Sec. 3.

We use open source software package Parallel Bilby, UltraNest [7, 47, 39, 11]
and the nested sampling Monte Carlo algorithm MLFriends [9, 10] to conduct
SPEs on $\bm{\theta}=\big{\\{}{\cal
M}^{(1)},\eta^{(1)},d_{L}^{(1)},t_{c}^{(1)}\big{\\}}$. For ${\cal M}$, $\eta$,
and $t_{c}$, we choose uniform priors. As for $d_{L}$, we first generate a
uniform prior for the comoving volume $V_{c}$ in the source frame, and then
convert it in the parameter space of $d_{L}$. Different prior choices may
affect PE results, but the difference between them can be ignored when
$\rho^{(1)}$ is sufficiently large [48].

When injecting mock GW data, $g(t)$, different noise realizations also affect
the PE results. Here we take $n=0$, which is the most convenient case for
generating data. More importantly, as can be seen from Eq. (9), the deviation
of the MLE result from the true value in this case is precisely the bias,
$\Sigma^{\alpha\beta}\big{(}h^{(1)}_{\beta}\,,h^{(2)}\big{)}$. It is not only
reasonable, but also necessary to use the SPE results when $n=0$ to verify the
validity of FM [33].

Figure 8: The contour plot of the parameters when $\Delta t_{c}=10.1\,{\rm
ms}$ in the Equal configuration. The maximum likelihood points in the full
Bayesian analysis and biases forecasted by the FM are labeled “$+$” and
“$\times$” respectively.

In Fig. 7, we show the marginalized distribution of each parameter in the
Equal configuration, with $\Delta t_{c}$ taken as those values corresponding
to the magenta reference lines in Fig. 2. The biases calculated using FM in
Sec. 3 are marked in the figure. The results of the other $5$ groups of
configurations are shown in D. For PE results with small biases, the
difference between MLE and the truth is almost the same as the biases
forecasted by the FM approximation. Furthermore, the contour plot of the
parameters when $\Delta t_{c}=10.1\,{\rm ms}$ in the Equal configuration is
shown in Fig. 8. In this case, all parameters have a large bias
($|B_{\alpha}|>1$), but the forecasted biases of FM still highly coincide with
the MLE of SPE.

##  5 Conclusion and Summary

With the continuous upgrading of GW detectors, more and more GW events will be
detected. For next generation ground-based detectors, a large number of GW
signals will overlap with each other [35, 40, 34, 38, 22]. In order to extract
physical information from these OSs more efficiently and accurately, it is
necessary to study the PE technology of OSs. However, before performing
complex PE on the whole OS, it is important to understand how the PE of a
specified signal in the OS is affected by other signals. This effect is
reflected in the bias between the SPE result and the true value of the
specified signal. If the overlapping effects on each SPE are small, then the
signals in the OS can be analyzed separately, and there is no need for the
time-consuming JPE. So far, many works have discussed the SPE biases of OSs,
covering various combinations of detectors, waveforms, and GW sources [35, 40,
34, 38, 22, 6, 27, 46, 23]. However, calculating lots of SPEs will cost
considerable computing resources, and it is difficult to study the systematic
dependence of bias on various parameters, resulting in our limited
understanding of their behaviors.

In this work, we use the FM method to investigate the origin of biases in SPE
and their dependence on GW parameters. Thanks to the fast computational speed
of the FM method [18, 48, 49], we can explore the parameter space in more
detail, and gain a deeper understanding of the biases. In Refs. [34, 6], the
FM method was considered, but not treated as the primary tool to analyze
biases. Here, with the help of FM, we establish an analytical expression of
SPE biases in Sec. 2. We mainly studied the dependence of bias on the mass,
luminosity distance, coalescence time, and phase, where OSs are generated by
two BBH systems in AdvLIGO. Since the direct cause of OS is the overlapping of
two signals in the time domain, the relationship between $\Delta t_{c}$ and
biases is particularly discussed in detail.

In Sec. 3, we find that there are large biases in the PE when the frequency
evolutions of the two signals overlap to a large extent. In the frequency
domain, this is equivalent to a constant phase difference between the two
signals in the corresponding frequency band. This conclusion can also be
intuitively obtained by considering the properties of matched filtering
technology, but in this work, we give a rigorous and quantitative
demonstration. For the whole parameter space, the existence of large biases is
common, but biases are more likely to occur when the masses of the two
signals, especially the chirp masses, are close. For most OSs, the biases
strongly depend on the merger order of the two BBHs. When the binary with
heavier component masses merges first, its frequency evolution is more likely
to overlap with the other one, leading to large biases. The biases always
reach their global extrema near $\Delta t_{c}\simeq 0$, and then oscillate and
tend to zero with the increase of $|\Delta t_{c}|$. For a specific parameter,
the behavior depends on how this parameter contributes to the waveform. The
zero points of the bias of the parameter that mainly contributes to the
frequency-domain phase are often the extreme points of the bias of the
parameter that mainly contributes to the amplitude.

We also used other methods to analyze the bias behaviors. In Sec. 4, we
calculate the correlation coefficient of the corresponding parameters of the
two signals, and compare their behaviors with biases. Using the explicit
expressions in Eqs. (3.2) and (22), we find that the behavior of correlation
coefficients is similar to the bias of the parameter that mainly contributes
to the amplitude. The numerical results verified this result (see Fig. 6). We
also conducted full Bayesian inference on some representative parameter
configurations to check the validity of the FM approximation. In conclusion,
using FM to forecast the biases beforehand provides a powerful guidance for a
full JPE or SPEs.

The main purpose of this work is to explain the origin of biases. For the
conciseness and universality of our conclusion, some secondary effects are
ignored. In the future, the results of this work should be further extended
before practical application. For parameters, we only selected the most
critical ones including the binary masses, luminosity distance, merger time,
and merger phase. When considering more parameters, such as the spins and the
sky location, the PE error will likely increase, and the absolute size of the
bias will also be affected. For detectors, we only considered OSs in AdvLIGO
in this work and expect the bias behaviors to be similar to those in XG
ground-based detectors, which is implied by recent works [6, 17, 22]. For long
signals, one needs to consider the rotation of the Earth. For space-based
detectors, the duration of each signal will be even longer, and biases in
these detectors may have some special properties. Besides, the utilization of
detector networks can also have a great impact on biases, depending on the
merger times of two signals in different detectors [38]. For waveforms, the
influence of waveform selection on bias can be studied. In addition, PE biases
can come from other effects, such as inaccurate waveform modeling [23] and
existence of glitches [19]. For the correlation between two signals, we
pointed out the similarity between correlation coefficients and biases more
quantitatively, extending the studies by [34] and [6]. Whether there is a
numerical relationship between them remains to be discussed. Finally, the
conclusions of this work can be probably extended to cases where more GW
signals overlap. In the FM approximation, the bias caused by different GW
signals can be linearly superimposed, so one will need to consider the biases
led by every other signal one by one.

It should be pointed out that, this work focuses on understanding the PE
behaviors from OSs, rather than seeking specific algorithms for analyzing OSs.
In real GW detection, we do not know the true values of each signal, and even
not know how many individual signals are present in the OS. The simple FM
method can provide some insights for PE, but cannot replace a focused
analysis. On the contrary, we expect that, the understanding of the mechanism
of biases in OS will help develop a faster and more efficient PE algorithm,
and serve in analyzing real OS events.

## Acknowledgements

This work was supported by the National Natural Science Foundation of China
(11975027, 12147177, 11991053, 11721303), the China Postdoctoral Science
Foundation (2021TQ0018), the National SKA Program of China (2020SKA0120300),
the Beijing Municipal Natural Science Foundation (1242018), the Max Planck
Partner Group Program funded by the Max Planck Society, and the High-
Performance Computing Platform of Peking University. J. Z. is supported by the
“LiYun” Postdoctoral Fellowship of Beijing Normal University, and C. L. is
supported by the China Scholarship Council (CSC).

##  Appendix A Waveforms and derivatives of waveforms

Figure 9: The frequency-domain phase, $\phi(f)$, and the normalized modulation
amplitudes, ${\cal A}^{c}_{\alpha}$ and ${\cal A}^{s}_{\alpha}$, of the
waveform. The number after each parameter is the maximum value of the
corresponding normalized modulation amplitude.

Figure 9 shows the frequency-domain phase of the IMRPhenomD waveform template,
as well as the derivatives of the waveform (expressed by ${\cal
A}^{c}_{\alpha}$ and ${\cal A}^{s}_{\alpha}$). Note that ${\cal A}^{s}_{3}=0$
and ${\cal A}^{c}_{4}=0$ are not drawn. The injected parameters are
$\big{\\{}m_{1}=30\,{\rm M_{\odot}},m_{2}=20\,{\rm M_{\odot}},d_{L}=500\,{\rm
Mpc},t_{c}=0\big{\\}}$. The frequency-domain phase grows rapidly in the low
frequency band ($\lesssim 25\,{\rm Hz}$), while in the high frequency band
($\gtrsim 100\,{\rm Hz}$) it grows slowly and approximately linearly. The
modulation amplitude typically reaches the maximum value within $30$–$50\,{\rm
Hz}$, consistent with the most sensitive frequency band of AdvLIGO. Also, for
${\cal M}$ and $\eta$, ${\cal A}^{s}_{\alpha}$ is at least two orders of
magnitude larger than ${\cal A}^{c}_{\alpha}$.

##  Appendix B Max-max and max-mean of reduced biases

Figure 10: Same as Fig. 1, but the parameters are taken from real GW events.

FIG. B1 (continued).

Figure 10 shows the dependence of the max-max value, max-mean average value,
and average asymmetry of $\bm{B}$ on $\big{(}m_{1}^{(2)},m_{2}^{(2)}\big{)}$.
The mass components $\big{(}m_{1}^{(1)},m_{2}^{(1)}\big{)}$ are taken as the
masses in the detector reference system of GW events GW150914, GW190602
175927, GW190814, and GW190924 021846\. Their behaviors are similar to those
in Fig. 1. When ${\cal M}^{(2)}\approx{\cal M}^{(1)}$, large biases occur, and
the average asymmetry is close to $0$.

##  Appendix C SPE biases behaviours in other internal parameter
configurations

Figure 11: Same as Fig. 2, but for the Asymmetric2, Random, Symmetric and
close, and Symmetric and not close configurations.

In the main text, we demonstrate that the biases come from the overlapping of
frequency evolution. After understanding the origin of biases, it is easy to
understand the global dependence of biases on intrinsic parameters in Sec. 3.3
as follows.

  1. (i)

For a given detector, its sensitive frequency band is fixed. Due to the
selection effect of the detector, the characteristic frequencies of the
detected signals are similar ($\sim 100\,\rm{Hz}$ for AdvLIGO). As long as the
merger time difference is appropriate, it is easy to make the corresponding
frequency evolutions overlap in the time domain. This explains why biased OSs
are inevitable.

  2. (ii)

For two GW signals with similar mass components, the frequency evolutions are
also similar. Therefore, they can easily overlap for a wide range of $\Delta
t_{c}$, and the stable stage of $\Delta\phi$ is longer, which is more likely
to produce large biases on average.

  3. (iii)

The frequency evolution $f(t)$ is strongly dependent on the masses of the GW
source. As long as the masses of two GW signals are not particularly close,
their evolution behaviors will be significantly different. Obviously, when the
signal with a lower frequency merges first, it is easier to overlap with the
frequency evolution of the signal with a higher frequency, which produces a
strong dependence of the biases on the merger order. On the contrary, if the
masses of the two signals are close, it does not matter which one merges
first, so $|\overline{\rm Asy}|$ is small. At the same time, according to
(ii), the biases in these cases are large in the average sense. This explains
the consistency between the regions where $|\overline{\rm Asy}|\sim 0$ and
those where the bias is large.

Without loss of generality, we selected three representative configurations in
the $\big{(}m_{1}^{(2)},m_{2}^{(2)}\big{)}$ parameter space, namely
Asymmetric2, Random, Symmetric and close, Symmetric and not close to further
validate these conclusions.

For the Asymmetric2 configuration, it is equivalent to exchanging the
positions of $h^{(1)}$and $h^{(2)}$, and $\Delta\phi$ is monotonically
decreasing with respect to $f$. Therefore, biased OSs will mainly occur when
$\Delta t_{c}<0$, which is consistent with the results in Fig. 11. For the
Random configuration, the source masses of $h^{(2)}$ are significantly larger
than that of $h^{(1)}$, so we expect that large biases mainly occur when
$\Delta t_{c}<0$, and the dependence of the bias on $\Delta t_{c}$ should be
similar to the Asymmetric2 configuration. The results in Fig. 2 verified this.
Further, for most $\big{(}m_{1}^{(2)},m_{2}^{(2)}\big{)}$ in the parameter
space, their $|\overline{\rm Asy}|$ is close to $1$, showing a strong
asymmetry of their bias behaviors. Therefore, Asymmetric and Asymmetric2
configurations can represent the dependence of biases on $\Delta t_{c}$ for
most OSs whose masses of two signals are not similar.

Next, we discuss the parameter configurations where $|\overline{\rm Asy}|\sim
0$. For the Symmetric and close configuration, we choose $m_{1}^{(2)}$ to be
slightly larger than $m_{1}^{(1)}$, and then adjust $m_{2}^{(2)}$ such that
$|\overline{\rm Asy}|$ is minimal. It can be seen that $\bm{B}$ is
approximately symmetrical about $\Delta t_{c}$. This is very similar to the
Equal configuration, except that all $B_{\alpha}(t)$ curves are slightly moved
to the left. Same as before, we take some $\Delta t_{c}$ to calculate the
frequency-domain phase, which is shown in Fig. 12. When $\Delta t_{c}=0$, we
find that the overall change of $\Delta\phi_{0}$ under this configuration is
not large, and it increases with $f$ first and then decreases. The change of
monotonicity leads to a long stable stage, i.e., a large bias. This is
significantly different from the Asymmetric configuration, where
$\Delta\phi_{0}$ increases monotonically with $f$, which requires additional
non-zero $\Delta t_{c}$ to generate a significantly stable stage.

Figure 12: Same as Fig. 3, but for the Symmetric and close configuration.

In the Symmetric and close configuration, $\Delta\phi_{0}$ does not change
much with $f$, which is why $\bm{B}$ behaves similarly to the Equal
configuration. When $\Delta t_{c}$ significantly deviates from $0$,
$|2\pi\Delta t_{c}f|$ will be much larger than $|\Delta\phi_{0}|$ for most
$f$, so $\Delta\phi=\Delta\phi_{0}-2\pi\Delta t_{c}f\approx-2\pi\Delta
t_{c}f$, which is exactly the phase difference in the Equal configuration. Of
course, on the basis of the relatively stable $\Delta\phi_{0}$, we can
slightly adjust $\Delta t_{c}$ to make $\Delta\phi$ more stable. When $\Delta
t_{c}=-0.7\,{\rm ms}$, $\Delta\phi$ is very close to $\pm\pi$ in a wide
frequency band, making the bias larger. In other words, when $\Delta
t_{c}=-0.7\,{\rm ms}$, the behavior of $\Delta\phi$ is more like a constant,
so the symmetry axis of $\bm{B}$ will move to $\Delta t_{c}=-0.7\,{\rm ms}$.
This explains the slight movement of the bias curves in the Symmetric and
close configuration.

Physically, the process of adjusting $m_{2}^{(2)}$ is equivalent to finding a
$h^{(2)}$ whose frequency evolution is most similar to $h^{(1)}$. As long as
the frequency evolution behavior of the two signals is similar, there will be
a long overlapping of the frequency evolutions no matter which binary merges
first, and the corresponding $|\overline{\rm Asy}|$ is small. Therefore, if
$m_{1}^{(2)}>m_{1}^{(1)}$, to minimize $|\overline{\rm Asy}|$ one must have
$m_{2}^{(2)}<m_{2}^{(1)}$. Otherwise, $h^{(2)}$ is generated by a heavier
binary, leading to similar results to that of the Asymmetric2 case.

Finally, we discuss the influence of $\eta$ on the bias or the bias behaviors
of binaries with extreme mass ratios. In general, the frequency evolutions are
mainly determined by the chirp mass $\cal M$. This is why we take
$\eta^{(2)}=\eta^{(1)}$ in the Asymmetric and Asymmetric2 configurations. In
the Symmetric and not close configuration, we fix $m_{1}^{(2)}=2m_{1}^{(1)}$
and adjust $m_{2}^{(2)}$ to minimize $|\overline{\rm Asy}|$. Obviously, to
obtain a similar frequency evolution, $m_{2}^{(2)}$ must be small.
Interestingly, at this time, although $M^{(2)}>M^{(1)}$, there is ${\cal
M}^{(2)}<{\cal M}^{(1)}$. This phenomenon was mentioned in Sec. 3.3, that is,
when $\big{(}m_{1}^{(1)},m_{2}^{(1)}\big{)}=(30,20)$, the region where
$|\overline{\rm Asy}|$ is small locates below the line where ${\cal
M}^{(2)}={\cal M}^{(1)}$. In Fig. 13, we draw $\Delta\phi(f)$, $D_{\alpha}$,
$R_{\alpha}$, whitened waveforms and frequency evolutions of the two signals
when $\Delta t_{c}=-9.8\,{\rm ms}$. Although ${\cal M}^{(2)}<{\cal M}^{(1)}$,
$m_{1}^{(2)}$ (also $M^{(2)}$) is too large. In general, the frequency of
$h^{(2)}$ is still lower than $h^{(1)}$. In other words, the frequency of a GW
signal with an extreme mass ratio is lower, so smaller chirp masses will be
needed to generate a similar frequency evolution to $h^{(1)}$. Therefore, in
Fig. 2, the region where $|\overline{\rm Asy}|\sim 0$ locates below the line
where ${\cal M}^{(2)}={\cal M}^{(1)}$. On the contrary, if $\eta$ is small,
such as GW190814, the region will be above (see Fig. 10).

Figure 13: Same as Fig. 5, but for the Symmetric and not close configuration
and $\Delta t_{c}=-9.8\,{\rm ms}$.

Back to the Symmetric and not close configuration, in fact, the biases are not
asymmetric about $\Delta t_{c}$, but just distributes on $\Delta t_{c}<0$ and
$\Delta t_{c}>0$ evenly. In this configuration, although $\Delta\phi_{0}$ also
has a stable stage similar to that in the Symmetric and close configuration,
the source masses and the frequency evolutions of the two signals, are too
different. So one cannot expect $\bm{B}$ to have an obvious symmetric axis
about $\Delta t_{c}$. However, in this case, we can still choose an
appropriate $\Delta t_{c}$ to increase the length of the stable stage and
cause a large bias. The use of $\Delta t_{c}=-9.8\,{\rm ms}$ is an example,
and its behavior will not be described repeatedly here.

##  Appendix D Full Bayesian SPE results in other internal parameter
configurations

Figure 14 shows the marginalized distribution of parameters in the SPE when
the internal parameter configurations are taken as Asymmetric, Asymmetric2,
Random, Symmetric and close, and Symmetric and not close. Same as the Equal
configuration, we find that the results forecasted by FM and the actual MLEs
are in good agreement.

Figure 14: Same as Fig. 7, but for other parameter configurations. The
configuration names are marked on the top.

FIG D1 (continued).

## References

  * [1] B.. Abbott  “GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs”  In _Phys. Rev. X_ 9.3, 2019, pp. 031040  DOI: 10.1103/PhysRevX.9.031040
  * [2] B.. Abbott  “Observation of Gravitational Waves from a Binary Black Hole Merger”  In _Phys. Rev. Lett._ 116.6, 2016, pp. 061102  DOI: 10.1103/PhysRevLett.116.061102
  * [3] R. Abbott  “GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run”  In _Phys. Rev. X_ 11, 2021, pp. 021053  DOI: 10.1103/PhysRevX.11.021053
  * [4] R. Abbott  “GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo during the Second Part of the Third Observing Run”  In _Phys. Rev. X_ 13.4, 2023, pp. 041039  DOI: 10.1103/PhysRevX.13.041039
  * [5] Pau Amaro-Seoane  “Laser Interferometer Space Antenna”  In _e-prints_ , 2017  DOI: 10.48550/arXiv.1702.00786
  * [6] Andrea Antonelli, Ollie Burke and Jonathan R. Gair  “Noisy neighbours: inference biases from overlapping gravitational-wave signals”  In _Mon. Not. Roy. Astron. Soc._ 507.4, 2021, pp. 5069–5086  DOI: 10.1093/mnras/stab2358
  * [7] Gregory Ashton  “BILBY: A user-friendly Bayesian inference library for gravitational-wave astronomy”  In _Astrophys. J. Suppl._ 241.2, 2019, pp. 27  DOI: 10.3847/1538-4365/ab06fc
  * [8] Emanuele Berti, Alessandra Buonanno and Clifford M. Will  “Estimating spinning binary parameters and testing alternative theories of gravity with LISA”  In _Phys. Rev. D_ 71, 2005, pp. 084025  DOI: 10.1103/PhysRevD.71.084025
  * [9] Johannes Buchner  “A statistical test for Nested Sampling algorithms”  In _Statistics and Computing_ 26.1-2, 2016, pp. 383–392  DOI: 10.1007/s11222-014-9512-y
  * [10] Johannes Buchner  “Collaborative Nested Sampling: Big Data versus Complex Physical Models”  In _Publ. Astron. Soc. Pac._ 131.1004, 2019, pp. 108005  DOI: 10.1088/1538-3873/aae7fc
  * [11] Johannes Buchner  “Nested sampling methods”  In _Statist. Surv._ 17, 2023, pp. 169  DOI: 10.1214/23-ss144
  * [12] Alessandra Buonanno et al.  “Comparison of post-Newtonian templates for compact binary inspiral signals in gravitational-wave detectors”  In _Phys. Rev. D_ 80, 2009, pp. 084043  DOI: 10.1103/PhysRevD.80.084043
  * [13] George Casella and Roger L Berger  “Statistical inference”  Cengage Learning, 2021
  * [14] Nelson Christensen and Renate Meyer  “Markov chain Monte Carlo methods for Bayesian gravitational radiation data analysis”  In _Phys. Rev. D_ 58, 1998, pp. 082001  DOI: 10.1103/PhysRevD.58.082001
  * [15] Nelson Christensen, Rejean J. Dupuis, Graham Woan and Renate Meyer  “A Metropolis-Hastings algorithm for extracting periodic gravitational wave signals from laser interferometric detector data”  In _Phys. Rev. D_ 70, 2004, pp. 022001  DOI: 10.1103/PhysRevD.70.022001
  * [16] Curt Cutler  “Angular resolution of the LISA gravitational wave detector”  In _Phys. Rev. D_ 57, 1998, pp. 7089–7102  DOI: 10.1103/PhysRevD.57.7089
  * [17] Yixuan Dang, Ziming Wang, Dicong Liang and Lijing Shao  “Impact of overlapping signals on parameterized post-Newtonian coefficients in tests of gravity”  In _arXiv e-prints_ , 2023, pp. arXiv:2311.16184  DOI: 10.48550/arXiv.2311.16184
  * [18] Lee S. Finn  “Detection, measurement and gravitational radiation”  In _Phys. Rev. D_ 46, 1992, pp. 5236–5249  DOI: 10.1103/PhysRevD.46.5236
  * [19] Sudarshan Ghonge et al.  “Characterizing the efficacy of methods to subtract terrestrial transient noise near gravitational wave events and the effects on parameter estimation”  In _arXiv e-prints_ , 2023, pp. arXiv:2311.09159  DOI: 10.48550/arXiv.2311.09159
  * [20] Yungui Gong, Jun Luo and Bin Wang  “Concepts and status of Chinese space gravitational wave detection projects”  In _Nature Astron._ 5.9, 2021, pp. 881–889  DOI: 10.1038/s41550-021-01480-3
  * [21] S. Hild  “Sensitivity Studies for Third-Generation Gravitational Wave Observatories”  In _Class. Quant. Grav._ 28, 2011, pp. 094013  DOI: 10.1088/0264-9381/28/9/094013
  * [22] Yoshiaki Himemoto, Atsushi Nishizawa and Atsushi Taruya  “Impacts of overlapping gravitational-wave signals on the parameter estimation: Toward the search for cosmological backgrounds”  In _Phys. Rev. D_ 104.4, 2021, pp. 044010  DOI: 10.1103/PhysRevD.104.044010
  * [23] Qian Hu and John Veitch  “Accumulating Errors in Tests of General Relativity with Gravitational Waves: Overlapping Signals and Inaccurate Waveforms”  In _Astrophys. J._ 945.2, 2023, pp. 103  DOI: 10.3847/1538-4357/acbc18
  * [24] Wen-Rui Hu and Yue-Liang Wu  “The Taiji Program in Space for gravitational wave physics and the nature of gravity”  In _Natl. Sci. Rev._ 4.5, 2017, pp. 685–686  DOI: 10.1093/nsr/nwx116
  * [25] Sascha Husa et al.  “Frequency-domain gravitational waves from nonprecessing black-hole binaries. I. New numerical waveforms and anatomy of the signal”  In _Phys. Rev. D_ 93.4, 2016, pp. 044006  DOI: 10.1103/PhysRevD.93.044006
  * [26] Soichiro Isoyama, Hiroyuki Nakano and Takashi Nakamura  “Multiband Gravitational-Wave Astronomy: Observing binary inspirals with a decihertz detector, B-DECIGO”  In _Prog. Theor. Exp. Phys._ 2018.7, 2018, pp. 073E01  DOI: 10.1093/ptep/pty078
  * [27] Justin Janquart et al.  “Parameter estimation methods for analyzing overlapping gravitational wave signals in the third-generation detector era”  In _e-prints_ , 2022  DOI: 10.48550/arXiv.2211.01304
  * [28] Vicky Kalogera  “The Next Generation Global Gravitational Wave Observatory: The Science Book”  In _e-prints_ , 2021  DOI: 10.48550/arXiv.2111.06990
  * [29] Steven M Kay  “Statistical signal processing: estimation theory”, 1993
  * [30] Sebastian Khan et al.  “Frequency-domain gravitational waves from nonprecessing black-hole binaries. II. A phenomenological model for the advanced detector era”  In _Phys. Rev. D_ 93.4, 2016, pp. 044007  DOI: 10.1103/PhysRevD.93.044007
  * [31] Chang Liu and Lijing Shao  “Neutron Star–Neutron Star and Neutron Star–Black Hole Mergers: Multiband Observations and Early Warnings”  In _Astrophys. J._ 926.2, 2022, pp. 158  DOI: 10.3847/1538-4357/ac3cbf
  * [32] Jun Luo  “TianQin: a space-borne gravitational wave detector”  In _Class. Quant. Grav._ 33.3, 2016, pp. 035010  DOI: 10.1088/0264-9381/33/3/035010
  * [33] Samaya Nissanke et al.  “Exploring short gamma-ray bursts as gravitational-wave standard sirens”  In _Astrophys. J._ 725, 2010, pp. 496–514  DOI: 10.1088/0004-637X/725/1/496
  * [34] Elia Pizzati, Surabhi Sachdev, Anuradha Gupta and Bangalore Sathyaprakash  “Toward inference of overlapping gravitational-wave signals”  In _Phys. Rev. D_ 105.10, 2022, pp. 104016  DOI: 10.1103/PhysRevD.105.104016
  * [35] Tania Regimbau and Scott A. Hughes  “Gravitational-wave confusion background from cosmological compact binaries: Implications for future terrestrial detectors”  In _Phys. Rev. D_ 79, 2009, pp. 062002  DOI: 10.1103/PhysRevD.79.062002
  * [36] David Reitze  “Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO”  In _Bull. Am. Astron. Soc._ 51.7, 2019, pp. 035  DOI: 10.48550/arXiv.1907.04833
  * [37] David Reitze  “The US Program in Ground-Based Gravitational Wave Science: Contribution from the LIGO Laboratory”  In _Bull. Am. Astron. Soc._ 51, 2019, pp. 141  DOI: 10.48550/arXiv.1903.04615
  * [38] Philip Relton and Vivien Raymond  “Parameter estimation bias from overlapping binary black hole events in second generation interferometers”  In _Phys. Rev. D_ 104.8, 2021, pp. 084039  DOI: 10.1103/PhysRevD.104.084039
  * [39] I.. Romero-Shaw  “Bayesian inference for compact binary coalescences with bilby: validation and application to the first LIGO–Virgo gravitational-wave transient catalogue”  In _Mon. Not. Roy. Astron. Soc._ 499.3, 2020, pp. 3295–3319  DOI: 10.1093/mnras/staa2850
  * [40] Anuradha Samajdar, Justin Janquart, Chris Van Den Broeck and Tim Dietrich  “Biases in parameter estimation from overlapping gravitational-wave signals in the third-generation detector era”  In _Phys. Rev. D_ 104.4, 2021, pp. 044003  DOI: 10.1103/PhysRevD.104.044003
  * [41] B. Sathyaprakash  “Scientific Objectives of Einstein Telescope” [Erratum: Class.Quant.Grav. 30, 079501 (2013)]  In _Class. Quant. Grav._ 29, 2012, pp. 124013  DOI: 10.1088/0264-9381/29/12/124013
  * [42] Sanjib Sharma  “Markov Chain Monte Carlo Methods for Bayesian Data Analysis in Astronomy”  In _Ann. Rev. Astron. Astrophys._ 55, 2017, pp. 213–259  DOI: 10.1146/annurev-astro-082214-122339
  * [43] Kevin J. Shuman and Neil J. Cornish  “Massive black hole binaries and where to find them with dual detector networks”  In _Phys. Rev. D_ 105.6, 2022, pp. 064055  DOI: 10.1103/PhysRevD.105.064055
  * [44] John Skilling  “Nested Sampling”  In _Bayesian Inference and Maximum Entropy Methods in Science and Engineering_ 735, American Institute of Physics Conference Series, 2004, pp. 395  DOI: 10.1063/1.1835238
  * [45] John Skilling  “Nested sampling for general Bayesian computation”  In _Bayesian Analysis_ 1.4, 2006, pp. 833–859  DOI: 10.1214/06-BA127
  * [46] Rory Smith  “Bayesian Inference for Gravitational Waves from Binary Neutron Star Mergers in Third Generation Observatories”  In _Phys. Rev. Lett._ 127.8, 2021, pp. 081102  DOI: 10.1103/PhysRevLett.127.081102
  * [47] Rory J.. Smith, Gregory Ashton, Avi Vajpeyi and Colm Talbot  “Massively parallel Bayesian inference for transient gravitational-wave astronomy”  In _Mon. Not. Roy. Astron. Soc._ 498.3, 2020, pp. 4492–4502  DOI: 10.1093/mnras/staa2483
  * [48] Michele Vallisneri  “Use and abuse of the Fisher information matrix in the assessment of gravitational-wave parameter-estimation prospects”  In _Phys. Rev. D_ 77, 2008, pp. 042001  DOI: 10.1103/PhysRevD.77.042001
  * [49] Ziming Wang, Chang Liu, Junjie Zhao and Lijing Shao  “Extending the Fisher Information Matrix in Gravitational-wave Data Analysis”  In _Astrophys. J._ 932, 2022, pp. 102  DOI: 10.3847/1538-4357/ac6b99
  * [50] Ziming Wang et al.  “Simultaneous bounds on the gravitational dipole radiation and varying gravitational constant from compact binary inspirals”  In _Phys. Lett. B_ 834, 2022, pp. 137416  DOI: 10.1016/j.physletb.2022.137416
  * [51] Junjie Zhao et al.  “Probing dipole radiation from binary neutron stars with ground-based laser-interferometer and atom-interferometer gravitational-wave observatories”  In _Phys. Rev. D_ 104.8, 2021, pp. 084008  DOI: 10.1103/PhysRevD.104.084008