# Using overlap of sky localization probability maps for filtering potentially
lensed pairs of gravitational-wave signals††thanks: Shared first authorship:
Henry W. Y. Wong, Lok W. L. Chan

Henry W. Y. Wong,1 Lok W. L. Chan,1 Isaac C. F. Wong,1 Rico K. L. Lo,2 and
Tjonnie G. F. Li1,3,4
1Department of Physics, The Chinese University of Hong Kong, Shatin, N. T.,
Hong Kong
2LIGO Laboratory, California Institute of Technology, Pasadena, CA 91125, US
3Institute for Theoretical Physics, KU Leuven, Celestijnenlaan 200D, B-3001
Leuven, Belgium
4Department of Electrical Engineering (ESAT), KU Leuven, Kasteelpark Arenberg
10, B-3001 Leuven, Belgium Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

When a gravitational wave passes by a massive object, its trajectory gets
curved which can result in multiple images separated by a time delay, a
phenomenon known as strong lensing. Although gravitational wave lensing has
not been observed yet, it is currently predicted that it will be possible to
detect with third-generation detectors. Moreover, the future detectors are
expected to significantly increase the gravitational-wave detection rate up to
$O(10^{5}-10^{6})$ $\mathrm{yr}^{-1}$ which could also increase the chance of
observing gravitational wave lensing. Therefore, it is necessary to develop
accurate and efficient tools to filter potentially lensed gravitational-wave
signal pairs. In the study, we propose three statistics to analyze the overlap
of two sky localization probability maps, which could provide a quick
preliminary lensing analysis of two gravitational-wave signals within 3
seconds. We simulate 200 lensed pairs of gravitational-wave signals for 5
signal-to-noise ratios to investigate the performance of the search
methodology. By setting up a threshold with a specific false positive rate
$\mathrm{FPR}=10^{-2}$ for the three overlap statistics, we conclude that it
is possible to filter out more than $99\%$ of the non-lensed event pairs while
keeping all lensed event pairs.

###### keywords:

Gravitational wave – Gravitational lensing – Skymap overlap

††pubyear: 2021††pagerange: Using overlap of sky localization probability maps
for filtering potentially lensed pairs of gravitational-wave signals††thanks:
Shared first authorship: Henry W. Y. Wong, Lok W. L. Chan–Using overlap of sky
localization probability maps for filtering potentially lensed pairs of
gravitational-wave signals††thanks: Shared first authorship: Henry W. Y. Wong,
Lok W. L. Chan

##  1 Introduction

When a gravitational wave (GW) passes by a massive object, such as a black
hole or a galaxy cluster, its path is bent towards the massive object
(Bartelmann, 2010), a phenomenon known as gravitational lensing. As a result,
the GW can take different paths to reach the observer and form multiple images
with different arrival times (Ohanian, 1974; Deguchi & Watson, 1986; Wang et
al., 1996; Nakamura, 1998; Takahashi & Nakamura, 2003). Consequently, it is
predicted that a number of GW signals with different amplitudes and arrival
times originating from the same source could be observed.

Since the detection of the first GW signal in 2015 (Abbott et al., 2016), 93
confident GW events have been detected by the Advanced LIGO (Aasi et al.,
2015) and the Advanced Virgo (Acernese et al., 2014a) in the first (O1), the
second (O2), and the third observing run (O3) (Abbott et al., 2019, 2021b,
2021a; The LIGO Scientific Collaboration et al., 2021a). Although there is no
compelling evidence in lensing of the observed GW events in O3a (The LIGO
Scientific Collaboration et al., 2021b), theoretical predictions suggest
roughly one strong GW lensing event per year at the Advanced LIGO-Advanced
Virgo-KAGRA design sensitivity (Wierda et al., 2021; Xu et al., 2021).
Therefore, it is expected that the first lensed GW event pair will be detected
in the coming years (Ng et al., 2018; Wierda et al., 2021).

There are different methods to detect strongly lensed GW event pairs (The LIGO
Scientific Collaboration et al., 2021b). One method (Haris et al., 2018;
Hannuksela et al., 2019) uses the overlap of the posterior distribution of the
source parameters to check for the consistency of the source properties
between events. The overlap of the two posterior distributions is expected to
be high when the two events come from the same source. Another method (Liu et
al., 2021; Lo & Hernandez, 2021) performs a full Bayesian analysis, namely the
joint-parameter estimation (joint-PE), on the detected events.

