# On the gravitational lensing interpretation of three gravitational wave detections in the mass gap by LIGO and Virgo

Matteo Bianconi, Graham P. Smith, Matt Nicholl, Dan Ryczanowski, Johan Richard, Mathilde Jauzac, Richard Massey, Andrew Robertson, Keren Sharon, Evan Ridley

<figure><img src="image_1.png" /></figure>

## To cite this version:

Matteo Bianconi, Graham P. Smith, Matt Nicholl, Dan Ryczanowski, Johan Richard, et al. On the gravitational lensing interpretation of three gravitational wave detections in the mass gap by LIGO and Virgo. Monthly Notices of the Royal Astronomical Society, 2023, 521 (3), pp.3421-3430. 10.1093/mnras/stad673. hal-03665331

HAL Id: hal-03665331

https://hal.science/hal-03665331

Submitted on 15 Feb 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the gravitational lensing interpretation of three gravitational wave detections in the mass gap by LIGO and Virgo

Matteo Bianconi<sup>1,★</sup> Graham P. Smith<sup>1</sup>, Matt Nicholl<sup>1,2</sup> Dan Ryczanowski<sup>1</sup>, Johan Richard<sup>1,3</sup>, Mathilde Jauzac<sup>4,5,6,7</sup> Richard Massey<sup>5</sup>, Andrew Robertson<sup>8</sup>, Keren Sharon<sup>9</sup> and Evan Ridley<sup>1</sup>

Accepted 2023 February 27. Received 2023 February 23; in original form 2022 April 28

## ABSTRACT

We search for gravitational wave (GW) events from LIGO-Virgo’s third run that may have been affected by gravitational lensing. Gravitational lensing delays the arrival of GWs, and alters their amplitude – thus biasing the inferred progenitor masses. This would provide a physically well-understood interpretation of GW detections in the ‘mass gap’ between neutron stars and black holes, as gravitationally lensed binary neutron star (BNS) mergers. We selected three GW detections in LIGO-Virgo’s third run for which the probability of at least one of the constituent compact objects being in the mass gap was reported as high with low latency – i.e. candidate lensed BNS mergers. Our observations of powerful strong lensing clusters located adjacent to the peak of their sky localization error maps reached a sensitivity AB $\simeq 25.5$ in the z’ band with the GMOS instruments on the Gemini telescopes, and detected no candidate lensed optical counterparts. We combine recent kilonova light-curve models with recent predictions of the lensed BNS population and the properties of the objects that we followed up to show that realistic optical counterparts were detectable in our observations. Further detailed analysis of two of the candidates suggests that they are a plausible pair of images of the same low-mass binary black hole merger, lensed by a local galaxy or small group of galaxies. This further underlines that access to accurate mass information with low latency would improve the efficiency of candidate lensed BNS selection.

**Key words:** gravitational lensing: strong – gravitational waves – galaxies: clusters: individual Abell 370, MACS J2135.2–0102, RX J2129.6+0005.

## 1 INTRODUCTION

Gravitational lensing is contributing to unveiling otherwise inaccessible regions of the Universe. Intervening mass along the line of sight of conventional observations acts to magnify the observed radiation, allowing the detection of background objects whose direct electromagnetic (EM) radiation would be otherwise too faint because of their distance. This enabled the detection of increasingly remote galaxy populations (e.g. Kneib et al. 2004; Bouwens et al. 2014) and even individual stars at high redshifts (e.g. Kelly et al. 2018; Welch et al. 2022). Similarly, gravitational lensing can play a key role in investigating distant populations of gravitational wave (GW) sources, reaching beyond the limited sensitivity of current detectors, in particular in the low-mass regime of binary neutron stars (BNSs; Smith et al. 2023, and references therein). Therefore, it is highly relevant to investigate the scope for lensing as it currently provides

the only window to low-mass GW progenitors beyond the local Universe.

The census of confirmed GW events has steadily increased since the first detection (Abbott et al. 2016) due to the continuous performance improvement of the network of instruments composing the LIGO, Virgo, and KAGRA collaboration (LVK). The landscape of compact object mergers responsible for the GW emission is being continuously populated, in particular of progenitor masses in the range $1 < m[\text{M}_\odot] < 100$. Interestingly, theoretical models of stellar evolution exclude the presence of compact objects in the mass interval $2 < m[\text{M}_\odot] < 5$, whose limits encompass the most massive neutron stars and the least massive black holes, respectively (Farr et al. 2011; Alsing, Silva & Berti 2018). On the other hand, core-collapse supernova models can enforce this mass gap or produce a smooth remnant mass distribution by changing physical assumptions on the onset of the supernova explosion (e.g Belczynski et al. 2012; Olejak et al. 2022). In their third run LVK detected candidate mass gap binary mergers hinting at different and varied formation channels responsible for the binary progenitors from the currently available

<footnote n="1"><sup>1</sup>School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK</footnote>

<footnote n="2"><sup>2</sup>Institute for Gravitational Wave Astronomy, University of Birmingham, Birmingham B15 2TT, UK</footnote>

<footnote n="3"><sup>3</sup>Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval, France</footnote>

<footnote n="4"><sup>4</sup>Centre for Extragalactic Astronomy, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="5"><sup>5</sup>Institute for Computational Cosmology, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="6"><sup>6</sup>Astrophysics Research Centre, University of KwaZulu-Natal, Westville Campus, Durban 4041, South Africa</footnote>

<footnote n="7"><sup>7</sup>School of Mathematics, Statistics & Computer Science, University of KwaZulu-Natal, Westville Campus, Durban 4041, South Africa</footnote>

<footnote n="8"><sup>8</sup>Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr, Pasadena, CA 91109, USA</footnote>

<footnote n="9"><sup>9</sup>Department of Astronomy, University of Michigan, 1085 S. University Ave, Ann Arbor, MI 48109, USA</footnote>

<footnote>*E-mail: mbianconi@star.sr.bham.ac.uk</footnote>
in literature (Gupta et al. 2020). Alternatively, gravitational lensing offers a physically well-understood mechanism, whose intervention would cause the incorrect estimate of the GW progenitor parameters.

LVK classifies a GW source as a binary comprising of neutron stars (BNS), black holes (BBH), or neutron star-black hole system, according to the inferred mass retrieved from modelling of the detected waveform. The amplitude of the GW strain signal that the detectors measure, $A$, is affected by gravitational lensing magnification by intervening mass along the line of sight and the luminosity distance to the source, such that from the inverse square law: $A \propto \sqrt{\mu/D}$ (Wang, Stebbins & Turner 1996). Here, $\mu$ is the magnification caused by gravitational lensing. The current LVK pipeline defaults to null gravitational lensing intervention, corresponding to magnification of $\mu = 1$. Therefore, if the GW source is actually strongly lensed (multiply imaged and thus $\mu \gtrsim 2-10$), LVK source distance posteriors will be biased low. Contextually, the inferred rest-frame mass $\tilde{m}$ of the compact objects responsible for the GW emission scales as $m(1+z) = \tilde{m}(1+z)$, where $m$ and $z$ are the true mass and redshift of the GW source, and $\tilde{m}$ and $\tilde{z}$ are the mass and redshift inferred with low latency assuming $\mu = 1$. Hence, the masses inferred assuming $\mu = 1$ require to be revised down if the source is lensed. This implies the released LVK event classification based on binary component mass could be biased (see also Smith et al. 2023 for additional details).

