Searching for Strong Gravitational Lenses

Abstract

Strong gravitational lenses provide unique laboratories for cosmological and astrophysical investigations, but they must first be discovered – a task that can be met with significant contamination by other astrophysical objects and asterisms. Here we review strong lens searches, covering various sources (quasars, galaxies, supernovae, FRBs, GRBs, and GWs), lenses (early- and late-type galaxies, groups, and clusters), datasets (imaging, spectra, and lightcurves), and wavelengths. We first present the physical characteristics of the lens and source populations, highlighting relevant details for constructing targeted searches. Search techniques are described based on the main lensing feature that is required for the technique to work, namely one of: (i) an associated magnification, (ii) multiple spatially-resolved images, (iii) multiple redshifts, or (iv) a non-zero time delay between images. To use the current lens samples for science, and for the design of future searches, we list several selection biases that exist due to these discovery techniques. We conclude by discussing the future of lens searches in upcoming surveys and the new population of lenses that will be discovered.

1 Introduction

The first measurement of the deflection of light by a gravitational field was the angular displacement of stars close to the Solar limb. The experiment, carried out during the 7-minute-long Solar eclipse visible on 29 May 1919, measured a small displacement of stars in the Taurus constellation, of the order of 1.75 arcseconds (Dyson et al. 1920), supporting the factor of two increase in the deflection angle prediction of general relativity over Newtonian gravity. Should the Sun be counted as the first lens “discovery”? Throughout this chapter we will be concerned with strong gravitational lensing, in which multiple images of the source are observed. Though the 1919 experiment would not count as the first strong lens, it is the first genuine observation of the gravitational lensing phenomenon. The first strong lens discovery would wait 60 years.

The relevant calculations of multiple imaging were first published by Chwolson (1924), but a flurry of work was sparked a decade later by Einstein’s 1936 publication on the subject (Einstein 1936). In fact, Einstein had made similar unpublished calculations as early as 1912, possibly to explain the intensity of a newly discovered Nova, Nova Geminorum 1912, but eventually concluded that the low probability, symmetrical lightcurve, and achromaticity of lensing were an incompatible explanation (Sauer 2008). As Einstein had only considered stars as both the lens and source, he correctly concluded that such image separations were unobservable. In 1937, Zwicky noted that galaxies or clusters are much more likely lenses to provide observable lens systems (Zwicky 1937). After the discovery of quasars in 1963, they were quickly recognised as a possible variable lensed source by Refsdal for measuring the Hubble constant, alongside supernovae (Refsdal 1964). Indeed, the first lens discovered in 1979 was a doubly imaged lensed quasar, QJ0957+561, a.k.a. “The Twin Quasar” (Walsh et al. 1979).

The 1980s saw a few more serendipitous discoveries, including the first quadruply imaged lensed quasars: PG1115+080 (Weymann et al. 1980) and Q2237+030 (a.k.a. the Einstein Cross, Huchra et al. 1985). The same decade saw the discovery of giant luminous arcs in the galaxy clusters Abell 370 and 2244–02 (Lynds and Petrosian 1986; Soucail et al. 1987). The arcs were soon spectroscopically confirmed to be behind the cluster (Soucail et al. 1988; Lynds and Petrosian 1989), confirming Zwicky’s prescient suggestion that clusters were likely gravitational lenses.

Systematic lens searches soon began, as high-resolution imaging with the Hubble Space Telescope (HST) and at radio wavelengths was used to look for multiple imaging of quasars. The first successful search in the optical was the doubly imaged quasar UM 673 (Surdej et al. 1987), quickly followed by the quadruply imaged quasar H1413+117 a.k.a. the “Cloverleaf” (Magain et al. 1988). In the radio, systematic searches found the first full Einstein ring (Hewitt et al. 1988) and multiply-imaged radio sources such as MG 2016+112 (Lawrence et al. 1984; Narasimha et al. 1984). By the early 2000s, the Cosmic Lens All-Sky Survey (CLASS; Myers et al. 2003) found 22 lenses within a sample of 13,500 radio sources. These searches are often called “source-selected”, because the sample pre-selection is made on the properties of the background, lensed source.

One difficulty with source-selected lens candidates resides in their confirmation as real lenses: the glare of the bright quasar images masks the often faint lensing galaxy. The task is particularly challenging with ground-based data, although image deconvolution techniques have been effective for deblending in optical imaging (e.g. Courbin et al. 1998; Burud et al. 1998) and spectroscopy (e.g. Eigenbrod et al. 2007). High-resolution data from HST or from ground-based Adaptive Optics (AO) remains the best way to confirm lenses. Source-selected lens samples are mainly limited to quasars, but there are ∼1000 times more galaxies than quasars and therefore about as many more galaxy-galaxy strong lenses than galaxy-quasar strong lenses. In the former, the lens galaxy light typically dominates, even taking into account the magnification of the source. Searches for galaxy-galaxy lenses have thus focused on a pre-selection of galaxies potentially acting as lenses and are therefore called “lens-selected”.

Cluster lenses are also typically lens-selected, as deeper imaging is obtained for known clusters to search for multiple images. The clusters are first found through their X-ray emission (e.g., Ebeling et al. 2001), or optical imaging surveys (e.g. Gladders and Yee 2005), and subsequent HST follow-up reveals, on average, at least 1 giant arc per cluster (Horesh et al. 2010) for the most massive (X-ray-bright) clusters, let alone many more compact sets of multiple images. While the discovery of multiple images in any cluster is very likely, it is a question of the depth and resolution of available imaging, as indeed deep observations of the most massive isolated galaxies can also often uncover faint background lensed sources (e.g., Collett and Smith 2020).

For both lens- and source- selected searches, wide-field imaging and spectroscopic surveys at various wavelengths are now the main catalyst for developing new searches, as the depth and quality of the data allows not only for many more lenses to be discovered, but also new types of lenses. The Sloan Digital Sky Survey (SDSS) was the first such survey, with targeted campaigns yielding new lensed galaxies from imaging (e.g., Belokurov et al. 2007) and spectra (e.g., Bolton et al. 2004), and lensed quasars from a combination of spectra and imaging (e.g., Pindor et al. 2004; Oguri et al. 2006). As modern surveys have ever increasing depth, area, and quality, both machine learning and citizen science have played a pivotal role in dealing with the enormous number of candidates from the increased depth and area (e.g., Marshall et al. 2016; Jacobs et al. 2019a; Rojas et al. 2021). Entirely unique datasets, such as the high-resolution catalogue of Gaia, have allowed for the efficient identification of bright lensed quasars across the whole sky (e.g., Stern et al. 2021; Lemon et al. 2023).

The aim of this chapter is to present a detailed review of the various search methods presented in the literature for finding strong gravitational lenses of all source and lens types, and all image separations. The above dichotomy between source- and lens-selected samples often represents only the first filtering of a survey’s catalogue; the subsequent selection algorithms are where most of the filtering happens. Data types, wavelengths, and source types are all possible ways to subdivide the many lens searches in the literature. We find it more helpful, however, to avoid repetition by broadly sorting the search methods into those targeting a common lensing feature, that is to say, either a magnification, multiple imaging, multiple redshifts, or an associated time delay.

This chapter is outlined as follows: in Sects. 2 and 3 we describe the properties of the lens and source populations respectively. In Sect. 4, we provide an overview of the literature in terms of selection methods following the categories described above. Possible biases and selection effects from these searches are described in Sect. 5, and in Sect. 6 we discuss how lens finding may evolve in the coming decade, in a context where key facilities such as Euclid, Rubin-LSST, and Roman, will be delivering actual data.

2 Deflector Classes

Gravitational lenses span a wide range of characteristics, with various brightnesses, image configurations, separations, lens types, source types, and environments. Figure 1 shows a gallery of gravitational lenses to demonstrate this variety of possibilities. Before describing the lens search techniques (Sect. 4), it is first prudent to describe the physical characteristics of the possible lens systems themselves. For simplicity, we split the descriptions between the deflector (this Section) and source populations (Sect. 3). We limit our discussion to macro-lenses, in which the multiple images of the source are separated by the order of arcseconds, and can be readily resolved by current observatories. Such separations are only possible for deflectors and sources at cosmological distances, and with masses above $10^{8}$ M_⊙ so the relevant lensing objects are isolated galaxies (early- and late-type), galaxy groups, and galaxy clusters.

Fig. 1

Various strong lens systems with quasar, galaxy, and supernova sources. Images are composed from multi-band HST imaging, a mix of HST and ground-based adaptive optics, or Very Long Baseline Interferometry (VLBI). Sources without optical, infrared, sub-mm, or radio imaging counterparts should not be forgotten, however no definite localised examples have yet been discovered. Credit: NASA/ESA/NRAO/AUI/NSF

For lens-selected searches, an initial catalogue of systems is based on finding examples of the expected lens type. For example, high-redshift galaxies lensed by galaxy clusters might be found through an all-sky selection of galaxy clusters. Below we describe the physical properties of these deflector classes to better understand the lens search techniques and their possible biases.

2.1 Galaxies

Galaxy families are often split by morphology; we can broadly categorise them into ellipticals (early-type), and spirals (late-type) via the Hubble sequence. Ellipticals are featureless quiescent galaxies with old stellar populations and no star formation, representing the final product of hierarchical galaxy formation. Spirals, on the other hand, show prominent star formation and coherent rotation, in contrast to pressure-supported ellipticals. Due to their larger numbers at high masses, ellipticals dominate over spirals as lensing galaxies (approximately 7 to 1), with the latter only contributing at image separations smaller than 2 arcseconds (Turner et al. 1984; Kochanek 1996). For more details on lensing by galaxies, we refer the reader to the Galaxy Lensing Chapter.

2.1.1 Early-Type Lenses

Elliptical galaxies, often found within groups and clusters, have ellipsoidal isophotes, are composed of old stellar populations, and lack significant star formation, dust, or gas. They are accordingly red, with the most massive examples termed Luminous Red Galaxies (LRGs). Due to both selection effects of brighter galaxies, and their increased lensing efficiency, massive LRGs constitute the majority of known high-confidence galaxy-scale lenses, and dominate the system’s flux (and therefore colours) in ground-based imaging surveys, as shown in Fig. 2. Ellipticals’ spectra contain clear absorption lines from old stars (Caii H and K, Mg-b, and Na-D), and are void of emission lines. A typical feature used to measure redshifts of ellipticals is the characteristic 4000Å break immediately redward of the Caii H and K lines. The light profiles of LRGs are often well-described by a single Sersic profile (with Sersic index $n\approx 4$), however there is evidence for a dichotomy of formation pathways for massive ellipticals, with one coreless population with extra central light and lower Sersic index fits (Kormendy et al. 2009).

Fig. 2

Ground-based grz images from the Legacy Survey (top) and HST-ACS WFC F814W (bottom) images of eight galaxy-galaxy strong lenses from the SLACS sample (Bolton et al. 2008). The images are ordered by Einstein radius from left to right. Each cutout is 5^′′ across

As ellipticals dominate the galaxy-scale lensing population, several lensing studies have constrained their mass distributions, revealing standard NFW halos, heavy initial mass functions (e.g. Sonnenfeld et al. 2015; Shajib et al. 2021a), and an alignment between total mass and light (e.g. Keeton et al. 1997; Shajib et al. 2019; Schmidt et al. 2022). Their total mass profiles are approximately isothermal, despite the relative non-isothermal contributions from baryons and dark matter (the so-called ‘bulge-halo conspiracy’, e.g., Treu and Koopmans 2004; Auger et al. 2010; Dutton and Treu 2014), however detailed individual lens studies prefer multiple baryonic components or departure from purely elliptical masses (e.g., Nightingale et al. 2019; Powell et al. 2022). The velocity-dispersion function (VDF) has been found to be consistent with no redshift evolution below $z\sim 0.5$ (Shu et al. 2012), however, the evolution beyond that is unclear, as large samples of spectroscopic measurements are currently limited to below $z\approx 1$. Complete gravitational lens samples, through their image separation statistics, are a sensitive probe of the elliptical VDF and its evolution (Oguri et al. 2012).

2.1.2 Late-Type Lenses

The mass distribution of spiral galaxies can be decomposed into their disk, bulge, and halo. The relative mass contributions of each component is unclear from kinematic data alone due to the so-called ‘disc-halo degeneracy’ (van Albada and Sancisi 1986). Gravitational lensing offers a unique method to break this degeneracy given the additional mass constraints from lensing, and constrain both the initial mass function and halo properties (Barnabè et al. 2012).

