of 31
D
raft version
1 D
ecember
2021
Typeset using L
A
T
E
X
twocolumn
style in AASTeX62
Search for lensing signatures in the gravitational-wave observations from the first half of LIGO-Virgo’s third observing run
T
he
LIGO S
cientific
C
ollaboration and the
V
irgo
C
ollaboration
(Dated: 13 May 2021)
ABSTRACT
We search for signatures of gravitational lensing in the gravitational-wave signals from compact binary coales-
cences detected by Advanced LIGO and Advanced Virgo during O3a, the first half of their third observing run.
We study: 1) the expected rate of lensing at current detector sensitivity and the implications of a non-observation
of strong lensing or a stochastic gravitational-wave background on the merger-rate density at high redshift; 2) how
the interpretation of individual high-mass events would change if they were found to be lensed; 3) the possibility
of multiple images due to strong lensing by galaxies or galaxy clusters; and 4) possible wave-optics e
ff
ects due to
point-mass microlenses. Several pairs of signals in the multiple-image analysis show similar parameters and, in
this sense, are nominally consistent with the strong lensing hypothesis. However, taking into account population
priors, selection e
ff
ects, and the prior odds against lensing, these events do not provide su
ffi
cient evidence for
lensing. Overall, we find no compelling evidence for lensing in the observed gravitational-wave signals from any
of these analyses.
1.
INTRODUCTION
Gravitational lensing occurs when a massive object bends
spacetime in a way that focuses light rays toward an observer
(see Bartelmann 2010 for a review). Lensing observations
are widespread in electromagnetic astrophysics and have been
used for, among other purposes, making a compelling case
for dark matter (Clowe et al. 2004; Markevitch et al. 2004),
discovering exoplanets (Bond et al. 2004), and uncovering
massive objects and structures that are too faint to be detected
directly (Coe et al. 2013).
Similarly to light, when gravitational waves (GWs) travel
near a galaxy or a galaxy cluster, their trajectories curve,
resulting in gravitational lensing (Ohanian 1974; Thorne
1982; Deguchi & Watson 1986; Wang et al. 1996; Nakamura
1998; Takahashi & Nakamura 2003). For massive lenses,
this changes the GW amplitude without a
ff
ecting the fre-
quency evolution (Wang et al. 1996; Dai & Venumadhav 2017;
Ezquiaga et al. 2021). Strong lensing, in particular, can also
produce multiple images observed at the GW detectors as re-
peated events separated by a time delay of minutes to months
for galaxies (Ng et al. 2018; Li et al. 2018; Oguri 2018), and
up to years for galaxy clusters (Smith et al. 2018, 2017, 2019;
Robertson et al. 2020; Ryczanowski et al. 2020). The detec-
tion of such strongly lensed GWs has been forecast within this
decade (Ng et al. 2018; Li et al. 2018; Oguri 2018), at design
sensitivity of Advanced LIGO and Advanced Virgo, assuming
that binary black holes (BBHs) trace the star-formation rate
density. In addition, if GWs propagate near smaller lenses
such as stars or compact objects, microlensing may induce ob-
servable beating patterns in the waveform (Deguchi & Watson
1986; Nakamura 1998; Takahashi & Nakamura 2003; Cao
et al. 2014; Jung & Shin 2019; Lai et al. 2018; Christian et al.
2018; Dai et al. 2018; Diego et al. 2019; Diego 2020; Pagano
et al. 2020; Cheung et al. 2021; Mishra et al. 2021). Indeed,
lensing can induce a plethora of e
ff
ects on GWs.
If observed, GW lensing could enable numerous scientific
pursuits, such as localization of merging black holes to sub-
arcsecond precision (Hannuksela et al. 2020), precision cos-
mography studies (Sereno et al. 2011; Liao et al. 2017; Cao
et al. 2019; Li et al. 2019b; Hannuksela et al. 2020), precise
tests of the speed of gravity (Baker & Trodden 2017; Fan et al.
2017), tests of the GW’s polarization content (Goyal et al.
2021), and detecting intermediate-mass or primordial black
holes (Lai et al. 2018; Diego 2020; Oguri & Takahashi 2020).
Here we perform a comprehensive lensing analysis of data
from the first half of the third LIGO–Virgo observing run,
called O3a for short, focusing on compact binary coalescence
(CBC) signals. We begin by outlining the expected rate of
strongly lensed events. Strong lensing is rare, but magnified
signals enable us to probe a larger comoving volume, thus
potentially giving us access to more sources (Dai et al. 2017;
Ng et al. 2018; Smith et al. 2018; Li et al. 2018; Oguri 2018;
Smith et al. 2017, 2019; Robertson et al. 2020; Ryczanowski
et al. 2020). We forecast the lensed event rates using stan-
dard lens and black hole population models (Sec. 3). These
expected rates are subject to some astrophysical uncertainty
but are vital to interpreting our search results in later sections.
The rate of lensing can also be inferred from the stochastic
GW background (SGWB) (Buscicchio et al. 2020a; Mukher-
jee et al. 2021a; Buscicchio et al. 2020b). Thus, we use the
non-observation of strong lensing and the stochastic back-
arXiv:2105.06384v3 [gr-qc] 30 Nov 2021
2
ground to constrain the BBH merger-rate density and the rate
of lensing at high redshifts.
In addition, lensing magnification biases the inferred GW
luminosity distance and source mass measurements, which
could lead to observations of apparently high-mass (or low-
mass, when de-magnified) binaries (Dai et al. 2017; Broad-
hurst et al. 2018; Oguri 2018; Hannuksela et al. 2019; Broad-
hurst et al. 2020a). Therefore, we analyze several LIGO–
Virgo detections with unusually high masses under the alterna-
tive interpretation that they are lensed signals from lower-mass
sources which have been magnified (Sec. 4).
We then move on to search for signatures of lensing-induced
multiple images, which should appear as repeated similar
signals, magnified and with waveform di
ff
erences determined
by the image type (Dai & Venumadhav 2017; Ezquiaga et al.
2021), separated in time by minutes to months (or even years).
Consequently, if an event pair is strongly lensed, we expect to
infer consistent parameters for both events (Haris et al. 2018;
Hannuksela et al. 2019).
We search for these multiple images by first comparing the
posterior overlap between pairs of events occurring during the
O3a period as reported in Abbott et al. (2021a) (Sec. 5.1). Af-
ter identifying a list of candidates from the posterior-overlap
analysis, we follow these up with more computationally ex-
pensive but more accurate joint-parameter estimation (PE)
procedures (Sec. 5.2). Next, we perform a targeted search for
previously undetected counterpart images of known events in
Sec. 5.3, images that could have fallen below the threshold of
previous wide-parameter space CBC searches (as discussed in
Li et al. 2019a; McIsaac et al. 2020; Dai et al. 2020). Finally,
we search for microlensing induced by point-mass lenses in
the intermediate and low mass range, including wave-optics
e
ff
ects (Sec. 6).
Several searches for GW lensing signatures have already
been performed in the data from the first two observing runs
O1 and O2 (Hannuksela et al. 2019; Li et al. 2019a; McIsaac
et al. 2020; Pang et al. 2020; Liu et al. 2021; Dai et al. 2020),
including strong lensing and microlensing e
ff
ects. We will dis-
cuss these previous studies in the appropriate sections. Given
the growing interest in GW lensing and the existing forecasts,
an analysis of the most recent GW observations for lensing
e
ff
ects is now timely.
Results of all analyses in this paper and associated data
products can be found in LIGO Scientific Collaboration and
Virgo Collaboration (2021). GW strain data and posterior
samples for all events from GWTC-2 are available (GWOSC
2020) from the Gravitational Wave Open Science Center (Ab-
bott et al. 2021b).
2.
DATA AND EVENTS CONSIDERED
The analyses presented here use data taken during O3a
by the Advanced LIGO (Aasi et al. 2015) and Advanced
Virgo (Acernese et al. 2015) detectors. O3a extended from
2019 April 1 to 2019 October 1. Various instrumental up-
grades have led to more sensitive data, with median binary
neutron star (BNS) inspiral ranges (Allen et al. 2012a) in-
creased by a factor of 1.64 in LIGO Hanford, 1.53 in LIGO
Livingston, and 1.73 in Virgo compared to O2 (Abbott et al.
2021a). The duty factor for at least one detector being on-
line was 97%; for any two detectors being online at the same
time it was 82%; and for all three detectors together it was
45%. Further details regarding instrument performance and
data quality for O3a are available in Abbott et al. (2021a);
Buikema et al. (2020).
The LIGO and Virgo detectors used a photon recoil based
calibration (Karki et al. 2016; Cahillane et al. 2017; Viets et al.
2018) resulting in a complex-valued, frequency-dependent
detector response with typical errors in magnitude of 7% and
4 degrees in phase (Sun et al. 2020; Acernese et al. 2018) in
the calibrated O3a strain data.
Transient noise sources, referred to as glitches, contaminate
the data and can a
ff
ect the confidence of candidate detections.
Times a
ff
ected by glitches are identified so that searches for
GW events can exclude (veto) these periods of poor data
quality (Abbott et al. 2016a, 2020a; Davis et al. 2021; Nguyen
et al. 2021; Fiori et al. 2020). In addition, several known noise
sources are subtracted from the data using information from
witness auxiliary sensors (Driggers et al. 2019; Davis et al.
2019).
Candidate events, including those reported in Abbott et al.
(2021a) and the new candidates found by the searches for sub-
threshold counterpart images in Sec. 5.3 of this paper, have
undergone a validation process to evaluate if instrumental
artifacts could a
ff
ect the analysis; this process is described in
detail in Sec. 5.5 of Davis et al. (2021). This process can also
identify data quality issues that need further mitigation for
individual events, such as subtraction of glitches (Cornish &
Littenberg 2015) and non-stationary noise couplings (Vajente
et al. 2020), before executing PE algorithms. See Table V
of Abbott et al. (2021a) for the list of events requiring such
mitigation.
The GWTC-2 catalog (Abbott et al. 2021a) contains 39
events from O3a (in addition to the 11 previous events from
O1 and O2) with a false-alarm rate (FAR) below two per
year, with an expected rate of false alarms from detector
noise less than 10% (Abbott et al. 2021a). We neglect the
potential contamination in this analysis. These events were
identified by three search pipelines: one minimally modeled
transient search
c
WB
(Klimenko et al. 2004, 2005, 2006, 2011,
2016) and the two matched-filter searches
G
st
LAL
(Sachdev
et al. 2019; Hanna et al. 2020; Messick et al. 2017) and
P
y
CBC
(Allen et al. 2012b; Allen 2005; Dal Canton et al.
2014; Usman et al. 2016; Nitz et al. 2017). Their param-
eters were estimated through Bayesian inference using the
3
LALI
nference
(Veitch et al. 2015) and
B
ilby
(Ashton et al.
2019; Smith et al. 2020; Romero-Shaw et al. 2020) pack-
ages. Both the matched-filter searches and PE use a variety
of CBC waveform models which generally combine knowl-
edge from post-Newtonian theory, the e
ff
ective-one-body for-
malism and numerical relativity (for general introductions
to these approaches, see Blanchet 2014; Damour & Nagar
2016; Palenzuela 2020; Schmidt 2020 and references therein).
The analyses in this paper rely on the same methods, and
the specific waveform models and analysis packages used are
described in each section.
Most of the 39 events from O3a are most probably
BBHs, while three (GW190425, GW190426
152155, and
GW190814) have component masses below 3
M
(Abbott
et al. 2020b, 2021a, 2020d), thus potentially containing a
neutron star. We consider these 39 events in most of the
analyses in this paper, except in the magnification analysis
(Sec. 4), which concerns only six of the more unusual events,
and the microlensing analysis (Sec. 6), which focuses on the
36 clear BBH events only.
Specifically, we use the following input data sets for each
analysis. The magnification analysis in Sec. 4 and posterior-
overlap analysis in Sec. 5.1 start from the Bayesian inference
posterior samples released with GWTC-2 (GWOSC 2020).
The joint-PE analyses in Sec. 5.2 and microlensing analysis
in Sec. 6 reanalyze the strain data in short segments around
the event times, available from the same data release, with
data selection and noise mitigation choices matching those
of the PE analyses in Abbott et al. (2021a). In addition, the
searches for sub-threshold counterpart images in Sec. 5.3
cover the whole O3a strain data set, using the same data
quality veto choices as in Abbott et al. (2021a) but a strain
data set consistent with the PE analyses: the final calibration
version of LIGO data (Sun et al. 2020) with additional noise
subtraction (Vajente et al. 2020).
3.
LENSING STATISTICS
In this section, we first forecast the number of detectable
strongly lensed events (Sec. 3.1). Then, we infer upper limits
on the rate of strongly lensed events using two di
ff
erent meth-
ods; the first uses only the non-detection of resolvable strongly
lensed BBH events (Sec. 3.2), while the second leverages addi-
tionally the non-observation of the SGWB (Sec. 3.3) (Callister
et al. 2020; Abbott et al. 2021c). Since the background would
originate from higher redshifts, this second method comple-
ments the first method.
Throughout this section, we model the mass distribution of
BBHs following the results for the
P
ower
L
aw +
P
eak model
of Abbott et al. (2021d). We consider two distinct models
of the BBH merger rate density. Model A brackets most of
the population synthesis results (Eldridge et al. 2019; Neijs-
sel et al. 2019; Boco et al. 2019; Santoliquido et al. 2021)
corresponding to Population I and II stars while Model B
assumes the Madau–Dickinson ansatz (Madau & Dickinson
2014) where the rate peaks at a particular redshift. For con-
sistency with previous analyses (e.g., Abbott et al. 2021c),
we take the Hubble constant from Planck 2015 observations
to be
H
0
=
67
.
9
km s
1
Mpc
1
(Ade et al. 2016). Detailed
discussion on both models and their respective parametriza-
tion is given in Appendix A. The obtained rates are subject to
uncertainty because of their dependence on the merger rate
density, which is model-dependent and only partially con-
strained. They are nevertheless vital to interpreting our search
results in later sections (see Sec. 5).
3.1.
Strong lensing rate
We predict the rate of lensing using the standard methods
outlined in the literature (Ng et al. 2018; Li et al. 2018; Oguri
2018; Xu et al. 2021; Mukherjee et al. 2021b; Wierda et al.
2021), at galaxy and galaxy-cluster lens mass scales. To
model the lens population, we need to choose a density profile
and a mass function. We adopt the singular isothermal sphere
(SIS) density profile for both galaxies and clusters. Moreover,
we use the velocity dispersion function (VDF) from the Sloan
Digital Sky Survey (SDSS; Choi et al. 2007) for galaxies and
the halo mass function from Tinker et al. (2008) for clusters
which have been used in other lensing studies as well (e.g.,
Oguri & Marshall 2010; Robertson et al. 2020). The SIS
profile can well describe galaxies. However, the mass distri-
bution of clusters tends to be more complicated. Nevertheless,
Robertson et al. (2020) have demonstrated that the SIS model
can reproduce the lensing rate predictions from a study of nu-
merically simulated cluster lenses. Thus, we adopt the same
model. Under the SIS model, we obtain two images with
di
ff
erent magnifications and arrival times. The rate of strong
lensing is
R
lens
=
d
N
(
M
h
,
z
l
)
d
M
h
d
D
c
d
z
l
R
m
(
z
m
)
1
+
z
m
d
V
c
d
z
m
σ
(
M
h
,
z
l
,
z
m
,ρ,ρ
c
)
×
p
(
ρ
|
z
m
) d
ρ
d
z
m
d
z
l
d
M
h
,
(1)
where d
N
(
M
h
,
z
l
)
/
d
M
h
is the di
ff
erential comoving number
density of lensing halos in a halo mass shell at lens redshift
z
l
,
D
c
and
V
c
are the comoving distance and volume, respectively,
at a given redshift,
R
m
(
z
m
) is the total comoving merger rate
density at redshift
z
m
, (1
+
z
m
) accounts for the cosmological
time dilation,
p
(
ρ
|
z
m
) is the distribution of signal-to-noise
ratio (SNR) at a given redshift, and
σ
is the lensing cross-
section (Appendix A). Throughout this section and in Sec. 3.3
we choose a network SNR threshold
ρ
c
=
8 as a point esti-
mator of the detectability of GW signals. We find it to be
consistent with the search results in Abbott et al. (2021a) and
in Sec. 5.3, and we estimate its impact to be subdominant with
respect to other source of uncertainties.
4
Table 1.
Expected fractional rates of observable lensed double events at current LIGO–Virgo sensitivity.
Merger Rate Density
Galaxies
Galaxy Clusters
Model
R
D
R
S
R
D
R
S
A
0
.
9–4
.
4
×
10
4
2
.
9–9
.
5
×
10
4
0
.
4–1
.
8
×
10
4
1
.
4–4
.
1
×
10
4
B
1
.
0–23
.
5
×
10
4
2
.
5–45
.
2
×
10
4
0
.
7–10
.
9
×
10
4
1
.
6–19
.
9
×
10
4
N
ote
—This table lists the relative rates of lensed double events expected to be observed by LIGO–Virgo at the current sensitivity where both
of the lensed events are detected (
R
D
) and only one of the lensed events is detected (
R
S
) above the SNR threshold. For Model A, the range
corresponds to the bracketing function (see Appendix A) and for Model B, the rates encompass a 90 per cent credible interval. We show the rate
of lensing by galaxies (
σ
vd
=
10–300
km s
1
) and galaxy clusters (
log
10
(
M
halo
/
M
)
14–16) separately. Besides their usage for forecasts, the
fraction of lensed events allows us to interpret the prior probability of the strong lensing hypothesis, which we require to identify lensed events
confidently.
In Table 1, we show our estimates of the relative rate of
lensing assuming di
ff
erent models (Models A and B) for the
merger rate density. The results are shown separately for
galaxy-scale (G) and cluster-scale (C) lenses. Furthermore,
these rates are calculated for events that are doubly lensed
and for two cases: when only a single event (i.e., the brighter
one) is detected (S), and when both of the doubly lensed
events are detected (D). The expected fractional rate of lens-
ing (lensed to unlensed rate), which will be necessary for the
multi-image analyses (Sec. 5), ranges from
O
(10
3
–10
4
), de-
pending on the merger rate density assumed. We estimate the
fractional rate of observed double (single) events for galaxy-
scale lenses in the range 0
.
9
4
.
4
×
10
4
(2
.
9
9
.
5
×
10
4
) when
using Model A for the merger rate density. Similarly, for
cluster-scale lenses, the fractional rate is estimated to be in the
range of 0
.
4
1
.
8
×
10
4
(1
.
4
4
.
1
×
10
4
) much rarer than the
rates at galaxy scales. These estimates suggest that observing
a lensed double image is unlikely at the current sensitivity
of the LIGO–Virgo network of detectors. Nevertheless, at
design sensitivity and with future upgrades, standard forecasts
suggest that the possibility of observing such events might
become significant (Ng et al. 2018; Li et al. 2018; Oguri 2018;
Xu et al. 2021; Mukherjee et al. 2021b; Wierda et al. 2021).
Our lensing rates are consistent with those predicted for sin-
gular isothermal ellipsoid models (e.g., Oguri 2018; Xu et al.
2021; Wierda et al. 2021). The main uncertainty in the rate
estimates derives from the uncertainties in the merger-rate
density at high redshift.
Depending on the specific distribution of lenses and the
source population, the time delays between images can change.
Models favoring galaxy lensing produce minutes to perhaps
months of time delay, while galaxy cluster lensing can pro-
duce time delays up to even years. However, the time delay
distribution for galaxy cluster lenses is more di
ffi
cult to model
accurately, owing to the more complex lensing morphology.
Since the merger rate density at high redshift is observation-
ally constrained only by the absence of the SGWB, these rates
are subject to uncertainty. Nevertheless, standard theoretical
0
2
4
6
8
z
10
1
10
2
10
3
10
4
10
5
R
m
(z) [Gpc
3
yr
1
]
No lensing
Lensing
GD
GS
CD
CS
Figure 1.
Merger rate density as a function of redshift based on the
GWTC-2 results without lensing constraints (blue) and with lens-
ing (red) included in LIGO–Virgo detections. We show the results
for galaxy-scale lenses (G) and cluster-scale lenses (C) separately.
Furthermore, S (or D) correspond to doubly lensed events where
single (or double) events are detected. Because lensed detections
occur at higher redshifts than unlensed events, their non-observation
can be used to constrain mergers at higher redshifts. The results
without lensing do not include constraints derived from the absence
of a SGWB.
models will still produce useful forecasts. We will later refer
to these rate estimates in the relevant sections (see Sec. 5).
3.2.
Implications from the non-observation of strongly lensed
events
Motivated by the absence of evidence for strong lensing
(Sec. 5), we assume that no strong lensing has occurred, in
order to constrain the merger rate density at high redshift.
We use the standard constraints on the merger rate density at
low redshift from the LIGO–Virgo population studies (Abbott
et al. 2021d). We assume the Madau–Dickinson form for
the merger rate density (Model B). This model’s free param-
eters include the local merger rate density, the merger rate
density peak, and the power-law slope. The non-observation
of lensing constrains the merger-rate density at high redshift,
which is unconstrained by the low-redshift observations alone
5
10
0
10
1
10
2
Magnification
μ
10
6
10
5
10
4
10
3
10
2
10
1
10
0
P
(
μ
)
Figure 2.
Cumulative fraction of lensed detectable BBH mergers
at any redshift with magnification greater than
μ
, constrained by
the non-observation of the SGWB. The solid line shows the value
obtained from the median BBH merger rate density posterior. The
shaded region corresponds to the 90% credible interval. Fewer than
1 in 10
3
events are expected to be lensed with magnification
μ >
2,
on average. Significantly higher magnifications (e.g.,
μ >
30) are
suppressed by a further factor of 10. The results here show the
probability of
observing
an event above a given magnification, which
includes the merger-rate density and magnification bias information.
(Fig. 1). These lensing constraints are complementary to the
current strictest high-redshift limits obtained through SGWB
non-observation (Abbott et al. 2021c).
3.3.
Constraints from stochastic background
We can also constrain the redshift evolution of the merger
rate density from the reported non-observation of the SGWB
from BBHs (Callister et al. 2020; Abbott et al. 2021c). This, in
turn, provides constraints on the relative abundance of distant
mergers, which are more likely to undergo lensing. Thus, the
non-observation of the SGWB can inform the estimate of the
probability of observing lensed BBH mergers (Buscicchio
et al. 2020a; Mukherjee et al. 2021a).
Following Buscicchio et al. (2020a), we forecast constraints
on the merger rate density in O3 using up-to-date constraints
on the mass distribution and redshift evolution of BBH merg-
ers obtained from the latest detections (Abbott et al. 2019a,b,
2021a,d), as well as those inferred from current upper limits
on the SGWB, given its non-observation (Abbott et al. 2021c).
While the measured parameters for each merger (redshifts,
source masses) are potentially biased by lensing, as discussed
in Sec. 4, we express all quantities as functions of non-biased
merger redshift
z
m
and chirp mass
M
(Buscicchio et al. 2020a)
for consistency with other sections. However, following Bus-
cicchio et al. (2020a), we do not assume as prior information
that lensing is not taking place. Instead, we include the mag-
nification bias self-consistently in the analysis, by imposing
population constraints in apparent masses and redshifts.
We model the di
ff
erential lensing probability following Dai
et al. (2017). The di
ff
erential merger rate in a redshift and
magnification shell is
d
2
R
d
z
m
d ln
μ
=
d
P
(
μ
|
z
m
)
d ln
μ
4
π
D
2
c
(
z
m
)
H
0
(
1
+
z
m
)
E
(
z
m
)
×
d
m
1
d
m
2
d
3
R
m
(
z
m
)
d
m
1
d
m
2
d
z
m
p
(
ρ>ρ
c
|
m
1
,
m
2
,
z
m
)
,
(2)
where d
3
R
m
(
z
m
)
/
d
m
1
d
m
2
d
z
m
is the di
ff
erential merger rate
density;
p
(
ρ>ρ
c
|
m
1
,
m
2
,
z
m
)
provides the probability of
observing mergers with source masses
m
1
,
m
2
, redshift
z
m
, and
magnified by a factor
μ
above a fixed network SNR threshold
ρ
c
=
8, integrated over the population distribution of source
parameters; the factor 4
π
D
2
c
(
z
m
)
/
[
H
0
(1
+
z
m
)
E
(
z
m
)] gives the
comoving volume of a redshift shell in an expanding Universe
(taking into account the redshifted rate definition with respect
to the source frame); and d
P
(
μ
|
z
m
)
/
(d
ln
μ
) is the lensing
probability. However, as noted by Dai et al. (2017), the di
ff
er-
ential magnification probability at 0
.
9
< μ <
1
.
1 and
z
m
<
2
is a
ff
ected by relative uncertainties up to 40%. We therefore
consider magnified detections only (
μ >
1), which are subject
to less uncertainty, and normalize our results accordingly. We
then integrate the di
ff
erential merger rate (Eq. 2) over redshift
and magnifications in
[
μ,μ
max
]
and divide it by the total rate
of magnified detections. By doing so, we obtain the cumu-
lative fraction of detected lensed events at any redshift with
magnifications larger than
μ
.
We show the result in Fig. 2. We find the observation of
lensed events to be unlikely, with the fractional rate at
μ >
2
being 4
.
9
+
1
.
7
1
.
3
×
10
4
. More significantly magnified events
are even more suppressed, with a rate of 3
.
5
+
0
.
6
0
.
4
×
10
5
at
μ >
30. These estimates suggest that most binary mergers that
we observe are not strongly lensed. However, as projected
in Buscicchio et al. (2020a); Mukherjee et al. (2021a), at
design sensitivity, the same probability will be enhanced, as a
widened horizon will probe the merger rate density deeper in
redshift.
Comparing the above predictions with the expected frac-
tional rates
R
S
of single lensed detections with Model B in
Table. 1, the predictions agree within a factor of 5 for the
relative rate of lensing. The di
ff
erences are due to a dif-
ferent underlying lens model and partly to the inclusion of
de-magnified events in Sec. 3.1.
4.
ANALYZING HIGH-MASS EVENTS
If a GW signal is strongly lensed, it will receive a magni-
fication
μ
defined such that the GW amplitude increases by
a factor
|
μ
|
1
/
2
relative to an unlensed signal. The luminosity
distance inferred from the GW observation will be degener-
ate with the magnification such that the inferred luminosity
distance
D
inferred
L
=
D
L
|
μ
|
.
(3)
6
Because of this degeneracy, lensing biases the inferred redshift
and thus the source masses. Consequently, the binary appears
to be closer than it truly is, and it appears to be more massive
than it truly is.
Broadhurst et al. (2018, 2020a,b) argued that some of the
relatively high-mass LIGO–Virgo events could be strongly
lensed GWs from the lower-mass stellar black hole popula-
tion observed in the electromagnetic bands. However, the
expected strong lensing rates and the current constraints on
the merger-rate density, based on the absence of a detectable
SGWB, disfavor this interpretation (Dai et al. 2017; Ng et al.
2018; Li et al. 2018; Oguri 2018; Hannuksela et al. 2019; Bus-
cicchio et al. 2020a,b) compared to the standard interpretation
of a genuine unlensed high-mass population (Abbott et al.
2019a; Roulet et al. 2020; Abbott et al. 2021d; Kimball et al.
2020). Hence, in the absence of more direct evidence, such as
identifying multiple images within LIGO–Virgo data (Sec. 5),
it is di
ffi
cult to support the lensing hypothesis purely based on
magnification considerations. Nevertheless, it is informative
to analyze the degree to which the lensed interpretation would
change our understanding of the observed sources.
Under the strong lensing hypothesis
H
SL
, the GW would
originate from a well-known, intrinsically lower-mass popu-
lation, and the LIGO–Virgo observations have been biased
by lensing. Using such a mass prior, we infer the required
magnification and corrected redshift and component masses
under
H
SL
. The posterior distribution of the parameters is
(Pang et al. 2020)
p
(
μ,θ
|
d
,
H
SL
)
p
(
d
|
θ
)
p
(
θ
|
μ,
H
SL
)
p
(
μ
|H
SL
)
,
(4)
where we distinguish the
apparent
parameters of the wave-
form received at the detector
θ
, which di
ff
er from the intrinsic
parameters
θ
due to bias by lensing magnification. There-
fore, we can compute the magnification posterior and other
parameters by simply re-weighting existing posteriors.
Studies along these lines were already done for the
GW190425 BNS event by Pang et al. (2020) and for the
GW190521 BBH event in Abbott et al. (2020e). Here we ex-
tend the approach to cover additional interesting O3a events,
focusing on two cases: (i) the (apparently) most massive
observed BBHs, and (ii) sources with an (apparent) heavy
neutron star component. In the BBH case, we take the prior
over component masses,
m
1
and
m
2
, and redshift,
z
of the
source
p
(
m
1
,
m
2
,
z
) from the power-law BBH population
model used in Abbott et al. (2019a) for O1 and O2 obser-
vations, with a mass power-law index
α
=
1, mass ratio
power-law index
β
q
=
0, and minimum component mass
m
min
=
5
M
, and assume an absence of BBHs above the pair
instability supernova (PISN) mass gap. As in the previous
GW190521 study (Abbott et al. 2020e), we consider two
di
ff
erent values to account for uncertainties on the edge of
the PISN gap,
m
max
=
(50
,
65)
M
. Such a simple model is
adequate for this analysis because our analysis results are
most sensitive to the mass cut (highest masses allowed by
the prior) and less sensitive to the specific shape of the mass
distribution. For events with an apparent heavy neutron star
component, we assume a Galactic BNS prior following a
total mass with a 2
.
69 M
mean and 0
.
12 M
standard devi-
ation (Farrow et al. 2019). In both cases, the magnification
could explain the apparent high mass of the events from the
LIGO–Virgo observations.
We assume that the redshift prior
p
(
z
)
τ
(
z
)d
V
c
/
d
z
, where
the optical depth of lensing by galaxies or galaxy clusters
τ
(
z
)
D
c
(
z
)
3
(Haris et al. 2018). The redshift dependence
of the optical depth is approximately the same for both
galaxies and galaxy clusters, while the overall scaling can
change (Fukugita & Turner 1991). We use the lensing prior
p
(
μ
|H
SL
)
μ
3
(Blandford & Narayan 1986) with a lower
limit
μ >
2 appropriate to strong lensing (Ng et al. 2018). This
prior is appropriate when we are in the high-magnification,
strong lensing limit, i.e. assuming that the observed masses
are highly biased. We do not consider weak lensing, which
does not produce multiple images and would require expanded
future GW data sets to study (Mukherjee et al. 2020a,b).
We analyze all O3a BBH events with the primary mass
above 50 M
at 90% probability using the Bayesian inference
posterior samples released with GWTC-2 (GWOSC 2020; Ab-
bott et al. 2021a). Moreover, we analyze GW190425, a high-
mass BNS (Abbott et al. 2020b), and GW190426
152155, a
low-significance potential neutron star–black hole (NSBH)
event (Abbott et al. 2021a) which was investigated as a pos-
sible lensed BNS event (Smith et al. 2019). We use the re-
sults for the
IMRP
henom
P
v
2 waveform (Hannam et al. 2014;
Boh
́
e et al. 2016) for most of the events. For GW190521,
where higher-order multipole moments are important to in-
clude in the analysis (Abbott et al. 2020e), we adopt the
NRS
ur
7
dq
4 waveform (Varma et al. 2019) results as in Ab-
bott et al. (2020f). Furthermore, for GW190425 (Abbott et al.
2020b), we use the
IMRP
henom
P
v
2
NRT
idal
(Dietrich et al.
2019) low-spin samples. Results are summarized in Table 2.
To interpret the heavy BBHs as lensed signals originat-
ing from the assumed lower-mass population, they should be
magnified at a moderate magnification
μ
10 at
z
1–2.
Depending on the lens model, this magnification may imply a
moderate chance of an observable multi-image counterpart as
events closer to the caustic curves experience more substantial
magnifications. Consequently, they often produce events with
similar magnification ratios and shorter time delays (compa-
rable magnifications and shorter time delays can be derived
from the lens’s symmetry, although if lensing by substructures
or microlenses is present, the magnifications between images
can di
ff
er even in the high-magnification limit). However, we
could not identify any multi-image counterparts for any of the
high-mass events in our multiple image search (Sec. 5).
7
Table 2.
Inferred properties of selected O3a events under the lensing magnification hypothesis.
Event name
m
1
[M
]
m
2
[M
]
z
μ
GW190425
1
.
42
+
0
.
16
0
.
12
1
.
27
+
0
.
12
0
.
15
0
.
3
+
0
.
1
0
.
1
68
+
163
44
GW190426
152155
1
.
89
+
0
.
40
0
.
55
0
.
90
+
0
.
25
0
.
40
1
.
3
+
0
.
5
0
.
2
497
+
452
272
Event name
m
50
1
(
m
65
1
) [M
]
m
50
2
(
m
65
2
) [M
]
z
50
(
z
65
)
μ
50
(
μ
65
)
GW190521
43
+
6
16
(55
+
9
22
)
36
+
10
15
(45
+
13
19
)
2
.
5
+
2
.
1
0
.
7
(1
.
8
+
1
.
7
0
.
5
)
13
+
55
8
(6
+
28
4
)
GW190602
175927
42
+
7
17
(48
+
14
19
)
31
+
13
16
(33
+
15
16
)
1
.
4
+
1
.
5
0
.
5
(1
.
1
+
1
.
4
0
.
4
)
10
+
65
7
(6
+
46
4
)
GW190706
222641
39
+
10
15
(42
+
17
17
)
29
+
12
13
(29
+
13
13
)
1
.
7
+
1
.
8
0
.
5
(1
.
6
+
1
.
7
0
.
6
)
5
+
26
3
(4
+
22
2
)
N
ote
—Under the hypothesis that the listed events are lensed signals from intrinsically lower-mass binary populations with
μ >
2, this table lists
the favored source masses, redshifts, and magnifications for the BNS and neutron star–black hole (NSBH) (top) and BBH (bottom) high-mass
events. For the BBHs, two sets of numbers are given for di
ff
erent assumptions about the edge of the pair instability supernova (PISN) mass gap
(a cut at 50 M
and 65 M
). For the BNSs, we presume that they originate from the Galactic BNS population. To interpret the heavy BBHs as
lensed signals originating from the assumed lower-mass population, they should be magnified at a moderate magnification
μ
∼O
(10) at
z
1 to
2. The BNS and NSBH events would require extreme magnifications.
The BNS and NSBH events, on the other hand, would
require extreme magnifications (68
+
163
44
and 497
+
452
272
, respec-
tively) to be consistent with the Galactic BNS distribution. At
these magnifications, we would expect the source to be close
to a caustic, and therefore it may be possible that the presence
of microlenses would produce observable e
ff
ects (Diego et al.
2019; Diego 2020; Pagano et al. 2020; Mishra et al. 2021).
Moreover, the event would likely be multiply imaged (Bland-
ford & Narayan 1986; Oguri 2018). A more detailed follow-
up study to quantify the likelihood of multiple images and
microlensing could produce more stringent evidence for the
lensing hypothesis for these events. We will briefly comment
on these events in the context of multi-image and microlensing
results in the sections that follow.
At this stage, we cannot set robust constraints on the lens-
ing hypothesis based on the magnification alone. Moreover,
as detailed in the following section, we have also not found
any other clear evidence to indicate that these GW events
are lensed. The prior lensing rate disfavors the lensing hy-
pothesis for most standard binary population and lens models,
as discussed in Sec. 3. However, if other BBH formation
channels exist that produce an extensive number of mergers
at high redshift, the lensing rates can change. In the future,
more quantitative constraints could be set by connecting the
inferred magnifications with lens modeling to make predic-
tions for the appearance of multiple images or microlensing
e
ff
ects.
5.
SEARCH FOR MULTIPLE IMAGES
In addition to magnification, strong lensing can produce
multiple images of a single astrophysical event. These multi-
ple images appear at the GW detectors as repeated events. The
images will di
ff
er in their arrival time and amplitude (Wang
et al. 1996; Haris et al. 2018; Hannuksela et al. 2019; Li
et al. 2019a; McIsaac et al. 2020). The sky location is the
same within the localization accuracy of GW detectors, given
that the typical angular separations are of the order of arcsec-
onds. Additionally, lensing can invert or Hilbert transform
the image (Dai & Venumadhav 2017; Ezquiaga et al. 2021),
introducing a frequency-independent phase shift. This trans-
formation depends on the image type, set by the lensing time
delay at the image position: Type-I, II, and III correspond to
a time-delay minimum, saddle point, and maximum, respec-
tively (Ezquiaga et al. 2021).
The multiply imaged waveforms
{
̃
h
L
j
}
of a single signal
̃
h
then satisfy (Dai & Venumadhav 2017; Ezquiaga et al. 2021)
̃
h
L
j
(
f
;
θ,μ
j
,
t
j
,
φ
j
)
=
|
μ
j
|
̃
h
(
f
;
θ,
t
j
) exp
(
i
sign(
f
)
φ
j
)
,
(5)
where
|
μ
j
|
is the lensing magnification experienced by the
image
j
and
φ
j
=
π
n
j
/
2 is the Morse phase, with index
n
j
=
0
,
1
,
2 for Type-I, II, and III images.
̃
h
(
f
;
θ,
t
j
) is the
original (unlensed) waveform before lensing, but evaluated as
arriving with a time delay
t
j
. The multi-image hypothesis
then states that most parameters measured from the di
ff
erent
lensed images of the same event are consistent.
The relative importance of di
ff
erent parameters for the
overall consistency under the multi-image hypothesis will
vary for di
ff
erent events. For example, the sky localization
match will have greater relevance for well-localized, high-
SNR events. Similarly, the overlap in measured chirp mass
(1
+
z
)
M
=
(1
+
z
) (
m
1
m
2
)
3
/
5
/
(
m
1
+
m
2
)
1
/
5
will be more signifi-
cant when the uncertainty in that parameter is lower, although
in this case the underlying astrophysical mass distribution will
play a key role. The similarities in other parameters such as
mass ratios or spins will be more important when they depart
from the more common astrophysical expectations. Evidence
of strong lensing could also be acquired with a single Type-II
8
(saddle point) image if the induced waveform distortions in
the presence of higher modes, precession, or eccentricity are
observed (Ezquiaga et al. 2021). Such evidence is unlikely
to be observed without next-generation detectors (Wang et al.
2021).
In this section, we perform three distinct but related anal-
yses. First, we test the lensed multi-image hypothesis by
analyzing, for all pairs of O3a events from GWTC-2, the
overlap of posterior distributions previously inferred for the
individual events. This allows us to set ranking statistics to
identify an initial set of candidates for lensed multiple images.
We perform a more detailed joint-PE analysis for these most
promising pairs, considering all potential correlations in the
full parameter space and the image type. This joint analysis
provides a more solid determination of the lensing probabil-
ity for a given GW pair. Finally, we search for additional
sub-threshold candidates that could be multiply imaged coun-
terparts to the previously considered events: some counterpart
images can have lower relative magnification compared with
the primary image and
/
or fall in times of worse detector sensi-
tivity or antenna patterns, and hence may not have passed the
detection threshold of the original broad searches. According
to the predictions of the expected lensing time delays and the
rate of galaxy and galaxy cluster lensing (Smith et al. 2018;
Oguri 2018; Dai et al. 2020), we expect it to be less likely
for counterpart images to the O3a events to be detected in
observing runs O1 or O2. Relative lensing rates for galaxies
and clusters are given in Table 1. Thus, we only search for
multiple images within O3a itself.
Previous studies have also searched for multiple images in
the O1–O2 catalog GWTC-1 (Hannuksela et al. 2019; Broad-
hurst et al. 2019; Li et al. 2019a; McIsaac et al. 2020; Dai
et al. 2020; Liu et al. 2021). The first search for GW lensing
signatures in O1 and O2 focused on the posterior overlap of
the masses, spins, binary orientation and sky positions (Han-
nuksela et al. 2019) and the consistency of time delays with
expectations for galaxy lenses, but found no conclusive ev-
idence of lensing. The search did uncover a candidate pair
GW170104–GW170814 with a relatively high Bayes factor
of
&
200. Still, this study disfavored the candidate due to its
long time delay and the low prior probability of lensing. In
parallel, Broadhurst et al. (2019) suggested that the candidate
pair GW170809–GW170814 could be lensed, but this claim is
disfavored by more comprehensive analyses (Hannuksela et al.
2019; Liu et al. 2021). Both Li et al. (2019a) and McIsaac et al.
(2020) performed searches for sub-threshold counterparts to
the GWTC-1 events, identifying some marginal candidates
but finding no conclusive evidence of lensing. More recently,
Dai et al. (2020) and Liu et al. (2021) searched for lensed
GW signals including the analysis of the lensing image type,
which can be described through the Morse phases,
φ
j
in
Eq.
(5)
. These analyses have revisited the pair GW170104–
GW170814 and demonstrated that the Morse phase is consis-
tent with the lensed expectation but would require Type-III
(time-delay maximum) images, which are rare from an obser-
vational standpoint. Dai et al. (2020) also pointed out that a
sub-threshold trigger, designated by them as GWC170620, is
also consistent with coming from the same source. However,
the required number and type of images for this lens system
make the interpretation unlikely given current astrophysical
expectations. Also, two same-day O3a event pairs (on 2019
May 21 and 2019 August 28) have already been considered
elsewhere, but were both ruled out due to vanishing localiza-
tion overlap (Singer et al. 2019; Abbott et al. 2020e).
5.1.
Posterior-overlap analysis
As a consequence of degeneracies in the measurements of
parameters, the lensing magnification can be absorbed into
the luminosity distance (Sec. 4), the time delay can be ab-
sorbed into the time of coalescence, and, when the radiation is
dominated by
`
=
|
m
|
=
2 multipole moments, the phase shifts
introduced by lensing (the Morse phases) can be absorbed into
the phase of coalescence. The multi-image hypothesis then
states that all other parameters except the arrival time, lumi-
nosity distance, and coalescence phase are the same between
lensed events, and thus there should be extensive overlap in
their posterior distributions, even if those have been inferred
without taking lensing into account.
Therefore, we use the consistency of GW signals detected
by LIGO and Virgo to identify potential lensed pairs. Follow-
ing Haris et al. (2018), we define a ranking statistic
B
overlap
to
distinguish candidate lensed pairs from unrelated signals:
B
overlap
=
d
Θ
p
(
Θ
|
d
1
)
p
(
Θ
|
d
2
)
p
(
Θ
)
,
(6)
where the parameters
Θ
include the redshifted masses (1
+
z
)
m
1
,
2
, the dimensionless spin magnitudes
χ
1
,
2
, the cosine of
spin tilt angles
θ
1
,
2
, the sky location (
α,
sin
δ
), and the cosine
of orbital inclination
θ
JN
, but they do not include the full 15-
dimensional set of parameters
Θ
to ensure the accuracy of
the kernel density estimators (KDEs) that we use to approxi-
mate the posterior distributions
p
(
Θ
|
d
1
,
2
) for each event when
evaluating Eq. (6). Here,
p
(
Θ
) denotes the prior on
Θ
.
The accuracy of the KDE approximation was demonstrated
in Haris et al. (2018) through receiver operating characteristic
curves with simulated lensed and unlensed BBH events. To
improve the accuracy further, we compute the sky localization
(
α
,
δ
) overlap separately from other parameters and combine it
with the overlap from the remaining parameters. Splitting the
two overlap computations is justified because the posterior
correlations of (
α
,
δ
) with other parameters are minimal.
We use posterior samples (GWOSC 2020) obtained us-
ing the
LALI
nference
software package (Veitch et al. 2015)
with the
IMRP
henom
P
v
2 waveform model (Hannam et al.