The LVK consortium has recently published a study on the impact of lensing on the GW detection rate and on searching for multiple images due to strong lensing within the events of the first half of O3, concluding against the occurrence of lensing (Abbott et al. 2021c). This study relied on the analysis of the GW data stream by LVK. We argue for the need of optical follow-up observations of candidate lensed GW sources in order to localize the source to a gravitational lens – i.e. sub-arcsecond accuracy. This is currently not achievable with the LVK data alone, despite the improvements in sky localization achieved with three detectors (Abbott et al. 2017a, b, c). This motivates concentrating on lensed BNS mergers because BNSs are now confirmed as being associated with EM counterpart, i.e. kilonovae (KNe, e.g. Abbott et al. 2017c). In contrast, BBH mergers are expected to have no or very faint counterpart that is beyond the reach of today’s telescopes (see also Graham et al. 2020), even when aided by gravitational lensing (e.g. Smith et al. 2019a, and references therein).

If on the one hand lensing masks intrinsic properties of the GW binary, on the other hand it can be used to search further along the recesses of the Universe for the EM counterpart of the GW event. Simple lensing arguments show the impact of magnification on the expected observed magnitude of the EM counterpart of a BNS. For example, consider the detection of a GW170817/AT2017gfo-like counterpart to a BNS merger located at $z \simeq 1$ that is lensed by a massive galaxy or cluster in the foreground. Within four observer-frame days from the GW event, such an object would have absolute magnitude of $M \lesssim -14$ (Arcavi 2018). If the EM counterpart to a gravitationally lensed BNS merger located at $z \simeq 1$ has a similar luminosity, then its apparent magnitude is given by $m \simeq -14 + 5\log(D/10\,\text{pc}) - 2.5\log(\mu)$, where $D$ is the intrinsic luminosity distance to the lensed GW, $\mu = (D/\tilde{D})^2$ is the gravitational magnification suffered by the source, and $\tilde{D}$ is the luminosity distance initially assigned to the source by LVK assuming $\mu = 1$. Interestingly, $m \simeq -14 + 5\log(\tilde{D}/10\,\text{pc})$, because the inverse square law cancels the gravitational magnification and k-corrections appear to be modest (Smith et al. 2023). Assuming that the GW is initially placed at $\tilde{D} = 600$ Mpc, this implies an apparent magnitude of $m \simeq 25$. Once a typical strong lensing magnification

for lensed BNS mergers of $\mu = 100$ (Smith et al. 2023) is taken into consideration, the effective magnitude that can be accessed reaches $m + 2.5\log(\mu)$, equalling $m \approx 30$. Therefore, lensing is crucial not only to boost the EM counterpart signal, but also necessary to access BNS mergers at redshift $z \gtrsim 1$.

This paper is organized as follows. In Section 2, we motivate our observational strategy and summarize the GW events from LVK’s third run that we have targeted in our follow-up campaign. In Section 3, we describe the data reduction and analysis. In Section 4, we introduce the assumptions and models used for assessing the impact of lensing on the detectability of BNS mergers and their EM counterparts, and present the models of KN light curves that we use to validate our follow-up observations. In Section 5, we discuss the physical interpretation of our observations and the candidates that we observed. We close by summarizing our results in Section 6, and by discussing the impact of posterior parameters on the observing strategy of future campaigns. We assume cosmology values presented in Planck Collaboration XIII (2016), with $h = 0.678$, $H_0 = 100\,\text{h km s}^{-1}\text{Mpc}^{-1}$, $\Omega_M = 0.309$, and $\Omega_\Lambda = 0.691$. All celestial coordinates are stated at the J2000 epoch, and all magnitudes are stated in the AB system.

## 2 CAMPAIGN RATIONALE

Our efforts are concentrated on identifying GW sources that may be affected by strong lensing. Therefore, we devised a selection criteria aiming at events that LVK initially classifies as having at least one progenitor component in the mass gap with a powerful strong lensing cluster located within its sky localization. Mass gap progenitors are compact object with masses in the range $2.5 < m[\text{M}_\odot] < 5$, which encompass the mass boundaries of the most massive neutron star (NS) and the lightest black hole (BH) known (LVK2021). Hence, any object within this mass range cannot be accounted for by the current models of stellar evolution. A natural solution to the detection of these object is invoking the impact of strong lensing, which would bias the recovered posterior masses. Recent studies from numerical simulations show that galaxies ($M_{200} < 10^{13}\text{M}_\odot$) and clusters ($M_{200} > 10^{13}\text{M}_\odot$) contribute roughly equally to the optical depth to strong lensing (Robertson et al. 2020). We focus our effort on known strong lensing clusters primarily because of the capabilities of current observational infrastructures. A search of all galaxy- and cluster-scale lenses would require $\ge 60\,\text{h}$ of observing with Subaru/HSC or CTIO/DECam to cover the full median sky localization error region of $\gtrsim 100\,\text{deg}^2$ to the necessary depth. This is stretching beyond the limit of observing capabilities at the current larger observational facilities, before the beginning of the Vera Rubin Observatory’s Legacy Survey of Space and Time. Hence, we focus on the most promising line of sight, i.e. cores of strong lensing clusters, which typically cover few square arcminutes and account for $\simeq 5$ per cent of the optical depth to strong lensing per GW sky localization (Robertson et al. 2020). We combine this with Smith et al.’s (2023) detailed predictions of the rate of lensed BNS detections in O3 ($0.02\,\text{yr}^{-1}$ for their baseline model) to estimate the probability of successfully confirming a candidate lensed BNS with detection of its EM counterpart with our strategy of $\approx 0.1$ per cent. We refer the interested reader to Robertson et al. (2020) and Smith et al. (2023) for full details of the underlying calculations.

### 2.1 Observing strategy and target selection