Joint-PE provides a solid statistical framework to calculate the probability
of a GW event pair being lensed images from the same source. The current
number of confident GW events is 93 (Abbott et al., 2019, 2021b, 2021a; The
LIGO Scientific Collaboration et al., 2021a), which gives 4278 pairs of GW
events. However, performing joint-PE on all event pairs is very time-
consuming. Moreover, for a proposed Einstein Telescope (ET) survey, the
predicted event rate of binary neutron star (BNS) mergers is as high as
$O(10^{5}-10^{6})$ $\mathrm{yr}^{-1}$ and $O(10^{4})$ $\mathrm{yr}^{-1}$ for
NS-BH mergers (Amaro-Seoane et al., 2009). This means that the detector
sensitivity and the number of GW signals in the future will be far higher and
will further increase the time for doing joint-PE on all possible event pairs.
For the purpose of sustainability, it is necessary to establish a low latency
filter to screen out event pairs that are definitely not lensed.

There are several fast algorithms for rapid identification of strongly lensed
GW events developed to address the sustainability issue. One of the algorithms
is aimed to search for strongly lensed multiple GW images - from the fact that
the signals come from the same source, it is possible to split the joint PE
into two easier PEs, that we can use the result from one event to perform the
inference of the another event (Janquart et al., 2021). Such methodology
enables relatively quick multiple-image analyses with high accuracy. Another
approach proposed by Goyal et at. uses a machine learning model to predict the
probability of lensing in event pairs from the time-frequency maps and
posterior distributions of GW events (Goyal et al., 2021). The computational
time of the lensing probability is estimated to be around 3 seconds with most
of the time spent on loading necessary files. Both of the mentioned works have
successfully reduced the computational expense of filtering potentially lensed
GW event pairs.

To establish a low latency filter, we propose the use of sky localization
probability maps (skymaps) to analyze the possibility of an event pair coming
from the same source. The deflection angle is usually a few hundred arc
seconds for black holes and galaxy clusters (Schneider et al., 1992), which is
negligibly small as it is far less than the uncertainty of the localization of
the source (Aasi et al., 2016). Therefore, the skymaps of the images from the
same source are expected to be highly overlapped. In other words, event pairs
with non-overlapping skymaps are definitely not lensed pairs. We use this fact
to analyze the possibility of whether two GW signals are a lensed pair and
filter potentially lensed pairs by comparing their skymaps. With the use of
overlap, we can reject a large fraction of event pairs that are definitely not
image pairs. However, one must note that having a high overlap in the skymaps
does not necessarily imply the events are a lensed pair. In contrast to the
machine learning approach in Goyal et al. (2021), we investigate using
measures of the overlap of skymaps without involving sophisticated learning
models, which can be easily incorporated into search pipelines such as the
subthreshold search (Li et al., 2021) to send alerts for performing follow-up
analyses. The goal of the low latency filter is to reduce time cost by rapidly
filtering out most of the non-lensed GW event pairs for further lensing
analysis. Therefore, we would like to quantitatively investigate the
efficiency of using skymap overlap to filter the non-lensed pairs while not
rejecting the lensed pairs. With a low latency skymap generated by a rapid sky
localization algorithm such as BAYESTAR (Singer & Price, 2016), a quick
preliminary lensing analysis on two GW signals can be done within a minute.

In this paper, we present three overlap statistics: normalized posterior
overlap, 90% credible region overlap and cross-HPD statistic, to
quantitatively investigate the efficiency of using skymap overlap to filter
the non-lensed pairs while not rejecting the lensed pairs. The details of the
statistics are described in Sec. 3. Using simulated GW events, we study the
distribution of skymap overlap of lensed GW event pairs under different
signal-to-noise ratios (SNRs). In order to filter potentially lensed event
pairs, a reasonable threshold of skymap overlap that can filter out non-lensed
event pairs while keeping all of the lensed event pairs should be set with a
specific false positive rate (FPR). Based on the results, we compare the
effectiveness of the three statistics and discuss the feasibility of using
skymap overlap to filter potentially lensed pairs of GW signals.

