# Mitigating the effect of population model uncertainty on strong lensing
Bayes factor using nonparametric methods

Damon H. T. Cheung1,2, Stefano Rinaldi3,4, Martina Toscani5 and Otto A.
Hannuksela1
1Department of Physics, The Chinese University of Hong Kong, Shatin, N.T.,
Hong Kong
2Physics Department, University of Michigan, Ann Arbor, MI 48109, USA
3Dipartimento di Fisica “E. Fermi”, Università di Pisa, I-56127 Pisa, Italy
4INFN, Sezione di Pisa, I-56127 Pisa, Italy
5Laboratoire des 2 Infinis - Toulouse (L2IT-IN2P3), Université de Toulouse,
CNRS, UPS, F-31062 Toulouse Cedex 9, France E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Strong lensing of gravitational waves can produce several detectable images as
repeated events in the upcoming observing runs, which can be detected with the
posterior overlap analysis (Bayes factor). The choice of the binary black hole
population plays an important role in the analysis as two gravitational-wave
events could be similar either because of lensing or astrophysical
coincidence. In this study, we investigate the biases induced by different
population models on the Bayes factor. We build up a mock catalogue of
gravitational-wave events following a benchmark population and reconstruct it
using both non-parametric and parametric methods. Using these reconstructions,
we compute the Bayes factor for lensed pair events by utilizing both models
and compare the results with a benchmark model. We show that the use of a non-
parametric population model gives a smaller bias than parametric population
models. Therefore, our study demonstrates the importance of choosing a
sufficiently agnostic population model for strong lensing analyses.

###### keywords:

methods: data analysis – methods: statistical – gravitational waves – stars:
black holes

††pubyear: 2023††pagerange: Mitigating the effect of population model
uncertainty on strong lensing Bayes factor using nonparametric
methods–References

##  1 Introduction

General relativity posits that massive objects bend spacetime in a way that
focuses the light rays that travel near them and causes gravitational lensing
(Schneider et al., 1992; Bartelmann, 2010). Through decades of technological
innovation, gravitational lensing has grown from an abstract theoretical idea
into a standard tool in the study of the cosmos. For example, it has
historically provided among the most stringent evidence for the existence of
dark matter (Clowe et al., 2004; Markevitch et al., 2004), aided in exoplanet
discoveries (Bond et al., 2004), and allowed scientists to measure the
expansion of the Universe (Treu & Marshall, 2016).

Similar to electromagnetic waves, gravitational waves (GWs) can also be
gravitationally lensed if they travel near massive astrophysical objects
(Ohanian, 1974; Wang et al., 1996; Nakamura & Deguchi, 1999; Deguchi & Watson,
1986). Such lensing of GWs could produce a plethora of observational
signatures. In the so-called geometrical optics limit,111The geometrical optic
limit is the situation in which the GW wavelength is significantly smaller
than the Schwarzschild radius of the lens. lensing of GWs is similar to
traditional lensing of light (Takahashi & Nakamura, 2003). In this limit,
lensing can

  1. (i)

magnify the GW (Pang et al., 2020), producing an apparently high-mass
distribution of binary black holes (BBHs) (Dai et al., 2016; Oguri, 2018;
Broadhurst et al., 2018; Cusin et al., 2021) that might be detectable if the
source is a neutron star merger (Pang et al., 2020),

  2. (ii)

invert the GW along one or both principal axes (Dai & Venumadhav, 2017;
Ezquiaga et al., 2020) that may be detectable as a discernible complex phase
shift of waveform with higher order modes (Dai & Venumadhav, 2017; Ezquiaga et
al., 2020; Wang et al., 2021a; Janquart et al., 2021b; Vijaykumar et al.,
2022a),

  3. (iii)

split the signal into multiple images detectable as repeated events in the
detectors (Haris et al., 2018; Hannuksela et al., 2019; Li et al., 2023;
Broadhurst et al., 2019; McIsaac et al., 2019; Liu et al., 2020; Dai et al.,
2020; Abbott et al., 2023a; Kim et al., 2022).

In the wave optics limit,222The wave optics limit is the situation in which
the GW wavelength is comparable to or larger than the Schwarzschild radius of
the lens. the GW can experience frequency-dependent deformations that might be
captured with template-based or machine-learning searches (Cao et al., 2014;
Lai et al., 2018; Christian et al., 2018; Dai et al., 2018; Singh et al.,
2018; Kim et al., 2020; Cheung et al., 2020; Yeung et al., 2021; Bulashenko &
Ubach, 2022; Seo et al., 2021; Tambalo et al., 2022). Thus far, there have
been a number of searches for these GW lensing effects in LIGO-Virgo-Kagra
(LVK) data (Hannuksela et al., 2019; Li et al., 2023; McIsaac et al., 2019;
Liu et al., 2020; Dai et al., 2020; Abbott et al., 2023a; Kim et al., 2022;
Abbott et al., 2021, 2023b), but there has been no widely accepted detection,
despite a recent reinterpretation of the LVK lensing search results (Diego et
al., 2021).

If lensed GW events are observed, they could enable many scientific pursuits,
including the detection of intermediate-mass and primordial BHs (Lai et al.,
2018; Diego, 2020; Basak et al., 2022), MACHOs (Basak et al., 2022), and a
better comprehension of the nature of dark matter (Basak et al., 2022; Oguri &
Takahashi, 2020, 2022; Wang et al., 2021b; Cao et al., 2022a). The effectively
expanded detector network that can detect lensed events and has excellent
timing accuracy may also enable the search for the associations of fast radio
bursts and GWs (Singh et al., 2023), precise tests of the cosmology (Jana et
al., 2022), speed of gravity (Collett & Bacon, 2016; Fan et al., 2017; Baker &
Trodden, 2016; Ezquiaga et al., 2022), polarization (Goyal et al., 2020), and
propagation (Mukherjee et al., 2019a, b; Chung & Li, 2021; Finke et al., 2021;
Iacovelli et al., 2022). In the wave optics limit, they may break the mass-
sheet degeneracy (Cremonese et al., 2021) and allow to measure the velocity of
the lens (Itoh et al., 2009). When strongly lensed, they might allow for sub-
arcsecond precision localization of merging BHs (Hannuksela et al., 2020;
Wempe et al., 2022), precision tests of cosmology (Sereno et al., 2010; Liao
et al., 2017; Hannuksela et al., 2020; Cao et al., 2022b), and combined
studies of lensed binaries in the GW and electromagnetic band (Smith et al.,
2022).

Strong lensing has been of particular interest in recent years. It has been
forecasted to become detectable at design detector sensitivity, assuming that
the current BBH detections originate from stellar merger remnants (Ng et al.,
2018; Oguri, 2018; Xu et al., 2021; Wierda et al., 2021; Smith et al.,
2022)333Observational limits on the rate of lensing have also been set,
agnostically, based on the non-detection of the stochastic GW background
(Buscicchio et al., 2020b, a; Mukherjee et al., 2020; Abbott et al., 2023a)..
Consequently, there has been a major effort to develop the necessary tools to
detect and analyse strong lensing (Haris et al., 2018; Dai et al., 2020; Liu
et al., 2020; Cheung et al., 2020; Hannuksela et al., 2020; Seo et al., 2021;
Lo & Magaña Hernandez, 2023; Janquart et al., 2021a, b; Goyal et al., 2021;
Janquart et al., 2022; Goyal et al., 2023).

The current strong lensing analysis methods all search for repeated events:
two or more GW events that appear to have the same frequency evolution and
originate from the same sky location, in line with the strong lensing
expectation. These waves would have identical detector-frame parameters and
sky location but differ in arrival times and luminosity due to magnification
after strong lensing and overall phases as they pass through caustics
(Takahashi & Nakamura, 2003; Dai & Venumadhav, 2017). Thus, in practice, a
strong lensing analysis inspects how similar/dissimilar the detector-frame
parameters and sky location are between two or more GW events. Alternatively,
one could state that a strong lensing search attempts to identify two or more
copies of binary compact object systems that appear to be near-identical.
However, in more technical terms, the analyses evaluate the strong lensing
Bayes factor – the evidence ratio between lensed and unlensed hypothesis. The
Bayes factor quantifies the support for an event that is strongly lensed
versus only appearing to be strongly lensed due to astrophysical coincidence
(see Fig. 1), and is therefore of utmost importance for strong lensing
detections.

Indeed, the probability of encountering strong lensing mimickers increases
quadratically as the number of detections increases, a challenge that could be
mitigated by higher signal-to-noise ratios or the inclusion of lensing time-
delay estimates which however necessitate the use of accurate lens population
models.

Another challenge associated with strong lensing detections is their reliance
on population modelling. In particular, the likelihood that two or more near-
identical binary compact objects appear in the GW detector through
astrophysical coincidence depends on the population model: if the Universe
produces many similar compact binaries, then the chance of them being an
astrophysical coincidence is high, and vice-versa. There has been a study on
the degree to which unlensed events mimic lensed ones because of the overlap
of parameters due to random coincidence (Çalışkan et al., 2022). However, the
statistic used does not include population terms that account for the
probability of coincidence. To make strong lensing detections, one must
consider the probability of astrophysical coincidence in the lensing
statistic. Therefore, one must have an accurate BBH population model. This
representation can, in principle, be inferred from LVK data directly: however,
one must be careful to use a model that is able to capture all the features
that might be present in the population. Currently, phenomenological
parametric population fits, motivated by physical arguments but not directly
linked to physical processes, are used, which might in principle capture the
characteristics of the source population, if the model is sophisticated
enough. However, there has not yet been a significant investigation to
understand to what extent the use of these parametric models could bias strong
lensing detections.

Thus, here we investigate the following two questions:

  1. (a)

Could the utilisation of an incorrect parametric population model lead to an
inaccurate strong lensing analysis?

  2. (b)

If so, can we overcome this challenge using non-parametric population models?

In particular, the use of different GW population models in the analysis will
affect the prediction and can result in a wrong classification (false
detection of lensing).

In this study, we build up a mock catalogue of GW events (benchmark population
model) that mimic the preferred model from (Abbott et al., 2023c), the
PowerLaw+Peak model (PP model, thereinafter), and generate the posterior
samples of the observed events. On top of this catalogue, we simulate also a
set of lensed pairs of events. Then, we calculate the Bayes factor of the
lensed pairs based on four population models:

  1. (i)

Benchmark population model that mimics PP model;

  2. (ii)

A Hierarchy of Dirichlet Process Gaussian Mixture Models ((H)DPGMM) - Non-
parametric population model inferred from the mock catalogue;

  3. (iii)

Power-law model - Parametric population model inferred from the mock
catalogue;

  4. (iv)

Uniform model - events are uniformly distributed across the parameter space.

We will then show how the Bayes factor of lensed pairs changes with these
different population models.