Our intent is to obtain new observations of strong lensing cluster cores soon after the detection of a GW event, to search for EM transients
<table><caption>**Table 1.** Summary of the salient posteriors and their 90 per cent confidence interval of the confirmed GW events considered in this study, from LVK data release assuming no lensing. From left to right, name of the confirmed GW event, trigger ID, mass gap probability released by LVK with low latency, individual masses and chirp mass of the binary, luminosity distance, redshift, network match-filtered SNR, and credible area of the sky localization from Abbott et al. (2021a) and LVK2021.</caption><thead><tr><th>Name</th><th>Trigger</th><th>p<sub>gap</sub></th><th>~m<sub>1</sub> [M<sub>☉</sub>]</th><th>~m<sub>2</sub> [M<sub>☉</sub>]</th><th>~M [M<sub>☉</sub>]</th><th>~D [Mpc]</th><th>~z</th><th>SNR</th><th>90 per cent Skymap area [deg<sup>2</sup>]</th></tr></thead><tbody><tr><td>GW 190930.133541</td><td>S190930s</td><td>&gt; 95 per cent</td><td>14.2<sup>+8.0</sup><sub>-4.0</sub></td><td>6.9<sup>+2.4</sup><sub>-2.1</sub></td><td>8.5<sup>+0.5</sup><sub>-0.4</sub></td><td>770<sup>+320</sup><sub>-320</sub></td><td>0.16<sup>+0.06</sup><sub>-0.06</sub></td><td>9.7<sup>+0.3</sup><sub>-0.5</sub></td><td>1600</td></tr><tr><td>GW 191216.213338</td><td>S191216ap</td><td>&gt; 99 per cent</td><td>12.1<sup>+4.6</sup><sub>-2.3</sub></td><td>7.7<sup>+1.6</sup><sub>-1.9</sub></td><td>8.33<sup>+0.22</sup><sub>-0.19</sub></td><td>340<sup>+120</sup><sub>-130</sub></td><td>0.07<sup>+0.02</sup><sub>-0.03</sub></td><td>18.6<sup>+0.2</sup><sub>-0.2</sub></td><td>490</td></tr><tr><td>GW 200115.042309</td><td>S200115j</td><td>&gt; 94 per cent</td><td>5.9<sup>+2.0</sup><sub>-2.5</sub></td><td>1.44<sup>+0.85</sup><sub>-0.29</sub></td><td>2.43<sup>+0.05</sup><sub>-0.07</sub></td><td>290<sup>+150</sup><sub>-100</sub></td><td>0.06<sup>+0.03</sup><sub>-0.02</sub></td><td>11.3<sup>+0.3</sup><sub>-0.5</sub></td><td>370</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Sky localization in spherical projection of the candidate GW events considered in this work, together with their distance probability assuming no lensing intervention, and strong lensing cluster candidates from the Smith et al. (2018) list. The maps and distance posteriors are the ones available at the time of observations. Left: spherical projection of the trigger S200115j skymap. The inset shows the location of the strong lensing cluster Abell 370, within the 90 per cent probability of the skymap. The lineplot highlights the probability distance of the event, with respect to the cluster redshift. Right: same as left panel, for trigger S190930s and S191216ap, and the selected clusters MACS J2135.2–0102 and RX J2129.6+0005. See Section 2.1 for additional details regarding the strong lensing cluster selection.</figcaption></figure>

by comparing them to archival data set and to communicate the finding within the transient community to extend the follow-up efforts. We require rapid observations because KNe fade quickly, therefore observations as soon as possible after the LVK detection are essential. Our previous attempts at detecting lensed EM counterparts to GW events (Smith et al. 2019a), and extensive follow-up campaign of the KN AT2017gfo associated with GW170817 (Villar et al. 2017) informed our choice of instrument, waveband, and observing time. As a result, we chose redder optical wavebands, in which the EM emission from the KN is more luminous and longer lasting (Arcavi 2018). In particular, we converged on 1-h integration time and z-band filter (central wavelength ≈900 nm) which allows to reach magnitude limits $m \approx 25$ depth using an 8-m class ground-based telescope (see discussion in Section 1).

Our list of galaxy clusters acting as gravitational telescopes comprises the sample of 130 strong lensing clusters discussed by Smith et al. (2018). Due to numerous HST and 8-m class ground-based telescope observing campaigns that have targeted X-ray luminous galaxy clusters (e.g. Kneib et al. 1996; Smith et al. 2005; Richard et al. 2010; Jauzac et al. 2016; Lagattuta et al. 2017; Mahler et al. 2018), these cluster lenses are all spectroscopically confirmed, with a well-constrained model of the cluster mass distribution. The choice of the most promising cluster to observe follows the identification of the closest object to the peak of the 2d probability distribution of each GW sky localization. Additional cluster characteristics, in particular Einstein radius size, are considered in the case of multiple clusters being available to observe.

## 2.2 Candidate events with possible EM counterpart

Hereafter we summarize the mass gap events detected by LVK's third run, and our strategy to select the most promising ones, and to observe the cluster targets within the localization maps of the GW events. These events, which are listed in Table 1, have been selected for having a high mass gap probability (>94 per cent) and a low false alarm rate at the time of announcement.

LVK announced the detection of six GWs with $p_{gap} > 90$ per cent during the period overlapping with our access to Gemini, from Spring 2019 to Spring 2020. We summarize these six briefly here, in chronological order. S190426c was assigned high probability of originating from a BNS or a mass gap event (GCN24250, Smith et al. 2019b), but its sky localization shared no overlap with our galaxy cluster list. S190924h was classified as mass gap event which we did not consider due to the sub-optimal position on the sky, resulting in estimated high airmass for our telescope pointings. This was coupled with delays in the release of the offline sky localization from LVK, which resulted in us prioritizing the subsequent event S190930s instead. S190930s, together with S191216ap and S200115j, were well suited for our EM counterpart search, and are the targets of our observations. Lastly, S200316bj was communicated after our share of observing time at Gemini had been used. The target decision tree is detailed in the following subsections for each GW event we have considered in this study. A visual representation of each GW sky localization at the time of observations, together with the most promising strong lensing cluster candidates, is presented in Fig. 1. The main properties of clusters that we have selected is presented in Table 2.
<table><caption>**Table 2.** Summary of the salient properties of the clusters Abell 370, MACS J2135.2–0102, and RX J2129.6+0005 located within the skymaps of the GW events listed in Table 1.</caption><thead><tr><th></th><th>Abell 370</th><th>MACS J2135.2–0102</th><th>RX J2129.6+0005</th></tr></thead><tbody><tr><td>Cluster redshift</td><td>0.375</td><td>0.33</td><td>0.235</td></tr><tr><td>Right ascension</td><td>02h39m52s9</td><td>21h35m15s192</td><td>21h29m39s6</td></tr><tr><td>Declination</td><td>-01°34′36′′5</td><td>-01°03′01′′70</td><td>00°05′21′′2</td></tr><tr><td>M<sub>500</sub> [10<sup>14</sup> M<sub>☉</sub>]</td><td>10.6<sup>a</sup></td><td>7.6<sup>b</sup></td><td>3.5<sup>c</sup></td></tr><tr><td>L<sub>X</sub><sup>bol</sup> [10<sup>44</sup> erg s<sup>-1</sup>]</td><td>11.1<sup>a</sup></td><td>4.1<sup>d</sup></td><td>21.1<sup>c</sup></td></tr><tr><td>Einstein radius [arcsec]</td><td>39 (z<sub>s</sub> = 2)<sup>e</sup></td><td>38 (z<sub>s</sub> = 2)<sup>e</sup></td><td>18 (z<sub>s</sub> = 2)<sup>e</sup></td></tr></tbody></table>

Notes. $^a$Morandi, Ettori & Moscardini (2007)

$^b$Giacintucci et al. (2017)

$^c$Okabe & Smith (2016)

$^d$Hlavacek-Larrondo et al. (2012)

$^e$Richard et al. (2010)

## 2.2.1 GW190930.133541