The efficiency (or cross-section) for lensing by these systems is strongly dependent on the inclination, as this significantly changes the projected mass density – a key quantity for lensing probabilities and configurations. Keeton and Kochanek (1998) show that the bulge component is necessary for strong lensing to occur in face-on systems, which would otherwise not surpass the surface mass density required for multiple imaging, and quad systems are expected almost exclusively from edge-on systems. A handful of serendipitous discoveries and a concerted search within SDSS spectra – SWELLS (Sloan WFC Edge-on Late-type Lens Survey, Treu et al. 2011, see also Sect. 4.4.1) – now provide ∼ 20 such examples. Figure 3 shows a selection of these lenses. The presence of dust is common in these galaxies, and the extinction of both point sources and extended sources must be accounted for during searches and modelling.

Fig. 3

Examples of lenses with late-type galaxy deflectors; the first four are from the SWELLS sample (Treu et al. 2011), and the last one is a serendipitously discovered lensed quasar (Castander et al. 2006)

2.2 Groups and Clusters

Galaxy groups are loosely defined as comprising several to ∼ 50 galaxies within $\sim 1$ Mpc, and are thought to be embedded within their dark matter halos. Galaxy clusters contain anywhere upwards from ∼50 members, and can span up to 10 Mpc. Equivalently, groups have masses $\sim 10^{13-14}$ M_⊙, and clusters have masses over $\sim 10^{14}$ M_⊙. Such clusters can lens many dozens of high-redshift sources, making them efficient tools for studying high-redshift galaxies individually or in a statistical sense. While these extra sources provide more constraints (e.g., 30 sources with a total of 90 multiple images in Abell 2744, Bergamini et al. 2023), the mass distributions are naturally also more complex than isolated galaxies, and even state-of-the-art models are limited to reproducing image positions of multiply sourced systems to within a few tenths of an arcsecond in the image plane (e.g., Acebron et al. 2022). As a result, the lensing configurations are more often non-standard (e.g., Orban de Xivry and Marshall 2009), often leading to larger source magnifications than isolated galaxies (e.g., Robertson et al. 2020), and the lensed images have wider image separations than produced by galaxy-scale lenses (e.g., More et al. 2012). Complex configurations are particularly common in disturbed clusters with ongoing mergers, reflective of the turbulent mass distributions. At cluster scales, the lack of baryonic cooling beyond a certain mass scale implies fundamentally different mass profiles relative to isolated galaxies. X-ray surface brightness and strong and weak lensing analyses show that the total density profiles of relaxed clusters are well-described by Navarro-Frenk-White (NFW) profiles, namely an inner density slope of $\alpha =1$ ($\rho (r)\propto r^{-\alpha }$; Schmidt and Allen 2007; Newman et al. 2013). This leads to less efficient lensing in clusters, and explains well the paucity of observed cluster-scale lensed quasars (Keeton 1998; Kochanek and White 2001). The associated shallower mass profiles lead to a predicted surplus of cusp geometries, in which only three bright images are apparent, and coupled with triaxiality, this can reduce the lensing efficiency by up to a factor of 4 (Oguri and Keeton 2004).

We briefly discuss the discovery methods for clusters, since this is often a significant step towards finding new giant arcs. The hot intracluster gas emits X-ray emission, which has been used to discover clusters in all-sky X-ray surveys, most notably the MAssive Cluster Survey (MACS, Ebeling et al. 2001). In optical imaging surveys, clusters are found as overdensities in the 4-dimensional space of position, colour, and magnitude, for example in the Red-Sequence Cluster Sequence (RCS, Gladders and Yee 2005) and the SDSS Giant Arcs Survey (SGAS, Sharon et al. 2020). The available datasets for these catalogues are key to the discovery of multiply imaged sources. The deepest datasets – such as those of the Frontier Fields (Lotz et al. 2017) and the Reionization Lensing Cluster Survey (RELICS, Coe et al. 2019) – reveal sources as intrinsically faint as 33 mag.

For a full description of lensing by clusters, we point the reader to the Cluster Lensing Chapter.

3 Lensed Source Classes

We briefly review the various sources of gravitational lenses, highlighting the physical characteristics that are important for their discovery and application to lens searching methods.

3.1 Galaxies

The deflector classes from Sect. 2 can naturally be the sources of lower redshift deflectors. However, ETG sources are intrinsically rare, with only 15 known examples of ‘Early-type Early-type Lenses’ (EELs, systems in which both the lens and source are ETGs, Oldham et al. 2017). While the mass function was more important when discussing the deflector population, the luminosity functions and source shapes are the key metrics for the source populations. We cover two general classes of common source galaxies: unobscured, blue star-forming galaxies and dust-obscured star-forming galaxies (or sub-mm galaxies).

3.1.1 Blue Star-Forming Galaxies

In the hierarchical picture of galaxy formation, galaxies are built up through mergers and accretion of smaller galaxies, accompanied by the onset of star formation. This star formation is seen to peak around ‘Cosmic Noon’ ($z\approx 2$), in which the merging galaxies are predominantly low-mass, late-type galaxies – irregulars and spirals. Such high-surface-brightness, compact galaxies dominate the galaxy population, and particularly the observed-frame optical luminosity function, and thus represent one of the most common lensed sources of gravitational lenses.

The spectrum of a star-forming galaxy composes a blue continuum with absorption lines from gas proximate to the star-forming regions, and any associated emission/nebular lines, such as [Oii] 3727 and [Oiii] 5007. These features are particularly useful for identification in spectra (see Sect. 4.4.1). At optical wavelengths, requiring the detection of multiple common emission lines typically limits the redshifts of such sources to $z\lesssim 1$. Beyond this, single emission lines can only suggest a redshift, however at $z\gtrsim 2$, Ly$\alpha$ enters the optical, and may be distinguished at sufficient signal-to-noise as its profile often has a sharp edge towards the blue, and a tail towards the red (Shu et al. 2016a).

The shapes of these galaxies compose multiple star-forming regions perhaps due to recent mergers, with several compact and diffuse components (Ritondale et al. 2019). For example, Fig. 4 shows the well-known lensed star-forming galaxy, the Cosmic Horseshoe (Belokurov et al. 2007), with the reconstructed source showing several typical clumps.

Fig. 4

Left: HST F475W/F606W/F814W colour image of the Cosmic Horseshoe; right: reconstructed source ($z=2.38$) from the F160W observation, showing four distinct star-forming regions. Reproduced from James et al. (2018)

3.1.2 Dusty Star-Forming Galaxies

Dusty star-forming galaxies (DSFGs), also known as submillimeter galaxies, are high redshift (median $z\sim 4$) galaxies with remarkably bright blackbody emission in the rest-frame far infrared, which is redshifted into the millimeter and submillimeter bands. These galaxies are some of the most active star-forming galaxies in the Universe, with star-formation rates exceeding thousands of solar masses per year (e.g., Casey et al. 2014). Their star-forming regions are obscured by dust, which, after being heated by UV radiation, are responsible for the bright far-infrared emission. In addition to copious amounts of dust, these galaxies contain vast reservoirs of molecular gas, which are commonly observed through their emission lines (e.g., Weiss et al. 2013). In contrast, these galaxies are notably difficult to detect in rest-frame optical or UV-bands, as such emission is absorbed and scattered by dust.

One of the unique properties of DSFGs is the fact that their (unlensed) luminosity function has a sharp cut-off at the bright end (e.g., Negrello et al. 2010). This simply means that (without lensing) there are large numbers of lower-flux galaxies and almost no extremely bright sources. Therefore if an unresolved extremely bright DSFG is observed, it is highly likely to be strongly lensed due to the associated magnification (see Sect. 4.2.1). Since most lensing galaxies are faint in the submillimeter, such lensed DSFGs do not suffer from lens-source deblending problems that are common at optical and infrared wavelengths.

3.2 Quasars

Supermassive black holes are thought to exist in the central potential of all massive galaxies, and can have luminosities of up to 10$^{15}L_{ \odot}$, due to the thermal emission associated with friction in the accretion disk. While in this actively accreting state, they are known as Active Galactic Nuclei (AGN), and the brightest high-redshift (typically $z>1$) examples are quasars, which offer a common source population for gravitational lenses, given their brightness and luminosity function peak at high-redshift.

Quasars emit over the whole electromagnetic spectrum, as shown in Fig. 5 in the case of 3C273, the apparently brightest quasar on the sky. We point the reader towards Beckmann and Shrader (2012) and Padovani et al. (2017) for explanations of the physical processes that cause the multi-wavelength emission of quasars. In brief, the accretion disk contributes a continuum power-law in the rest-frame UV to infrared, and surrounding broad and narrow line regions of ionised gas add broad and narrow lines attributed to specific ions. Strong radio emission is found in about 10% of quasars, and is attributed to the presence of powerful jets that can extend out to many times the size of the host galaxy. A variety of observed quasar ‘types’ are possible depending on viewing angle: when viewed face-on, sources have prominent blue continua and broad emission lines (type-I); edge-on viewing leads to an obscured accretion disk and narrow lines (type-II). This continuum of viewing angles leads to a variety of observed quasar properties, and is further divided by levels of dust in the host galaxy, the presence of radio jets, and absorption due to outflows.

Fig. 5

Multi-wavelength observations of the 3C273 quasar at $z = 0.158$ (data taken from Türler et al. 1999)

A typical scale for an accretion disk is several light days, thus quasars are unresolved at cosmological distances, and appear as point sources often outshining their host galaxies. To select these systems without contamination from stars, several wavelength-specific selection techniques have been used: UV-excess selection (Richards et al. 2009), unique infrared colours (Stern et al. 2012a), and radio and X-ray detections. Quasars are also seen to be significantly variable at all wavelengths on timescales above a week. This can be used as a selection method, with the variability amplitude being found to positively correlate with bluer rest-frame wavelengths and smaller black hole masses (MacLeod et al. 2010). We note that UV-excess selections are not reliable for higher quasar redshifts, as the Lyman alpha dropout appears in the UV filters, leading to colours similar to stars. Equivalently for infrared colour selections, the local minimum at $1~\upmu \text{m}$ is redshifted into redder filters, increasing stellar contamination.

3.3 Transients

Transients cover all variable astrophysical phenomena which eventually fade away on human timescales. While supernovae have been known for millennia, several new classes of transient have only been discovered in recent years, and our current understanding suggests that they should exist at high-redshift offering sources for gravitational lenses. The discovery of the first definite lensed examples is only a matter of time as observatories open up new volumes and resolving capabilities. For complete reviews of lensed transients, we direct the reader to Oguri (2019a) and Liao et al. (2022).

3.3.1 Supernovae

Supernovae are intrinsically rare events, constituting the final moments in the evolution of only the most massive stars or specific binary star systems. The former are known as Type Ia supernovae and the latter are core-collapse supernovae. A wide array of supernova subtypes exist based on the presence, or lack thereof, of various narrow absorption lines and exact lightcurve shape. All supernova lightcurves are smoothly varying with a maximum brightness at 15-25 days after explosion, followed by a slow fade, and typical supernovae rates are roughly once per galaxy per century. Refsdal (1964) first suggested supernovae as sources in gravitational lens systems, highlighting that their time variations allowed precise determination of a time delay, and in turn, a measurement of the Hubble constant.

For detection of such transients, the Zwicky Transient Facility has played a pivotal role by surveying the entire sky every two nights down to a magnitude of $r\sim 20$ thanks to its mosaic camera, capable of imaging 47 square degrees of sky at once. From 2016 to 2023, over 100,000¹ transients consistent with supernovae have been reported by ZTF and various cadenced imaging surveys. However only 10% of these have spectra, often limited to those brighter than ∼18.5 mag.

The standard-candle nature of Type Ia supernovae makes them valuable not only for cosmological probes (see, e.g., Goobar and Leibundgut 2011, for a review), but also as a method for detecting an associated magnification from strong lensing. Only a spectroscopic redshift, classification, and peak magnitude are needed (see Sect. 4.2.1). We refer to the reader to the Lensed Supernova Chapter for a full description.

3.3.2 Fast Radio Bursts

Fast Radio Bursts (FRBs) are a new class of radio transients at cosmological distances whose durations range from nano- to milli-seconds. Around 100 FRBs have been discovered, and their population-inferred redshift peaks at $z = 0.4$ (Shin et al. 2022). The individual redshifts of most bursts are unknown since they are unlocalised, though a handful have been localised to $z>0.5$ galaxies (e.g., Ryder et al. 2022). The physical origin of FRBs is unknown, however they are divided into repeating and non-repeating sources. We point readers to Petroff et al. (2019) and Cordes and Chatterjee (2019) for detailed reviews.