This paper is structured as follows: in Section 2, we review the strong
lensing mechanism of GW events and posterior overlap analysis using Bayes
factor. In Section 3, we give the details of the population models we use in
our analysis, both in the simulation and in the recovery. Section 4
illustrates the performance of the various population models in lensing
detections. In Section 5, we discuss the results and conclude, arguing that
the use of non-parametric population models is suitable for strong lensing
detections, while the use of parametric population models may lead to biases.

Figure 1:  A strong lens can produce multiple images of a binary black hole
(BBH). Although we can not discern the images in the sky with a ground-based
gravitational-wave detector due to poor angular resolution, we can detect them
as ”repeated events” – two or more gravitational-wave copies that only differ
in their amplitudes, phases, and overall arrival times, but are otherwise
identical. Indeed, to detect strong lensing, we look for multiple copies of
the same event. However, due to finite detector accuracy, a non-lensed BBH may
also appear to be similar to other events by astrophysical coincidence that
mimics lensing, as indicated in the diagram. The probability of finding such
lensing mimickers depends on the underlying population of BBHs. Correctly
modelling the BBH population on lensed detections is therefore vital to
detections.

##  2 Gravitational-wave lensing

If a BBH is strongly lensed by a galaxy or a galaxy cluster, GWs would split
into multiple “images”, observable as repeated near-identical444Differing only
in amplitude, arrival time, and overall phase. copies of the same binary. If,
in addition, the GW encounters smaller lenses (such as stars, primordial BHs,
globular clusters, or dark matter subhalos) along its travel path, they can
induce frequency-dependent beating patterns and wave optics modulations (often
referred to as GW microlensing and millilensing). In the limit of weak
lensing, the GW would merely be magnified, and the resulting weakly lensed GW
would have a different amplitude. In this work, we focus on strong lensing of
binary stellar remnants by galaxies and galaxy clusters, and neglect the
effect of micro-, milli- and wave optics lensing.

In the strong lensing scenario, that is, when the binary is sufficiently
aligned with a galaxy or a galaxy cluster such that strong lensing occurs, the
lens produces several images of the BBH. More specifically, the strong lens
produces $2N+1$ images of the source (usually 2 or 4 visible images for a
galaxy lens); each image would be related to the others by some magnification
$\mu$, lensing time delay $t_{d}$, and Morse phase $n$, such that the $i$th GW
strain (Takahashi & Nakamura, 2003)

$h_{L}^{i}(f)=|\mu_{i}|^{1/2}e^{i(2\pi ft_{i}^{d}-\pi
n_{i}\text{sign}(f))}h(f;\theta)\,,$ (1)

where $h(f;\theta)$ is the GW strain in the absence of strong GW lensing,
$\theta:=\\{m_{1},m_{2},\vec{s}_{1},\vec{s}_{2},z_{s},\theta_{J_{N}},\alpha,\delta,\psi,t^{\rm
u}\\}$,555$m_{1,2}$ denote the component masses, $\vec{s}_{1,2}$ the BH spins,
$z_{s}$ the source redshift, $\theta_{J_{N}}$ the inclination between total
angular momentum $J$ and the direction of propagation, $\alpha$ and $\delta$
the right ascension and declination, $\psi$ the polarization angle, and
$t^{\rm u}$ the arrival time of the emitted GW at the detector are the BBH
source parameters. The Morse phase $n_{i}$ equals $0$, $1/2$, and $1$ for
minima, saddle, and maxima points of the lensing time-delays and can be solved
when the lens geometry is known.

Note that although it is common to use a strong lensing model to evaluate the
image properties $\nu_{i}$ (magnifications, time delays, and Morse phases), it
is typically more convenient for us to work with the image properties
directly. Even more convenient is working with so-called effective parameters.
In particular, absorbing the magnification into an ”effective” luminosity
distance of the $i$th image,

$d_{L}^{\rm eff,i}=d_{L}/|\mu_{i}|^{1/2}\,,$ (2)

and the lensing time-delay into the effective arrival time of the $i$th image

$t_{i}=t_{0}+t_{i}^{d}\,,$ (3)

is more convenient from an analysis perspective. By using effective
parameters, one can model the two lensing effects using a general-relativistic
waveform without modifying anything but the luminosity distance and arrival
time terms. I.e., we can absorb the lensing modification and time delay by
simply replacing all luminosity distance and arrival time variables with their
effective counterparts. Such a procedure is attractive because it does not
require changing existing analysis infrastructure.666Similar to many prior
works, here we ignore the Morse phase term, which could be absorbed in the
coalescence phase when the GW consists primarily of (2,2) modes; the inclusion
of the Morse phase would merely add a phase term into the analysis that might
allow for a slightly improved detection criterion (Dai & Venumadhav, 2017;
Ezquiaga et al., 2020; Vijaykumar et al., 2022b).

###  2.1 Lens modelling

With the form of the strongly lensed GW strain (Eq. 1) known, one can model
strongly lensed GWs, given that we know the image properties. However, it is
still important to ask the question of how one obtains the image properties in
the first place. These can be solved from the lens equation777We refer those
interested in relating the image parameters $\nu$ to the lens model to
standard lensing textbooks such as (Congdon & Keeton, 2018)., given a strong
lens model. It was realized early on that relatively simple _smooth_ lens mass
models888E.g., power-law ellipsoidal profiles with shear. fit the strong
galaxy lens flux observations on large scales quite well (Mao & Schneider,
1997; Metcalf & Madau, 2001), while modelling substructure within the lens and
along the line-of-sight can be captured relatively well with a Gaussian
uncertainty term to the magnification (e.g. Wempe et al., 2022, and references
therein). Such galaxy lens models could be used for data-analysis purposes.
However, for the purposes of mock data and benchmarking, simpler lens models
are often favored, as they can be used to produce data relatively quickly
without complex lens morphology and population modelling. Since we are
primarily interested in the source population modelling and not the lens
modelling, we opt for this latter approach.

In particular, we assume that the _a priori_ distribution of magnifications
for each GW image, $\mu$, follows a power-law distribution (Narayan &
Bartelmann, 1996)

$p(\mu)\propto\mu^{-3}\,,$ (4)

which is appears to be correct near the caustic (i.e., when the magnification
is high, $\mu\gtrsim 10$). Based on simulations and forecasts for the current
and near future detector network, most detected GW events are relatively
highly magnified (Dai et al., 2016; Ng et al., 2017; Xu et al., 2021; Wierda
et al., 2021; Smith et al., 2022), which justifies the high-magnification
limit.

While the magnification distribution of any given image follows the
distribution above, the magnifications between the images are correlated.
Specifically, in most strong lensing configurations where the signal is
heavily magnified, the source is close to a caustic, and due to near-symmetry
typical strong lenses, the magnifications of both images are comparable
($\mu_{1}\sim\mu_{2}$, More & More 2022)999For the singular isothermal sphere
(SIS) model, $\mu_{1}=\mu_{2}$ in the limit $\mu_{1}\rightarrow\infty$. A
similar result can be derived for the singular isothermal ellipsoid (SIE)
models.. However, the magnifications are rarely precisely the same, except in
highly simplified lens models. Assuming that the magnification ratios were
exactly known would produce highly optimistic results. Thus, to correct for
lens modeling uncertainties, including ellipticities, shears, and
substructure, one might assume, as a simple approximation, that the
distribution of the second magnification follows

$\displaystyle
p(\mu_{2}|\mu_{1})\propto\mathcal{N}(\mu_{1},0.15\times\mu_{1})\,,$ (5)

where $\mu_{1}$ and $\mu_{2}$ are the magnification for the first image and
second image respectively, $\mathcal{N}(\mu,\sigma)$ denotes a normal
distribution with mean $\mu$ and standard deviation $\sigma$. Such a spread in
the magnification induced by microlensing and the substructure is conservative
for high-signal events when neglecting the non-linear effect and is used in
many prior studies (Hannuksela et al., 2020; Wempe et al., 2022), and it can
generally be regarded as better than presuming that the magnifications were
fixed, for example, by a singular isothermal sphere lens model.

###  2.2 Strong lensing Bayes factor

The quantity we will use to decide whether a pair of events is likely to be
lensed or not is the strong lensing Bayes factor
$\mathcal{B}_{{}^{U}}^{{}^{L}}$. The strong lensing hypothesis states that the
set of events originates from the same astrophysical event. Therefore, each of
the detector-frame parameters will be the same under this hypothesis (Haris et
al., 2018). The unlensed hypothesis states that the set of signals was
produced by unrelated astrophysical events (in this case, compact binary
coalescences).

In principle, the strong lensing Bayes factor should include all detector-
frame binary parameters and lens parameters (Liu et al., 2020; Lo & Magaña
Hernandez, 2023; Janquart et al., 2021a). Evaluating the Bayes factor would
require constructing and estimating the joint likelihood between the set of
waveforms. However, such an evaluation is highly demanding from a
computational point of view and thus impossible for us to use at high
throughput when performing an analysis at the population level.

Therefore, we opt to use a well-developed and widely used technique, called
_posterior overlap analysis_ , to estimate a “subset” of the Bayes factor,
used for example in (Haris et al., 2018). In this study, we focus on lensing
that produces 2 images with the subset of the binary parameters
$\bm{\theta}^{\prime}$ are $m_{1}^{z},m_{2}^{z}$, which correspond to the
detector-frame primary mass and secondary mass are unaffected by the
gravitational lensing. We marginalise over $\bm{\theta}^{\prime}$ to get the
following strong lensing Bayes factor (see Appendix B for the full
derivation):

$\displaystyle\mathcal{B}_{{}^{U}}^{{}^{L}}$
$\displaystyle=\frac{p(\bm{d}|\mathcal{H}_{L})}{p(\bm{d}|\mathcal{H}_{U})}=\frac{p(d_{1},d_{2}|\mathcal{H}_{L})}{p(d_{1}|\mathcal{H}_{U})p(d_{2}|\mathcal{H}_{U})},$
(6)

where $p(d_{1},d_{2}|\mathcal{H}_{L})$ is the evidence of having data
$d_{1},d_{2}$ under lensed hypothesis and $p(d_{i}|\mathcal{H}_{U})$ is the
evidence of having data $d_{i}$ under unlensed hypothesis. The Bayes factor
(Eq. 6) can then be expanded in binary parameter space as

$\displaystyle\mathcal{B}_{{}^{U}}^{{}^{L}}$
$\displaystyle=\dfrac{\int\dfrac{p_{\mathrm{obs}}(\bm{\theta}^{\prime})p(\bm{\theta}^{\prime}|d_{1})p(\bm{\theta}^{\prime}|d_{2})}{p_{\mathrm{PE}}^{2}(\bm{\theta}^{\prime})}d\bm{\theta}^{\prime}}{\int\dfrac{p_{\mathrm{obs}}(\bm{\theta}^{\prime})p(\bm{\theta}^{\prime}|d_{1})}{p_{\mathrm{PE}}(\bm{\theta}^{\prime})}d\bm{\theta}^{\prime}\int\dfrac{p_{\mathrm{obs}}(\bm{\theta}^{\prime})p(\bm{\theta}^{\prime}|d_{2})}{p_{\mathrm{PE}}(\bm{\theta}^{\prime})}d\bm{\theta}^{\prime}},$
(7)