The event (announced as S190930s) was detected by both LIGO L1 and H1 on 2019 September 30 at 13:35:41 UTC with an initial false alarm rate of 1/10.52 yr (GCN Circular 25871). The BAYESTAR pipeline (Singer & Price 2016) assigned 50 per cent probability of the sky localization spanning 706 deg², mass gap probability >95 per cent, and distance 752 ± 224 Mpc. We located four strong lensing clusters from our list within the 90 per cent skymap, including MACS J2135.2–0102 (redshift z = 0.33). This cluster lays on the contour enclosing p = 0.177 of the localization probability density and is the closest cluster to the peak of the skymap, corresponding to a region subtending 196.92 deg², and has an Einstein radius θ_E = 38 arcsec, the largest among the four clusters. Subsequently, the 50 per cent skymap and distance were updated to 536 deg² and 709 ± 191 Mpc by the LALINFERENCE pipeline (Veitch et al. 2015), in which only two strong lensing clusters were contained. MACS J2135.2–0102 lays on the contour enclosing p = 0.825 per cent corresponding to a region subtending Area = 1556.31 deg². This map was released prior to our observations, in which we have triggered Gemini-North observations totalling 1 h in the z' band (program ID:GN-2019B-Q-205). The GW event has been subsequently confirmed by several offline pipelines and included in the LIGO/Virgo source catalogue (Abbott et al. 2021b). The final skymap (50 per cent covering 569 deg², distance 770 ± 320 Mpc, source masses m₁ = 14.2M☉, m₂ = 6.9M☉) includes MACS J2135.2–0102 on the contour enclosing p = 0.680 per cent corresponding to a region subtending Area = 924 deg². The false alarm rate has been updated to be of ≈ 1/30 yr. Interestingly, the only other cluster consistently selected by our search scheme was RX J2129.6+0005, across the refinement process of the sky localization maps (see following section).

## 2.2.2 GW191216_213338

The event (announced as S191216ap) was detected on 2019 December 16 at 21:33:38 UTC by the online pipeline (BAYESTAR) analysing the data stream from LIGO H1 and Virgo. The event presented a false alarm rate of ≈1 in 10¹⁵ yr, a 50 per cent skymap encompassing 85 deg², mass gap probability >99 per cent, and distance 324 ± 78 Mpc. Our finder routine identified two strong lensing clusters within the 90 per cent probability of the event peak, of which RX J2129.6+0005 (redshift z = 0.23, Einstein radius θ_E = 17 arcsec) was closest to the peak and located within p = 0.717 corresponding to an area = 158.29 deg². An updated skymap from the LALinference pipeline showed the 50 per cent probability

to 68 deg², false alarm rate of ≈ 3 in 10¹⁵ yr, and a distance of 376 ± 70 Mpc. The cluster RX J2129.6+0005 was the only cluster left from our search at p = 0.876 corresponding to an area = 226.05 deg². Furthermore, at $\tilde{D} \simeq 300$ Mpc, this detection is close to the peak of the predicted population of lensed BNS mergers in LVK's third run (Smith et al. 2023). Following the GW trigger announcement, two GCN circulars confirmed the detection of a contextual neutrino by IceCube (GCN 26460) and a sub-threshold gamma-ray source by HAWC (GCN 26472) from a sky position offset from RX J2129.6+0005 (HAWC Collaboration 2019; IceCube Collaboration 2019). Therefore, we applied for additional Director Discretionary observing time at UKIRT to image the location of the HAWC detection, as detailed in Section 3.2.

## 2.2.3 GW200115_042309

The event (announced as S200115j) was detected on 2020 January 15 at 04:23:09.742 UTC with false alarm rate of ≈ 1 in 10³ yr and a mass gap probability > 94 per cent. The initial skymap produced by the online pipeline BAYESTAR was updated three times within a few hours of the initial event release. As suggested in the LIGO/Virgo GCN 26759, the third skymap release was the preferred choice presenting a 50 per cent sky probability of 186 deg². The mass gap probability was updated to >99 per cent and distance 331 ± 97 Mpc. Within this skymap, we identified three strong lensing clusters including Abell 370 (z = 0.38, Einstein radius: 45 arcsec). The cluster is the closest to the sky localization peak and located at p = 0.163 corresponding to an area = 31.53 deg². Due to the large Einstein radius of the cluster, together with the plethora of ancillary archival data which include deep multiband HST imaging, serving as reference image for our science case, we triggered a 1-h observation with Gemini South in the z' band (program ID: GS-2020A-Q-136). A subsequent skymap update from the LALinference resulted in a 50 per cent sky probability of 153 deg², distance 340 ± 79 Mpc, and false alarm rate ≈ 1 in 1500 yr, in which Abell 370 was the only cluster identified as located at p = 0.595, corresponding to an area = 222.76 deg².

# 3 DATA COLLECTION AND ANALYSIS

Our choice to use the Gemini telescopes is motivated primarily by the 5.5 sq. arcmin field of view of the GMOS instrument, which is large enough to cover the typical angular extent of the strong lensing region of the clusters in our sample, together with the capability of accessing both northern and southern sky. In addition, we used
<table><caption>**Table 3.** Follow-up observations of strong lensing clusters within sky localizations of the GW 190930-133541, 191216-213338, and 200115-042309.</caption><thead><tr><th>Visit</th><th>Start of observation (UTC)</th><th>Airmass<sup>a</sup></th><th>Integration time (ks)</th><th>Seeing<sup>b</sup> (arcsec)</th><th>Sensitivity<sup>c</sup></th></tr></thead><tbody><tr><td>1</td><td>2019 Oct 7, 07: 31: 44</td><td>GMOS-N observations of MACS J2135.2-0102<br>1.12<br>2.7</td><td></td><td>0.68</td><td>25.5</td></tr><tr><td>1</td><td>2019 Dec 19, 04: 47: 55</td><td>GMOS-N observations of RX J2129.6+0005<br>1.61<br>2.7</td><td></td><td>0.78</td><td>25.5</td></tr><tr><td>1</td><td>2020 Jan 16, 00: 58: 52</td><td>GMOS-S observations of Abell 370<br>1.27<br>2.7</td><td></td><td>0.76</td><td>25.4</td></tr></tbody></table>

Notes. $^a$The airmass at the mid-point of the observation.

$^b$Mean full width at half maximum of point sources in the reduced data.

$^c$5σ point source sensitivity within a photometric aperture of diameter 2 arcsec, estimated from the magnitude at which the photometric uncertainty is 0.2 mag.

Director Discretionary Time at the UKIRT telescope to access the WFCAM instrument, whose field of view is well matched to the size of the ICECUBE localization uncertainty of the neutrino detection discussed above.

## 3.1 Gemini observations

Our target of opportunity observing programmes at the Gemini Observatory covered the duration of LVK's third run, from 2019 April to 2020 May (Program IDs GN-2019B-Q-205, GN-2020A-Q-139, and GS-2020A-Q-136; PI: Bianconi). Each of these programs allowed for a single 1-h visit with the GMOS instrument on the Gemini-North and South telescopes. Table 3 summarizes the conditions and integration time of the individual observations. We have applied the same data reduction pipeline to each cluster observation, as described hereafter. Individual GMOS exposures were de-biased, dark-subtracted, flat-fielded, and de-fringed using Gemini DRAGONS Python package<fnref n="1" />, to produce a single science frame per exposure comprising the mosaiced individual chips. The individual exposure frames are stacked to produce a single frame per visit, after masking bad pixels. The full width at half-maximum (FWHM) of point sources in the reduced frames is consistently sub-arcsecond in all our observations (see Table 3). The central cut-outs of the GMOS images are presented in Fig. 2.