3.3.3 Gamma Ray Bursts

Gamma Ray Bursts (GRBs) are brief highly energetic explosions that appear at gamma ray frequencies and display afterglows at all wavelengths. They are divided into long (70%) and short (30%) GRBs, separated by durations around ∼2 seconds (Kouveliotou et al. 1993). Long GRBs are associated with supernovae, and short GRBs with neutron star mergers.

Since GRBs are particularly bright and occur at cosmological redshifts, they are expected to be gravitationally lensed (Oguri 2019b). Of order $10^{4}$ GRBs have been observed, one third of those with the Fermi Gamma-ray Space Telescope. Were the redshift distribution for GRBs the same as that for quasars, one would expect several of these to have been gravitationally lensed. The search techniques and potential discoveries are presented in Sect. 4.5.2.

3.3.4 Gravitational Waves

Since the first detection of gravitational waves (GW) from the merger of binary black holes in 2015 by the Laser Interferometer Gravitational-wave Observatory (LIGO), nearly 100 GW events from binary systems have been detected. These binaries comprise either a pair of black holes (up to 100$M_{ \odot}$), neutron stars ($\sim 2$M_⊙), or a neutron star and a black hole. The current ground-based GW observatories, LIGO–Virgo–Kagra, are sensitive to GW signals in the frequency range $10-10^{3}$ Hz and the angular accuracy of the sky localisation contours spans an area of $\mathcal{O}(100)$ deg². Supermassive black hole mergers (up to $10^{10} M_{\odot}$) will emit in the nano-Hz frequency range, in which Pulsar Timing Arrays (PTAs) are sensitive (Nanograv Collaboration et al. 2023; Xu et al. 2023). Gravitational lensing of GWs is speculated to be discovered within the next few years with much theory already laid down a few decades ago (e.g. Deguchi and Watson 1986; Nakamura 1998).

Binary black holes (BBHs) are expected to be the first lensed GW detections owing to their cosmological distances. The exquisite precision on the time delay between multiply lensed BBH will make them an attractive probe of cosmological parameters such as the Hubble constant (e.g. Liao et al. 2017). This is feasible via their distance measurements when combined with the complementary electromagnetic (EM) data constraints, for instance, optical redshifts, accurate lensed image positions and lensing configuration. If the binary has a direct EM counterpart as is expected in the case of neutron stars (e.g. GW170817, the first BNS), lensing can shed light on the physical processes right before the merger through an early warning study (Magare et al. 2023).

4 Selection Methods

One single search method for all lenses is impossible, given the wide variety of lens types, source types, image configurations, colours, available surveys, etc. Therefore lens searches are tailored to target specific lenses in specific datasets. Search techniques can be fully described by: (i) the type of lens being searched for, (ii) the data being used, and (iii) the algorithm that prioritises or selects lens candidates. These algorithms often target one of four unique features of gravitational lenses: source magnification, multiple well-separated images, multiple redshifts, or time delays. We divide this section by these lensing attributes. For example, a discovery method that searches for multiple blue arcs around a red galaxy requires only the second feature, multiple imaging, and is independent of the other three. We remind the reader, however, that modern selection techniques can rely on a combination of these methods. We describe serendipitous discoveries in their own subsection, and citizen science searches within the multiple imaging section. Machine learning techniques are given their own section though principally rely on the resolved-images feature of lensing. A summary of the various discovery methods, divided by source type, are listed in Table 1 with approximate numbers of known systems.

Table 1 Summary of discovery methods and numbers of known lens systems. Text in italics implies that the method is a prospective idea for a search

4.1 Serendipitous Discoveries

Nearly all of the first dozen lenses were serendipitous discoveries, and even now lenses continue to be found fortuitously, and merit single publications. We do not aim to describe all these discoveries but highlight some of the most interesting and important examples.

The first known lens, Q0956+561, was part of a follow-up programme for quasars (at that time defined as radio sources) from the Jodrell Bank radio survey, to identify optical counterparts, and determine redshifts (Walsh et al. 1979). They obtained spectra of two blue point sources that might be associated to the radio emission, separated from one another by 6 arcseconds. The spectra were close to identical (at $z=1.405$) and soon after the lensing galaxy was discovered (Young et al. 1980). Within a year of this first lens discovery, a known quasar, PG1115+080, was being followed up for high resolution spectroscopy, when two further stellar objects were observed on the TV guiding monitor within ∼3 arcseconds. Weymann et al. (1980) proposed this as a lensed triple (at $z=1.722$), and later observations showed that the brightest image A was in fact a close pair of images (Hege et al. 1981). A few more years would pass before a third gravitational lens was discovered, also by serendipity. Q2237+0305, also known as the Einstein Cross, was discovered by Huchra et al. (1985), during a spectroscopic survey of bright galaxies. The spectrum showed a quasar at $z=1.7$ at the nucleus of a $z=0.039$ spiral galaxy, and later imaging revealed the four individual quasar images (Yee 1988). Lenses that are outliers in terms of their extreme properties, such as the very low redshift lens of Q2237+0305, are, perhaps unsurprisingly, often found serendipitously. For example, one of the brightest known lensed quasars, APM08279+5255, at $z=3.87$, was a serendipitous discovery within a survey of distant cool carbon stars (Irwin et al. 1998), which share similar infrared colours to high-redshift quasars. Also, the lensed quasar with the lowest source redshift, RXJ1131−1231, was discovered serendipitously while measuring polarisations of unlensed quasars (Sluse et al. 2005). Even in the last few years, bright quadruply lensed quasars are discovered as a byproduct of visual inspection – the two diamonds of Lucey et al. (2018) were identified when removing contaminants for the target catalogue of the Taipan Galaxy Survey.

Many lensed galaxies have also been discovered unexpectedly, including many of the earliest examples. The giant arc in the galaxy cluster Abell 370 was noticed on the CCD images during a spectroscopic campaign of the galaxy members (Soucail et al. 1987). Originally it was thought to be due to tidal interactions or star formation, but soon after, spectroscopy placed the source at $z=0.59$, i.e. lensed by the $z=0.374$ cluster. Galaxy-scale complete Einstein rings are striking and thus relatively easy to spot even within a wide field of data, and are thus regular serendipitous discoveries in the literature (e.g., Tanaka et al. 2016; Bettinelli et al. 2016), and similarly for bright galaxy-scale arcs (e.g., Allam et al. 2007).

4.2 Magnification Selection

A general property of strong gravitational lensing is the magnification theorem: at least one image is as bright as the source (Schneider 1984). Therefore a magnification $>1$ is expected in any strong lens system, implying that the brightest sources will contain an enhanced fraction of lenses. Some systems will show anomalously bright sources that would disobey limits imposed by physics if not lensed. We discuss these two consequences as lens search techniques below.

4.2.1 The Brightest Objects

Simply obtaining high-resolution imaging of the brightest quasars was an early strategic search for lensed quasars above $z\approx 1$ (for example with HST; Bahcall et al. 1992; Maoz et al. 1992), and is still a potentially fruitful method for the brightest high-redshift quasars (McGreer 2013; Fan et al. 2019). The bright quasar catalogues of the Hamburg Quasar and Hamburg/ESO surveys have continued to provide small-separation lensed quasars even long after their initial identification (e.g. Wisotzki et al. 1996; Blackburne et al. 2008). This increased fraction of lensed sources at the brightest end of the observed luminosity function is even stronger for sources with particularly steep luminosity functions – namely very few bright sources and many faint sources. This is apparent at submillimeter wavelengths with dusty star-forming galaxies (DSFGs) (Coppin et al. 2006). Negrello et al. (2010) took advantage of this by selecting bright sources within the 550 deg² of the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS), providing both a large enough area and high enough resolution to select five new lensed dusty star-forming galaxies without contamination from unlensed DSFGs, but still requiring low-redshift interlopers to be removed through colour cuts. Extensions of this technique have been used to build sizeable samples of strongly lensed DSFGs from the Herschel and South Pole Telescope surveys (e.g., Vieira et al. 2013; Wardlow et al. 2013), and to carry out the selection at even brighter magnitude limits across wider areas, using the Planck all-sky survey (e.g., Cañameras et al. 2015).

4.2.2 Outliers from Physical Relations

The associated magnification of strong lensing can be used to identify outliers from non-linear physical relations. The brightest sub-mm galaxies (Sect. 4.2.1) can be thought of as outliers from an underlying true luminosity function, however in this section we highlight objects with other physical properties or relations that are disobeyed when magnification changes the observed brightness.

The X-ray – UV luminosity relation for quasars has a clear flattening at higher luminosities. As outlined in Stern and Walton (2020), if a source is strongly lensed, the UV and X-ray luminosities are increased by the same factor, taking them above the known relation, suggestive of lensing, as they demonstrated for the system MG1131+0456 (see Fig. 6). It has also been used to argue against strong lensing for high-redshift lensed quasar candidates (Connor et al. 2021): the physical size of the proximity zone of a quasar at high-redshift should be related to the intrinsic luminosity of the ionising quasar source. However, strong gravitational lensing will change the inferred luminosity, while the inferred size remains the same, and a discrepancy between the size-luminosity relation could support the lensing hypothesis. Davies et al. (2020) showed that this technique can successfully recover the known $z=6.5$ lensed quasar, UHSJ0439+1634, while it conclusively rules out lensing for the candidate SDSSJ0100+2802 (see Fig. 6).

Fig. 6

Examples of lens discovery as outliers from physical relations, which can otherwise be explained by an associated magnification from lensing. Left: MG1131+0456 UV - X-ray luminosity before and after the magnification against the observed relation, reproduced from Stern and Walton (2020). Right: proximity zone size for the lensed quasar UHSJ0439+1634 appearing as an outlier among mock systems with similar apparent luminosities, reproduced from Davies et al. (2020)

Another physical limit is that of accretion rates onto SMBHs, i.e. an Eddington ratio of unity. This has been used to propose an associated magnification to a quasar, in the case of J2329–0522 at $z=4.85$ (Yue et al. 2023). In this case, however, no multiple imaging has been observed with HST data, and thus has been dubbed an intermediately-lensed system, where the weak lensing magnification is still significant.

The magnification associated with lensing is most apparent when the source has a specific absolute luminosity. This is the case with type Ia supernovae – they are standard candles, all having similar peak absolute brightnesses. This feature is notably used to probe cosmological parameters, as the redshift and luminosity distances can both be measured. Assuming a cosmology, simply measuring the redshift of a type Ia supernova allows the prediction of the observed peak brightness. In the case of gravitational lensing, the additional magnification will lead to an observed peak brightness brighter than predicted. This is how the first galaxy-scale lensed supernova, iPTF16geu, was discovered (Goobar et al. 2017). By obtaining spectroscopic redshifts and classifications of supernovae discovered in the intermediate Palomar Transient Factory (iPTF), iPTF16geu – at z = 0.409 – was found to be 4 magnitudes (a factor of 50) too bright. High-resolution follow-up imaging clearly revealed the 4 multiple images, with a fraction of this extra magnification associated to microlensing (e.g., More et al. 2017). The discovery spectra also showed clear absorption signatures associated to the lensing galaxy, presenting a further sign of strong lensing. A second example, SN Zwicky, has also been found in similar fashion (Goobar et al. 2023), with an Einstein radius below 0.2^′′, demonstrating that this selection technique can uniquely probe the full mass range of lensing galaxies. Even without a spectroscopic redshift, Quimby et al. (2014a) show how irregular SN Ia colours relative to observed magnitude can be used to infer a higher redshift and hence associated magnification.

Requiring a spectroscopic redshift creates a strong dependence on magnitude, as the brightest supernovae candidates are followed up preferentially. For current imaging surveys – such as the Zwicky Transient Facility – very few lensed supernova are expected to be brighter than this magnitude cutoff, due to the cross-section to lensing being small at low redshift. Future surveys, such as the LSST, will discover several hundreds of lensed supernovae. However, many of these cases will not require spectra for confirmation, since the LSST imaging quality will be able to resolve individual images. (Goldstein and Nugent 2017) proposed a method that relies neither on resolving the multiple images, nor a spectroscopic redshift, nor a classification of the supernova: by considering supernovae near ellipticals to be hosted by them, one can calculate the apparent magnitude if the supernova were type Ia using the photometric redshift of the elliptical. Since the vast majority of other supernovae are fainter than type Ia, a surplus in observed flux should then be associated to lensing magnification of a source supernova (not necessarily even type Ia when the magnification is large).

4.3 Multiple-Imaging Selection

