A fast and precise methodology to search for and analyse strongly lensed gravitational-wave events

Gravitational waves, like light, can be gravitationally lensed by massive astrophysical objects such as galaxies and galaxy clusters. Strong gravitational-wave lensing, forecasted at a reasonable rate in ground-based gravitational-wave detectors such as Advanced LIGO, Advanced Virgo, and KAGRA, produces multiple images separated in time by minutes to months. These images appear as repeated events in the detectors: gravitational-wave pairs, triplets, or quadruplets with identical frequency evolution originating from the same sky location. To search for these images, we need to, in principle, analyze all viable combinations of individual events present in the gravitational-wave catalogs. An increasingly pressing problem is that the number of candidate pairs that we need to analyse grows rapidly with the increasing number of single-event detections. At design sensitivity, one may have as many as $\mathcal O(10^5)$ event pairs to consider. To meet the ever-increasing computational requirements, we develop a fast and precise Bayesian methodology to analyse strongly lensed event pairs, enabling future searches. The methodology works by replacing the prior used in the analysis of one strongly lensed gravitational-wave image by the posterior of another image; the computation is then further sped up by a pre-computed lookup table. We demonstrate how the methodology can be applied to any number of lensed images, enabling fast studies of strongly lensed quadruplets.


INTRODUCTION
Gravitational waves (GWs) result from cataclysmic events that distort the fabric of space and time. The Advanced LIGO (Aasi et al. 2015) and Advanced Virgo (Acernese et al. 2015) detectors have found dozens of GW signals emitted by the mergers of binary neutron stars or black holes (Abbott et al. 2019;Abbott et al. 2020). Advanced LIGO and Virgo are scheduled to be upgraded further. Meanwhile, Japan's KAGRA (Somiya 2012;Aso et al. 2013;Akutsu et al. 2019;Akutsu et al. 2020) has become online. As the sensitivities of these instruments improve, many novel avenues in research could become observationally feasible. One such avenue could be the study of GW lensing.
When GWs travel from their source to the Earth, they can be gravitationally lensed by intervening objects (galaxies, galaxy clusters, and other compact objects) (Ohanian 1974;Deguchi & Watson 1986;Wang et al. 1996;Nakamura 1998;Takahashi & Nakamura 2003;Ezquiaga et al. 2021;Oguri 2018;Liu et al. 2021;Lo & Magana Hernandez 2021). Lensing could produce several observable effects on GW signals detectable by ground-based GW detectors. A small fraction of binary black hole mergers will be strongly lensed by intervening galaxies Ng et al. 2018;Li et al. 2018;Oguri 2018) or galaxy clusters (Smith et al. 2018(Smith et al. , 2017Smith et al. 2019;Robertson et al. ★ E-mail: j.janquart@uu.nl 2020; Ryczanowski et al. 2020), producing multiple images detectable as repeated events (Wang et al. 1996;Haris et al. 2018). The images can be magnified by the lens, resulting in a biased luminosity distance measurement Ng et al. 2018;Pang et al. 2020), or inverted along one or both principal axes, resulting in an overall shift in the GW phase Ezquiaga et al. 2021), and they arrive a few minutes up to years apart. However, because the GW wavelength is typically much smaller than the lens size, lensing does not affect the signal's frequency content, referred to as the geometrical optics limit (Takahashi & Nakamura 2003). Thus, GWs strongly lensed by galaxies or galaxy clusters appear as repeated events with identical frequency evolution and sky location but separated in time. If smaller compact objects (stars, black holes) intervene with the GW, microlensing can also occur. In that case, the geometrical optics approximation can break down, and frequency-dependent "beating patterns" to the waveform appear (Takahashi & Nakamura 2003;Cao et al. 2014;Lai et al. 2018;Christian et al. 2018;Jit Singh et al. 2018;Hannuksela et al. 2019;Meena & Bagla 2020;Pagano et al. 2020;Cheung et al. 2021;Kim et al. 2020).
If identified, lensed GWs could give rise to several exciting possibilities. For example, they might enable us to locate merging black holes at a sub-arcsecond precision by matching the GW image properties with the properties of the lenses discovered in electromagnetic surveys when quadruple images are available (Hannuksela et al. 2020). Similar localisation might be possible with galaxy cluster lensing even when only a signal double is found, owing to the rarity of the cluster lenses (Smith et al. 2017;Ryczanowski et al. 2020) (see also Refs. (Sereno et al. 2011;Yu et al. 2020)). When accompanied by an electromagnetic counterpart, they could enable precision cosmography studies owing to the sub-millisecond lensing time-delay measurements granted by GW observations (Sereno et al. 2011;Liao et al. 2017;Cao et al. 2019;Li et al. 2019b). Additionally, lensed GWs might enable improved GW propagation tests by comparing the lensing time delay of strongly lensed waves with their transient electromagnetic counterparts to measure the speed of gravity relative to the speed of light (Baker & Trodden 2017;Fan et al. 2017). Moreover, strongly lensed events allow us to detect the same event multiple times at different detector orientations, effectively multiplying the number of detectors by the number of images to arrive at an enlarged (synthetic) detector network; this could be exploited to probe the full GW polarisation content, including alternative polarizations (Goyal et al. 2021). Another prospective avenue is detecting intermediate-mass and primordial black holes through microlensing observations (Lai et al. 2018;Jung & Shin 2019;Diego 2020;Oguri & Takahashi 2020).
The principal idea to search for strong lensing is to locate event pairs (and from there, possibly triplets and quadruplets) with similar detector-frame parameters arriving at the detector from the same sky location. To assess if an event is strongly lensed, we compare the likelihood that an event pair is lensed against the likelihood that they were produced by astrophysical coincidence. Two parameterestimation-based approaches exist to do this. The first one is the posterior-overlap methodology, which performs a Gaussian kernel density estimation (KDE)-based fit on the single event posterior density functions. It then tests if any given event pair is consistent with lensing by assessing the consistency of the posteriors (Haris et al. 2018). This approach allows for rapid tests on large quantities of data. The second approach, joint parameter estimation analysis, abandons the KDE-based fits in favour of sampling the full joint likelihood, enabling improved accuracy (Liu et al. 2021;Lo & Magana Hernandez 2021).
An increasingly pressing problem is the rising computational demand of these strong lensing analyses. As an order-of-magnitude estimate, at design sensitivity, the number of observed GW events is likely to reach O (10 3 ) (Ng et al. 2018),which results in ∼ 5 × 10 5 event pairs. Even after dropping event pairs which show significant mismatch in the inferred parameters, there would be a large number of event pairs that need to be analysed. To address the rising computational demand, we develop a new method: A fast and precise way to analyse strongly lensed GW event pairs by using the posterior of one event of the event pairs as a prior for the analysis of the other event in an importance sampling procedure. We further accelerate the methodology with a fast likelihood computation. The methodology fills the gap between posterior-overlap (Haris et al. 2018) and joint parameter estimation (Liu et al. 2021;Lo & Magana Hernandez 2021) methodologies in terms of speed and precision, enabling fast but still relatively accurate strong lensing analyses.
We summarise the relevant strong lensing theory and how it manifests itself in the detectors in Sec. 2. We describe the Bayesian framework and the approximations involved in our methodology to study strongly lensed events in Secs. 3, 4, and 5. We provide an example of a typical use case by analysing a non-lensed and a lensed pair in Sec. 6. We detail how one can use our approach to analyse multiple images in Sec. 7 and demonstrate an example application in Sec. 8. A summary, conclusions, future outlook are given in Sec. 9.