This paper is organized as follows: In Sec. 2 and 3, we establish the
methodology of the search by introducing posterior distribution and the three
statistics mentioned. In Sec. 4, we investigate the feasibility of using
skymap overlap for filtering potentially lensed pairs of GW signals with the
receiver operating characteristic (ROC) curves. In Sec. 5, we conclude our
findings on the feasibility of skymap overlap in this paper and its potential
work in the future.

##  2 Skymap

###  2.1 Posterior Distribution

To generate a skymap, we use the observed data to estimate the source location
of a GW projected on the celestial sphere under a fixed waveform model. This
process relies on parameter estimation, which allows us to quantify the
uncertainty of the GW location on the celestial sphere. We first define a
hypothesis $H$ that the GW signal fits the fixed waveform model. Then for each
GW source parameter, we set up a prior probability distribution $p(\theta|H)$,
which represents the a priori knowledge of the model parameters before
observing the data. Using Bayes’ theorem, the prior probability distribution
$p(\theta|H)$ and the data $d$ of the GW signals are used to compute the
posterior distribution $p(\theta|d,H)$ of the waveform parameters (Veitch et
al., 2015) defined as

$p(\theta|d,H)\equiv\frac{p(\theta|H)p(d|\theta,H)}{p(d|H)}\,.$ (1)

###  2.2 Sky Localization Probability Map

A waveform model typically has many parameters that we can denote as
$\theta=\\{\theta_{1},\theta_{2},...,\theta_{N}\\}$. The posterior
distribution $p(\theta|\textit{d,H})$ represents the joint posterior
distribution of all parameters. To determine the posterior distribution of a
specific parameter, we marginalize the nuisance parameters, where the
posterior distribution of a specific parameter $p(\theta_{1}|d,H)$ is obtained
by

$p(\theta_{1}|d,H)\equiv\int\mathrm{d}\theta_{2}...\mathrm{d}\theta_{N}p(\theta|d,H)\,.$
(2)

After obtaining the posterior distribution of the right ascension and the
declination of the source location, which is referred as skymap in this paper,
we compute the three overlap statistics of a skymap pair described in the
following section.

##  3 Overlap Statistics

###  3.1 Posterior Overlap

After performing the parameter estimation on two GW events separately, we
obtain two posterior distributions $p(\hat{\Omega})$ and $q(\hat{\Omega})$ of
their skymap coordinate $\hat{\Omega}$. We determine whether a GW signal pair
is potentially lensed by the inner product of $p(\hat{\Omega})$ and
$q(\hat{\Omega})$, which quantifies the overlap between two probability
densities. The posterior overlap $\mathcal{D}_{\rm PO}$ between two skymaps is
defined similarly to Haris et al. (2018) as

$\mathcal{D}_{\rm PO}\equiv\int
p(\hat{\Omega}^{\prime})q(\hat{\Omega}^{\prime})\mathrm{d}\Omega^{\prime}$ (3)

which is not normalized. In order for the posterior overlap of different GW
signal pairs to be comparable to each other, we normalize the posterior
overlap with skymap of the two GW events, which is defined as

$\mathcal{D}_{\rm NPO}\equiv\frac{\int
p(\hat{\Omega}^{\prime})q(\hat{\Omega}^{\prime})\mathrm{d}\Omega^{\prime}}{\sqrt{\int
p(\hat{\Omega}^{\prime})p(\hat{\Omega}^{\prime})\mathrm{d}\Omega^{\prime}}\sqrt{\int
q(\hat{\Omega}^{\prime})q(\hat{\Omega}^{\prime})\mathrm{d}\Omega^{\prime}}}\in[0,1]\,.$
(4)

###  3.2 90% Credible Region Overlap

For a GW signal with an extremely large signal-to-noise ratio (SNR), the
posterior distribution of the source location could be narrowly peaked. This
may lead to an extremely small value of posterior overlap and normalized
posterior overlap with another GW signal. As a result, the two GW signals may
be recognized as a non-lensed pair and filtered out. Therefore, we propose to
use a 90% credible region overlap.