In this section, we discuss the techniques that identify lenses as systems containing spatially separated components as potential images of a background source. Typically this requires the analysis of optical and infrared imaging data (i.e., pixels), but in some cases catalogued detections can be used directly. In most cases, the starting data for a search are derived from wide-field ground-based imaging surveys, and higher resolution follow-up data are needed to not only confirm the lens nature, but also for most scientific analyses. Figure 7 shows the difference in imaging quality between discovery data and subsequent higher-resolution imaging for four quadruply imaged lensed quasars.

Fig. 7

Top: Discovery images of four lensed quasars – from left to right: DESJ0029–3814, ATLASJ0259–1608, DESJ0405–3308, and ATLASJ2344–3056. The pixel sizes of 0.26^′′ (DES) and 0.21^′′ (ATLAS) are typical of current ground-based surveys. Middle: Magellan images obtained with pixel sizes of 0.07^′′ (left three) and 0.11^′′ (right) in excellent seeing. Bottom: HST F814W imaging

4.3.1 Citizen Science

Citizen science is one of the novel techniques that has emerged in the last two decades as a systematic means of doing science in Astronomy. Typically, it involves citizens from all over the world collecting and/or analysing data to collectively achieve a scientific goal. Repetitive tasks requiring special skills, such as pattern recognition, and that are difficult to carry out with computer algorithms are deemed suitable for citizen science. One of the first citizen science projects, Galaxy Zoo, was successfully used to study the morphology of galaxies (Lintott et al. 2008). Another major citizen science initiative – The Zooniverse (Simpson et al. 2014) – was setup to support many more citizen science projects, both within and outside of astronomy, including Planet Hunters, Bat Detectives, and Weather Rescue. Searching for gravitational lenses in modern wide-field, deep, multi-band imaging datasets is well suited to citizen science.

Space Warps (Marshall et al. 2016) was the first citizen science project designed specifically to discover strong lenses from astronomical imaging surveys. The classification interface has tools to aid visual analysis, a spotter’s guide which displays typical examples of lenses and non-lenses and the data stream is interspersed with training images unknown to the users. On the backend, a Bayesian pipeline combines classifications from multiple users to assign probabilities for an image to contain a lens. The end product is not only a high-confidence catalogue of lens candidates with calibrated probabilities, but also a vetted sample of contaminants which can be used by the community to improve their lens search algorithms, and for specific studies of spiral galaxies, red star-forming galaxies or ring galaxies.

Space Warps has to-date conducted lens searches in VISTA-CFHT Stripe 82 Survey (Geach et al. 2015), the CFHT Legacy Survey (More et al. 2016a) and the Hyper Suprime-Cam Survey (Sonnenfeld et al. 2020). A unique feature of Space Warps is the supervised learning approach implemented in the analysis of the image classification wherein realistic training images are used. As a result, the final sample of candidate lenses can be well-characterised in terms of completeness and purity which is useful for the statistical analyses of lenses. Recently, a byproduct of the Zooniverse Hubble Asteroid Hunter project has been to identify strong lenses in archival HST observations, identifying ∼200 new high-quality candidates (Garvin et al. 2022). Citizen science will remain an important way to test automated lens discovery algorithms in future surveys, such as LSST and Euclid.

4.3.2 Visual Inspection

Before citizen science projects, the brute force approach of inspecting all available imaging to look for lenses was only possible with small area surveys. In particular HST has been extensively searched thanks to its depth and data quality, in the Hubble Deep Field (Hogg et al. 1996; Barkana et al. 1999), GOODS/ACS fields (Fassnacht et al. 2004), the COSMOS fields (Faure et al. 2008; Jackson 2008), the fields of known lenses (Fassnacht et al. 2006; Newton et al. 2009), and various other fields (e.g., Moustakas et al. 2007; More et al. 2011). High-resolution imaging can also be obtained at radio wavelengths: the Cosmic Lens All-Sky Survey (CLASS) imaged over 13,500 radio sources at a resolution of 0.2^′′ with the Very Large Array (VLA) at 8.4 GHz (Myers et al. 2003; Browne et al. 2003). Together with the ∼2400 images from the Jodrell Bank VLA Astrometic Survey (JVAS, King et al. 1999), the sample was inspected for systems with multiple components, and promising candidates were followed-up for higher-resolution imaging with the Multi-Element Radio Linked Interferometer Network (MERLIN) and Very Long Baseline Array (VLBA), each at 5GHz. These efforts led to the discovery of 22 new lensed quasars. Similar investigations have been carried out with the VLA in the Southern hemisphere, which uncovered 4 further radio-loud lensed quasars (e.g., Winn et al. 2002), and again more recently with the VLBA using a different parent sample (Spingola et al. 2019).

In lower-resolution imaging, the lensing nature is often not obvious, except for the widest separation systems. The 14,000 square degrees of SDSS imaging fuelled a variety of such searches for wider-separation lensed galaxies (typically with images separated by >3^′′ from the lensing galaxy), through visual inspection of red galaxies with nearby blue catalogued detections (e.g., Shin et al. 2008; Kubo et al. 2009; Diehl et al. 2009; Stark et al. 2013). Before visual inspection, images can be pre-processed to remove the bright lensing galaxies that often hide faint background arcs, either by simply subtracting a rotated version of the system (Anguita et al. 2012), or through Principal Component Analysis (Joseph et al. 2014; Paraficz et al. 2016). Similarly, giant arc searches have been performed by visual inspection of pre-identified clusters in both HST (Sand et al. 2005), SDSS (e.g., Estrada et al. 2007; Hennawi et al. 2008), the Dark Energy Survey (DES, Diehl et al. 2017; O’Donnell et al. 2022), and at group-to-cluster scales in the Hyper Suprime Cam Subaru Strategic Program Survey (HSC-SSP, Jaelani et al. 2020).

However, for those lenses with smaller separations, comparable to the Point Spread Function (PSF) width of the survey, contaminant systems resembling lenses quickly become overwhelming and outnumber the real lenses. For lensed quasars, these contaminants are most often quasars projected close to foreground stars, or compact star-forming galaxies. Thus, such searches are designed at specifically removing these contaminants.

In the remainder of this section, we will discuss techniques that apply some algorithms for reducing the number of candidates to a small fraction of the input catalogue. However, these searches still require a final step of visual inspection of cutouts, sometimes many tens of thousands. Biases from individual inspector’s subjectivity can naturally lead to biases towards certain types of lenses (see Sect. 5), and it is common practice to have multiple authors inspect the same stamps, or inject true positives into the sample to assess this bias statistically. Rojas et al. (2023) investigated the effects of visual classification by using a labelled training set of 1500 mock and real images of lenses and non-lenses. 55 classifiers with varying levels of experience in lens discovery graded the systems through visual inspection of their colour images. As expected, arcs with low signal-to-noise or Einstein radii less than or equal to the PSF width are rarely recovered. They found that substantial variations are found between individual classifiers which can be mitigated by combining scores from 6 or more individual classifiers.

4.3.3 Catalogue Searches

Wide-field imaging and spectroscopic surveys are now available across various wavelengths, meaning that all-sky classification and cataloguing of galaxies and quasars is possible down to a few arcseconds resolution. Simple queries and self-cross-matches within these catalogues can be used to look for nearby sources with similar colours as lens candidates. One of the first dedicated catalogue-based lens searches of this type was the SDSS Quasar Lens Search (SQLS), which obtained deeper imaging and/or long-slit spectroscopy of spectroscopically confirmed quasars from SDSS which were either extended or had multiple nearby components in the SDSS imaging catalogue. Colour cuts on companions reduced the many contaminants, and the multi-year campaign discovered 49 new lensed quasars (Oguri et al. 2006; Inada et al. 2012). These searches have been extended by including additional infra-red imaging for candidate selection (Jackson et al. 2012), or simply by application to newer spectroscopic samples (More et al. 2016b).

Lensed quasars can also be discovered as outliers from the unlensed quasar population in their catalogued magnitudes, due to the additional contribution from the lens galaxy, which leads to an apparent near-infrared flux excess. Ofek et al. (2007) applied this colour cut on SDSS spectroscopic quasars with red $g-H$ colours using the Two Micron All Sky Survey (2MASS), successfully identifying a new double quasar.

Even without spectra, quasars can be selected photometrically thanks to their unique infrared colours with the Wide-field Infrared Survey Explorer all-sky survey (WISE, Wright et al. 2010; Stern et al. 2012a), and various radio and X-ray surveys (see, e.g., Flesch 2021). Searches using WISE photometric quasar candidates coupled with resolved or extended detections in ground-based imaging data such as DES and VST-ATLAS have yielded a handful of new lenses (e.g., Agnello et al. 2015; Ostrovski et al. 2017), but typically with significant contamination from quasar-star projections and star-forming galaxies which mimic the infrared colours of quasars (e.g., Williams et al. 2018). WISE alone is limited due to its resolution of ∼5^′′, resulting in the photometry of any quasar-star projection becoming blended and appearing very similar to that of a lensed quasar. Significant visual inspection is often required as a final step in these processes (Rusu et al. 2019; Dawes et al. 2023; He et al. 2023).

The high level of contamination in photometric searches was partially addressed by Gaia, which provided an all-sky, space-based-resolution map of bright point sources (Gaia Collaboration et al. 2016). For known lens systems, the multiple lensed images were often all catalogued in Gaia, and contaminant systems could be removed with various catalogue parameters (e.g., Lemon et al. 2018). Cross-matches between the aforementioned photometric quasar catalogues and Gaia yielded several dozen new lenses (Agnello et al. 2018; Spiniello et al. 2018; Desira et al. 2022), aided by the proper motion and parallax parameters and improved completeness in the second Gaia data release. For those systems with only one image bright enough to be detected or simply with undetected images, the single Gaia detection was useful when combined with ground-based imaging. By looking for astrometric offsets between the resolved Gaia detection of one image and the blended weighted centroid of the ground-based image, coupled with a difference in the overall photometry, Lemon et al. (2017) discovered several candidate lensed quasars which were otherwise only seen as single point sources in ground-based imaging. This is similar to the method employed by Jackson and Browne (2007), who looked for offsets between radio positions from FIRST and optical SDSS galaxy positions. Methods making use of Gaia data alone and its improved completeness and various catalogued parameters has unveiled a potentially large population of new lenses and dual quasar candidates, thanks to high-resolution HST follow-up imaging (Shen et al. 2021; Chen et al. 2022b; Mannucci et al. 2022). However, contamination from quasar-star projections still dominates, as detections with close neighbours either do not have recorded proper motions or parallax, or these are consistent with zero for most faint stars. Simultaneous use of Gaia and pixel-fitting of other imaging surveys has overcome this remaining contamination (Sect. 4.3.4). Quadruply imaged point-like sources occupy specific regions in the position/flux-ratio space, such that any experienced observer can easily recognize lens systems from spurious star/galaxy asterisms by only looking at the image positions and their relative magnifications. Indeed, depending on the position of the source compared to the tangential caustics, one usually ends up with three type of configurations: (i) cusps, when the sources stands near the intersection of two caustics, (ii) folds, when close to a single caustic and (iii) crosses, when the source lies close to the centre of the caustics.

The position of four postulated images can be forward modelled by a simple lensing potential (e.g., Falor and Schechter 2022), such as the Singular Isothermal Ellipsoid Potential (SIEP), which often recovers image positions for most known lensed quasars, and requires seven free parameters for the fit (2 for the source position, 2 for the galaxy position, and 3 for galaxy parameters, i.e., Einstein radius and shape). Witt (1996) demonstrated that the galaxy, source, and images of a SIEP must all lie on a hyperbola with asymptotes along the major and minor axes of the potential. This reduces the dimensionality of the model to three degrees of freedom, while the constraints are simply the positions of the four images. Wynne and Schechter (2018) also shows that the images of an SIEP must lie on an ellipse centred on the source and axes parallel to the hyperbola’s asymptotes. They subsequently suggest a two-dimensional minimisation (source position along the hyperbola and potential axis ratio) to find the best fit model, from which one can report a figure of merit which determines the disagreement from a SIEP model. They report that for a certain threshold of this figure of merit, 98% of random quartets can be removed while only removing 20% of known quad quasars – namely those with peculiar potentials such as double lensing galaxies. Schechter and Wynne (2019) present an even faster determination of the best-fit model by insisting that the two brightest images fall on a shifted copy of the Witt hyperbola, ensuring always a four-image solution and no need for solving a minimization problem.

The application to lens finding simply requires four positions, which promises an extremely quick lens discovery technique in the catalogues of high-resolution imaging surveys. Indeed, this method recovers the majority of known lensed quasars from Gaia, suggesting that there are no bright missing quads. Machine learning has also been applied to this problem, and is discussed in Sect. 4.6.2.