where $p(\bm{\theta}^{\prime}|d_{1})$ and $p(\bm{\theta}^{\prime}|d_{2})$ are
the shared posterior parameters (all other detector-frame parameters except
effective luminosity distance, arrival times, and phase) of both images,
$p_{\mathrm{obs}}$ is the observed population probability and
$p_{\mathrm{PE}}$ is the simple prior used in GW parameter estimation. As a
result, under the lensed hypothesis, all the images come from the same
astrophysical source. In the meantime, under the unlensed hypothesis, they
come from two different astrophysical sources, resulting in different overall
prior volumes (one population term $p_{\mathrm{obs}}$ in the nominator and two
in the denominator). Overall the Bayes factor is proportional to the prior
volume (inversely proportional to the population prior), meaning that events
with similar parameters are more likely to be lensed if the population
probability is low.

Interpreting the results of the Bayes factor requires some care, because
marginalising over some of the parameters will reduce the discriminatory power
of the Bayes factor; a moderate Bayes factor could become more significant if
a full joint analysis was performed. This would primarily affect our ability
to confidently detect strong lensing. However, in this study, we will not
focus on detectability. Instead, we focus on the effect of the population
model on the Bayes factor, which depends only on the parameters that enter the
BBH population model, masses which we include in the Bayes factor. Because of
this, the Bayes factors should be interpreted as effectively rank-ordering of
the candidates, and not as an absolute indicator of detection (for a
detection, computation of the posterior Odds would be required). To obtain the
Bayes factors, we need both the posterior samples from a parameter estimation
run and the BBH population probability reconstructed with selection effects
(Sec. 3).

####  2.2.1 Posterior generation

As discussed above, the posterior overlap analysis consists of comparing the
posterior probability distribution for the parameters of two strong lensing
candidates (two different images of the same event). To this end, we need to
generate posterior distributions for the various events.

The posterior probability density for a particular parameter depends on two
quantities: firstly, the true value of the considered parameter; secondly, the
specific noise realisation that comes with the GW itself. In particular, in
the absence of noise, and assuming a Gaussian probability density function,
the posterior probability density function would peak at the true value.
However, due to noise fluctuations, the peak is typically shifted. More
specifically, the distribution of the peaks follows a normal distribution
centred around the mean value (e.g. Çalışkan et al., 2023).

In this limit, for a generic parameter $\bm{\theta}^{\prime}$ and the strain
data $d_{i}=h_{i}+n_{i}$ around the $i$th candidate strong GW lensing image,
we can write

$\displaystyle
p(\bm{\theta}^{\prime}|d_{i})=f(\bm{\theta}^{\prime}|\hat{\bm{\theta}}^{\prime},n_{i})\,,$
(8)

where $\hat{\bm{\theta}}^{\prime}$ is the set of true parameters and
$f(\bm{\theta}^{\prime}|\hat{\bm{\theta}}^{\prime},n_{i})$ is a Gaussian
distribution with a mean value drawn from a normal distribution centered
around the true value. If we consider the case in which two different signals
are two images of the same gravitational event and choosing
$\bm{\theta}^{\prime}$ among the parameters which are not affected by the
gravitational lensing (all detector-frame parameters except effective
luminosity distances, arrival times, and overall phases), the differences
among the two posterior distributions are due to the different noise
realisations:

$\displaystyle p(\bm{\theta}^{\prime}|d_{1})$
$\displaystyle=f(\bm{\theta}^{\prime}|\hat{\bm{\theta}}^{\prime},n_{1})$ (9)
$\displaystyle p(\bm{\theta}^{\prime}|d_{2})$
$\displaystyle=f(\bm{\theta}^{\prime}|\hat{\bm{\theta}}^{\prime},n_{2})\,.$
(10)

On the other hand, if $\bm{\theta}^{\prime}$ is affected by the gravitational
lensing, the detector will see the two events as if they had two different
true values for the same parameter. Therefore,

$\displaystyle p(\bm{\theta}^{\prime}|d_{1})$
$\displaystyle=f(\bm{\theta}^{\prime}|\hat{\bm{\theta}}^{\prime}_{1},n_{1})$
(11) $\displaystyle p(\bm{\theta}^{\prime}|d_{2})$
$\displaystyle=f(\bm{\theta}^{\prime}|\hat{\bm{\theta}}^{\prime}_{2},n_{2})\,.$
(12)

In this case, the two distributions are not expected to overlap. Since the
frequency evolution depends on the detector-frame masses
$m_{i}^{z}=(1+z)m_{i}^{\mathrm{s}}$, where $m_{1}^{\mathrm{s}}$ denotes the
source-frame masses, and the evolution is not affected by the gravitational
lensing, the detector-frame masses are expected to be the same between
different GW images and unbiased by strong lensing.

In this study, we only consider detector-frame masses components for posterior
overlap analysis. Instead of generating the $m_{1}^{z},m_{2}^{z}$ posteriors,
we generate $m_{1}^{z},q$ posteriors for lensed event pairs according to the
recipe in Appendix A for more efficient Bayes factor computation.

##  3 Population

We aim to investigate if careless use of a parametric population model or a
clearly incorrect model could bias strong lensing detections and if these
biases can be mitigated by the use of flexible non-parametric models. To
perform the study, we will use four population models:

  1. (i)

a benchmark model that corresponds to the population model used in the mock
catalogue, the correct model;

  2. (ii)

a non-parametric model, which is used to reconstruct the underlying model
encoded in the available data without assumptions.

  3. (iii)

a physically motivated population model that differs from the benchmark model,
the _wrong model_ ;

  4. (iv)

a uniform population model;

In the following, we will give the details of the parametric models as well as
the general idea of the non-parametric model. For the scope of this paper, the
relevant astrophysical properties of the BH population are the BH source-frame
mass (both primary and secondary) and the source redshift.

###  3.1 Benchmark model

In this Section, we describe the population used to generate the mock
catalogue and thus serve as a benchmark model in this work. For the benchmark
model, the primary and secondary masses are taken from the same distribution,
which presents an explicit redshift dependence:

$p(m^{\mathrm{s}}_{1,2},z)=p(m^{\mathrm{s}}_{1,2}|z)p(z)\,.$ (13)

The redshift evolution function, $p(z)$, is motivated by analytical fits to
the star-formation rate density, similarly to (Madau & Dickinson, 2014),

$\displaystyle
p(z)=1+(1+z_{p})^{-\alpha_{z}-\beta_{z}}\frac{(1+z)^{\alpha_{z}}}{1+\quantity(\frac{1+z}{1+z_{p}})^{\alpha_{z}+\beta_{z}}}\,,$
(14)

where $\alpha_{\rm z}=2.7$, $\beta_{\rm z}=2.9$ and $z_{\rm{max}}=1.9$.

On the other hand, the source-frame mass distribution $p(m^{\mathrm{s}}|z)$ is
modeled after the PP model: in this model, the power law accounts for the
initial mass function (IMF), as in (Salpeter, 1955), while the Gaussian
distribution models a pileup around $45M_{\odot}$ that is expected to be due
to certain processes – e.g., the pulsational pair-instability supernova
(PPISN) process – that some classes of stars (e.g., stars with helium core
masses in the range of $\sim 30-50M_{\odot}$) may undergo during their life
(Woosley, 2017). The main difference between the mass distribution we adopted
in this work and the PP model used in (Abbott et al., 2023c) is that we allow
for a redshift dependence in the mass function. In particular,
$p(m^{\mathrm{s}}|z)$ reads

$P(m_{\rm s}|z)=(1-w(z))\times\mathcal{C}m_{\rm s}^{-2}+\\\
+w(z)\mathcal{N}(m_{\rm s}|f_{m}(z))\,,$ (15)

where $\mathcal{C}$ is the normalisation constant of the mass distribution,
$w(z)$ is the redshift-dependent weight,

$\displaystyle w(z)=w_{0}\times\quantity(1+\frac{z}{z_{\rm max}})\,,$ (16)

with $w_{0}\leq 0.12$, and $f_{m}(z)$ is the standard deviation of the
Gaussian distribution

$\displaystyle f_{m}(z)=m_{0}+m_{1}\frac{z}{z_{\mathrm{max}}}\,,$ (17)

with parameters $m_{0}=45$ and $m_{1}=20$.

###  3.2 Non-parametric population model

In this subsection, we will briefly introduce (H)DPGMM, the non-parametric
model we use as an agnostic population model. We refer the interested reader
to (Rinaldi & Del Pozzo, 2022a) and references therein, where the model is
introduced, and to (Rinaldi & Del Pozzo, 2022b), where figaro – the
implementation we used to infer the population distribution – is presented for
more details.

(H)DPGMM is a model based on the Dirichlet process Gaussian mixture model (or
DPGMM) that is designed to infer probability densities given a GW catalogue.
The idea behind the DPGMM is to approximate an arbitrary probability density
with an infinite weighted sum of multivariate Gaussian distributions:

$f(x)\sim\sum^{\infty}_{i}w_{i}N(x|\mu_{i},\sigma^{2}_{i})\,.$ (18)

This model has a countably infinite number of parameters, namely the weights,
the means and the variances. These parameters can be inferred using samples
drawn from the distribution $f(x)$.

We apply this framework to the problem of inferring the primary mass and the
mass ratio of the primary and secondary mass from GW observations. The outcome
of this model is a phenomenological distribution that describes the observed
events.

###  3.3 Wrong model

The _wrong_ population model is a physically inspired population model that
doesn’t properly capture the features which are present in the true underlying
distribution of the data. We included this model to show that biases might
arise from wrong population assumptions.

In particular, the source-frame primary mass distribution and the mass ratio
$q=m_{2}/m_{1}$ follow power-laws:

$\displaystyle p(m_{1})\propto m_{1}^{X}\,,$ (19) $\displaystyle p(q)\propto
q^{Y}\,,$ (20)

where the power-law indices $X,Y$ can be inferred from the observed catalogue.
The redshift distribution is taken to be the same as Eq.[14].

###  3.4 Uniform model

The uniform model prescripts a uniform distribution on all the involved event
parameters. This population model is used in literature – for instance in
(Diego et al., 2021) – to express agnosticism towards the population
distribution and to keep it as conservative as possible. We believe, however,
that such a choice is incorrect and leads to bias in the Bayes factor
computation: in a nutshell, the posterior overlap analysis compares the
probability of observing two events that _may look the same_ by chance with
the probability of having it caused by gravitational lensing. By assuming a
flat population, one is neglecting the possibility of having regions in the
parameter space that are more likely to be observed than others. Such choice,
therefore, artificially enhances the value of the Bayes factor, leading to a
significantly higher chance of misclassifying an unlensed event as lensed. In
order to demonstrate the dramatic effects that this incorrect choice can have,
we included this model in our analysis.