Instead of finding the overlap probability of the source position on the sky,
we directly compute the overlapping area of the 90% credible regions. The
overlapped area is further normalized by the minimum area between two 90%
credible regions. The minimum value is used to keep a high 90% credible region
overlap when the difference between two event SNRs is high for lensed pairs.
Therefore, it reduces the chance of lensed pairs being filtered out. The 90%
credible region overlap $\mathcal{D}_{\rm 90\%CR}(p,q)$ of two GW signals is
computed by masking their skymaps, which is defined as

$\mathcal{D}_{\rm 90\%CR}(p,q)=\frac{\displaystyle\int\mathbbm{1}_{\rm
90\%CR}\left[p(\hat{\hat{\Omega}}^{\prime})\right]\mathbbm{1}_{\rm
90\%CR}\left[q(\hat{\Omega}^{\prime})\right]\mathrm{d}\Omega^{\prime}}{{\rm
min}(\displaystyle\int\mathbbm{1}_{\rm
90\%CR}\left[p(\hat{\Omega}^{\prime})\right]\mathrm{d}\Omega^{\prime},\displaystyle\int\mathbbm{1}_{\rm
90\%CR}\left[q(\hat{\Omega}^{\prime})\right]\mathrm{d}\Omega^{\prime})},$ (5)

where $\mathbbm{1}_{\rm 90\%CR}\left[p(\hat{\Omega})\right]$ and
$\mathbbm{1}_{\rm 90\%CR}\left[q(\hat{\Omega})\right]$ are the indicator
function. This indicator function is 1 when the sky location is in the 90%
credible region of the corresponding skymap, and is 0 otherwise.

###  3.3 Cross-HPD Statistic

We introduce another statistic called the cross-highest posterior density
(HPD). The purpose of this statistic is similar to that of the 90% credible
region overlap. We use the maximum a posteriori (MAP), which is the mode of
the posterior distribution, as an alternative of the 90% confidence level to
decide the credible region used for computing the overlap of the skymap
between two GW signals. We first compute the probability of the MAP of one
skymap on another skymap, which is referred as the search probability
$\mathcal{H}$. It is defined as

$\mathcal{H}(p,q)=\int\mathbbm{1}\left[q(\hat{\Omega}^{\prime})>q(\hat{\Omega}_{p\textrm{,
MAP}})\right]q(\hat{\Omega}^{\prime})d\hat{\Omega}^{\prime},$ (6)

where $\mathbbm{1}$ is the indicator function, $q(\hat{\Omega})$ is the
posterior distribution of sky position of an event, and
$\hat{\Omega}_{p\textrm{,MAP}}$ is the MAP estimate of $\hat{\Omega}_{p}$,
also known as the coordinates of the point at highest probability density.

$\mathcal{H}(p,q)$ quantifies the overlap of the skymap between the event
pairs. The lower the value of $\mathcal{H}(p,q)$, the higher the overlap
between the skymaps. The cross-HPD statistic $\mathcal{D}_{\rm HPD}(p,q)$ is
defined as

$\mathcal{D}_{\rm
HPD}(p,q)\equiv\max\left\\{1-\mathcal{H}(p,q),1-\mathcal{H}(q,p)\right\\},$
(7)

which is symmetric in $p$ and $q$.

##  4 Results & Discussion

We simulate 200 lensed pairs of GW signals for 5 different SNRs to investigate
the performance of the search methodology. We first generate the simulated
data of 200 lensed pairs of GW signals with PyCBC (Nitz et al., 2021). Primary
mass follows a power-law profile with slope -2.35, while the mass ratio and
sky position are uniformly distributed. With the simulated data, we re-scale
all the events to 5 different SNRs by adjusting the luminosity distance of
each event. In total, there are 1,000 event pairs in the injection campaign.
For details, see Sec. 4.1. We use a Bayesian inference library Bilby (Ashton
et al., 2019; Romero-Shaw et al., 2020) to perform parameter estimation on the
GW events. After obtaining the posterior distribution of the sky location of
each event, the skymaps of all events are generated using ligo.skymap (Singer
et al., 2016). We pair up the skymaps for every possible combination and
separate them into lensed-pair group and non-lensed-pair group. The three
overlap statistics of all event pairs are then computed. Details of the
overlap statistics can be found in Sec. 3. Finally, we plot the receiver
operating characteristic (ROC) curve to compare the performance of the overlap
statistics. Details of the ROC curve can be found in Sec. 4.3.