## 3.2 UKIRT observations

The GW trigger S191216ap was followed by a neutrino detection obtained by IceCube, whose sky localization is consistent with that of the GW event (IceCube Collaboration 2019). Subsequently, a sub-threshold gamma-ray detection was announced by the HAWC collaboration, with a sky localization consistent with both LIGO and IceCube (HAWC Collaboration 2019), but offset with respect to the position of the cluster RXJ 2129.6+0005, which had been targeted by our Gemini observations. Therefore, we identified a circle of radius 0.3 deg encompassing 68 per cent probability of the HAWC detection, and obtained Director Discretionary time observations (Program ID: U/19b/D05) with the WFCAM instrument on UKIRT through the z-band covering 0.75 sq. deg. We performed a 1-h long observation centred at RA 21:32:00 Dec. +05:13:48 on 2019 December 20, characterized by airmass ≈1.5, FWHM ≈ 1.4 arcsec, and a 5σ magnitude limit $m_z \approx 21.5$. A second 1-h visit was completed on the following day, but yielded lower magnitude limits with respect to the previous night, hence rendering the search for a

fading transient difficult. From the comparison between the Epoch 1 and Pan-STARRS1 archival data, we identified two candidate transients not associated with any known or candidate gravitational lens, one of which is located ≈7 arcsec from an edge-on galaxy at RA: 21:32:45.97 Dec:+5:19:57.0 (GCN 26605, Smith et al. 2019c). Further analysis highlighted the probable bogus nature of these candidates, due to detector-related artefacts.

## 3.3 Comparison with archival HST observations

The search for candidate lensed EM counterparts was performed by means of visual identification of new sources in proximity to the strong lensing regions of the clusters considered here. Archival HST WFC3-IR images through filters F105W and F110W were used as reference for the search of new transients in the GMOS frames. This search was performed by several participant of the collaboration promptly after the collection of the new data. The depth reached by the archival HST imaging is listed hereafter for each cluster at the 5σ level. Stacked HST images of Abell 370, MACS J2135.2-0102, and RX J2129.6+0005 reach $m \approx 29$ (Proposal ID: 14038, PI: Lotz), $m \approx 25.5$ (Proposal ID: 12166, PI: Ebeling, filter F110W), and $m \approx 27$ (Proposal ID: 12457, PI: Postman), respectively.

## 3.4 Photometric calibration

The angular extent and depth of the GMOS frames meant that there was no overlap between unsaturated bright stars in our field of view and measured in all-sky surveys. Therefore, we benchmarked the photometric calibration by measuring the $(g' - z')$ colours using archival SDSS catalogues and comparing it with model colours for massive early-type galaxies at the clusters' redshift. These model colours were computed using the EZGAL code<fnref n="2" />, and considering a single stellar population that formed at high redshift and evolved passively to the relevant cluster redshifts based on the Bruzual & Charlot (2003) populations. Formation redshift and the metallicity do not impact significantly the predicted colours. We obtained a consistent 5σ sensitivity level across each cluster observations reaching $m \approx 25.5$ in the z' band.

# 4 LENSING AND EM RADIATION MODELS

The following sections introduce the main aspects of the theoretical framework used to predict and interpret the signatures of EM lensed counterparts in our observations. We summarize the methods used to compute BNS rates as a function of mass and redshift, the delay time

<footnote n="1">¹https://dragons.readthedocs.io/projects/gmosimg-drtutorial/en/stable/index.html</footnote>

<footnote n="2">²www.baryons.org/ezgal</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 2. Central cut-outs of the Gemini GMOS z'-band images of Abell 370 (left panel), RX J2129.6+0005 (middle panel), and MACS 2135.2-0102 (right panel). North is up and east is to the left. Critical curves are plotted as red lines, and mark the location of infinite magnification for source redshifts of z = 1 and z = 2.</figcaption></figure>

expected for a GW event as a function of lens properties, and state-of-the-art model light curves for KNe. We refer the reader to Smith et al. (2023) and Nicholl et al. (2021) for a complete description of the models considered here.

## 4.1 GW progenitor model

The baseline model assumes a rate of BNS mergers as follows:

$$R(z, m) = R_0 g(z) f(m), \quad (1)$$

where $R_0$ is the comoving merger rate density of BNS mergers in the local universe (LVK2021). The functions $g(z)$ and $f(m)$ regulate the redshift evolution of GW sources, following the typical cosmic star formation rate density (Madau & Dickinson 2014), and the BNS mass function, respectively. The latter takes the shape of a top-hat function within the mass interval $1 < m$ [M$_{\odot}$] < 2.5, following LVK2021.

The number of lensed GWs per year arriving at Earth can be obtained by integrating the rate of BNS mergers $R(z, m)$ along redshift. In doing so, the differential source-plane optical depth is included, $\tau_{\mu}^{S}$, which describes the fraction of the source plane that is magnified via lensing. A lower limit on redshift $z_0$ is imposed and corresponds to the distance $D_0$ for a magnification yielding multiple images, i.e. $\mu = (D_0/\bar{D})^2 = 2$. Finally, the sensitivity of the GW detectors is taken into account (Martynov et al. 2016; Chen et al. 2021) in order to capture those lensed GWs that arrive at Earth that are detectable by current instruments.

The lensing model also permits predictions of the time delay between the arrival of multiply lensed GWs. Qualitatively, the number of images a source generates depends on the structure of the gravitational potential of the lens, and the alignment between source and lens caustics, i.e. theoretical surfaces corresponding to infinite magnification. Smith et al. (2023) showed that the pseudo catastrophe of the singular isothermal lens and the fold catastrophe bracket the range of time delays measured to date for lensed quasars, and thus form a solid basis for predicting time delays for other lensed transients. We therefore adopt the following convenient scaling relations for lenses with an isothermal slope:

$$\frac{\Delta t_{\text{SIS}}}{92 \text{ d}} = \left[ \frac{\theta_{\text{E}}}{\text{arcsec}} \right]^2 \left[ \frac{\mu_{\text{p}}}{4} \right]^{-1} \left[ \frac{D}{3.3 \text{ Gpc}} \right]; \quad (2)$$

$$\frac{\Delta t_{\text{fold}}}{3.9 \text{ d}} = \left[ \frac{\theta_{\text{E}}}{\text{arcsec}} \right]^2 \left[ \frac{\mu_{\text{p}}}{4} \right]^{-3} \left[ \frac{D}{3.3 \text{ Gpc}} \right], \quad (3)$$

where $\theta_{\text{E}}$ is the Einstein radius, $\mu_{\text{p}}$ is the combined lens magnification encountered by the two images, $D = D_{\text{L}}^{C} D_{\text{S}}^{C} / D_{\text{LS}}^{C}$, where $D_{\text{L}}^{C}$, $D_{\text{S}}^{C}$,