###  3.5 Selection effects

GWs are detected using the so-called _online search pipelines_ , such as MBTA
(Aubin et al., 2021), pyCBC (Usman et al., 2016), GstLAL (Sachdev et al.,
2019), and cWB (Klimenko et al., 2016): these pipelines release triggers,
using the signal-to-noise ratio (SNR) $\rho$ as detection statistics. Whereas
the SNR exceeds a threshold value $\rho_{\mathrm{th}}$, a new event candidate
is identified. The SNR is not the sole detection statistics employed to
determine whether a stretch of data contains an event, it is accompanied with
the false alarm rate (FAR) and the probability of the astrophysical origin of
the event $p_{\mathrm{astro}}$ (Farr et al., 2015; Kapadia et al., 2020).

In principle, all these quantities play a role in the promotion of an event
candidate to confident detection: however, due to the fact that computing the
FAR and $p_{\mathrm{astro}}$ for an event is a computationally expensive task,
in what follows we will consider SNR alone.

Each event will produce a different SNR in the detector depending on two
factors: its parameters (both intrinsic and extrinsic) and the noise
realisation surrounding the event. Due to the stochastic nature of the noise,
the SNR of a specific event will be a realisation of a stochastic process,
whose expected value is determined by the event’s parameters $\bm{\theta}$,
$\hat{\rho}(\bm{\theta})$.

The detection probability of an event follows the probability of observing
$\rho$ greater than $\rho_{\mathrm{th}}$. According to (Allen et al., 2012),
$\rho$ is normally distributed with unitary variance:

$\rho\sim N(\rho|\hat{\rho}(\bm{\theta}),1)\,.$ (21)

Therefore, $p_{\mathrm{det}}$ for a set of parameters $\bm{\theta}$ is given
by

$p_{\mathrm{det}}(\bm{\theta})=1-\int_{0}^{\rho_{\mathrm{th}}}N(\rho|\hat{\rho}(\bm{\theta}),1)\differential\rho\,,$
(22)

where $p_{\mathrm{det}}$ is the _selection function_. Concerning
$\hat{\rho}(\bm{\theta})$, we relied on the python code observation-bias-
gw101010Publicly available at https://github.com/dveske/observation-bias-gw
(Veske et al., 2021), which provides a script to compute the expected SNR
given the primary and secondary mass and the redshift of the source. The
angular parameters of the source (sky position, inclination, and polarisation)
are marginalised out while computing the expected SNR.

###  3.6 Mock data catalogue

We simulated an observed GW catalogue composed of 1658 events using the
benchmark population model with the selection effect. Then we sample
$m_{1}^{z},q$ and $z$ posteriors for each event according to the recipe
Appendix A. The number of events included in the catalogue is a middle ground
between the number of events that would be needed to have a precise
reconstruction with the non-parametric model111111Preliminary studies suggest
that we would need tens of thousands of events..Once the mock posteriors are
generated, we perform a hierarchical population inference to get the (H)DPGMM
and the power-law model.

In order to compute the Bayes factor for lensed pairs, we need the posteriors
for lensed GW pairs. We repeat the previous process to simulate another
astrophysical catalogue. Then for each astrophysical source, we generate two
magnification factors using Eq. 4, 5 for two images. We apply Eq. 2 to
calculate the effective luminosity distance for each image. As a result, the
images from the same source have the same $m_{1}^{z}$ and $q$ but different
luminosity distance. After that, we sample $m_{1}^{z},q$, and $z$ posteriors
for each image separately.

##  4 Results

###  4.1 Population inference

Figure 2:  Population models for observed distribution. Moving from left to
right, the columns correspond to the benchmark model, the inferred (H)DPGMM
model, the inferred power-law model, and the uniform model, respectively. The
first row represents the marginalised distribution of redshifted primary mass,
while the second row represents the marginalised distribution of mass ratio.

We perform the population inference to obtain
$p_{\mathrm{obs}}(m_{1}^{z},q,z)$ for the (H)DPGMM using the mock observed
catalogue prepared in Sec. 3.6. We then marginalise over $z$ to get
$p_{\mathrm{obs}}(m_{1}^{z},q)$. For the power-law model, we infer the power-
law indices using the mock observed catalogue to obtain
$p_{\mathrm{obs}}(m_{1},q,z)$. The power-law indices for $m_{1},q$ are
$0.76\pm 0.02,-0.28\pm 0.04$ with one standard deviation. Then, we convert it
to be in detector-frame and marginalise over $z$ to get
$p_{\mathrm{obs}}(m_{1}^{z},q)$. Fig. 2 shows the marginalised $m_{1}^{z}$ and
$q$ observed distribution for all population models.

Now, we have 4 population models for the observed distribution:

  1. (i)

a benchmark model with the high-dimensional Gaussian peak, which provides the
most optimal inference;

  2. (ii)

an inferred (H)DPGMM, which makes fewer model assumptions and could be used in
future lensing analyses;

  3. (iii)

an inferred power-law model, which assumes the primary mass and mass ratio
follow power-law relation when we do not have enough data from observation;

  4. (iv)

a uniform prior, which assumes all the parameters follow a uniform
distribution.

We use these four models to evaluate the Bayes factor in favour of strong
lensing, with the purpose of investigating the population model’s impact on
lensed detections.

###  4.2 Population model effects on Bayes factor

Figure 3:  The percentage error of lensed pair Bayes factors
$\mathcal{B}_{{}^{U}}^{{}^{L}}$ on $m_{1}^{z}$ posterior overlap using
different marginalised $m_{1}^{z}$ observed population models. The top panel
shows the marginalised $m_{1}^{z}$ observed distribution from the benchmark
model (black), (H)DPGMM (blue) with 68% and 95% confidence intervals, and the
power-law model (orange). The bottom panel shows the percentage difference of
the Bayes factors using (H)DPGMM (blue), power-law model (orange) compare to
the benchmark model with 1 standard deviation. The Bayes factors using
(H)DPGMM give $\sim 10\%$ error which is smaller than using the power-law
model for most of the lensed pairs.

We first evaluate the Bayes factor using only $m_{1}^{z}$ posteriors of lensed
pairs and marginalised $m_{1}^{z}$ observed population model.

We randomly select 4 event pairs for each $20M_{\odot}$ interval from the
lensed catalogue to perform the Bayes factor calculation for better
illustration. In the upper panel of Fig. 3, we show the comparison of
benchmark, (H)DPGMM, and power-law models. The (H)DPGMM agrees with the
benchmark model within the 68% confidence interval while the power-law model
does not. As a result, the percentage errors for lensed pair Bayes factors
between using (H)DPGMM and benchmark model (blue) are smaller than those
between using the power-law model and benchmark model (orange)(Fig. 3 lower
panel). The Bayes factors calculated with (H)DPGMM give $\pm 10\%$ error with
1 standard deviation while those using the power-law model give a larger
percentage error for most of the lensed pairs. Within the range of
$10M_{\odot}<m_{1}^{z}<100M_{\odot}$, the power-law model exhibits a lower
probability compared to the benchmark model. It favors the probability of
having lensed pair when two observed events have similar event posteriors. As
a result, the Bayes factors are overestimated, and the opposite occurs for
$m_{1}^{z}>100M_{\odot}$ as the power-law model has a higher probability
compared to the benchmark model. The result shows how the population model
affects the Bayes factor. A parametric model that fails to capture the
population’s structures will bias the Bayes factor.

Figure 4:  The percentage error of lensed pair Bayes factors
$\mathcal{B}_{{}^{U}}^{{}^{L}}$ on $m_{1}^{z},q$ posterior overlap using
different observed population models. The top panel shows the marginalised
$m_{1}^{z}$ observed distribution from the benchmark model (black), (H)DPGMM
(blue) with 68% and 95% confidence intervals, and the power-law model
(orange). The bottom panel shows the percentage difference of the Bayes
factors using (H)DPGMM (blue), power-law model (orange) compare to the
benchmark model with 1 standard deviation. The Bayes factors using (H)DPGMM
give $\sim 20-40\%$ error which is smaller than using the power-law model for
most of the lensed pairs.

Next, we evaluate the Bayes factor using both $m_{1}^{z}$ and $q$ posteriors
of lensed pairs to see the error propagation in higher dimensions. As shown in
Fig. 4, the percentage error for using (H)DPGMM increased from around $\pm
10\%$ to around $\pm 20-40\%$ while half of the lensed pair Bayes factors
using power-law model give $>100\%$ error. In Fig. 2, the (H)DPGMM uncertainty
on marginalised $q$ distribution is greater than $m_{1}^{z}$ distribution,
increasing the percentage error when calculating the Bayes factor using both
$m_{1}^{z}$ and $q$ posteriors. The power-law model, on the other hand, gives
greater inaccuracy for some lensed pairs because it cannot capture all
features of the joint $(m_{1}^{z},q)$ distribution.

###  4.3 Prior assumptions on Bayes factor

Figure 5:  The percentage error of lensed pair Bayes factors
$\mathcal{B}_{{}^{U}}^{{}^{L}}$ on $m_{1}^{z}$ posterior overlap using the
uniform model with different ranges compared to the benchmark observed
population model. The top panel shows the marginalised $m_{1}^{z}$ observed
distribution from the benchmark model (black), the uniform model in
$(5,140)M_{\odot}$ (green solid line), and the uniform model in
$(5,300)M_{\odot}$ (green dashed line). The bottom panel shows the percentage
difference of the Bayes factors using a uniform model in $(5,140)M_{\odot}$
(green circle), a uniform model in $(5,300)M_{\odot}$ (green diamond) compare
to using the benchmark model with 1 standard deviation. Using a uniform model
bias the Bayes factor.