###  4.1 Injection Campaign

Number of lensed event pairs | 200
---|---
Waveform approximant |  IMRPhenomPv2 (Husa et al., 2016; Bohé et al., 2017)
Spin | No
Re-scaled SNR | 8, 18.5, 29, 39.5, 50
Detector network | HLV at design sensitivity (Acernese et al., 2014b; Buikema et al., 2020)
Primary mass distribution | Power-law profile with slope -2.35
Mass ratio distribution | Uniform
Sky location distribution | Uniform
Table 1: The parameters of the waveform simulated and the set-up of the
injection campaign. 200 event pairs are generated. For simplicity, the spins
of all events are set to zero. All events are rescaled to 5 different SNRs by
adjusting the luminosity distance of each event. Finally, there are 2,000
events in total.

To study the distribution of skymap overlap of lensed-pair candidates under
different SNRs and the efficiency of overlap statistics on filtering
potentially lensed pairs of GW signals, we set up an injection campaign with
200 lensed GW event pairs. The localization uncertainty of a GW signal depends
on its $\mathrm{SNR}$, therefore we are interested in the effect of different
$\mathrm{SNRs}$ on the skymaps of lensed GW signal pairs. The uncertainty of
the localization is smaller when the $\mathrm{SNR}$ is higher. Therefore, we
expect the overlap to be generally lower for GW event pairs with a higher
$\mathrm{SNR}$ than those with a lower $\mathrm{SNR}$. We use 5 different
$\mathrm{SNR}$s to study the difference in the filtering threshold for
different SNR pairs by skymap overlap.

The largest $\mathrm{SNR}$ of the event pairs is chosen to be 50, it is a
sufficiently high $\mathrm{SNR}$ upper limit due to the rarity of detecting a
signal with such $\mathrm{SNR}$. After setting up the upper limit and the
lower limit of the event $\mathrm{SNR}$, we uniformly divide the range between
8 and 50 to obtain 5 different $\mathrm{SNR}$s and rescale the luminosity
distances to attain the desired SNR.

After performing the parameter estimation, we generate the skymaps for all
events. The prior of chirp mass is set to be uniform from 4 $M_{\odot}$ to 150
$M_{\odot}$ and the two component masses are constrained between 5 $M_{\odot}$
to 150 $M_{\odot}$. The prior of geocent time is chosen to be a delta function
that peaks at the injected parameter. The sampler used for the parameter
estimation is dynesty (Speagle, 2020). Finally, the three overlap statistics
of each event pair are computed. The lensed event pairs consist of 1000
samples which are 200 lensed events pairs for each of the 5 SNRs. For the non-
lensed event pairs, the overlap of all the possible event pairs from the
original 2,000 events is computed, which results in about 700,000 samples. By
plotting the ROC curve of each statistic, we study the performance of the
three statistics on filtering lensed candidate event pairs.

###  4.2 Coverage of Samples

SNR | Mean area of 90% credible region ($\text{deg}^{2}$)
---|---
8 | 186.8
18.5 | 34.8
29 | 12.0
39.5 | 6.4
50 | 4.1
Table 2: Mean of the 90% credible region area of the skymaps of 200 events
simulated using HLV network for 5 SNRs. Details of the injection campaign
refer to Table 1. Figure 1: The contours of 90% credible region of 200 skymaps
with $\mathrm{SNR}=8$ are shown in the figure. Details of the injection
campaign are referred to Table 1. The figure shows the sky coverage of the
injections with $\mathrm{SNR}=8$.

Table 2 shows the mean area of 90% credible region of the skymaps in the
injection campaign. The size of the 90% credible region relates to the
uncertainty of the localization. For GW signals with a smaller $\mathrm{SNR}$,
the uncertainty of localization is larger. Therefore, the mean area of 90%
credible region for events with $\mathrm{SNR}=8$ is significantly larger than
that of the events with $\mathrm{SNR}=50$.

Figure 1 illustrates the 90% credible region of events with $\mathrm{SNR}=8$.
The contours of their 90% credible regions are demonstrated in different
colors for distinguishing different contours. The 200 events come from
different parts of the sky and their 90% credible regions covered about half
of the sky. Note that there are overlaps of 90% credible regions of non-image
pairs. Those non-image event pairs with relatively high overlap values are the
false positives. The coverage ensures that there is a sufficient number of
lensed event pairs in the injection campaign to study the distribution of
skymap overlap for non-image pairs.