We conclude our review of catalogue-based methods with cluster lensed quasars. Typically, giant arcs in clusters can be difficult to detect due to their low signal-to-noise and require pixel-level search algorithms (see Sect. 4.3.4). Quasar sources, however, are often bright enough to be available in ground-based imaging catalogues. Shu et al. (2019) cross-matched their Gaia AGN catalogue to known bright cluster galaxies, and discovered a new 21^′′separation doubly lensed quasar. Similarly Yantovski-Barth et al. (2023) created a catalogue of clusters from the Legacy Surveys imaging, and found 5 new high-grade lens candidates (of which 2 were known) by searching for multiple DESI quasar targets within the Einstein radius of their clusters.

4.3.4 Pixel-Level Searches

For lensed quasar discovery, the astrometric lensing information is entirely in two (doubles) or four (quads) relative image positions (often without a lensing galaxy visible). Doubles are often indistinguishable from common contaminant quasar-star or star-star pairs, and thus careful pixel-level analysis that looks for colour similarity through the use of multiple survey datasets remains the leading discovery methodology. For lensed galaxy discovery however, the state-of-the-art is machine learning. We defer the discussion of these searches to their own section (see Sect. 4.6), but here describe the earlier techniques that were applied to look for arcs.

A major limitation of catalogue searches in ground-based imaging datasets is related to deblending – for any lens system with images separated by less than the size of the PSF, the separate images and lensing galaxy risk being catalogued as a single object. Other methods, such as identifying multiple redshifts in spectra (e.g., Bolton et al. 2004, see Sect. 4.4.1), offer ways to overcome this limitation, however it is still possible to detect the multiple images by analysing the pixels themselves. Schechter et al. (2017) fit multiple components to the multi-band pixels of VST-ATLAS for photometric quasar candidates, prioritising systems containing components with similar quasar-like colours, discovering four quasar pairs and three lensed quasars. They remove star-forming galaxies, a common source of contamination for WISE-selected photometric quasars, as systems with components incompatible with a given image’s PSF. This technique has also been extended to DES pixels (Anguita et al. 2018). Chan et al. (2022) have applied the same principle to the 0.6^′′ median-seeing survey of the Ultraviolet Near Infrared Optical Northern Survey (UNIONS), identifying nearby point sources with similar $u-r$ colours via a convolution of the Laplacian of the point spread function.

Given the precise positions of the Gaia detections, Lemon et al. (2019) analysed the pixels of WISE, to extract the photometry at the known positions of the individual point sources from Gaia. Simple colour cuts on the two components were able to remove 80% of previously identified quasar/star projections and retain 99% of lensed quasars. Applying this selection, alongside inspection of Pan-STARRS and DECaLS imaging, long-slit spectroscopic follow-up for 208 candidates was obtained with 2m- to 4m-class telescopes (ISIS on the WHT, ALFOSC on the NOT, and EFOSC2 on the NTT), with a success rate (quasar pairs or lenses) of 88% (Lemon et al. 2019, 2023). Gaia’s limiting magnitude is ∼21, while ground-based imaging surveys are now reaching beyond this, and should resolve many fainter lensed quasars.

Pixel-level analysis for compatibility with lens configurations has also been used for final candidate selection. Typical features in lensed galaxies at both galaxy- and cluster-scales are arcs, which can be targeted through algorithms built to identify such elongated and curved features. Several ‘arc-enchancing’ algorithms have been published, which simply require detection of significant structures in the resulting images: Horesh et al. (2005) targeted giant arcs in clusters by applying Sextractor multiple times, repeatedly removing non-elliptical sources from the segmentation map; Lenzen et al. (2004) applied anisotropic diffusion to enhance arcs along their direction while keeping them the same thickness and reducing noise; and Alard (2006) defined a local elongation map based on second order moments within a few times the effective PSF size. The latter was successfully applied to the Canada-France-Hawaii Telescope Legacy Survey (CFHTLS) by More et al. (2012). Seidel and Bartelmann (2007) similarly use moments to determine local elongation but in cells of a grid, which are subsequently grouped with nearby cells sharing similar local elongation directions, resulting in groups of arc-like structures which can have curvature, and do not require a subsequent detection algorithm. This arcfinder has since been extended to require colour information for the arcs (Maturi et al. 2014). Gavazzi et al. (2014) presented and applied their RingFinder algorithm to CFHTLS; by analysing PSF-matched $g-\alpha i$ images of galaxies (i.e., scaled $i$-band images subtracted from $g$-band images to remove the galaxy light), RingFinder finds connected tangential pixels between 0.5^′′ and 2.7^′′ from the galaxy, requiring at least two detections or one elongated detection. Lee (2017) use a circular Hough transform to identify overdensities of light lying on rings in SDSS LRG cutouts, recovering the four known Einstein rings in the sample, with 40% purity.

Requiring a more strict compatibility with a source lensed by a simple mass model has been applied at the full pixel-level light distribution (Marshall et al. 2009) and on extracted image positions (Chitah, Chan et al. 2015; Sonnenfeld et al. 2018; Chan et al. 2020), which results in high purities, but at the risk of missing systems with more complex lensing potentials. YattaLens combines arc detection and lens modelling for lens discovery (Sonnenfeld et al. 2018). Arc candidates are first found as tangential detections in galaxy-subtracted $g$-band images of the HSC-SSP survey; subsequently the pixels are fit with a lens model (Sersic for lensing galaxy + lensed Sersic for arcs) and with a non-lens model consisting of only Sersic profiles. The systems with a better fit from a lens model were kept as candidates. The algorithm has been further applied to find new lenses in more recent data releases of HSP-SSP (Sonnenfeld et al. 2020; Wong et al. 2022).

4.3.5 Variability

In recent years, cadenced wide-field optical imaging surveys have opened up the possibility of finding variables and transients by comparing different epochs of imaging, via a method called difference imaging. If the source quasars are sufficiently variable, lensed quasars appear as unique sources in such difference images, since they consist of multiple crowded variable point sources, unlike most variable objects, as first proposed by Kochanek et al. (2006). They predicted that the dominant source of extended variable objects as seen by examining combined difference images would be lensed quasars, with some small contamination from binary quasars, quasar/variable star projections and variable star pairs. Figure 8 shows such an extended stacked variability map of a known lensed quasar. Lacki et al. (2009) first demonstrated the efficacy of this technique on 20,536 candidates in the SDSS Supernova Survey region, finding only eight extended variables recovering the only known lensed quasar in the sample. Kostrzewa-Rutkowska et al. (2018) searched 670 square degrees of OGLE (Optical Gravitational Lensing Experiment; an optical imaging survey of the Magellanic Clouds region with cadence of ∼3 days) for red WISE sources, i.e. photometric quasar candidates, with at least two variable components within 6 arcseconds yielding 63 candidates, of which 3 showed similar lightcurves. Two were spectroscopically confirmed as quasar-star projections, while OGLEJ0218-7335 was confirmed as a new lensed quasar. Within this area and to the single epoch depth of OGLE, approximately 10 lensed quasars are expected. Lemon et al. (2020) used component modelling of DES epoch data to measure variability properties of lensed quasars and quasar-star projections, further demonstrating variability as a powerful selection and prioritisation tool. They also uncovered a variability bias associated with finding lensed quasars through variability alone, further discussed in Sect. 5. More recently attempts have been made to implement such a variability search in HSC by Chao et al. (2021), however the implementation of the difference imaging and selection criteria can limit the purity and efficiency of such a method, and many artefacts can lead to large numbers of false positives. Dux et al. (2023b) use the multiple epochs of Pan-STARRS and apply difference imaging on close Gaia pairs with quasar-like infrared colours in WISE; they present spectroscopic follow-up with only 25% contamination from stellar systems.

Fig. 8

Fig. 3 from Kochanek et al. (2006) showing the field of RXJ1131-1231 (left) and the variability map built from several epochs (right). The lens system is clearly the only extended variable source in the field

We note that such a search could be applied to any lensed transient in which the image separation is greater than the resolving ability of the imaging survey, for example with the discovery of lensed supernovae in the optical, however, this technique would typically be complemented with other methods to ensure earlier detection. Finally, we note that as long as the time delay between images is much shorter than the survey length, each image should have approximately similar variability (if ignoring microlensing). Indeed a non-zero time delay is not necessary for this technique to work, in contrast to the discovery methods outlined in Sect. 4.5.

4.4 Multiple-Redshift Selection

In this section, we discuss the techniques that identify lenses from spectra exhibiting features at multiple redshifts. The Einstein cross was the first spectroscopic lens discovery (Huchra et al. 1985), thanks to the serendipitous detection of broad emission lines in a galaxy spectrum during a redshift survey of nearby galaxies. Despite no clear signs of the multiple images, the presence of a bright $z=1.7$ quasar in the spectrum of a $z=0.0396$ bright spiral galaxy led to the likely conclusion of strong gravitational lensing. In most cases, fibre spectra have led to the identification of lens candidates in exactly the same way: identifying a higher redshift galaxy in the spectra of lower redshift galaxies. This technique can be extended to resolved spectroscopy at a smaller scale, i.e. Integral Field Units (IFUs), while overcoming the drawbacks of fibre spectra.

4.4.1 Fibre Spectra

In the early 2000s, the arrival of the first SDSS data releases provided hundreds of thousands of high-quality spectra for targeted searches of background source features. Bolton et al. (2004) searched 50,000 LRG spectra for higher redshift oxygen and hydrogen nebular emission lines, finding 50 high-quality candidates which became the basis of the Sloan Lens ACS Survey (SLACS) – an HST snapshot imaging survey for strongly lensed galaxies. Figure 9 shows an example of such a lens discovered by SLACS. The efficiency of requiring multiple emission lines, well-fit LRG spectra, and larger Einstein radii led to a final sample of 70 definite strong lens systems (Bolton et al. 2008). This selection has since been repeated many times, with new datasets (Holwerda et al. 2015; Talbot et al. 2021) or machine learning selection (Li et al. 2019), and targeting specific lenses and/or source types: edge-on late-type lenses (Treu et al. 2011), higher-redshift (Brownstein et al. 2012) and lower-mass lenses (Shu et al. 2015), and higher-redshift sources (Shu et al. 2016a; Cao et al. 2020) or early-type sources (Auger et al. 2011; Oldham et al. 2017). Quasars acting as lenses have been found in a similar way (Courbin et al. 2012; Meyer et al. 2019), and examples of lensed quasars with galaxy-quasar blended SDSS spectra have been found (Lemon et al. 2019).

Fig. 9

SDSSJ1205+4910 discovered by Bolton et al. (2008). Left: an SDSS spectrum of a $z=0.215$ LRG shows additional emission lines not fit by the model template, but that of a higher-redshift source at $z=0.481$; middle: SDSS grz imaging with fibre overlaid; right: HST ACS F555W, F814W colour composite shows clear signs of lensing

Due to increased fibre numbers, telescope sizes, and instrumentation improvements, next-generation spectroscopic surveys, such as the Dark Energy Spectroscopic Instrument (DESI), will obtain spectra for tens of millions of galaxies, allowing for the discovery of lensed galaxies and quasars with fainter, higher-redshift lensing galaxies.

4.4.2 IFU Spectra

While fibre spectra can efficiently provide spectra of many galaxies, they provide only integrated light within the typically 3^′′-wide fibres. This leads to two major drawbacks for lens searches: (i) if the lensed images lie outside of the fibre, the object will not be selected, and (ii) high-resolution follow-up imaging is almost always required to rule out a wide separation projection without evidence of strong lensing. On the other hand, imaging-only searches (Sect. 4.3) can immediately confirm the lensing nature, but follow-up spectroscopy is always required to convert lensing values into physical quantities. IFUs overcome these drawbacks, by simultaneously providing both spatial and spectral information. It is unsurprising that IFUs with wide fields of view have thus been used for targeted lens searches.

Targeted IFU observations have been used to locate the rare low-redshift lenses – particularly powerful probes of measuring the stellar mass-to-light ratio and hence the Initial Mass Function (IMF) due to probing the lens at smaller radii (where dark matter can be neglected, and thus lensing probes only the stars). The SINFONI Nearby Elliptical Lens Locator Survey (SNELLS, Smith et al. 2015) surveyed 27 galaxies selected to have the largest velocity dispersions (and thus largest lensing probability), with observations in the near-infrared to optimise the chances of observing lensed H$\alpha$, discovering two new lenses. This technique has been repeated with SDSS-MANGA (Smith 2017; Talbot et al. 2018), and MUSE with new targets and archival observations (Collier et al. 2018, 2020). Even those systems with only singly imaged nearby background objects provide competitive constraints on the IMF of nearby galaxies. Future surveys, such as Hector, will provide targeted narrow-field IFU observations of over 50,000 galaxies, and will thus offer a chance of creating larger low-redshift lens samples. The ever-increasing area of archival wide-field IFU observations also promises further lens discoveries of various source and lens types (e.g., Galbany et al. 2018).