However, what about the effect of the population model on existing detections?
(Diego et al., 2021) advocated that the Bayes factor evaluated using a uniform
prior is less biased than a $\mathcal{B}_{{}^{U}}^{{}^{L}}$ utilizing a prior
informed by population analyses. In Fig. 5, we consider uniform $m_{1}^{z}$
prior compared to the benchmark observed model with different ranges. The
uniform prior overestimates the $\mathcal{B}_{{}^{U}}^{{}^{L}}$ in the low
mass region since it has no peak feature compared to the benchmark model. In
the high mass region, it either overestimates or underestimates the
$\mathcal{B}_{{}^{U}}^{{}^{L}}$ depending on the chosen $m_{1}^{z}$ range. The
result demonstrates that assuming a uniform prior can significantly bias the
$\mathcal{B}_{{}^{U}}^{{}^{L}}$ depending on the range of the uniform prior.
The population probability strongly depends on the volume of the parameter
space, a larger parameter space gives a smaller population probability for a
uniform prior. For our benchmark model, the probability remains about the same
as long as we set the upper limit of $m_{1}^{z}$ to be larger than
$140M_{\odot}$ since most of the event is concentrated in the region
$m_{1}^{z}<140M_{\odot}$. This makes the Bayes factor using a uniform prior
arbitrary depending on the range we set and gives inaccurate calculation. On
the other hand, the result may be counter-intuitive because a uniform prior is
typically assumed to yield conservative results. However, the strong Bayes
factor contrasts the probability that two events are strongly lensed against
astrophysical coincidence, which depends sensitively on the prior knowledge of
the BBH population. Indeed, a uniform population model would produce very
dissimilar BBHs and therefore highly dissimilar waveforms. Thus, under this
model, it would be predicted that overlap in two waveforms (and thus the
probability of lensing mimickers due to astrophysical coincidence) is
unlikely, giving false credence to the strong lensing hypothesis. In contrast,
a more realistic population model would place greater weight on the
astrophysical coincidence. Therefore, using a uniform mass prior does not
better discriminate strongly lensed events than models that presume more about
the population of binaries and lenses. Instead, using a uniform mass prior can
easily result in a number of false detections of strong lensing (See Section
4.4).

Other than assuming a uniform population prior, if we assume the prior used in
parameter estimation to be the same as the observed population (i.e.,
$p_{\mathrm{PE}}=p_{\mathrm{obs}}$), the Bayes factor (Eq. 7) will reduce to

$\displaystyle{\mathcal{B}_{{}^{U}}^{{}^{L}}}^{*}=\int\frac{p(\bm{\theta}|d_{1})p(\bm{\theta}|d_{2})}{p_{\mathrm{obs}}(\bm{\theta})}d\bm{\theta},$
(23)

where ${\mathcal{B}_{{}^{U}}^{{}^{L}}}^{*}$ denotes the Bayes factor with such
an assumption. The Bayes factor (Eq. 23) is often used in various GW lensing
literature but it does not correctly account for the population probability.
We compare the Bayes factor difference using (Eq. 7 and Eq. 23) with the
benchmark population model.

Figure 6:  The percentage error of lensed pair Bayes factors on $m_{1}^{z},q$
posterior overlap when assuming prior used in parameter estimation to be the
same as the observed population ($p_{\mathrm{PE}}=p_{\mathrm{obs}}$) versus
without such an assumption (${\mathcal{B}_{{}^{U}}^{{}^{L}}}^{*}$ versus
$\mathcal{B}_{{}^{U}}^{{}^{L}}$). Bayes factors by assuming
$p_{\mathrm{PE}}=p_{\mathrm{obs}}$ give stronger bias with smaller
$\mathcal{B}_{{}^{U}}^{{}^{L}}$ magnitude.

Fig. 6 shows the percentage difference between the Bayes factors with and
without assuming $p_{\mathrm{PE}}=p_{\mathrm{obs}}$, which denote by
${\mathcal{B}_{{}^{U}}^{{}^{L}}}^{*}$ and $\mathcal{B}_{{}^{U}}^{{}^{L}}$
respectively. The ${\mathcal{B}_{{}^{U}}^{{}^{L}}}^{*}$ is biased stronger
with a smaller $\mathcal{B}_{{}^{U}}^{{}^{L}}$. As the Bayes factor is
inversely proportional to $p_{\mathrm{obs}}$, smaller $p_{\mathrm{obs}}$ gives
a larger Bayes factor and a larger discrepancy between
$\mathcal{B}_{{}^{U}}^{{}^{L}}$ and ${\mathcal{B}_{{}^{U}}^{{}^{L}}}^{*}$. It
is critical to properly assess the effect of the population model in low-
probability regions since two events with similar parameters in such regions
are very likely to be lensed pair than astrophysical coincidence.

###  4.4 Detection performance

In order to understand how the bias in the Bayes factor affects the strong
lensing detection which is essentially a classification problem using the
Bayes factor, we plot the receiver operating characteristic (ROC) curves. We
use 200 lensed pairs and 5000 unlensed pair to generate the ROC curves by
showing the true positive rate (TPR) and false positive rate (FPR), where TPR
is the portion of lensed pairs classified correctly as lensed detection and
FPR is the portion of unlensed pairs classified wrongly as lensed detection
using Bayes factor.

Figure 7:  Receiver operating characteristic (ROC) curves of strong lensing
classification using benchmark population model (black), (H)DPGMM (blue),
power-law model (orange), and uniform model (green). The ROC curve using
(H)DPGMM gives the most comparable performance to the benchmark population
model.

In Fig. 7, the ROC curve using the benchmark model (black) shows the best
performance on strong lensing detection. The ROC curve using the (H)DPGMM
(blue) gives the most comparable performance to the benchmark ROC curve when
compared to the ROC curves using the power-law (orange) and uniform model
(green). This again demonstrates non-parametric population model is more
suitable for strong lensing detections.

Figure 8:  Receiver operating characteristic (ROC) curves of strong lensing
classification when assuming prior used in parameter estimation to be the same
as the population ($p_{\mathrm{PE}}=p_{\mathrm{obs}}$) versus without such an
assumption (${\mathcal{B}_{{}^{U}}^{{}^{L}}}^{*}$ versus
$\mathcal{B}_{{}^{U}}^{{}^{L}}$). ROC curve without such assumption gives
better performance for the classification.

Fig. 8 shows the ROC curve comparison on utilising/not utilising the
assumption $p_{\mathrm{PE}}=p_{\mathrm{obs}}$. The ROC performance is better
if such an assumption is not used, implying that Eq. 7 should be used for
strong lensing detection.

Indeed, we advocate for the use of the corrected Bayes factor formula for the
posterior overlap analyses (Eq. 7). We note that the posterior overlap Bayes
factor is not used in the joint parameter estimation Bayes factor Liu et al.
(2020); Lo & Magaña Hernandez (2023); Janquart et al. (2021a) and so the
discussion on the Bayes factor formulation here does not apply to them.

##  5 Discussion and Conclusions

The strong lensing analysis relies on accurate modelling of the population of
BBHs and strong lenses: an incorrect population model could lead to either
inappropriate or missed strong lensing detections.

Here we have investigated the use of non-parametric population models to
detect strong lensed events. Unlike parametric modelizations based on
phenomenological fits, the non-parametric ones can capture more complicated
formation scenarios without assuming a functional form a priori. We find that
the use of a non-parametric population model as inferred from mock GW
catalogue is suitable for strong lensing analyses: while the Bayes factor will
end up with broader error bounds than a parametric model that appropriately
captures the population, the Bayes factor is less biased by prior assumptions.

Nevertheless, the convergence of the non-parametric (H)DPGMM observed
population model depends on the underlying population. Given that the non-
parametric model is informed by the data alone and not by the physics behind
it, the number of detected events required to properly capture all the
features of the underlying distribution depends on the distribution itself –
the more complex the distribution, the more events are needed. Nonetheless,
some preliminary studies that we run on a plausible model based on the 3D
$m1-q-z$ model used in this work show that we need order of thousands of
events.

Furthermore, (Diego et al., 2021) suggests the analysis with no population
prior assumption provides better merger-rate estimation and hence less biased
results than the analysis with a PP model in (The LIGO Scientific
Collaboration et al., 2021). Here, we have shown that a uniform population
prior is not an appropriate choice for strong lensing analyses. It strongly
biases the result of our analysis on the parameter overlap without considering
the merger rate and time-delay distribution. The uniform prior flattens out
the probability and hence gives optimistic estimates of the strong lensing
Bayes factor for most of the lensed events, which can lead to more false
detections of strong lensing as shown in Fig. 7.

In addition, the selection effect would not take place in the reconstruction
of the population from observed data, which means, we should use the observed
population instead of the astrophysical population. In Eq. 7, we include the
population term accounts for the astrophysical coincidence, however, we can
not see all the astrophysical coincidences because of the selection effect.
Therefore, instead of using an astrophysical population, we should use the
observed population when computing the Bayes factor. The Bayes factor is used
to rank the possibility of lensing among the possible pair seen by our
detector instead of all the possible pairs in the Universe. Therefore, by
using the observed population, the Bayes factor for unlensed pair with similar
parameters will be penalised by accounting for the probability of coincidence
that’s seen by the detector. This saves lots of effect on reconstructing the
astrophysical population from the observed population.

We have employed simplifying assumptions in our modelling of the BBH
population, parameter inference, and the strong lensing Bayes factor
evaluation. Specifically, our parameter inference considers a simplified,
parametrised model consistent with current observations, to avoid excessive
computational cost. Furthermore, our strong lensing Bayes factor is evaluated
using the overlap in mass and mass-ratio using the posterior overlap method of
Haris et al. 2018, because a full joint parameter estimation evaluating the
overlap in the full 15-dimensional BBH parameter space (see (e.g. Janquart et
al., 2021a)) could not be used in conjunction with the parametrized model that
we use for parameter inference. A full analysis employing joint parameter
estimation would likely yield higher Bayes factors and would therefore be a
better discriminator for whether an event is strongly lensed or not. However,
with the current methodologies, an extensive joint parameter study is
computationally too exhaustive to be applied to all possible pairs of lensed
candidates found in a GW catalogue. Furthermore, the posterior overlap
analysis, equipped with the non-parametric population model, can be used to
reduce the computational burden of joint PE by pre-filtering the pairs that
make a good lensed candidate. However, we note that including more dimensions
in the Bayes factor computation would exasperate the problem: an incorrect
population model would lead to even greater discrepancies. While these
simplified models are adequate for studying the applicability of non-
parametric models to strong lensing analyses, future studies utilising strong
lensing studies, especially on real data should perform joint parameter
estimation, which can better quantify if an event is finally strongly lensed
or not. Indeed, we hope that our work motivates the inclusion of non-
parametric models as a core part of strong GW lensing research.

##  6 Acknowledgements

The authors are grateful to Walter Del Pozzo and Nicola Tamanini for the
useful discussions. The authors are thankful to Shasvath J. Kapadia for useful
comments on this manuscript. MT gratefully acknowledges support from the
CNRS/IN2P3 Computing Center (Lyon - France) for providing computing and data-
processing resources needed for this work. OAH acknowledges support by grants
from the Research Grants Council of Hong Kong (Project No. CUHK 2130822 and
4443086), and the Direct Grant for Research from the Research Committee of The
Chinese University of Hong Kong.

##  Appendix A Mock posterior

Here we present a simple mock posterior generation for GW events in our mock
catalogue. For each GW event with true event parameter $(m_{1}^{z},q,z)$, we
want to generate $(m_{1}^{z,\mathrm{obs}},q^{\mathrm{obs}},z^{\mathrm{obs}})$
samples that mimic the inference result from parameter estimation using the
strain data.