###  4.3 Performance of Statistics

With the overlap statistics of all the event pairs, the performance of the
statistics is shown in a ROC curve. In a ROC curve, which shows the
performance of a classification model at all classification thresholds by
plotting the true positive rate (TPR) against the false positive rate (FPR).

The y-axis of the ROC curves is TPR, which is defined as

$\mathrm{TPR}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}},$ (8)

where TP is the number of lensed event pairs which overlap value is higher
than the filtering threshold, while FN is the number of lensed event pairs
which overlap values are smaller than the filtering threshold.

The x-axis of the ROC curves is FPR, which is defined as

$\mathrm{FPR}=\frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{TN}},$ (9)

where FP is the number of non-lensed event pairs which overlap value is higher
than the filtering threshold, while TN is the number of non-lensed event pairs
which overlap value is smaller than the filtering threshold.

Given a specific FPR or a specific filtering threshold, a higher TPR indicates
a better performance of the statistic because less lensed event pairs are
recognized as non-lensed event pairs by the statistic. We aim to set up a
reasonable threshold with $\mathrm{TPR}\approx 1$ to filter non-lensed event
pairs while keeping most of the lensed event pairs for further analysis. Ten
ROC curves are shown in Figure 2 and Figure 3. The five ROC curves in Figure 2
indicate the performance of the three statistics for event pairs with one of
the event SNRs is small ($\mathrm{SNR}=8$), while the five ROC curves in
Figure 3 indicate the performance of the three statistics for event pairs with
one of the event SNRs is large ($\mathrm{SNR}=50$).

(a) $\mathrm{SNR}=8$ vs $\mathrm{SNR}=8$

(b) $\mathrm{SNR}=8$ vs $\mathrm{SNR}=18.5$

(c) $\mathrm{SNR}=8$ vs $\mathrm{SNR}=29$

(d) $\mathrm{SNR}=8$ vs $\mathrm{SNR}=39.5$

(e) $\mathrm{SNR}=8$ vs $\mathrm{SNR}=50$

Figure 2: Receiver operating characteristic (ROC) curves of event pairs with
one event $\mathrm{SNR}=8$. The black solid line $x=y$ represents the ROC
curve of a random classifier. ROC curves indicate the performance of the
overlap statistics. The classifier is better if it can achieve a higher true
positive rate given a false positive rate. The three statistics do not show a
large difference in the region of FPR between $10^{-2}$ and 1.

Figure 2 shows the ROC curves of the three overlap statistics for event pairs
with one event SNR$=8$ while another is one of the 5 SNRs. In this case, the
TPR of the three statistics shows nearly no difference when the FPR is fixed
to a value between $10^{-2}$ and 1. The difference in the performance of the
statistics is noticeable only when we fix FPR to a value smaller than around
$10^{-2}$.

For a relatively small FPR value ($\mathrm{{FPR}<10^{-2}}$), 90% credible
region overlap shows the worst overall performance as it shows the most
significant drop at a lower FPR. However, all three statistics show a
significant drop when the FPR is lower than $10^{-2}$. Within this range, the
losses of lensed event pairs are significant for all three statistics.
Therefore, the threshold should be fixed higher than $10^{-2}$ if we aim to
keep most of the lensed event pairs. Thus, the performance of the three
statistics is similar for a reasonable choice of threshold that can filter out
non-lensed event pairs while keeping most of the lensed event pairs. The
straight line of the TPR drop for 90% credible region overlap in
$\mathrm{FPR}\lessapprox 10^{-3}$ is due to a lack of samples. The data points
of the ROC curve are generated by using a threshold to separate the input data
into FP, TP, FN, and TN.

(a) SNR=8 vs SNR=50

(b) SNR=18.5 vs SNR=50

(c) SNR=29 vs SNR=50

(d) SNR=39.5 vs SNR=50

(e) SNR=50 vs SNR=50

Figure 3: Receiver operating characteristic (ROC) curves of event pairs with
one event $\mathrm{SNR}=50$. The black solid line $x=y$ represents the ROC
curve of a random classifier. The classifier is better if it can achieve a
higher true positive rate given a false positive rate. The three statistics do
not show a large difference in the region of FPR between $10^{-2}$ and 1.