4.5 Temporal Selection

A unique feature of gravitational lenses is their associated time delay. When the source is variable, this time delay can be measured. However, it can also be used to search for lenses with variable sources. While there is overlap with the idea of using variability of nearby objects (as outlined in Sect. 4.3.5), we focus on lens searches that require a non-zero time delay. For the following techniques, the length of the dataset must be at least the length of the time delay, and the exact success rates will depend strongly on the sampling, signal-to-noise, and intrinsic source lightcurve shape. Typically, variations with large absolute second time derivatives are most useful for detecting a time delay.

4.5.1 Lightcurves

Pindor (2005) showed that resolved lightcurves can significantly improve a search for lensed quasars by targeting systems with well-detected time delays through the use of a dispersion statistic between pairs of lightcurves. Even earlier, Geiger and Schneider (1996) proposed using autocorrelations of blended lightcurves to measure time delays in radio monitoring of lensed quasars, bypassing the need for expensive interferometric monitoring to resolve the multiple images. While their aim was not discovery, the principle has been recently revisited as a way to find lensed quasars in blended optical lightcurves, fueled by the promise of many hundreds of thousands of well-sampled quasar lightcurves in upcoming large imaging surveys. After subtracting a best-fit damped random walk (DRW) model, Shu et al. (2021) found robust features in the autocorrelation function at the expected true time delay, with true positive rates of up to 50% while keeping the false positive rate below 10%. Application on real COSMOGRAIL lightcurves returned ∼ 20-25% of the systems. Bag et al. (2022) propose a method of reconstructing lightcurves and identifying lenses as those with minimum fluctuations at non-zero time delays. Denissenya et al. (2022) have applied a similar technique to blended lensed supernovae using freeform templates. By considering the best-fit models for one, two, or four images, and using the Akaike Information Criterion, they find they can always recover the number of images for systems with modest magnification ratios and time delays longer than 10 days. Krone-Martins et al. (2019) note that the blended lightcurve is less stochastic than that of an unlensed quasar, and thus has lower entropy. Applying this low-entropy selection to quasar lightcurves in the Catalina Real-Time Transient Survey with multiple Gaia detections led to the successful discovery of new doubly imaged lensed quasars.

4.5.2 Transients

Observations of GWs, GRBs, and FRBs often provide uniquely well-sampled time baselines that would allow for an exact time delay measurement if a source were strongly lensed into multiple images. The principle for detection of a lensed event is straightforward: for each image, all the parameters should be identical (except for the arrival time, luminosity distance, and coalescence phase). Therefore, a search for signals with similar parameters may highlight pairs of lensed images from the same source. For GWs, Abbott et al. (2021) explore a ranking statistic for such a posterior overlap, and a joint posterior estimation between pairs of events. In most of the (circular, non-spinning) binary systems, we expect lensing to only change the amplitude of the signal and introduce a constant phase shift between a pair of lensed events. Thus, some methods simply search for similarity in the shape of the GW signals in the time and/or frequency dimensions as a means of identifying lensed counterparts. Accounting for overlap of sky localisation regions, population priors, selection effects, and the odds against lensing, no convincing lens pairs have yet been found. We point the reader to the Gravitational Wave Chapter for a full description of lensed GW detection methods.

For GRBs, 6 pairs have been identified as possible image pairs of a single event. A summary of relevant properties is given in Table 2. A common feature of these candidate events is their short time delays. This has led the studies’ authors to posit that their lenses are black holes in the range $10^{4}$ $M_{\odot} < M < 10^{6}$ $M_{\odot}$.

Table 2 We list the main properties of paired GRBs: times between the pairs of events in seconds, the observed flux ratios (namely the measured luminosity of the first pulse divided by the second pulse), total magnification (summing both images) for a point mass lensing model, and the inferred lensing mass (including a factor of $1+z$, where $z$ is the unknown lens redshift)

FRB searches are again similar, but are constrained to detect time delays below the duration of the data dumped by detection telescopes (see Sect. 5), corresponding to lens masses of $\sim 10^{4} M_{\odot}$. Incoherent FRB lensing searches use detailed studies of individual lightcurve morphologies to distinguish between a doubly-imaged FRB and the more-common scenario of repeat emission from the same source (Liao et al. 2020; Krochek and Kovetz 2022; Zhou et al. 2022) – approximately $5-10\%$ of FRB sources emit repeat bursts. These searches probe minimum time delays comparable to the pulse width. Another class of searches overcomes the degeneracy between the pulse morphology and gravitational lensing by relying on the compact emission region of the FRB and detecting the interference pattern between the two images in the time domain (Kader et al. 2022; Leung et al. 2022). These coherent searches are limited by the fundamental time resolution of the telescope, though radio scattering by inhomogeneous gas in foreground haloes remains an obstacle to detection. These complementary search methods probe different timescales and therefore different lens mass ranges, and upper limits on the fraction of dark matter comprised of compact lenses have been published, though detections of lensing have not yet been claimed.

4.6 Machine Learning

Machine learning is the field of artificial intelligence devoted to the development and study of methods that are able to ‘learn’. From a user-provided sample of training examples, these algorithms either:

• build arbitrarily complex relations between input attributes and numerical or discrete output values that maximise a given score measure in supervised learning methods; these models are then used on observations to perform predictions

• group observations according to the similarity of their input attributes according to a given similarity measurement in unsupervised learning methods.

Machine learning dates back to the 1950s. Examples include teaching machines to play the game of Checkers (Samuel 1959) and Rosenblatt’s perceptron – a linear neural network capable of recognising letters. In the search for strong gravitational lenses, supervised machine learning methods are commonly used as two-class classification problems (lens vs non-lens) with output values representing the probability of the input being a lens. These classification techniques can be viewed as automatically drawing frontiers in a multidimensional input space of attributes in order to isolate objects of two intrinsically different types. Machine learning’s first successful applications for lens finding were as late as 2017 (Petrillo et al. 2017; Jacobs et al. 2017; Lanusse et al. 2018), partly due to the advent of wide-field, deep surveys, but also to the requirement of large computational resources associated with the processing of astronomical images. A comparison of lens-finding algorithms, including conventional methods (e.g., arc finders and visual inspection), and machine learning methods, such as support vector machines (SVM, e.g., Hartley et al. 2017) and convolutional neural networks (CNNs), was presented in the ‘strong gravitational lens-finding challenge’ (Metcalf et al. 2019). Here, we will focus on the main machine learning techniques that have been applied to lens searches: CNNs and decision trees, however the field is quickly evolving, with new techniques constantly being developed and applied to new datasets. We direct the reader to Bishop (2006) for a more in-depth introduction to machine learning, and to Huertas-Company and Lanusse (2023) for an application of machine learning to astronomical survey data (which includes an overview of its application to lens finding and lens modelling).

4.6.1 CNNs

Since the early efforts by Dieleman et al. (2015) and Huertas-Company et al. (2015), deep learning algorithms have become state-of-the-art methods for galaxy morphological classifications based on pixels. These techniques similarly offer a solution to extract the pixel-level information contained in strong gravitational lenses, especially for systems with extended, spatially-resolved lensed features. The machine-learning algorithms best optimised for image analysis are supervised CNNs (LeCun et al. 1998). CNNs are able to capture the morphological patterns in imaging data by learning the coefficients of convolutional kernels and creating a range of two-dimensional feature maps. In practice, between the input layer of astronomical images, and the output layer of corresponding labels, CNN architectures comprise a range of convolutional, pooling, and fully connected hidden layers. To differentiate the network from a simple linear regression, non-linear activation functions, such as the hyperbolic tangent function, are inserted between each one of these layers.

Like all supervised learning methods, CNNs rely on a ground-truth dataset comprising images and their associated labels. During training, CNNs adjust their convolutional kernel weights by minimising a loss function, which encodes the difference between the ground truth and predicted labels. After each pass through the entire training set (refereed to as an epoch), the model is evaluated on a validation set to evaluate its performance on independent data. The evolution of training and validation losses are monitored, and the best CNN model with good generalisation to the validation set corresponds to the epoch with lowest validation loss. Information learned by the CNN is eventually stored in the feature maps. The size of labeled datasets required for training and validating CNNs for image classification vary between 10⁴ to 10⁶ images depending on the number of classes, image size and complexity, as well as network depth. Since the number of known strong lenses is smaller, it is necessary to create mock training sets.

The success of pioneering searches (e.g., Jacobs et al. 2017; Petrillo et al. 2017) and the results from the strong-lens classification challenge (Metcalf et al. 2019) have demonstrated the efficiency of supervised CNNs trained on strong-lens simulations to identify galaxy-galaxy configurations. While improving the CNN architectures only leads to minor gains (e.g., Schaefer et al. 2018), generating highly-realistic lens simulations and training sets that account for the complexity of the actual survey datasets is the main ingredient to reach optimal classification performances (e.g., Lanusse et al. 2018). Building on these early studies, CNNs have been applied to a large range of single and multiband imaging surveys in the optical and near-infrared to select strong-lens candidates, including DES (Jacobs et al. 2019a,b; Rojas et al. 2021), KiDS (Petrillo et al. 2019; Li et al. 2020, 2021), PanSTARRS (Cañameras et al. 2020), HSC (see Fig. 10; Cañameras et al. 2021; Shu et al. 2022), the DESI Legacy Imaging Surveys (Huang et al. 2020, 2021; Storfer et al. 2022), UNIONS-CFIS (Savary et al. 2022), and HST fields (Pourrahmani et al. 2018). In all cases, a visual inspection stage appears to be necessary for validating the CNN candidates and for increasing the sample purity by removing non-lens contaminants such as groups, spirals, and especially ring galaxies (see Rojas et al. 2021).

Fig. 10

Colour composite images ($10^{\prime \prime} \times 10^{ \prime \prime}$) of 105 high-quality strong lens candidates identified from the Hyper Suprime-Cam Subaru Strategic Program by Shu et al. (2022) using CNNs

Each of these searches has found tens to hundreds of high-quality candidates with apparent strong lensing features and, while the grading schemes vary between different teams²², $\sim 10^{4}$ new probable or definite galaxy-galaxy strong-lens candidates have been discovered with CNNs. Because confirming these candidates requires spectroscopic follow-up, ideally in combination with high-resolution imaging, the process remains highly time-consuming. Among the recent spectroscopic confirmation campaigns, Tran et al. (2022) have conducted a systematic follow-up of high-quality candidates from DES and DeCALS. The 88% success rate they obtain provides evidence for the high purity of CNN-selected galaxy-scale strong-lens candidate samples. The use of CNNs for finding group- and cluster-scale lenses has been much more limited, mainly due to the higher difficulty in simulating the diversity of this population.

CNN-based strong lens classification studies follow the supervised approach, with binary classification (lens or not lens). Multi-class models have been shown not to offer a significant benefit in lens classification (Teimoorinia et al. 2020). To the zero-th order, a supervised classifier can be regarded as a similarity checker. As a result, if the training set is biased in any specific way, decisions from such a classifier will be biased in the same way. For example, Shu et al. (2022) trained two strong-lens classifiers constructed from the same architecture with separate training sets that have different distributions in the lens galaxy redshift. Applied to the same set of strong-lens systems, it was found that, for both classifiers, the performance in lens galaxy redshift ranges that were under-represented in the training set was significantly worse than in lens galaxy redshift ranges that were well sampled. As illustrated with this example, the complex selection functions from supervised learning algorithms can fortunately be quantified with extensive tests based on real survey data (see also Cañameras et al. 2021). On the one hand, biases resulting from the design of training sets complicate the development of a self-consistent CNN classifier able to recover the variety of strong-lens morphologies, including the rare configurations (Wilde et al. 2022), while minimising confusion with non-lens contaminants. On the other hand, this sample bias can actually be exploited to design classifiers that target specific types of strong-lens systems. Jacobs (2022) quantified the sensitivity of lens-finding CNNs to various input properties in the training sample, finding the most important variables to include realistically are colour and PSF. Herle et al. (2023) report that CNNs preferentially select systems with both larger Einstein radii and larger, more concentrated sources.