First, due to the fact that noise is included in the real strain data, the
maximum likelihood from parameter estimation is not referring to the true
event parameter. In addition, the noise introduced measurement uncertainties
in the parameter estimation. To include these effects, we first draw a sample
$\hat{x}$ from the true value and then sample the posterior from a
distribution characterized by $(\hat{x},\sigma_{x})$, where $\sigma_{x}$ is
the parameter that takes the measurement uncertainty into account.

First, consider detector-frame primary mass $m_{1}^{z}$, we draw a sample
$\hat{m}_{1}^{z}$ by

$\log\hat{m}_{1}^{z}=\mathcal{N}(\log m_{1}^{z},\sigma_{m}),$ (24)

where $\mathcal{N}(\mu,\sigma)$ denotes a normal distribution with mean $\mu$
and standard deviation $\sigma$. Then, we sample the posteriors from the same
distribution but shifted the center to $\log\hat{m}_{1}^{z}$ as

$\log m_{1}^{z,\mathrm{obs}}=\mathcal{N}(\log\hat{m}_{1}^{z},\sigma_{m}),$
(25)

where $m_{1}^{z,\mathrm{obs}}$ is the detector-frame primary mass posteriors.

Then, we consider the mass-ratio component. We first draw a variable $y$ by

$y=\mathcal{N}(\log\frac{q}{1-q},\sigma_{q}),$ (26)

then we transform it back to a mass-ratio by

$\hat{q}=\frac{e^{y}}{1+e^{y}}.$ (27)

Then, we can sample the mass-ratio posteriors $q^{\mathrm{obs}}$ by replacing
$q$ to $\hat{q}$ in Eq. 26 and 27.

For redshift, we use the same recipe as $m_{1}^{z}$ which draws a sample
$\hat{z}$ from

$\log\hat{z}=\mathcal{N}(\log z,\sigma_{z}),$ (28)

and then we sample the redshift posterior $z^{\mathrm{obs}}$ by

$\log z^{\mathrm{obs}}=\mathcal{N}(\log\hat{z},\sigma_{z}).$ (29)

The measurement uncertainty parameters are chosen to be
$\sigma_{m}=0.088,\sigma_{q}=1.0671,\sigma_{z}=0.3$ respectively. The choices
of $\sigma_{m},\sigma_{q},\sigma_{q}$ yield measurement uncertainties that
match the similar width of credible intervals in O1 and O2 data (Vitale et
al., 2017). By following Eq. [24-29], we sample the posterior for each GW
event in the mock catalogue. The posterior samples represent the likelihood of
the observed GW event parameters.

##  Appendix B Odds Ratio derivation

Here, we present the derivation of $\mathcal{B}_{{}^{U}}^{{}^{L}}$ in Eq. 7
following the formalism in Lo & Magaña Hernandez (2023). To determine whether
a set of gravitational-wave events are lensed or not, we make use of the
quantity called the odds ratio

$\mathcal{O}_{{}^{U}}^{{}^{L}}=\frac{p(\mathcal{H}_{L}|\bm{d})}{p(\mathcal{H}_{U}|\bm{d})}=\frac{p(\mathcal{H}_{L})}{p(\mathcal{H}_{U})}\frac{p(\bm{d}|\mathcal{H}_{L})}{p(\bm{d}|\mathcal{H}_{U})}=\mathcal{P}_{{}^{U}}^{{}^{L}}\mathcal{B}_{{}^{U}}^{{}^{L}},$
(30)

where $\bm{d}=\\{d_{1},\ldots d_{N}\\}$ is the data set contains $N$ signals,
$\mathcal{P}_{{}^{U}}^{{}^{L}}$ is the prior odds and
$\mathcal{B}_{{}^{U}}^{{}^{L}}$ is the Bayes factor. And lensed and unlensed
hypotheses are denoted as:

  * •

$\mathcal{H}_{L}$: The data set $\bm{d}$ contains lensed signals from a single
binary black hole merger event with parameters
$\bm{\theta}^{(1)}=...=\bm{\theta}^{(N)}$.

  * •

$\mathcal{H}_{U}$: The data set $\bm{d}$ contains signals from independent
binary black hole merger events with parameters $\bm{\theta}^{(i)}$.

Here, we focus on the derivation of the Bayes factor. First, we focus on
$p(\bm{d}|\mathcal{H}_{L})$ in $\mathcal{B}_{{}^{U}}^{{}^{L}}$. We can expand
the probability in event parameter $\bm{\theta}$ space. In addition, we can
separate $\bm{\theta}$ into two independent subsets:
$\bm{\theta}_{\mathrm{U}}$ and $\bm{\theta}_{\mathrm{eff}}$, where

  * •

$\bm{\theta}_{\mathrm{U}}$: parameters that do not affect by the lensing
mechanism, such as masses, and spins.

  * •

$\bm{\theta}_{\mathrm{eff}}$: parameters that will be affected by the lensing
mechanism such as luminosity distance, arrival time, and phase difference.

Then we have:

$p(\bm{d}|\mathcal{H}_{L})=\int
p(\bm{d}|\bm{\theta}_{\mathrm{U}},\bm{\theta}_{\mathrm{eff}})p(\bm{\theta}_{\mathrm{U}}|\mathcal{H}_{L})p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{L})d\bm{\theta}_{\mathrm{U}}d\bm{\theta}_{\mathrm{eff}},$
(31)

where $p(\bm{\theta}_{\mathrm{U}}|\mathcal{H}_{L})$and
$p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{L})$ are the priors used for
$\bm{\theta}_{\mathrm{U}}$and $\bm{\theta}_{\mathrm{eff}}$ under lensed
hypothesis. Then, we can rewrite it with the following trick:

$\displaystyle p(\bm{d}|\mathcal{H}_{L})$ $\displaystyle=\int
p(\bm{d}|\bm{\theta}_{\mathrm{U}},\bm{\theta}_{\mathrm{eff}})p(\bm{\theta}_{\mathrm{U}}|\mathcal{H}_{L})p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{L})d\bm{\theta}_{\mathrm{U}}d\bm{\theta}_{\mathrm{eff}}$
(32) $\displaystyle=\int
p(\bm{d}|\bm{\theta}_{\mathrm{U}},\bm{\theta}_{\mathrm{eff}})p(\bm{\theta}_{\mathrm{U}}|\mathcal{H}_{L})p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{L})\dfrac{p(\bm{\theta}_{\mathrm{U}}|\mathcal{H}_{0})p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{0})}{p(\bm{\theta}_{\mathrm{U}}|\mathcal{H}_{0})p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{0})}d\bm{\theta}_{\mathrm{U}}d\bm{\theta}_{\mathrm{eff}}$
$\displaystyle=\int
p(\bm{d}|\bm{\theta}_{\mathrm{U}},\bm{\theta}_{\mathrm{eff}})p(\bm{\theta}_{\mathrm{U}}|\mathcal{H}_{0})p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{0})\dfrac{p(\bm{\theta}_{\mathrm{U}}|\mathcal{H}_{L})p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{L})}{p(\bm{\theta}_{\mathrm{U}}|\mathcal{H}_{0})p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{0})}d\bm{\theta}_{\mathrm{U}}d\bm{\theta}_{\mathrm{eff}},$

where $p(\bm{\theta}_{\mathrm{U}}|\mathcal{H}_{0})$and
$p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{0})$ are the priors used for
$\bm{\theta}_{\mathrm{U}}$and $\bm{\theta}_{\mathrm{eff}}$ in parameter
estimation. Here, we can use Bayes’ theorem to write the product of the
posterior and the evidence as

$p(\bm{\theta}_{\mathrm{U}},\bm{\theta}_{\mathrm{eff}}|\bm{d})p(\bm{d}|\mathcal{H}_{0})=p(\bm{d}|\bm{\theta}_{\mathrm{U}},\bm{\theta}_{\mathrm{eff}})p(\bm{\theta}_{\mathrm{U}}|\mathcal{H}_{0})p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{0}).$
(33)

We apply Eq.[33] on Eq.[32] to replace the first two terms in the integral and
it gives:

$p(\bm{d}|\mathcal{H}_{L})=p(\bm{d}|\mathcal{H}_{0})\int
p(\bm{\theta}_{\mathrm{U}},\bm{\theta}_{\mathrm{eff}}|\bm{d})\dfrac{p(\bm{\theta}_{\mathrm{U}}|\mathcal{H}_{L})p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{L})}{p(\bm{\theta}_{\mathrm{U}}|\mathcal{H}_{0})p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{0})}d\bm{\theta}_{\mathrm{U}}d\bm{\theta}_{\mathrm{eff}}.$
(34)

Notice that we assumed the signals are not overlapping and the parameter
estimations are not correlated for different events. Therefore, we can write
it as:

$\displaystyle p(\bm{d}|\mathcal{H}_{L})$
$\displaystyle=p(\bm{d}|\mathcal{H}_{0})\int\dfrac{p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{L})\prod^{N}_{i}p(\bm{\theta}_{\mathrm{U}}^{(i)}|\mathcal{H}_{L})p(\bm{\theta}_{\mathrm{U}}^{(i)},\bm{\theta}_{\mathrm{eff}}^{(i)}|d_{i})}{\prod_{i}p(\bm{\theta}_{\mathrm{U}}^{(i)}|\mathcal{H}_{0})p(\bm{\theta}_{\mathrm{eff}}^{(i)}|\mathcal{H}_{0})}d\bm{\theta}_{\mathrm{U}}^{(1)}\ldots
d\bm{\theta}_{\mathrm{U}}^{(N)}d\bm{\theta}_{\mathrm{eff}},$ (35)

And we expect $\bm{\theta}_{\mathrm{U}}^{(i)}$ to be the same across $N$
events under the lensing hypothesis. Therefore, we have $(N-1)$ Dirac-delta
distribution to reduce the dimension of the integral:

$\displaystyle p(\bm{d}|\mathcal{H}_{L})$
$\displaystyle=p(\bm{d}|\mathcal{H}_{0})\int\dfrac{p(\bm{\theta}_{\mathrm{U}}^{(N)}|\mathcal{H}_{L})p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{L})\prod^{N}_{i}p(\bm{\theta}_{\mathrm{U}}^{(i)},\bm{\theta}_{\mathrm{eff}}^{(i)}|d_{i})}{\prod_{i}p(\bm{\theta}_{\mathrm{U}}^{(i)}|\mathcal{H}_{0})p(\bm{\theta}_{\mathrm{eff}}^{(i)}|\mathcal{H}_{0})}d\bm{\theta}_{\mathrm{U}}^{(N)}d\bm{\theta}_{\mathrm{eff}},\,\bm{\theta}^{(1)}=...=\bm{\theta}^{(N)}.$
(36)