and $D_{\text{LS}}^{C}$ are the comoving distances from the observer to the lens, from the observer to the source, and from the lens to the source, respectively. Assuming values of lens redshift $z_{\text{L}} = 0.5$, which corresponds to the peak of optical depth to strong lensing (Robertson et al. 2020), and source redshift $z_{\text{S}} = 1.6$, which corresponds to the predicted peak of the true distances of the lensed BNS population (Smith et al. 2023), yields $D = 3.3 \text{ Gpc}$. We present an application of these estimates in Section 5.2, when we consider lenses for which the density profile matches that of an SIS at the Einstein radius.

## 4.2 Kilonova light curves

In parallel with our observing campaign, new models for KN light curves were being developed following the extensive follow-up efforts on KN AT2017gfo associated with GW170817. We take advantage of this progress to revisit the sensitivity of our observations. In particular, we use KN models from Nicholl et al. (2021) to predict the time evolution of the apparent magnitude of the EM counterpart to GW events. These models take as input parameters properties of GW mergers that can be recovered directly from the detected waveform, i.e. chirp mass, binary mass ratio, and inclination angle, as well as the equation-of-state-dependent parameters of tidal deformability and maximum NS mass.

We give a brief overview of Nicholl et al.'s (2021) models and refer the reader to their paper for full details. The EM emission from a KN event is constructed using five main ingredients. The first component describes the mass ejecta due to dynamical forces during the binary merger and is modelled using a highly opaque, equatorial component (labelled 'red' and responsible for the production of heavy elements). The second bluer, polar ejecta is caused by shocks at the binary contact surface and is characterized by a ≈20 times lower opacity compared to the red ejecta. The relative contribution of red and blue ejecta is most sensitive to the mass ratio of the binary system. The third component describes the mass ejecta that follows the binary merger, resulting from the interaction of neutrino winds with diffuse merger remnants, which are typically more massive than the pre-merger ejecta, and characterize by an opacity that lies between the red and blue dynamical ejecta. In summary, the geometry of the tri-phase ejecta is approximated by a sphere dominated by the red component for angles $|\theta| < 45$ deg from the orbital plane, and the lower opacity components carving a biconical polar cap orthogonal to the orbital plane. A fourth model component accounts for enhanced emission of blue ejecta due to magnetically driven winds. The total output luminosity driven by r-process decay is proportional to time
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Distributions of luminosity distance and mass for BNS models, and for detected GW events in O3, following Smith et al. (2023). The filled points mark the median posterior mass-distances of the individual binary components that LVK infer assuming $\mu = 1$ for the O3 detections. The red points mark the events classified as mass gap with low latency by LVK. The large red points mark the individual component masses of the progenitors of the GW events considered in this study. The dashed contours show the intrinsic mass-distance distribution expected for the lensed BNS population detected by LIGO-like instruments, while the solid contours show the inferred mass-distance distribution (assuming no lensing) for the same population. In each case the thicker (inner) and thinner (outer) contours encircle 50 and 90 per cent of the predicted magnified population, respectively. The grey horizontal bands show the mass range for typical NSs and stellar BHs.</figcaption></figure>

$t^{-1.3}$, and is computed considering each ejecta component mass, together with its velocity and opacity, which is assumed constant. An additional fifth component allows for the onset of a gamma-ray burst jet which can be responsible for additional shock-heating of the ejecta.

# 5 RESULTS

We have continued the deep search for EM counterparts to candidate lensed GW events started by Smith et al. (2019a). Our campaign targeting the appearance of new transient phenomena within the central, strong lensing, regions of three galaxy clusters yielded no significant candidates. Due to the inclusion of prescriptions for BNS lensing and the state-of-the-art KN light-curve models we are able to further inform our strategy and findings, and discuss the implications of lensing on the posterior parameters from LVK.

Fig. 3 shows the true and lensed distribution of mass and distance for the model BNS population discussed in Section 4.1. We note that lensed BNS are predicted to be located primarily within the boundaries of the mass gap interval, confirming this as a viable mechanism responsible for the detection of binary components with masses inferred to be in the range $2.5 < m [M_\odot] < 5$. The median of the posteriors of mass and distance for individual binary components from confirmed GW events from LVK's third run are marked as filled circles, and include the events considered in this study in red (see Section 2.2). We note that the final mass posteriors released from LVK for the mass gap events considered here have values outside of

the mass gap. This is due to the different offline pipeline deployed by LVK, with respect to the ones used for the initial detection and classification of the GW trigger. This strengthens the case for an improved analysis of the GW events with low latency from LVK, which would yield robust mass posteriors informing the selection of the most promising events to follow-up. Nevertheless, we note that the predicted lensed population of BNS extend beyond the mass gap, and that the mass-distance posteriors of several progenitors of the events considered here are consistent with the lensing scenario (see overlap between red points and solid contours in Fig. 3).

Lensing allows us to reconcile the mass posterior from LVK to that expected for the GW sources considered here by adjusting the true redshift of the merger. We adopt a true chirp mass $M = (m_1m_2)^{3/5}(m_1 + m_2)^{-1/5} = 1.18 M_\odot$, which corresponds to the best fit describing the KN AT2017gfo observed in conjunction with GW170817 (Nicholl et al. 2021). Following $M(1+z) = \bar{M}(1+z)$, the true redshifts of GW 190930_133541, GW191216_213338, and GW 200115_042309 would be $z = [6.0, 6.0, 0.9]$, respectively. We note that this would yield extreme magnification ($\mu > 3000$) for GW 190930_133541 and GW191216_213338, compared to the predicted values for O3, whose distribution peaks at $\mu \approx 1000$ (Smith et al. 2023). Therefore, we also consider a higher intrinsic chirp mass $M[M_\odot] = 2$ for these two events, which translates into individual binary masses $m_1 = m_2 = 2.3 M_\odot$ when considering equal mass ratio merger, and a true redshift $z \approx 3.0$. This increase in true intrinsic chirp mass is within the range of NS masses inferred by LVK2021.

## 5.1 Expected light curves

Models for light curves require the input of parameters describing the BNS merger and remnants. As delineated above, we use the best-fitting parameters from the AT2017gfo counterpart to GW170817 to inform our choice of parameters. This is so as to produce an empirically motivated set of light curves. Specifically, we adopt an intrinsic chirp mass $M = 1.18 M_\odot$, mass ratio $q = 0.92$, a shocked cocoon opening angle $c = 24$ deg, blue ejecta enhanced by a factor 1.6, post-merger ejecta of 12 per cent of the total disc mass, and a viewing angle of $\theta = 32$ deg, to our reference light-curve model. We produce two additional models: one with the same parameters but a larger chirp mass $M = 2.0 M_\odot$; and a 'conservative' model with typical NS binary parameters $M = 1.18 M_\odot$, $q = 0.9$, viewed at $\theta = 60$ deg, with no shock cooling or blue ejecta enhancement.