Figure 3 shows the ROC curves of the three overlap statistics for event pairs
with one event $\mathrm{{SNR}=50}$ while another is one of the 5 SNRs. In this
case, the TPR of the three statistics shows nearly no difference when the FPR
is fixed to a value between $10^{-2}$ and 1. The difference in the performance
of the statistics is noticeable only when we fix FPR to a value smaller than
around $10^{-3}$.

For a relatively small FPR value ($\mathrm{{FPR}<10^{-3}}$), 90% credible
region overlap shows the worst overall performance. However, again, choosing a
threshold with FPR lower than around $10^{-3}$ will lead to a significant loss
in the number of lensed event pairs. Therefore, the filtering threshold should
be fixed in an FPR that is higher than $10^{-3}$ in order to keep all of the
lensed event pairs. Thus, the performance of the three statistics is similar
for a reasonable choice of filtering threshold that can filter out non-lensed
event pairs while keeping all of the lensed event pairs. Again, the straight
line of the TPR drop for 90% credible region overlap in
$\mathrm{{FPR}\lessapprox 10^{-3}}$ is due to the lack of samples.

For the purpose of filtering potentially lensed candidate event pairs, the
three skymap overlap statistics demonstrated a similar performance while
fixing FPR between $10^{-2}$ and 1. In the mock data study, we can set up a
specific $\mathrm{FPR}=10^{-2}$ for event pairs in all SNR to filter out 99%
of the non-lensed candidate pairs, which demonstrated that a large fraction of
non-image pairs could be rapidly filtered away using the skymap overlap.

##  5 Conclusions

In this paper, we have demonstrated our proposed search methodology using
skymap overlap to filter potentially lensed pairs of GW signals. By setting up
a threshold with a specific $\mathrm{FPR}=10^{-2}$ for the three overlap
statistics, we can filter out more than 99% of the non-lensed event pairs
while keeping most of the lensed event pairs, which shows the feasibility of
using skymap overlap to filter potentially lensed event pairs of GW signals.

Since the computation of the three statistics for each event pair takes less
than 3 seconds, a new trigger could be constructed to compute the skymap
overlap of all existing events within a short period of time using the low
latency skymaps generated by BAYESTAR (Singer & Price, 2016) or other rapid
sky localization algorithms.

## Acknowledgements

ICFW and TGFL are partially supported by grants from the Research Grants
Council of the Hong Kong (Project No. 24304317 and 14306419) and Research
Committee of the Chinese University of Hong Kong. RKLL and TGFL would also
like to gratefully acknowledge the support from the Croucher Foundation in
Hong Kong. The authors acknowledge the use of the IUCAA LDG cluster Sarathi
for the computational/numerical work. The authors are also grateful for the
computational resources provided by the Chinese University of Hong Kong.