where $p(\bm{d}|\mathcal{H}_{0})$ is the evidence from Bilby, the first two
terms of nominator in $\langle...\rangle$ are the population prior, followed
by the event posteriors from Bilby. And the denominator is the Bilby prior.

We can do the same on the unlensed hypothesis,

$\displaystyle p(\bm{d}|\mathcal{H}_{U})$ $\displaystyle=\int
p(\bm{d}|\bm{\theta})p(\bm{\theta}|\mathcal{H}_{U})d\bm{\theta}=\int
p(\bm{d}|\bm{\theta})p(\bm{\theta}|\mathcal{H}_{U})\dfrac{p(\bm{\theta}|\mathcal{H}_{0})}{p(\bm{\theta}|\mathcal{H}_{0})}d\bm{\theta}=p(\bm{d}|\mathcal{H}_{0})\int
p(\bm{\theta}|\bm{d})\dfrac{p(\bm{\theta}|\mathcal{H}_{U})}{p(\bm{\theta}|\mathcal{H}_{0})}d\bm{\theta}$
(37) $\displaystyle=p(\bm{d}|\mathcal{H}_{0})\prod_{i}\int
p(\bm{\theta}^{(i)}|d_{i})\dfrac{p(\bm{\theta}^{(i)}|\mathcal{H}_{U})}{p(\bm{\theta}^{(i)}|\mathcal{H}_{0})}d\bm{\theta}^{(i)}.$

In conclusion:

$\displaystyle\mathcal{B}_{{}^{U}}^{{}^{L}}$
$\displaystyle=\frac{p(\bm{d}|\mathcal{H}_{L})}{p(\bm{d}|\mathcal{H}_{U})}=\dfrac{\int\dfrac{p(\bm{\theta}_{\mathrm{U}}^{(N)}|\mathcal{H}_{L})p(\bm{\theta}_{\mathrm{eff}}|\mathcal{H}_{L})\prod^{N}_{i}p(\bm{\theta}_{\mathrm{U}}^{(i)},\bm{\theta}_{\mathrm{eff}}^{(i)}|d_{i})}{\prod_{i}p(\bm{\theta}_{\mathrm{U}}^{(i)}|\mathcal{H}_{0})p(\bm{\theta}_{\mathrm{eff}}^{(i)}|\mathcal{H}_{0})}d\bm{\theta}_{\mathrm{U}}^{(N)}d\bm{\theta}_{\mathrm{eff}}}{\prod_{i}\int
p(\bm{\theta}^{(i)}|d_{i})\dfrac{p(\bm{\theta}^{(i)}|\mathcal{H}_{U})}{p(\bm{\theta}^{(i)}|\mathcal{H}_{0})}d\bm{\theta}^{(i)}}.$
(38)

For our case, we assume only two images are produced. In addition, only masses
components are considered (no $\bm{\theta}_{\mathrm{eff}}$).

$\displaystyle\mathcal{B}_{{}^{U}}^{{}^{L}}$
$\displaystyle=\dfrac{\int\dfrac{p(\bm{\theta}|\mathcal{H}_{L})p(\bm{\theta}|d_{1})p(\bm{\theta}|d_{2})}{p(\bm{\theta}|\mathcal{H}_{0})p(\bm{\theta}|\mathcal{H}_{0})}d\bm{\theta}}{\int
p(\bm{\theta}^{(1)}|d_{1})\dfrac{p(\bm{\theta}^{(1)}|\mathcal{H}_{U})}{p(\bm{\theta}^{(1)}|\mathcal{H}_{0})}d\bm{\theta}^{(1)}\int
p(\bm{\theta}^{(2)}|d_{2})\dfrac{p(\bm{\theta}^{(2)}|\mathcal{H}_{U})}{p(\bm{\theta}^{(2)}|\mathcal{H}_{0})}d\bm{\theta}^{(2)}}$
(39)

Also, we assume the population prior under $\mathcal{H}_{L}$ and
$\mathcal{H}_{U}$ to be the same since the probability of strong lensing is
very low. We denote the astrophysical population prior to be $p_{\mathrm{A}}$
and the simple prior used in Bilby to be $p_{\mathrm{PE}}$. Then the Bayes
factor can be expressed as:

$\displaystyle\mathcal{B}_{{}^{U}}^{{}^{L}}$
$\displaystyle=\dfrac{\int\dfrac{p_{\mathrm{A}}(\bm{\theta}^{\prime})p(\bm{\theta}^{\prime}|d_{1})p(\bm{\theta}^{\prime}|d_{2})}{p_{\mathrm{PE}}^{2}(\bm{\theta}^{\prime})}d\bm{\theta}^{\prime}}{\int\dfrac{p_{\mathrm{A}}(\bm{\theta}^{\prime})p(\bm{\theta}^{\prime}|d_{1})}{p_{\mathrm{PE}}(\bm{\theta}^{\prime})}d\bm{\theta}^{\prime}\int\dfrac{p_{\mathrm{A}}(\bm{\theta}^{\prime})p(\bm{\theta}^{\prime}|d_{2})}{p_{\mathrm{PE}}(\bm{\theta}^{\prime})}d\bm{\theta}^{\prime}},$
(40)

where $p(\bm{\theta}^{\prime}|d_{1})$ and $p(\bm{\theta}^{\prime}|d_{2})$ are
the shared posterior parameters of images 1 and 2 respectively. If we use a
simple uniform prior in the parameter estimation, the formula is simplified to

$\displaystyle\mathcal{B}_{{}^{U}}^{{}^{L}}$ $\displaystyle=\dfrac{\int
p_{\mathrm{A}}(\bm{\theta}^{\prime})p(\bm{\theta}^{\prime}|d_{1})p(\bm{\theta}^{\prime}|d_{2})(\bm{\theta}^{\prime})d\bm{\theta}^{\prime}}{\int
p_{\mathrm{A}}(\bm{\theta}^{\prime})p(\bm{\theta}^{\prime}|d_{1})d\bm{\theta}^{\prime}\int
p_{\mathrm{A}}(\bm{\theta}^{\prime})p(\bm{\theta}^{\prime}|d_{2})d\bm{\theta}^{\prime}}.$
(41)