Fig. 4 shows the different model KN light curves for each of the three GW events considered in this study. Overall, we note that GW170817-like models predict a detectable transient with our follow-up strategy, together with highlighting the necessity of quick activation given the rapid fading of the transient magnitude within $\simeq 4$ d of the merger. Both GW 190930_133541 and GW 191216_213338 present short spans of visibility, due to the high true redshift required to reconcile the posterior chirp mass to values compatible to that of BNS. Hence, the high-magnification results in short, high-luminosity peaks associated with rest-frame ultraviolet emission that rapidly fades in the observer frame notwithstanding the cosmological time dilation. Light curves for GW200115_042309 allow visibility up to 8 d with the current magnitude limits from our campaign, benefitting from relatively low true redshift ensuring that the observations would have probed the rest-frame optical emission around 500 nm. A common feature of all the light curves presented here is that their duration, defined as the time required to fade by a factor 2 in flux, is below 2 d. This makes them the fastest fading transients and suggest short duration as key signature of lensed KNe (for additional details see Smith et al. 2023).
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Predicted evolution of z'-band magnitudes for KN emission following the models by Nicholl et al. (2021) for each of the GW events considered here. The reference model, labelled KN170817, is plotted as yellow lines and uses the best-fitting parameters to the KN observed in conjunction with GW 170817 (Nicholl et al. 2021). In addition, we test light curves using the reference model parameters and $M = 2 M_\odot$, plotted as dot-dashed lines. A conservative model is included for comparison as teal lines. For this, we assume a merger with chirp mass $M = 1.18 M_\odot$, mass ratio $q = 0.9$, viewing angle $\theta = 60$ deg, and no additional blue ejecta from magnetic winds, or shock cooling emission. The magnitude values are boosted according to the ratio between the LVK posterior distances, which assume no lensing, and a range of redshifts $z_{\text{true}} \in [0.9, 3.0, 6.0]$, marked with dashed, solid, and dot-dashed line, respectively, and assumed as the true location of the binary merger. The grey area marks magnitudes beyond the reach of our observations.</figcaption></figure>

We note that the GW170817-like model is brighter than the conservative model by 2–3 mag, owing to the combined effects of larger blue ejecta mass, favourable viewing angle, and cooling emission from the shocked cocoon. These two models therefore provide optimistic (but plausible) and more pessimistic scenarios for EM counterparts. When considering the implication of varying the true redshift of the merger, we note that the steep decrease of the blue component of the KN, which is responsible for the luminosity peak, is longer lasting in the observer frame than the rest frame due to cosmological time dilation. This effect delays the onset of the light-curve red component, and accentuates the fading of the observed magnitudes. In this respect, the red component is observable for the closer event considered here GW 200115_042309 when assuming a true redshift $z = 0.9$. Regarding GW 190930_133541 and GW 191216_213338, using $M = 2 M_\odot$ results in a decrease of the overall luminosity due to the significant reduction of merger ejecta, due to the rapid collapse of the massive merger product to a BH (Nicholl et al. 2021).

## 5.2 Multiple image interpretation

We now turn to a more detailed discussion of GW190930_133541 and GW191216_213338. The overlap between the skymaps associated with these detections, in which cluster RX J2129.6+0005 is included, motivated us to discuss the hypothesis that they are the manifestations of a single, gravitationally lensed event. In addition, the median values of merger mass ratio $q$, which is lensing-invariant, are consistent between the two events, but are typically not well constrained (Abbott et al. 2021b). While this does not provide conclusive evidence of strong lensing, we utilize the difference between the arrival time of these two GWs (≈78 d) to explore the hypothesis that they are strongly lensed images of the same source.

The ratio between these events' peak posterior distances is $\tilde{D}(\text{GW190930})/\tilde{D}(\text{GW191216}) \approx 2$, which implies a magnification ratio of $\mu(\text{GW191216})/\mu(\text{GW190930}) \approx 4$. This is firmly in the regime of low-magnification strong lensing, in which isothermal galaxy-scale lenses are more efficient than cluster-scale lenses (Smith

et al. 2023). We therefore adopt the following expression under the assumption that the putative lens is close to isothermal:

$$ \mu_{\pm} = 1 \pm \theta_{E} \beta^{-1}, \qquad (4) $$

where $\beta$ is the source position and the respective images are denoted ‘+’ and ‘−’. This enables us to estimate the individual magnifications suffered by the two images as $\mu_+ \approx 2.67$ and $|\mu_-| \approx 0.67$, and thus a combined magnification of $\mu_p \approx 3.3$. This shows that GW 190930_133541 and 191216_213338 are de-magnified and magnified, respectively. The true distance of the source can be rewritten as

$$ D = \tilde{D}_{+} \mu_{+}^{1/2} = \tilde{D}_{-} |\mu_{-}|^{1/2}, \qquad (5) $$

which yields $D \approx 600$ Mpc which corresponds to $z \approx 0.13$. Given the low redshift, the GW posteriors are only marginally affected. Then substituting $\mu_p \approx 3.3$, $\Delta t = 78$ d, and $D \approx D \approx 600$ Mpc in to equation (2), we obtain a $\theta_E \approx 2$ arcsec. The following relation for an isothermal lens

$$ \theta_E = \frac{4\pi\sigma^2}{c^2} \frac{D_{LS}}{D_S} \qquad (6) $$

then allows us to estimate the velocity dispersion of the lens as $\sigma = 370$ km s⁻¹, assuming that $D_{LS} = D_S/2$. This value is typical of small group-sized haloes that are commonly inhabited by a massive early-type galaxy – i.e. a fairly typical low magnification strong lens configuration. Repeating this exercise to characterize the properties of a lens producing more than two images using equation (3) yields an Einstein radius is $\theta_E = 7.8$ arcsec and $\sigma = 735$ km s⁻¹, which is typical of a cluster-scale lens, and thus less plausible given the low magnifications at play here (Smith et al. 2023).

## 6 CONCLUSIONS

The steady progress in the detection of GWs revealed growing inconsistency in the understanding of the mass regimes accessible to compact objects. The detection of mass gap objects during LVK's third run highlights this tension. The intervention of gravitational
lensing offers a physically well-understood solution to the detection of these GW events, without invoking modifications to the current models of stellar evolution. To test this hypothesis, we have selected promising mass gap events detected during LVK's third run, and performed a search for EM counterparts to explore if they were affected by strong lensing. In particular, we have selected three GW events as containing at least one component whose mass falls within the mass gap $2.5 < m [M_\odot] < 5$, and within hours of their detection, selected and observed one galaxy cluster within each event localization map using Gemini-GMOS. The data obtained reach magnitude $m_z \approx 25.5$ at $5\sigma$, and we detect no clear transient within the data limit, when comparing the newly obtained observations to deep HST archival imaging. This is the continuation of the deep pilot search for EM transients to GW events started by Smith et al. (2019a) in anticipation of the upcoming opportunities afforded by the Vera Rubin Observatory.

Due to state-of-the-art light-curve models for KN-like transients, we are able to confirm that EM counterparts to lensed GW events would have been detectable due to lensing magnification using 8-m class ground-based telescope within days from the binary mergers. Lensing arguments highlight that posterior values of chirp masses up to $\mathcal{M} \approx 9M_\odot$ can be compatible with that of BNS mergers. Furthermore, we tested the scenario in which GW 190930_133541 and 191216_213338 are multiple images of a single, multiply imaged gravitationally lensed event. We compared the arrival time difference between the detection of these two events with theoretical predictions of arrival time differences and found that lensing arguments allow for both events to be interpreted as strongly lensed images of an individual GW source multiply imaged by a modest group-scale lens. Given the low magnifications involved, the true masses of the compact objects involved in the merger would be consistent with those inferred by LVK under their assumption that lensing was not at play.