Overall, these limitations are encouraging the development of complementary data-driven approaches. Among these efforts, Cheng et al. (2020) have conducted unsupervised deep learning classification of Euclid-like images, by combining a convolutional autoencoder with a Gaussian mixture model to cluster the extracted morphological features without the need for a training set, and Stein et al. (2021) have explored the capabilities of self-supervised representation learning using ground-based images. Keerthi Vasan et al. (2023) have demonstrated the use of semi-supervised models in which the training sets are augmented with a Generative Adversarial Network (GAN), discovering ∼20 promising lens candidates in the deep, small-area survey, the Deep Lens Survey. Thuruthipilly et al. (2022) demonstrated that self-attention-based encoders – namely input-dependent weighting of features – can outperform conventional CNNs in finding lenses in mock space- and ground-based imaging. These novel methods are promising, particularly for identifying exotic lens configurations that are difficult to include in training sets.

Perhaps surprisingly, CNNs have very rarely been used to discover lensed quasars. This is likely due to current samples having faint lensing galaxies, leading to confusion between quasar-star projections and double lensed quasars when using optical imaging data alone. However, quadruply imaged quasars are uniquely identifiable when their image separation is large enough. Indeed Akhazhanov et al. (2022) have shown that CNNs can reach 80% classifcation performance on mixtures of mocks and known systems, and recently such a CNN-based search has discovered a new lens in Pan-STARRS gri data alone (2023 a). Andika et al. (2023) find reduced contamination in real data when combining CNNs and vision transformers, and discover new candidate lensed quasars in HSC, as well as lensed galaxies as a fortuitous by-product.

4.6.2 Decision Trees

Decision trees are one of the oldest structures for building supervised learning models. They can be thought of as a sequence of if-else conditions, that visually take the form of a tree. An input passes through the tree and finishes with a classification. While passing through the tree, various attributes of the input decide which subsequent branches to take. The main advantage of decision trees is that they are easy to interpret, though they often provide poor performance given their sensitivity to the training set that is used. They are said to have a high variance. In order to mitigate against this effect, one can instead build $N$ training sets, compute $N$ trees and average all their predictions. This approach, known as Bootstrap aggregation, is used in practice in the Random Forest method along with a random selection of the input attributes. Another variant includes Extremely Randomised Trees (ERT), in which both the input attributes and observation splits are randomly selected, and many such trees are combined to provide predictions.

As such, decision tree algorithms perform well in classifying extragalactic sources, in particular quasars (e.g. Ball et al. 2006). Working at the catalogue level, they allow the user to explore large datasets with little human intervention and affordable computing time, thus selecting candidates with less stringent pre-selection criteria, maximising the precision (recovery rate), and minimising the stellar contamination. This was used as a way to remove stellar contamination from the lensed quasar search of Khramtsov et al. (2019).

As already mentioned in Sect. 4.3.3, quadruply-imaged quasars occupy specific regions of image-position and flux-ratio space. Large astrometric catalogues of individual quasar positions thus constitute an opportunity to detect lensed quasars through the use of decision trees, i.e. without having to rely on complexities of processing astronomical images. Gaia provides an all-sky high-resolution catalogue of quasars, ideally suited to the search of strongly lensed quasars. Using an ERT on Gaia DR2 data, Delchambre et al. (2019) achieved an identification rate of quadruply-imaged quasars that is better than 90% along with a misclassification rate of stellar asterisms below 1%. Nevertheless, it should be noted that the miss of a single image leads to a misidentification rate between 10 and 20%, depending on the lensed image that is unobserved, due to the dimensionality reduction and associated loss of constraints.

5 Selection Biases

Having described various lens searches, we now turn to the possible biases from these searches. Lenses are often used to probe the source or lens population, and require mass modelling. The priors that apply to these models should not be uniform because of the way that the lens system was selected. For example, magnitude-limited lens samples are biased towards higher source magnifications, and thus the sources do not lie uniformly in the source plane, but rather they preferentially lie closer to the high-magnification caustics. This might not strongly bias results when the data are strongly constraining, however population level studies can fall foul to neglecting selection effects. It is important for anyone using lenses to keep in mind the various biases of existing samples of lenses. Furthermore, understanding the limitations of previous selections might initiate novel searches aimed at removing any harmful biases from these samples.

Mandelbaum et al. (2009) provide an in-depth exploration of the selection effects of lensed point sources using a mock generation pipeline with mass models composed of both stellar and dark matter components. Sonnenfeld et al. (2023) extend this investigation to extended sources paying particular attention to stellar and dark matter masses, lens half-mass radius, and source size and surface brightness. Below, we briefly summarise several biases that exist in strong lens samples and explain their origin.

• Magnification bias: For a magnitude-limited sample, all of the intrinsically faint sources in the sample must have large magnifications. This leads to biases towards sources near caustics (by definition the highest magnification regions), and thus with image configurations corresponding to cusps and folds. For the same reason, lenses with larger caustics are more prevalent, which are those with more elliptical potentials or with multiple lensing galaxies (see also the quad ellipticity-shear bias this causes, described later).

• Redshift bias: Spectroscopic searches are limited to requiring multiple emission lines within the spectral range. The majority of spectroscopic searches are performed in the optical, and thus are limited to only one strong emission line above $z\approx 1$. Furthermore, photometric selection of sources, such as quasars, has a redshift-dependent completeness. Above $z=2.7$, the UV-excess used to identify quasars is no longer present due to the $u$-band dropout of Lyman alpha, so more stringent constraints are used to remove stellar contamination, leading to less complete quasar discovery. A redshift bias can also manifest through a bias towards brighter lensing galaxies: as lensed quasar candidates with brighter lensing galaxies (as seen in existing optical imaging) are preferentially selected for limited spectroscopic confirmation, those systems with higher redshift (and thus fainter) lenses are not confirmed. Lower redshift sources naturally have lower redshift lensing galaxies, and are thus brighter and more likely to be confirmed. This is demonstrated in the known lensed quasar sample in Fig. 11.

  Fig. 11

  

  Left: Distribution of image separations against source redshifts for known lensed quasars (red) satisfying certain ‘discoverable’ criteria, overlaid on an updated mock catalogue of Oguri and Marshall (2010). There is a clear lack of known lenses (regions shaded blue) with smaller separations and higher redshifts. Right: examples of quasar pairs with similar spectra from Lemon et al. (2023) that may represent these missing systems, lacking only a lens galaxy detection

• Einstein radius bias: Many search techniques rely on resolving the multiple images of a background source, however if these images are unresolved at the catalogue and/or pixel level then the lens cannot be selected in the original parent sample. This leads to a lower limit on the lensing galaxy Einstein radii. Conversely, fibre-selection of lens candidates (see Sect. 4.4.1) can miss those systems in which the brightest images lie outside the fibres, i.e. when the Einstein radii are larger than the fibre radius.

• Model bias: The mass models that are used to create mock lenses do not capture the complexity of real lens systems, and thus both humans and machines trained on such mocks are at risk of missing real lenses with non-standard mass distributions.

• Microlensing bias: Lens masses with Einstein radii comparable to the source size, such as individual stars and quasar accretion disks, can cause considerable magnifications and demagnifications of the background source. This effect depends on the convergence, internal shear, and stellar mass density at the image positions, and is more dramatic for saddle points (Schechter and Wambsganss 2002). These discrepant flux ratios can manifest as ‘missing’ images in shallow imaging (see Fig. 12), which can lead to systems being missed by various magnitude-dependent selection techniques.

  Fig. 12

  

  $i$-band images of two quadruply imaged lensed quasars: SDSSJ0924+0219 shows one image being severely demagnified by microlensing over two decades, and PG1115+080 during an epoch without strong microlensing. The top two images were taken with MagIC, and the bottom two with IMACS. We note that PG1115+080 has a larger Einstein radius, and thus a larger cutout

• Variability bias: Lemon et al. (2020) uncovered a variability bias associated with finding lensed quasars through variability alone. This bias is due to a well-known correlation between variability amplitude and quasar luminosity (see, e.g., MacLeod et al. 2012), resulting in the lenses with the faintest sources having the most variable images. Quads have, on average, higher magnifications than doubles, and so for the image brightness, quads will vary more since they probe fainter sources.

• Survey length bias: Oguri et al. (2003) pointed out that cadenced surveys for detecting lensed supernovae will miss images and lens systems as the duration of the survey approaches the time delay, leading to biases towards both more symmetric and smaller Einstein radius systems.

• Time-delay selection function Similar to the Survey length bias above, constraints on observations – both spatially or temporally – prohibit the discovery of lensed transients found in the time domain through the reappearance of extra images. For example, FRBs are of millisecond duration, and many widefield instruments (e.g., CHIME/FRB) are incapable of slewing. This imposes a “time-delay selection function” (Connor and Ravi 2022), which selects for strongly-lensed FRB systems whose time delays are either (i) less than the amount of time spent by one sky position within the beam pattern, or (ii) an integer multiple of the sidereal day.

• Caustic area bias: Baldwin and Schechter (2021) explore the Malmquist-like bias of quadruple image lensed quasars having larger caustic areas. Models that do not account for this selection effect typically have caustic areas that are biased low, and as such the predicted time delays are also biased low, leading to an important bias in the Hubble constant from time-delay cosmography of the order of ∼1%.

• Quad ellipticity-shear bias: Quadruple image lenses are biased towards more flattened lensing potentials, due to more elliptical lensing galaxy masses and/or larger external shears. This effect is less strong for extended sources than point sources, as sources can often overlap a large fraction of the background caustics. Conversely, very circular lenses have smaller caustic areas and are thus less likely to be discovered as bright quads, though this might be compensated for by the higher magnifications associated with near-circular lenses.

6 Prospects for Future Lens Searches

As emphasised throughout this Chapter, strong gravitational lensing is a rare phenomenon that spans a very broad range in lens and source properties. It is thus important that very different techniques are used to discover strong lenses. Our playground for this exercise resides in the many multi-wavelength, time-domain and even multi-messenger datasets and surveys that are ongoing or planned in the very near future. We outline which future surveys and observatories will provide the drive for finding orders of magnitudes more lenses. The section is divided by source type, though there is naturally significant overlap between surveys providing the relevant data for finding each type of lens.

6.1 Lensed Galaxies

The ESA-NASA Euclid survey will image 15,000 square degrees of extragalactic sky in one wide optical filter (550–920 nm) and three near-IR filters: Y, J, and H (Laureijs et al. 2011; Euclid Collaboration et al. 2022). The first strength of Euclid will be its spatial resolution of 0.18^′′, giving access to the most compact systems with Einstein radii down to 0.3-0.5^′′. The sharp Euclid PSF increases the contrast between the smooth lens light and the thin arc-like images of the background source. While the broad band-pass of the Euclid optical filter will optimise the detection of lensed sources with strong continuum emission, it may be less effective at detecting lensed sources with flux dominated by strong emission lines. Other future space-based surveys like the Roman space telescope or the Chinese Space Station Telescope (CSST) will complement the capabilities of Euclid through increased depth and number of filters.

Collett (2015) create mock Euclid lens observations to estimate the number of discoverable systems, assuming at least partially resolved sources, magnifications above 3, and sufficient signal-to-noise. If the lens can be well subtracted, 170,000 lenses should exist in Euclid with a peak Einstein radius around 0.5^′′. For comparison, an existing survey like DES is predicted to contain 2400 lenses with peak Einstein radii around 1.25^′′.

The Rubin-LSST survey will image the whole Southern sky in 6 optical filters every few days. This adds the time and colour dimensions to Euclid and will surpass the depth of Euclid after a few months of observation but, as a ground-based survey, will be limited by atmospheric turbulence. It will play a major role in the colour pre-selection of targets, providing photometric redshifts, and in the lens-source deblending using machine learning or morpho-spectral techniques (e.g. Joseph et al. 2016; Melchior et al. 2018). Collett (2015) also predict the number of lenses to be discovered in LSST data alone: 120,000 in the optimal stack, and if the lensing galaxy must be removed through a blue-red band difference imaging technique, then 62,000. The peak of the Einstein radii is around 1^′′, and source redshifts around $z=2$.

The depth and coverage of both LSST and Euclid will expand galaxy cluster catalogues to higher redshifts, through simple conventional searches of overdensities in catalogued nearby galaxies with similar colours and magnitudes (e.g., Gladders and Yee 2005). The ongoing eROSITA X-ray survey will also play a role in finding ∼ 100,000 X-ray-bright clusters out to $z>1$. The discovery of giant arcs in many of these clusters will be possible with the discovery datasets themselves.