## References

  * Aasi et al. (2015) Aasi J., et al., 2015, Classical and quantum gravity, 32, 074001
  * Aasi et al. (2016) Aasi J., et al., 2016, Living Reviews in Relativity, 19
  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2019) Abbott B. P., et al., 2019, Phys. Rev. X, 9, 031040
  * Abbott et al. (2021a) Abbott R., et al., 2021a, arXiv preprint arXiv:2108.01045
  * Abbott et al. (2021b) Abbott R., et al., 2021b, Physical Review X, 11
  * Acernese et al. (2014a) Acernese F., et al., 2014a, ] 10.1088/0264-9381/32/2/024001, 32, 024001
  * Acernese et al. (2014b) Acernese F., et al., 2014b, Classical and Quantum Gravity, 32, 024001
  * Amaro-Seoane et al. (2009) Amaro-Seoane P., et al., 2009, Codified Document ET-030-09
  * Ashton et al. (2019) Ashton G., et al., 2019, The Astrophysical Journal Supplement Series, 241, 27
  * Bartelmann (2010) Bartelmann M., 2010, Classical and Quantum Gravity, 27, 233001
  * Bohé et al. (2017) Bohé A., et al., 2017, Phys. Rev. D, 95, 044028
  * Buikema et al. (2020) Buikema A., et al., 2020, Phys. Rev. D, 102, 062003
  * Deguchi & Watson (1986) Deguchi S., Watson W. D., 1986, Phys. Rev. D, 34, 1708
  * Goyal et al. (2021) Goyal S., D. H., Kapadia S. J., Ajith P., 2021, Rapid Identification of Strongly Lensed Gravitational-Wave Events with Machine Learning (arXiv:2106.12466)
  * Hannuksela et al. (2019) Hannuksela O., Haris K., Ng K., Kumar S., Mehta A., Keitel D., Li T., Ajith P., 2019, The Astrophysical Journal Letters, 874, L2
  * Haris et al. (2018) Haris K., Mehta A. K., Kumar S., Venumadhav T., Ajith P., 2018, arXiv preprint arXiv:1807.07062
  * Husa et al. (2016) Husa S., Khan S., Hannam M., Pürrer M., Ohme F., Forteza X. J., Bohé A., 2016, Physical Review D, 93, 044006
  * Janquart et al. (2021) Janquart J., Hannuksela O. A., Haris K., Van Den Broeck C., 2021, Monthly Notices of the Royal Astronomical Society, 506, 5430–5438
  * Li et al. (2021) Li A. K. Y., Lo R. K. L., Sachdev S., Chan C. L., Lin E. T., Li T. G. F., Weinstein A. J., 2021, Finding diamonds in the rough: Targeted Sub-threshold Search for Strongly-lensed Gravitational-wave Events (arXiv:1904.06020)
  * Liu et al. (2021) Liu X., Hernandez I. M., Creighton J., 2021, The Astrophysical Journal, 908, 97
  * Lo & Hernandez (2021) Lo R. K., Hernandez I. M., 2021, arXiv preprint arXiv:2104.09339
  * Nakamura (1998) Nakamura T. T., 1998, Phys. Rev. Lett., 80, 1138
  * Ng et al. (2018) Ng K. K., Wong K. W., Broadhurst T., Li T. G., 2018, Physical Review D, 97, 023012
  * Nitz et al. (2021) Nitz A., et al., 2021, gwastro/pycbc:, doi:10.5281/zenodo.5347736, https://doi.org/10.5281/zenodo.5347736
  * Ohanian (1974) Ohanian H. C., 1974, International Journal of Theoretical Physics, 9, 425
  * Romero-Shaw et al. (2020) Romero-Shaw I., et al., 2020, Monthly Notices of the Royal Astronomical Society, 499, 3295
  * Schneider et al. (1992) Schneider P., Ehlers J., Falco E. E., 1992, Gravitational Lenses, doi:10.1007/978-3-662-03758-4.
  * Singer & Price (2016) Singer L. P., Price L. R., 2016, Physical Review D, 93, 024013
  * Singer et al. (2016) Singer L. P., et al., 2016, The Astrophysical Journal Supplement Series, 226, 10
  * Speagle (2020) Speagle J. S., 2020, Monthly Notices of the Royal Astronomical Society, 493, 3132
  * Takahashi & Nakamura (2003) Takahashi R., Nakamura T., 2003, The Astrophysical Journal, 595, 1039
  * The LIGO Scientific Collaboration et al. (2021a) The LIGO Scientific Collaboration et al., 2021a, GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run (arXiv:2111.03606)
  * The LIGO Scientific Collaboration et al. (2021b) The LIGO Scientific Collaboration et al., 2021b, Search for lensing signatures in the gravitational-wave observations from the first half of LIGO-Virgo’s third observing run (arXiv:2105.06384)
  * Veitch et al. (2015) Veitch J., et al., 2015, Physical Review D, 91, 042003
  * Wang et al. (1996) Wang Y., Stebbins A., Turner E. L., 1996, Phys. Rev. Lett., 77, 2875
  * Wierda et al. (2021) Wierda A. R. A. C., Wempe E., Hannuksela O. A., Koopmans L. V. E., Broeck C. V. D., 2021, Beyond the detector horizon: Forecasting gravitational-wave strong lensing (arXiv:2106.06303)
  * Xu et al. (2021) Xu F., Ezquiaga J. M., Holz D. E., 2021, Please repeat: Strong lensing of gravitational waves as a probe of compact binary and galaxy populations (arXiv:2105.14390)