Looking to the future, significant progress in the search for lensed BNS will arrive with the increase in GW detector sensitivity prospected for LIGO's fifth run, with predicted lensed rates around one per year (Smith et al. 2023). Our candidate lensed BNS strategy is effective in selecting the most promising GW events affected by lensing, and offers physically well-motivated and sensitive observations of any false positives. Additional progress can be achieved by having access to robust event posteriors on mass and mass ratios with low latency, i.e. as soon as possible after the detection of a GW event, and by increasing the census of lines of sight with high lensing probability (Ryczanowski et al. 2023). This is particularly important in light of LVK's upcoming fourth run (early 2023), which features increased sensitivity, and for which the Vera Rubin Observatory will not be accessible yet. This is essential to select more efficiently the most promising candidates, inform the search for lensing intervention, and focus on the observational follow-up efforts.

## ACKNOWLEDGEMENTS

MB and GPS acknowledge support from STFC grant numbers ST/N021702/1 and ST/S006141/1. DR and ER acknowledge STFC studentships. GPS acknowledges support from The Royal Society and the Leverhulme Trust. MN is supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 948381) and by a Fellowship from the Alan Turing Institute. We thank the staff and Directors of the Gemini and UKIRT Observatories and the Cambridge Astronomy Survey Unit for their superb support

of our observing programmes. UKIRT is owned by the University of Hawaii (UH) and operated by the UH Institute for Astronomy. When (some of) the data reported here were obtained, the operations were enabled through the cooperation of the East Asian Observatory. Based on observations obtained at the international Gemini Observatory, a program of NSF's NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation on behalf of the Gemini Observatory partnership: the National Science Foundation (USA), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea).

## DATA AVAILABILITY

The Gemini/GMOS data set used in this work is available at the Gemini Observatory Archive. The UKIRT/WFCAM data set is available upon request to the authors.

## REFERENCES

* Abbott B. P. et al., 2016, Phys. Rev. Lett., 116, 061102

* Abbott B. P. et al., 2017a, Phys. Rev. Lett., 118, 121101

* Abbott B. P. et al., 2017b, Phys. Rev. Lett., 119, 141101

* Abbott B. P. et al., 2017c, ApJ, 848, L12

* Abbott R. et al., 2021a, preprint (arXiv:2108.01045)

* Abbott R. et al., 2021b, Phys. Rev. X, 11, 021053

* Abbott R. et al., 2021c, ApJ, 923, 14

* Alsing J., Silva H. O., Berti E., 2018, MNRAS, 478, 1377

* Arcavi I., 2018, ApJ, 855, L23

* Belczynski K., Wiktorowicz G., Fryer C. L., Holz D. E., Kalogera V., 2012, ApJ, 757, 91

* Bouwens R. J. et al., 2014, ApJ, 795, 126

* Bruzual G., Charlot S., 2003, MNRAS, 344, 1000

* Chen H.-Y., Holz D. E., Miller J., Evans M., Vitale S., Creighton J., 2021, Class. Quantum Gravity, 38, 055010

* Farr W. M., Sravan N., Cantrell A., Kreidberg L., Bailyn C. D., Mandel I., Kalogera V., 2011, ApJ, 741, 103

* Giacintucci S., Markevitch M., Cassano R., Venturi T., Clarke T. E., Brunetti G., 2017, ApJ, 841, 71

* Graham M. J. et al., 2020, Phys. Rev. Lett., 124, 251102

* Gupta A., Gerosa D., Arun K. G., Berti E., Farr W. M., Sathyaprakash B. S., 2020, Phys. Rev. D, 101, 103036

* HAWC Collaboration, 2019, GRB Coordinates Network, 26472, 1

* Hlavacek-Larrondo J., Fabian A. C., Edge A. C., Ebeling H., Sanders J. S., Hogan M. T., Taylor G. B., 2012, MNRAS, 421, 1360

* IceCube Collaboration, 2019, GRB Coordinates Network, 26460, 1

* Jauzac M. et al., 2016, MNRAS, 463, 3876

* Kelly P. L. et al., 2018, Nat. Astron., 2, 334

* Kneib J. P., Ellis R. S., Smail I., Couch W. J., Sharples R. M., 1996, ApJ, 471, 643

* Kneib J.-P., Ellis R. S., Santos M. R., Richard J., 2004, ApJ, 607, 697

* Lagattuta D. J. et al., 2017, MNRAS, 469, 3946

* Madau P., Dickinson M., 2014, ARA&A, 52, 415

* Mahler G. et al., 2018, MNRAS, 473, 663

* Martynov D. V. et al., 2016, Phys. Rev. D, 93, 112004

* Morandi A., Ettori S., Moscardini L., 2007, MNRAS, 379, 518

* Nicholl M., Margalit B., Schmidt P., Smith G. P., Ridley E. J., Nuttall J., 2021, MNRAS, 505, 3016

* Okabe N., Smith G. P., 2016, MNRAS, 461, 3794

* Olejak A., Fryer C. L., Belczynski K., Baibhav V., 2022, MNRAS, 516, 2252

* Planck Collaboration XIII, 2016, A&A, 594, A13
Richard J. et al., 2010, MNRAS, 404, 325

Robertson A., Smith G. P., Massey R., Eke V., Jauzac M., Bianconi M., Ryczanowski D., 2020, MNRAS, 495, 3727

Ryczanowski D., Smith G. P., Bianconi M., McGee S., Robertson A., Massey R., Jauzac M., 2023, MNRAS, 520, 2547

Singer L. P., Price L. R., 2016, Phys. Rev. D, 93, 024013

Smith G. P., Kneib J.-P., Smail I., Mazzotta P., Ebeling H., Czoske O., 2005, MNRAS, 359, 417

Smith G. P., Jauzac M., Veitch J., Farr W. M., Massey R., Richard J., 2018, MNRAS, 475, 3823

Smith G. P. et al., 2019a, MNRAS, 485, 5180

Smith G. P., Bianconi M., Massey R., Robertson A., 2019b, GRB Coordinates Network, 24250, 1

Smith G. P., Nicholl M., Sharon K., Bianconi M., Varricatt W. P., Benigni S., Ridley E. J., Gravitationally Wave Gravitational Lensed Hunters, 2019c, GRB Coordinates Network, 26605, 1

Smith G. P. et al., 2023, MNRAS, 520, 702

The LIGO Scientific Collaboration, 2021, preprint (arXiv:2111.03634)

Veitch J. et al., 2015, Phys. Rev. D, 91, 042003

Villar V. A. et al., 2017, ApJ, 851, L21

Wang Y., Stebbins A., Turner E. L., 1996, Phys. Rev. Lett., 77, 2875

Welch B. et al., 2022, Nature, 603, 815

This paper has been typeset from a *TEX/LATEX* file prepared by the author.