STRONGLY LENSED GRAVITATIONAL WAVES
Strong lensing splits the GWs into multiple potentially observables images, which are categorised by their image type. Type-I images correspond to the minimum of the Fermat potential and leave the overall shape of the GWs unaffected. Type-II images are saddle points of the potential and Hilbert-transform the GWs. Finally, Type-III images correspond to the maximum of the Fermat potential and invert the original image and thus the waveform. Such images are typically captured by an overall phase shift of the waveform, referred to as the Morse phase Ezquiaga et al. 2021). Besides the potential image inversions and Hilbert transforms, the images can also be magnified.
The lens is often a galaxy or galaxy cluster. Galaxy lensing typically produces two or four Type-I/II bright images (although a fifth, Type-III bright central image can be observed in rare cases when the density slope of the galaxy is very shallow (Collett & Bacon 2016;Dahle et al. 2013;Collett et al. 2017)) separated by a few minutes to months (Ng et al. 2018;Haris et al. 2018).
Compared to galaxy lenses, galaxy cluster lenses can have much more complex lens morphologies, and therefore can produce a much richer spectrum of images, separated by up to years (Smith et al. 2018(Smith et al. , 2017Smith et al. 2019;Robertson et al. 2020;Ryczanowski et al. 2020). Irrespective of the lens configuration, each lensed image will arrive at the detector as a Type-I, Type-II, or Type-III image, separated by times , and have their intensities magnified by factors . Thus, the search approach is usually phenomenological and agnostic to the specific lens configuration that produced the images.
However, we note that the lensed time delays and image types can be forecasted statistically, at least for galaxy lenses (Oguri 2018;Haris et al. 2018). Such statistical information would improve the discriminatory power of our searches, and we expect to pursue this in future work.
Consequently, in the frequency domain, the GW waveform of the ℎ image is modified such that it arrives with a time delay , experiences an overall magnification due to the focusing by lensing, 1 and can be inverted/Hilbert-transformed, an effect captured in an overall complex "Morse factor" (phase shift) Ezquiaga et al. 2021): whereh ( ; ) is the frequency-domain waveform in the absence of lensing with the usual set of binary parameters , and ℎ ( ; , , , ) is the lensed frequency-domain waveform; for an illustration of a lensed waveform, see Fig. 1. Note that the image parameters enter the GW waveform in an overall multiplicative factor, which can easily be decoupled from the rest of the waveform. For brevity, we will denote the individual image properties (magnifications, time delays, Morse factors) by The lens magnifies the two waves differently, and that is why their amplitudes are different. One of the images has also been inverted by a lens, and thus there is a relative phase difference between the waveforms. Normally, the two images would arrive at different times. However, here the two waveforms have been time-shifted to superpose them.
. An overview of the notations used in this work can be found in Table 1.
We can absorb the 1/2 into an observed luminosity distance obs,j and the time-delay due to lensing into an observed coalescence time obs,j , such that where and are the true luminosity distance and coalescence time.