## References

  * Abbott et al. (2021) Abbott R., et al., 2021, SoftwareX, 13, 100658
  * Abbott et al. (2023b) Abbott R., et al., 2023b, Open data from the third observing run of LIGO, Virgo, KAGRA and GEO (arXiv:2302.03676)
  * Abbott et al. (2023a) Abbott R., et al., 2023a
  * Abbott et al. (2023c) Abbott R., et al., 2023c, Phys. Rev. X, 13, 011048
  * Allen et al. (2012) Allen B., Anderson W. G., Brady P. R., Brown D. A., Creighton J. D. E., 2012, Phys. Rev. D, 85, 122006
  * Aubin et al. (2021) Aubin F., et al., 2021, Class. Quant. Grav., 38, 095004
  * Baker & Trodden (2016) Baker T., Trodden M., 2016, Physical Review D, 95
  * Bartelmann (2010) Bartelmann M., 2010, Class. Quant. Grav., 27, 233001
  * Basak et al. (2022) Basak S., Ganguly A., Haris K., Kapadia S., Mehta A. K., Ajith P., 2022, The Astrophysical Journal Letters, 926, L28
  * Bond et al. (2004) Bond I. A., et al., 2004, The Astrophysical Journal, 606, L155
  * Broadhurst et al. (2018) Broadhurst T., Diego J. M., Smoot III G., 2018, doi:10.48550/arxiv.1802.05273
  * Broadhurst et al. (2019) Broadhurst T., Diego J. M., Smoot III G. F., 2019, doi:10.48550/arxiv.1901.03190
  * Bulashenko & Ubach (2022) Bulashenko O., Ubach H., 2022, Journal of Cosmology and Astroparticle Physics, 2022, 022
  * Buscicchio et al. (2020a) Buscicchio R., Moore C. J., Pratten G., Schmidt P., Vecchio A., 2020a, Physical Review D, 102
  * Buscicchio et al. (2020b) Buscicchio R., Moore C. J., Pratten G., Schmidt P., Bianconi M., Vecchio A., 2020b, Physical Review Letters, 125, 141102
  * Cao et al. (2014) Cao Z., Li L.-F., Wang Y., 2014, Physical Review D, 90, 62003
  * Cao et al. (2022a) Cao S., Qi J., Cao Z., Biesiada M., Cheng W., Zhu Z.-H., 2022a, Astronomy and Astrophysics, 659
  * Cao et al. (2022b) Cao M.-D., Zheng J., Qi J.-Z., Zhang X., Zhu Z.-H., 2022b, The Astrophysical Journal, 934, 108
  * Cheung et al. (2020) Cheung M. H. Y., Gais J., Hannuksela O. A., Li T. G. F., 2020, Monthly Notices of the Royal Astronomical Society, 503, 3326
  * Christian et al. (2018) Christian P., Vitale S., Loeb A., 2018, Physical Review D, 98, 103022
  * Chung & Li (2021) Chung A. K. W., Li T. G., 2021, Physical Review D, 104, 124060
  * Clowe et al. (2004) Clowe D., Gonzalez A., Markevitch M., 2004, ApJ, 604, 596
  * Collett & Bacon (2016) Collett T. E., Bacon D., 2016, Physical Review Letters, 118
  * Congdon & Keeton (2018) Congdon A. B., Keeton C. R., 2018, Principles of Gravitational Lensing. Springer International Publishing, doi:10.1007/978-3-030-02122-1
  * Cremonese et al. (2021) Cremonese P., Ezquiaga J. M., Salzano V., 2021, Physical Review D, 104
  * Cusin et al. (2021) Cusin G., Durrer R., Dvorkin I., 2021, Universe 2022, Vol. 8, Page 19, 8, 19
  * Dai & Venumadhav (2017) Dai L., Venumadhav T., 2017
  * Dai et al. (2016) Dai L., Venumadhav T., Sigurdson K., 2016, Physical Review D, 95
  * Dai et al. (2018) Dai L., Li S.-S., Zackay B., Mao S., Lu Y., 2018, Physical Review D, 98
  * Dai et al. (2020) Dai L., Zackay B., Venumadhav T., Roulet J., Zaldarriaga M., 2020
  * Deguchi & Watson (1986) Deguchi S., Watson W. D., 1986, Phys. Rev. D, 34, 1708
  * Diego (2020) Diego J. M., 2020, Phys. Rev. D, 101
  * Diego et al. (2021) Diego J. M., Broadhurst T., Smoot G., 2021, Physical Review D, 104
  * Ezquiaga et al. (2020) Ezquiaga J. M., Holz D. E., Hu W., Lagos M., Wald R. M., 2020, Physical Review D, 103
  * Ezquiaga et al. (2022) Ezquiaga J. M., Hu W., Lagos M., Lin M.-X., Xu F., 2022, Journal of Cosmology and Astroparticle Physics, 2022, 016
  * Fan et al. (2017) Fan X. L., Liao K., Biesiada M., Piórkowska-Kurpas A., Zhu Z. H., 2017, Physical Review Letters, 118, 091102
  * Farr et al. (2015) Farr W. M., Gair J. R., Mandel I., Cutler C., 2015, Phys. Rev. D, 91, 023005
  * Finke et al. (2021) Finke A., Foffa S., Iacovelli F., Maggiore M., Mancarella M., 2021, Physical Review D, 104
  * Goyal et al. (2020) Goyal S., Haris K., Mehta A. K., Ajith P., 2020, Physical Review D, 103
  * Goyal et al. (2021) Goyal S., D. H., Kapadia S. J., Ajith P., 2021, Physical Review D, 104
  * Goyal et al. (2023) Goyal S., Kapadia S., Cudell J.-R., Li A. K. Y., Chan J. C. L., 2023, A rapid method for preliminary identification of subthreshold strongly lensed counterparts to superthreshold gravitational-wave events (arXiv:2306.04397)
  * Hannuksela et al. (2019) Hannuksela O. A., Haris K., Ng K. K. Y., Kumar S., Mehta A. K., Keitel D., Li T. G. F., Ajith P., 2019, The Astrophysical Journal, 874, L2
  * Hannuksela et al. (2020) Hannuksela O. A., Collett T. E., Çalışkan M., Li T. G. F., 2020, Monthly Notices of the Royal Astronomical Society, 498, 3395
  * Haris et al. (2018) Haris K., Mehta A. K., Kumar S., Venumadhav T., Ajith P., 2018
  * Iacovelli et al. (2022) Iacovelli F., Finke A., Foffa S., Maggiore M., Mancarella M., 2022
  * Itoh et al. (2009) Itoh Y., Futamase T., Hattori M., 2009, Physical Review D - Particles, Fields, Gravitation and Cosmology, 80
  * Jana et al. (2022) Jana S., Kapadia S. J., Venumadhav T., Ajith P., 2022, Cosmography using strongly lensed gravitational waves from binary black holes (arXiv:2211.12212)
  * Janquart et al. (2021a) Janquart J., Hannuksela O. A., Haris K., Van Den Broeck C., 2021a, Monthly Notices of the Royal Astronomical Society, 506, 5430
  * Janquart et al. (2021b) Janquart J., Seo E., Hannuksela O. A., Li T. G. F., Broeck C. V. D., 2021b, The Astrophysical Journal Letters, 923, L1
  * Janquart et al. (2022) Janquart J., More A., Broeck C. V. D., 2022, doi:10.48550/arxiv.2205.11499
  * Kapadia et al. (2020) Kapadia S. J., et al., 2020, Classical and Quantum Gravity, 37, 045007
  * Kim et al. (2020) Kim K., Lee J., Yuen R. S. H., Hannuksela O. A., Li T. G. F., 2020, The Astrophysical Journal, 915, 119
  * Kim et al. (2022) Kim K., Lee J., Hannuksela O. A., Li T. G. F., 2022, doi:10.48550/arxiv.2206.08234
  * Klimenko et al. (2016) Klimenko S., et al., 2016, Phys. Rev. D, 93, 042004
  * Lai et al. (2018) Lai K.-H., Hannuksela O. A., Herrera-Martín A., Diego J. M., Broadhurst T., Li T. G. F., 2018, Physical Review D, 98
  * Li et al. (2023) Li A. K. Y., Lo R. K. L., Sachdev S., Chan J. C. L., Lin E. T., Li T. G. F., Weinstein A. J., 2023, Phys. Rev. D, 107, 123014
  * Liao et al. (2017) Liao K., Fan X. L., Ding X., Biesiada M., Zhu Z. H., 2017, Nature Communications, 8, 1
  * Liu et al. (2020) Liu X., Hernandez I. M., Creighton J., 2020, The Astrophysical Journal, 908, 97
  * Lo & Magaña Hernandez (2023) Lo R. K. L., Magaña Hernandez I., 2023, Phys. Rev. D, 107, 123015
  * Madau & Dickinson (2014) Madau P., Dickinson M., 2014, Annual Review of Astronomy and Astrophysics, 52, 415
  * Mao & Schneider (1997) Mao S., Schneider P., 1997, Monthly Notices of the Royal Astronomical Society, 295, 587
  * Markevitch et al. (2004) Markevitch M., Gonzalez A. H., Clowe D., Vikhlinin A., Forman W., Jones C., Murray S., Tucker W., 2004, ApJ, 606, 819
  * McIsaac et al. (2019) McIsaac C., Keitel D., Collett T., Harry I., Mozzon S., Edy O., Bacon D., 2019, Physical Review D, 102
  * Metcalf & Madau (2001) Metcalf R. B., Madau P., 2001, The Astrophysical Journal, 563, 9
  * More & More (2022) More A., More S., 2022, Monthly Notices of the Royal Astronomical Society, 515, 1044
  * Mukherjee et al. (2019a) Mukherjee S., Wandelt B. D., Silk J., 2019a, Physical Review D, 101
  * Mukherjee et al. (2019b) Mukherjee S., Wandelt B. D., Silk J., 2019b, Monthly Notices of the Royal Astronomical Society, 494, 1956
  * Mukherjee et al. (2020) Mukherjee S., Broadhurst T., Diego J. M., Silk J., Smoot G. F., 2020, Monthly Notices of the Royal Astronomical Society, 501, 2451
  * Nakamura & Deguchi (1999) Nakamura T. T., Deguchi S., 1999, Progress of Theoretical Physics Supplement, 133, 137
  * Narayan & Bartelmann (1996) Narayan R., Bartelmann M., 1996, doi:10.48550/arxiv.astro-ph/9606001
  * Ng et al. (2017) Ng K. K. Y., Wong K. W. K., Li T. G. F., Broadhurst T., 2017, Physical Review D, 97
  * Ng et al. (2018) Ng K. K. Y., Wong K. W. K., Broadhurst T., Li T. G. F., 2018, Phys. Rev. D, 97
  * Oguri (2018) Oguri M., 2018, Monthly Notices of the Royal Astronomical Society, 480, 3842
  * Oguri & Takahashi (2020) Oguri M., Takahashi R., 2020, The Astrophysical Journal, 901, 58
  * Oguri & Takahashi (2022) Oguri M., Takahashi R., 2022, Phys. Rev. D, 106, 043532
  * Ohanian (1974) Ohanian H. C., 1974, International Journal of Theoretical Physics, 9, 425
  * Pang et al. (2020) Pang P. T. H., Hannuksela O. A., Dietrich T., Pagano G., Harry I. W., 2020, Monthly Notices of the Royal Astronomical Society, 495, 3740
  * Rinaldi & Del Pozzo (2022a) Rinaldi S., Del Pozzo W., 2022a, MNRAS, 509, 5454
  * Rinaldi & Del Pozzo (2022b) Rinaldi S., Del Pozzo W., 2022b, MNRAS, 517, L5
  * Sachdev et al. (2019) Sachdev S., et al., 2019, arXiv e-prints
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Schneider et al. (1992) Schneider P., Ehlers J., Falco E. E., 1992, Gravitational Lenses. Astronomy and Astrophysics Library, Springer Berlin Heidelberg, Berlin, Heidelberg, doi:10.1007/978-3-662-03758-4, http://link.springer.com/10.1007/978-3-662-03758-4
  * Seo et al. (2021) Seo E., Hannuksela O. A., Li T. G. F., 2021
  * Sereno et al. (2010) Sereno M., Sesana A., Bleuler A., Jetzer P., Volonteri M., Begelman M. C., 2010, Physical Review Letters, 105, 251101
  * Singh et al. (2018) Singh A. J., Li I. S. C., Hannuksela O. A., Li T. G. F., Kim K., 2018, American Journal of Undergraduate Research, 16, 5
  * Singh et al. (2023) Singh M. K., Kapadia S. J., Basak S., Ajith P., Tendulkar S. P., 2023, Déjà-vu et Déjà-entendu: Associating fast radio bursts with compact binary mergers via gravitational lensing (arXiv:2304.02879)
  * Smith et al. (2022) Smith G. P., et al., 2022
  * Takahashi & Nakamura (2003) Takahashi R., Nakamura T., 2003, The Astrophysical Journal, 595, 1039
  * Tambalo et al. (2022) Tambalo G., Zumalacárregui M., Dai L., Cheung M. H.-Y., 2022, doi:10.48550/arxiv.2210.05658
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration et al., 2021, The Astrophysical Journal, 923, 14
  * Treu & Marshall (2016) Treu T., Marshall P. J., 2016, The Astronomy and Astrophysics Review, 24
  * Usman et al. (2016) Usman S. A., et al., 2016, Class. Quant. Grav., 33, 215004
  * Veske et al. (2021) Veske D., Bartos I., Márka Z., Márka S., 2021, ApJ, 922, 258
  * Vijaykumar et al. (2022a) Vijaykumar A., Mehta A. K., Ganguly A., 2022a, Detection and parameter estimation challenges of Type-II lensed binary black hole signals (arXiv:2202.06334)
  * Vijaykumar et al. (2022b) Vijaykumar A., Mehta A. K., Ganguly A., 2022b
  * Vitale et al. (2017) Vitale S., Lynch R., Raymond V., Sturani R., Veitch J., Graff P., 2017, Phys. Rev. D, 95
  * Wang et al. (1996) Wang Y., Stebbins A., Turner E. L., 1996, Phys. Rev. Lett., 77, 2875
  * Wang et al. (2021a) Wang Y., Lo R. K. L., Li A. K. Y., Chen Y., 2021a, Physical Review D, 103
  * Wang et al. (2021b) Wang J.-S., Herrera-Martín A., Hu Y.-M., 2021b, Physical Review D, 104
  * Wempe et al. (2022) Wempe E., Koopmans L. V. E., Wierda A. R. A. C., Hannuksela O. A., Broeck C. v. d., 2022
  * Wierda et al. (2021) Wierda A. R. A. C., Wempe E., Hannuksela O. A., Koopmans L. V. E., Broeck C. V. D., 2021, The Astrophysical Journal, 921, 154
  * Woosley (2017) Woosley S. E., 2017, ApJ, 836, 244
  * Xu et al. (2021) Xu F., Ezquiaga J. M., Holz D. E., 2021, The Astrophysical Journal, 929, 9
  * Yeung et al. (2021) Yeung S. M. C., Cheung M. H. Y., Gais J. A. J., Hannuksela O. A., Li T. G. F., 2021
  * Çalışkan et al. (2022) Çalışkan M., María Ezquiaga J., Hannuksela O. A., Holz D. E., 2022, arXiv e-prints, p. arXiv:2201.04619
  * Çalışkan et al. (2023) Çalışkan M., Ezquiaga J. M., Hannuksela O. A., Holz D. E., 2023, Phys. Rev. D, 107, 063023