At sub-mm wavelengths, the models that correctly reproduce the dominance of lensed systems in sub-mm galaxies (SMGs) at flux densities above 100 mJy, predict that $>1$% of systems down to 1 mJy should be lensed – 10 lensed SMGs per square degree – in agreement with the number of lenses found in a common JWST-SCUBA-2 imaging field (Pearson et al. 2023). Discovery and confirmation of these systems is not possible at optical wavelengths, but could be possible in the wide-field infrared imaging of Euclid. JWST, ALMA, and the ngVLA (next-generation Very Large Array, to be fully operational by 2034, Selina et al. 2018) will be key to demonstrating the efficiency of these searches, and providing higher quality datasets to study the source and lens populations.

6.2 Lensed Quasars

Lensed quasar searches have relied predominantly on source selection, i.e. starting from a catalogue of known quasars and high-confidence quasar candidates. These have come from spectroscopic datasets, such as the 2dF Quasar Redshift Survey (3 × 10⁴), SDSS/BOSS (5 × 10⁵), as well as photometric selections due to UV-excess or typical red colours in the infrared (e.g. 10⁶ from Secrest et al. 2015). These previous searches will soon be extended to fainter magnitudes, as an order of magnitude more quasars are expected from upcoming spectroscopic surveys in the next 5 years, such as the Dark Energy Spectroscopic Instrument (DESI; 4 × 10⁶) and the 4-metre Multi-Object Spectroscopic Telescope (4MOST) extragalactic surveys (3 × 10⁶). Pushing to even deeper magnitude limits, the Maunakea Spectroscopic Explorer (MSE) is planned on an upgraded version of the Canada France Hawaii Telescope (CFHT) with a 11 m mirror, and the capability of simultaneously obtaining 4000 spectra. Photometric selection of $z<2.7$ quasars will also be extended thanks to the $u$-band imaging of the Rubin-LSST survey in the South.

However, at fainter magnitudes, lensed quasar searches can no longer be source-selected at optical or infrared wavelengths. As demonstrated in Fig. 13, for lensed quasars with total image brightnesses fainter than 21, the galaxy light begins to dominate. In such cases, the selection must rely on galaxies, as spectroscopic and photometric quasar catalogues will miss lensed quasars due to the lens light contamination. Lens-selected searches are already the standard for lensed galaxy searches, and should readily be applicable to quasar sources as well; for example, a SLACS-style selection to search for quasar emission lines in galaxy spectra (Sect. 4.4.1). In the X-ray and radio, there is often no contamination from the lens, offering a promising initial source-selected catalogue containing lensed quasars. The ongoing eROSITA X-ray survey will provide a deep, full-sky map at 2-10 keV at 15^′′ resolution, expected to detect 3 × 10⁶ AGN. Cross-matching optical catalogues to X-ray sources has seen success for lensed quasar discovery (e.g., between HSC and Chandra, Jaelani et al. 2021), and even already between eROSITA and Gaia (Tubín-Arenas et al. 2023).

Fig. 13

Distribution of the lens-to-quasar flux ratio against total quasar magnitude for the mock lensed quasars of Oguri and Marshall (2010). The lensing galaxy population is naturally the same at all source brightnesses, so the galaxy flux (and hence colours) will start to dominate the total system flux at the fainter quasar source magnitudes expected to be discovered in future surveys like LSST and Euclid

Another promising selection uncontaminated by the brightness of the lensing galaxy is through extended variability (see Sect. 4.3.5). Briefly, through multiple epochs of the same field in the same filter, all lensed quasars and quasar pairs may be found through difference imaging. LSST will provide the ideal dataset for this search: an approximate cadence of 3 days over 20,000 square degrees in 6 filters. By stacking epochs over the 10 year nominal lifetime, variability should be detectable down to $r\sim 26.5$, or deeper when using multiple bands for detection. Taak and Treu (2023) predict ∼1000 lensed quasars to have variability detectable within LSST.

We can also expect a push towards the discovery of smaller separation lensed quasars. High-resolution imaging of a handful of radio-loud or optically bright quasars has led to the smallest separation known lenses: only four systems confirmed below 0.4^′′. While Gaia offers a 0.2^′′-resolution map of point sources across the whole sky, cataloguing issues in existing data releases at small separations and in crowded local fields leads to the loss of several detections. If such issues are resolved in future releases, or the 1D single-epoch windows are provided, 100s of predicted, bright small-separation lenses might be able to be confirmed within the survey data themselves. At fainter magnitudes, Euclid’s space-based resolution and wide area (see Sect. 6.1) will yield many new lens systems not only from analysis of known quasars and quasar candidates, but will also detect the lensing galaxies of many existing unclassified quasar pairs, for which ground-based imaging is too shallow to detect any lenses.

At radio wavelengths, the International LOFAR Telescope (ILT) will observe 15,000 square degrees at ∼150 MHz with 0.35^′′ resolution and a 1$\sigma$ sensitivity of 90$\mu$Jy. Rezaei et al. (2022) show that CNNs can recover 95% of lensed radio sources with Einstein radii above 0.5^′′ and flux densities above 2 mJy, while removing nearly all contamination from unlensed double-lobed sources. At similar resolution, a deeper all-sky survey with the Square Kilometre Array (SKA) could uncover 10,000 new radio-loud lenses, with sources mixed between lensed radio jets and lensed starbursts, with source redshifts potentially at $z>10$ (McKean et al. 2015). Though the exact surveys of SKA are not yet confirmed, observations are planned to begin around 2028/2029.

6.3 Lensed Supernovae

The planned Legacy Survey of Space and Time (LSST), to be carried out with the Simonyi Telescope of the Vera C. Rubin Observatory, promises an unprecedented dataset for finding lensed supernovae, thanks to its high-cadence, wide-field, and deep imaging. As discussed in Sect. 4 there are both catalogue-based schemes and pixel-based schemes for discovering lensed supernovae. The simplest of catalogue-based schemes are based on the assumption that two or more images of the lensed supernovae are recorded as transients. Given a conservative estimate of a detection limit of 22.6 mag in the $i$-band, an effective survey length of 2.5 years, and image separations above 0.5^′′, Oguri and Marshall (2010) predict a total of 130 galaxy-scale lensed supernovae from LSST. Unresolved catalogue schemes based on strong magnifications of type Ia supernovae from measuring redshifts either photometrically or spectroscopically (see Sect. 4.2.2) are also possible, and Goldstein et al. (2019) predict 380 lensed supernovae in LSST (10 years) with such a method, broadly agreeing with the 180 lensed type Ia supernova estimate of Sainz de Murieta et al. (2023). Wojtak et al. (2019) consider both discovery methods, finding a much higher rate of 340 lensed supernovae per year, with ∼2/3 of these discoverable through their magnification alone. Combining multiple epoch observations and a variety of selection methods, such as CNN classification (e.g., Kodi Ramanah et al. 2022), will be key to discovery of these systems. The delivered image quality, cadence, observing strategy, and access to follow-up facilities will all be important factors for the true discovery rate (Huber et al. 2019a).

As LSST will also provide deep existing imaging of any candidate, and, in the absence of lens confirmation through spectroscopy, deconvolution schemes (e.g., Akhaury et al. 2022) may be used to detect a lensed host galaxy before possible further images arrive. In a similar fashion one might look for a leading image that was too faint to catalogue – “pre-covery” – as the source transient might fall near a fold caustic or due to extrinsic microlensing, similar to the case of iPTF16geu (see Fig. 14).

Fig. 14

Left: HST observation of the lensed supernova iPTF16geu; middle: observation after the supernova has faded; right: subtraction of the constant components, showing just the supernova images. Reproduced from Dhawan et al. (2020)

The number of known lensed galaxies, detected through large spectroscopic surveys or through arc detection in imaging (see Sect. 6.1), will increase by an order of magnitude before the end of the lifetime of LSST. Cross-matching known transients to lens candidates (Magee et al. 2023; Sheu et al. 2023) and monitoring of known lenses (Craig et al. 2021), will become a fruitful method for lensed supernova discovery, as it already has been for monitoring of cluster lenses (e.g., Kelly et al. 2015; Rodney et al. 2021). JWST will become a key tool for lensed supernova discovery in galaxy clusters given its unprecedented depth and resolution and thus ability to reach higher redshift fainter supernovae. Indeed, Frye et al. (2023) already report three images of a new lensed supernova in galaxy cluster G165. Based on star formation rates of source galaxies, they suggest as many as 1 more lensed supernova to be expected from monitoring of that one cluster alone. Shu et al. (2018) predict over 10 lensed supernovae occur per year in a sample of 128 known lensed galaxies from various SDSS lensed galaxy searches.

Such monitoring of clusters with the HST has also recently revealed individual lensed stars and star clusters (Miralda-Escude 1991; Kelly et al. 2018) thanks to the magnifications (both macro and micro) near critical curves of factors of several thousand. JWST is already delivering new systems (e.g., Chen et al. 2022a), and will be the most powerful observatory for the discovery of such systems for the foreseen future.

6.4 Lensed GWs

As LIGO–Virgo–Kagra (LVK) currently enter their fourth observing run (O4) at design sensitivity, detecting lensed GWs from stellar-mass black hole (BH) mergers becomes a possibility, and very likely in the fifth observing run (O5) sensitivity (estimated to start in 2027). Predictions for lensed GWs suggest rates around 0.1 – 5 lensed BH-BH mergers per year in O4, and 1 – 50 per year in O5, with the range caused by our uncertainty on underlying compact merger rates (e.g., Li et al. 2018; Ng et al. 2018). Lensed NS-NS events are unlikely to be detected, even during O5 (Smith et al. 2023). Most time delays will be shorter than one month, and ∼30% of lensed events are expected to be quadruply imaged (Li et al. 2018).

Several next-generation GW observatories (so-called third-generation) are currently planned and will be sensitive to lensed GWs of various origins. Many of these observatories will also have superior angular resolutions of 1^′′or better, and extend the detection volume out to high-redshift, offering the chances of detecting optical counterparts of lensed NS-NS mergers. The Einstein telescope will have sensitivity to NS-NS events out to $z=2$, and thus even higher redshift for lensed events. Piórkowska et al. (2013) predict up to several tens of detections per year, and an order of magnitude more lensed BH-BH mergers (Biesiada et al. 2014).

The Laser Interferometer Space Antenna (LISA) is a planned (2030s) space-based set of spacecraft that will be separated by 8 light-seconds, and will provide sensitivity to larger volumes, at lower frequencies, and accordingly higher total merger masses. Optimistic models suggest up to 4 well-detected strongly lensed events within the 5-year nominal mission (Sereno et al. 2010). Between the LISA-LIGO frequency gap, DECIGO, a planned space-based mission, will operate at deci-hertz frequencies. The principle mission is to detect primordial GWs expected from inflation, however it will also detect Intermediate Mass Black Hole mergers and their lensed counterparts (Piórkowska-Kurpas et al. 2021).

At even lower frequencies, we expect to see GW signals from supermassive BH mergers. Pulsar Timing Arrays (PTAs) have recently reported strong evidence for a stochastic background consistent with such a population of BH mergers (Nanograv Collaboration et al. 2023). Khusid et al. (2022) predict 10-30 lensed binaries might be detected using PTAs with the Square Kilometre Array. Tantalising evidence at optical wavelengths for such a system already exists in a known lensed quasar (Millon et al. 2022). Figure 15 summarises the frequencies and astrophysical sources of next-generation GW observatories. The sky localisation uncertainties depend on the signal-to-noise of the signal but typical localisation areas for the observatories (in square degrees) are ∼1 (LISA, Kocsis et al. 2008), ∼10 (advanced LIGO), and ∼40 (SKA and iPTA, Sesana and Vecchio 2010; Finn and Lommen 2010).

Fig. 15

Gravitational wave spectrum. It shows the frequencies at which various GW observatories (current and upcoming) are going to operate and correspondingly astrophysical sources. Figure is reproduced from K. Islo’s Thesis. Credit: C. Moore, R. Cole, C. Berry and C. F. Mingarelli

6.5 Lensed FRBs

Perhaps the most promising path to detecting galaxy- or cluster-lensed FRBs involves a combination of microlensing searches and high-precision FRB localisations to their host galaxies. In this scenario, optical follow-up can confirm the lensed host galaxy. This will be enabled by widefield FRB surveys such as CHIME/FRB Outriggers (Leung et al. 2021; Mena-Parra et al. 2022), CHORD, the DSA-2000, or future coherent all-sky monitors such as BURSTT (Lin et al. 2022a), which will overcome the time-delay selection bias by extending the detection field-of-view down to the horizon.