THE STRONG LENSING HYPOTHESIS
Under the lensed hypothesis H , when a GW from an image enters an interferometer, the data stream where ( ) is the detector noise contribution, and ℎ ( ; , ) is the lensed GW as seen by the detectors. For two observed lensed images, the two different data streams containing a lensed signal are connected through the binary parameters , such that the joint evidence (neglecting selection effects) (Liu et al. 2021;Lo & Magana Hernandez 2021) where ( 1 | , 1 ) and ( 2 | , 2 ) are the individual likelihoods (Veitch & Vecchio 2010) and ( , 1 , 2 ) is the prior. Under the usual non-lensed hypothesis H , the data streams take the form where are the binary parameters of the ℎ waveform, and ℎ ( ; ) is the waveform projected onto the detector frame. Under the unlensed hypothesis, the data streams contain unrelated signals, such that the joint evidence (Liu et al. 2021;Lo & Magana Hernandez 2021) ( (7) where ( 1 ) and ( 2 ) are the usual priors. Therefore, to test the lensed hypothesis, we adopt the ratio of evidences, or "coherence ratio" which describes how (dis)similar the two signals are. Note that here we do not call C a Bayes factor because it does not include selection effects (Lo & Magana Hernandez 2021) and it is sensitive to the binary black hole population prior. 2

THE CONDITIONED EVIDENCE
In this work, instead of evaluating the full joint evidence above, we evaluate the conditioned evidence and the individual evidences to obtain the coherence ratio. This way we can accelerate the computation of the conditioned evidence using importance sampling and a lookup table, as will be explained below. First, we can rewrite the joint evidence as (see Appendix A) Defining the relative magnification between images rel , the relative time delay between images 21 , and the relative Morse factor between images 21 = 2 − 1 , the conditioned evidence where = { rel , 21 , 21 } are the relative image properties and are the effective parameters which absorb the lensing magnification into an observed luminosity distance and the lensing time delay into an observed coalescence time (see Eq. (2)) and includes also the Morse factor of the first image 1 . That is, the usual prior in the evaluation of the second event has been replaced by the posterior of the first event ( | 1 ) and the first event is fully described by . The second event can be related to the first event through the difference in the image properties obs,2 = objs,1 + 21 , In terms of the conditioned evidence, the coherence ratio Thus, instead of evaluating the full joint evidence, we can evaluate the conditioned evidence and the individual evidences to obtain the coherence ratio.
The main desirable attribute of the conditioned evidence is that the integral converges faster than a regular joint parameter estimation run, as the posterior of the first event ( | 1 ), which replaces the usual prior in the conditioned evidence, is already concentrated around the relevant, shared detector-frame parameters. Indeed, the integral has reduced partially into an importance sampling problem over the shared waveform parameters.

EVALUATING THE CONDITIONED EVIDENCE
There are, in principle, several ways to solve the conditioned evidence. Here we solve the conditioned evidence is by re-writing the conditioned evidence in terms of a "marginalised" likelihood where the likelihood of the second event averaged over the posterior samples of the first event. 3 Evaluating the conditioned evidence instead of the full joint evidence is desirable when it comes to speed, as the posterior ( | 1 ) is already concentrated around the relevant waveform parameters, which accelerates the convergence of the integral. Furthermore, by recycling the posterior samples of the first event 4 , we can evaluate the marginalised likelihood without having to generate new trial GW waveforms, through a lookup table. The lookup table further accelerates the computation (see Appendix B for details).
Once we have the conditioned evidence, we have access to the posterior distributions of the lensing parameters ( | 1 , 2 ). We can then reweigh the posterior to obtain (Appendix A) This way, we can recover the joint parameters given by the combined observation; see Fig. 3 for an illustration. We note that Eqs. 10 and 13 are exact relationships. However, in evaluating them, we have limited the number of samples drawn from the posterior of the first event, to make the construction of the lookup table for each sample feasible memory-wise.
We have implemented our methodology as a module to the B GW parameter estimation tool (Ashton et al. 2019), and named it G (Gravitational-wave analysis Of Lensed and Unlensed waveform Models). Since we apply the Morse factor to the overall waveform as in Eq. (1), we can use IMRP P 2 (Khan et al. 2016b) (a precessing-spin waveform model), IMRP D (Khan et al. 2016a) (an aligned-spin waveform model) and other waveforms with precession or higher order modes in our analysis. In our trial analyses, a single evaluation of the conditioned evidence 3 Choosing 1 to be the event with a better-constrained posterior of the two events will ensure faster convergence of the marginalized likelihood. 4 To be able to use the posteriors of the first event, the analysis of the first image has been adapted to incorporate Morse factor as a discrete parameter takes less than O (1) CPU hour 5 , adopting the P M sampler (Feroz et al. 2009(Feroz et al. , 2019.

EXAMPLE ANALYSIS
As a practical example, we inject the signal from a spinning, precessing binary black hole merger generated with IMRP P 2, with parameters listed in the second column of Table 2, into synthetic stationary, Gaussian noise for a network of 3 detectors (LIGO-Livingston, LIGO-Hanford, and Virgo) at design sensitivity (Barsotti et al. 2021;Acernese et al. 2015). This event has a network SNR of ∼ 23. We then inject the event's lensed counterpart image, with relative magnification between the two events rel = 2, relative Morse factor 21 = 0.5, and relative time delay 21 = 14 hr. Throughout this analysis, we use a uniform prior for the relative magnification ( rel ∈ [0.01, 20]), the time-delay ( 21 ∈ [ 21 − 0.1, 21 + 0.1] s), the chirp mass (M ∈ [10, 100] , with M = ( 1 2 ) 3/5 /( 1 + 2 ) 1/5 where 1,2 are the component masses), and mass ratio ( ∈ [0.1, 1], with = 2 / 1 ); the spin distribution is isotropic. The prior on the Morse factor is a discrete uniform distribution over the three possible values ( 1 ∈ {0, 0.5, 1}) and the prior in Morse factor difference is a discrete uniform distribution over the four possible values ( 21 ∈ {0, 0.5, 1, 1.5}) 6 . These priors do not include results from lens modelling of the inferred astrophysical population of binary black holes.
To analyse this lensed pair, we perform four nested sampling runs. Firstly, we analyse the two injections under the non-lensed hypothesis. Secondly, we estimate the parameters of one of the events under the lensing hypothesis. Thirdly, we obtain the conditioned evidence ( 2 | 1 , H ), by sampling the second event's likelihood based on the earlier lensed parameter estimation, thereby also obtaining the relative image properties (see Fig. 2).
Combining the four runs, we obtain the coherence ratio C . In our example lensed simulation, we find log = 23.6, correctly consistent with lensing. We then inject two events that are unrelated (see Table 2 for the parameters) and repeat the analysis. In this case, the coherence ratio log C = −14, not consistent with lensing. We can also combine information from the two lensed images to better constrain the binary parameters. In particular, we can use the posterior of the lensed parameters that we obtain from the combined run to reweigh the posterior samples of the first run as in Eq. (15). The most notable impact is on the sky localisation, where the 90% confidence sky area improves by about a factor of two in our example case (Fig. 3). This is particularly important for the strong 5 This is the total time to compute the evidence for the lensed image with a Intel(R) Core(TM) i7-9750H CPU @ 2.60GHz processor. 6 The negative value 21 = −0.5 is equivalent to the transformation 21 = 1.5, and the negative value 21 = −1 is equivalent to 21 = 1. Thus, we do not consider them.

Parameter
Value lensed event Value unlensed event In the lensing framework, the distance of the event is the apparent one, as both images are in fact affected by the lensing parameters. lensing science case, as an improved sky localisation might help narrow down the number of possible host galaxies when combining the GW information with electromagnetic data ).
As an additional example, we analyse a sub-threshold trigger (a signal hidden in the noise background). In a targeted sub-threshold search, one uses a template bank to cover the source parameters posteriors recovered from the primary super-threshold event (Li et al. 2019a;McIsaac et al. 2020) -such searches may uncover many additional candidates, which would need to be analysed. We assume that a super-threshold event has already been observed (with  identical parameters to the event described in the earlier example). The lensing parameters for this event are rel = 25 7 , 2 = 1, and 21 = 16 hr, leading to a network SNR of ∼ 5.5, which is below the value typically required for detection. The resulting posteriors are shown in Fig. 4. We recover injected values, but (as expected) the relative magnification measurement is less accurate than it is for typical super-threshold events. The coherence ratio log C = 9.3, consistent with lensing. . Posterior distribution of the magnification, and time delay (recentered at zero) between two strongly lensed gravitational-wave images, where one of the events is sub-threshold (with an SNR of ∼ 5.5). The simulated relative magnification rel = 25, relative Morse factor 21 = 0.5, and relative time delay 21 = 57600 s = 16 hr. The parameters are well recovered and the Morse factor difference is also uniquely recovered in this scenario.

MULTIPLE IMAGE ANALYSES
In practice, strong lensing will often produce more than two images (Schneider et al. 1992b). Our approach enables one to study multiple images within a Bayesian framework in a computationally tractable way. The methodology detailed below is applicable to any number of images produced by strong lensing (see also (Lo & Magana Hernandez 2021) for a derivation in the context of full joint parameter estimation). Under the hypothesis that individual GW events are lensed, their detector-frame parameters are related to one another, while, under the unlensed hypothesis, they are not. Accounting for multiple images, the coherence ratio with the numerator and the denominator The conditioned evidence for the ℎ image can then be solved similarly to the case with two images (see Eq.(13)): where represents the lensing parameters for the ℎ image, i.e the parameters linking the first observed image to the one under consideration. In this expression, the marginalized likelihood is Indeed, to analyse a set of images, one uses the posterior of the first image as a prior for the analysis of the second image; then, the posterior of the combined first two events for the analysis of the third image; ad infinitum until all images have been analysed.

QUADRUPLE IMAGE ANALYSIS: SKY LOCALISATION
Let us analyse an example quadruplet of lensed images. We assume that the first and second images have the same parameters as in Sec. 6. We inject two more lensed signals, with relative magnifications of 4 and 5, time delays of 16 hours and 21 hours, and Morse factors 3 = 0 and 4 = 1, respectively. We begin by analysing the first two images. We then use the reweighed samples obtained from the joint analysis of the first two images to analyse the third image. As a consequence, we retrieve the lensing parameters ( 3 | 1 , 2 , 3 , H ) and the conditioned evidence ( 3 | 1 , 2 , H ) for the third image. We then reweigh the posterior samples from the second run with the results from third, obtaining ( , 3 | 1 , 2 , 3 ). Using those reweighed samples, we analyse the fourth image similarly, obtaining ( , 4 | 1 , 2 , 3 , 4 ).
A particularly noteworthy improvement is in the sky localization, which we show in Fig. 5. The initial 90% sky area of ∼ 20 deg 2 of the first image is reduced to a final area of ∼ 2 deg 2 when accounting for the four images. Such an improvement is particularly important for studies involving lensed host galaxy localisation, which rely on an accurate sky map estimate . The entire analysis was performed in around 12 CPU hours 8 .

CONCLUSIONS
This work introduces an approximate joint parameter estimation methodology, which allows for fast and precise strong lensing analyses. We have demonstrated its use in the analysis of simulated lensed and unlensed events.
Such a methodology enables us to analyse lensed events at a relatively high accuracy with low computational cost. Furthermore, it enables quick multiple-image analyses. The combination of speed and precision allowed by our method will likely become crucial in the future when we expect the number of detected individual events, each of which could, in principle, be a lensed image, to rise rapidly. Future work may focus on more realistic population priors and selection effects (Lo & Magana Hernandez 2021). We anticipate this methodology to play a complementary role to the existing posterior-overlap (Haris et al. 2018) and joint parameter estimation (Liu et al. 2021;Lo & Magana Hernandez 2021) methodologies, where its role would be to perform strong lensing estimates and multiple-image analyses in an accelerated fashionsituated between the two existing methodologies in terms of speed and precision. A three-tier analysis may be possible in the future, where we first analyse the strongly lensed events with a posterioroverlap method, after which we analyse a reduced set of events with our methodology, and finally, the best candidate(s) could be passed to the joint parameter estimation. where,