Nine lensed quasars and quasar pairs discovered through spatially-extended variability in Pan-STARRS

We present the proof-of-concept of a method to find strongly lensed quasars using their spatially-extended photometric variability through difference imaging in cadenced imaging survey data. We apply the method to Pan-STARRS, starting with an initial selection of 14 107 Gaia multiplets with quasar-like infrared colours from WISE. We identify 229 candidates showing notable spatially-extended variability during the Pan-STARRS survey period. These include 20 known lenses, alongside an additional 12 promising candidates for which we obtain long-slit spectroscopy follow-up. This process results in the confirmation of four doubly lensed quasars, four unclassified quasar pairs and one projected quasar pair. Only three are pairs of stars or quasar+star projections, the false positive rate is thereby 25%. The lenses have separations between 0.81"and 1.24"and source redshifts between z = 1.47 and z = 2.46. Three of the unclassified quasar pairs are promising dual quasars candidates with separations ranging from 6.6 to 9.3 kpc. We expect that this technique will be a particularly efficient way to select lensed variables in the upcoming Rubin-LSST, which will be crucial given the expected limitations for spectroscopic follow-up


Introduction
Gravitationally lensed quasars span a vast range of astrophysical and cosmological applications.Among them, time-delay cosmography has the potential to address reliably the so-called Hubble tension (e.g.Verde et al. 2019).The method relies on the photometric variability of quasars, which enables the measurement of the time delays between the lensed images, a quantity that is directly -and inversely -proportional to the Hubble constant (Refsdal 1964).Currently, only seven lensed quasars have been used (Wong et al. 2020;Shajib et al. 2020), in part because of the paucity of known lensed quasars with sufficient variability.This scarcity has triggered a variety of searches, in recent wide-field datasets, for lensed quasars and the number of systems has doubled in the last few years alone (e.g.Chan et al. 2022;Lemon et al. 2023).These searches have also generated useful by-products, such as dual quasars, i.e. physically associated quasars that are thought to be the progenitors of supermassive binary black holes responsible for the nano-hertz gravitational wave background (e.g.Mingarelli 2019), and a crucial element for understanding galaxy evolution (e.g.Kormendy et al. 2009).On the other hand, as with any search for rare cosmic objects, contaminants often outnumber the targeted objects by several orders of magnitude.Among the most common contaminants are quasar-star projections and star-forming galaxies (e.g.Treu et al. 2018).Such contaminants typically feature only one or no variable point sources, so looking for variability in more than one point source -hereafter extended variability -can be highly effective in removing contaminants, as originally suggested by Kochanek et al. (2006).Lacki et al. (2009) perform such a search with difference imaging in the SDSS Supernovae Survey region, recovering known lensed quasars with a false positive rate of ∼ 1 in 4000.Kostrzewa-Rutkowska et al. (2018) present the first variabilityselected lensed quasar discovery.This was not done through difference imaging (as in the present work), but through similar variability features in their light curves, using the 3-day-cadence 700 square degrees field of the Optical Gravitational Lensing Experiment (Udalski et al. 1992, OGLE).Their initial selection required quasar-like mid-infrared colours, and two candidates were confirmed as quasar-star projections.Chao et al. (2020) apply a difference-imaging-based search to the Hyper Suprime Cam (HSC) Transient Survey data, designed to effectively identify quadruply-lensed quasars.However, application to pre-selected variables only reduces the candidate list by ∼20 per cent, highlighting the need for additional visual inspection and further analysis (Chao et al. 2021).Lemon et al. (2020)  the less variable point source component of each system is indeed an effective way to recover all known lenses in the considered footprint while removing 94 percent of contaminants.
The Pan-STARRS 3π survey (Chambers et al. 2016;Waters et al. 2020) offers an unprecedented dataset for lensed quasar searches, having imaged 3/4 of the sky on average 7.6 times (Magnier et al. 2013) in each of 5 filters between 2009 and 2012, with further observing epochs between 2012 and 2014.Furthermore, Oguri & Marshall (2010) (hereafter OM10) predicts just under 2000 lensed quasars should be found in the Pan-STARRS dataset, to be compared with only ∼300 known systems as of today. 1n this paper, we present a difference imaging algorithm and its application to lensed quasar candidate identification in Pan-STARRS.We demonstrate the algorithm's efficacy in identifying lensed quasar candidates, discuss the by-products of our search, and explore potential applications for future surveys.We first describe the details of the difference imaging method used herein in Section 2. Next, we describe the available datasets and selection of candidates in Section 3. In Section 4, we present follow-up imaging and spectroscopy and a discussion on individual systems.We draw conclusions in Section 5.
In all this work we assume a flat ΛCMD cosmology with H 0 = 70 km s −1 Mpc −1 and Ω m = 0.3.

Difference imaging pipeline
Detecting photometric variations of astronomical objects implies obtaining images at different times and therefore of different quality.The seeing conditions and sky level in particular are bound to vary from one image to the next.In order to detect an intrinsic difference in the total luminosity of a given object, the simplest strategy is to subtract a reference image from all other images in the set.In an ideal situation, the only signal rising above the noise in such a difference-image would be that of the sought intrinsic variability.This will naturally never be the case because of the varying sky seeing and sky brightness.Hence a proper procedure needs to be devised to account for the different point spread functions (PSF) and overall brightnesses of the images.
In this section, we first introduce the main idea of difference imaging and introduce variability maps, which are the stacks of difference images used to assess extended variability.Lastly, we benchmark the method on real data.We make the resulting pipeline available online as a Python package.2 .

Basics of difference imaging
Given a 2-dimensional image I(x), and a reference image R(x) of the same size, we adopt the following simple model that blurs R(x) so it matches the resolution of I(x): where the reference image R is convolved with the kernel K that minimizes the following χ 2 : with n the total number of pixels in the image.Usually the reference image is the one with the best seeing such that K can be a low-pass filter as detailed in Appendix A.1.The constant term B in Eq. ( 1) accounts for the background level of the sky, and is assumed to be constant across the field of view.Additional terms (e.g.polynomials in the image coordinates) can potentially be added to better fit any significant light gradient, however the considered fields are small enough to make this unnecessary.Next, the kernel is decomposed onto a set of basis functions whose weights can be tuned to minimize Eq. ( 2).Choosing Gaussian basis functions yields the approach taken by Alard & Lupton (1998) for the analysis of Galactic microlensing light curves.We adopt instead the more robust (albeit slower) approach of Bramich (2008) and we set the kernel to be a sum of delta distributions, namely a pixelated kernel.With such a kernel, the convolution (1) reads With R i j the j th pixel of the reference image in the neighborhood centered on the i th pixel, K j are the values of the pixelated kernel, and N is the extent of the kernel.Now we can cast this last equation as a least squares problem by requesting that the blurred reference should be equal to the target image.
Each column of the matrix corresponds to a translation of the reference image, and each line corresponds to a pixel of the reference image.Eliminating bad pixels from the problem is thereby very easy, we just need to remove the corresponding line from the matrix equation.We also further refine the problem by dividing both the reference image and the target image with their respective noisemaps before building the matrix.Finally, we can use the closed-form solution of Eq. ( 4) to obtain both the kernel values K i and background B which minimize the χ 2 of Eq. (2).This is much faster than using an iterative minimization algorithm.

Variability maps
Provided the difference images at all observing epochs have been perfectly computed, they can be combined into a variability map that we define as follows: where R blurred→i is the blurred version of the reference image such that it matches the PSF and background of the i th target image.In words, the variability map is the sum of the absolute difference images weighted by the corresponding noise maps.Note that we also adapt the resolution of the noisemap of the reference image, see Appendix A.3.The overall difference-imaging process is illustrated using Pan-STARRS data in Fig. 1.

Tests with Pan-STARRS data
We perform two types of tests of the method on real data.In the present case it is simply a low-pass filter coupled to a small translation accounting for the shift between the target and reference images.
quasars in Lemon et al. (2019) that fall in the Pan-STARRS footprint, we find that 133 show extended variability.This is encouraging performance, even though 34 of the variability maps were not reliable due to the high brightness of the targets themselves.In these cases, the transition kernel K is not properly constrained from the fainter stars in the field, causing a poor subtraction of the PSFs of the candidate and leaving non-PSF-like residuals.This experiment allows us to set the threshold of reliability of our method at a Gaia G-mag of 17, given the magnitude distribution of PSF stars available in practice.Second, we examine the case of apparent quasar+star pairs, which are the most common contaminants to lens searches.Fortunately, the fraction of variable quasars is much higher than the fraction of variable stars, as the majority of quasars tend to vary by more than 0.1 mag over two years (Cimatti et al. 1993).To benchmark the correct rejection of stars, we build a small selection of quasar+star systems by cross-matching the MILLIQUAS (Flesch 2021) and Gaia EDR3 (Gaia Collaboration et al. 2021) catalogues and imposing additional restrictions: we request a large enough proper motion significance (PMSIG, see Lemon et al. 2018 for the definition) of the companion to make sure we are not including another quasar, and a small Gaia G-magnitude difference between the two.Ideally, the resulting systems should all produce variability with at most one PSF visible.This is achieved for 164 of the 170 test-systems.Among the remaining 6, 3 do not have sufficient Pan-STARRS coverage to perform difference imaging, while for the other 3 the difference imaging performs poorly as there are too few reference stars available even when using large cutouts.In such cases the variability maps wrongly show variability.This goes to show that the method is able to efficiently reject quasar+star pairs, although care must still be taken to assess the quality of the difference imaging.

Candidate Selection
To maximize the chances of lens discovery we first build a catalogue of sources that are likely to contain a quasar.To do this, we combine the Gaia (Brown et al. 2018) and WISE (Wright et al. 2010) catalogues.We start with all the WISE detections that are red enough, that is the WISE color fulfils: This criterion is stricter than that used by Lemon et al. (2019), but helps reducing the number of candidates to a manageable number 999000999, page 3 of 11 for this proof of concept at the cost of losing bluer candidates.Then, among the red WISE detections we select only those that have a first Gaia detection within 3 , and a second Gaia detection between 0.8 and 3 away from the first one.All the detections with a PMSIG larger than 3 are discarded.Also eliminated are all the detections within 15 • of the galactic plane to reduce the likelihood of selecting stars.Given the above criteria, our final selection consists of 14 591 objects.
Next, we make use of all the available Pan-STARRS cutouts in the r and i band of the objects selected above.The choice of the r band is motivated by the fact that quasars are more variable in bluer bands (e.g.MacLeod et al. 2010), but the g band typically comes with harder-to-model PSFs.The i band is also included to patch the regions where the r band is not available.The analysis could not be performed for 484 candidates due to insufficient data, as fewer than two epochs were available for these.Consequently, a total of 14 107 variability maps are generated using our differenceimaging method.Out of these, 9 692 maps are discarded due to the absence of any discernible signal above the background noise level.This leaves 4 415 maps requiring visual inspection.We exclude 3 300 more that show a single PSF-like structure in the variability maps but no extended variability.For the remaining maps displaying extended variation, we make sure that the stars in the rest of the field are indeed properly subtracted.This is not the case for 720 more maps, in which the difference-imaging procedure fails.This is usually caused by a lack of bright enough stars in field, leading to a poor constraint of the transition kernel K, or other defects: see Appendix A.6.Lastly, we inspect the Pan-STARRS gri color stamps of each candidate, giving priority to those showing consistent color across PSFs.We exclude 166 more objects which evidently do not look like lenses.Ultimately we identify a total of 229 variable objects, including 20 known lenses.From this subset, we further narrow down our selection to 14 particularly promising candidates for spectroscopic follow-up.

Spectroscopy
Eight systems were observed with grism #13 with the ESO Faint Object Spectrograph and Camera 2 (EFOSC2) on the ESO New Technology Telescope (NTT), over two runs in October 2020 and January 2021 3 , providing a dispersion of 2.77Å per pixel.Each object received between 900 and 1 200 seconds of exposure time.Four other systems were observed with grism #4 and the Alhambra Faint Object Spectrograph and Camera (ALFOSC) on the Nordic Optical Telescope (NOT) in April 2021 4 , providing a dispersion of 3.3Å per pixel.Each object received between 900 and 1 200 seconds of exposure time.
3 Timo Anguita, 106.218K.001and 106.218K.002. 4Cameron Lemon, The general spectral reduction is carried out exactly as described in Lemon et al. (2023).The key steps are as follows.Each 2D long-slit spectrum image is bias-subtracted, and cosmic rays are masked using a Laplacian threshold method.The threshold is adjusted based on the brightness of the objects and a visual inspection of the cosmic ray mask.Next, the sky background is subtracted at each wavelength by determining the median value within the pixels surrounding the object trace, providing a skysubtracted, debiased, cosmic-ray-masked 2D long-slit spectrum image.A corresponding 2D noise map is then generated, combining a Poisson noise map derived from the detector gain with the sky noise.Finally, a wavelength solution is established by fitting lines from a pre-observation arc exposure (HeNe or CuAr), and linearly shifted for each exposure to fit several major sky background emission lines.At each wavelength, a model of the 1D spatial profile of the two components is generated as two Moffat profiles separated by the Gaia separation, and integrated perpendicularly in the slit.The Moffat parameters as a function of wavelength are then determined by fitting this model to several bins across wavelength.Finally the flux of each component can be derived at each wavelength by using the best-fit Moffat parameters at each wavelength.

Imaging
The Pan-STARRS data are limited in depth when trying to detect any potential lens galaxy between the quasar pairs.Given the narrow separations and faint magnitudes of our candidates, we obtained deeper images of the 7 Southern double-type candidates in the R c band using the Wide Field Imager (WFI) on the Max-Planck-Institut für Astrophysik (MPIA) 2.2-meter telescope (2p2), in nights of excellent seeing between March and June 2022.We obtained a total exposure time of 4 × 320 seconds per target.
Due to the complex PSF in these new images, we build a PSF model oversampled by a factor of 3, using up to 6 stars in each field.Next, for each object we perform a joint deconvolution, fitting the positions and amplitudes of two point sources as well as a Starlet-regularized pixelated background.Because galaxies are known to be sparse in Starlet space, we expect the lensing galaxy to naturally emerge in the background, if detectable in the dataset.We use STARRED (Michalewicz et al. 2023) to perform both the PSF modelling and the image deconvolution.The result is shown in Fig. 4.
When deeper Legacy Survey data or Hyper Suprime-Cam Subaru Strategic Program (HSC, Aihara et al. 2022) data are available, we perform a similar analysis but we replace the pixelated background with an analytical Sersic (Sérsic 1963) 1).At this depth and filter, the lensing galaxy is detectable only in the two brightest objects.
vastly superior signal-to-noise.We make the code developed for this analysis publicly available as well5 .

Discussion: quasar pairs
Figure 2 shows the Pan-STARRS gri stacks and variability maps of all our double-type candidates for which we have follow-up imaging, ordered by right ascension.The variability of PS J0125−1012 clearly suffers from of an artifact in one of the cutouts, but we still select it given the general look of the Pan-STARRS gri cutout.We show the acquired spectra in Fig. 3.Among the 12 double-type candidates, 8 are spectroscopic quasar pairs (similar spectra in both components), 1 is a projected quasar pair, 2 are quasar+star pairs, and 1 is a pair of stars.Below we tentatively classify the spectroscopic quasar pairs with comments referring to the spectra of Fig. 3, the 2p2 images of Fig. 4, and the available survey data of Fig. 5.

PS J0125−1012 −
The NTT/EFOSC2 resolved spectrum shows two quasars at redshift z = 1.22 with very similar continua and broad emission lines.We jointly fit the iz Legacy Survey and rzy Pan-STARRS bands.The data is well modelled with only two point sources, and the resulting image separation is 1.12 , compatible with the Gaia measurement.Although the reduced χ 2 of the fit can be slightly improved by allowing for an additional Sersic profile, the improvement is not significant and the point sources are pushed further apart, breaking the compatibility with the Gaia separation.Considering the similar 50 OM10 mock doubles, i.e. with source redshifts within ∆z = 0.1 and image separations within 0.1 , their faintest lens galaxy magnitudes are i ∼ 21.1.This should be detectable at the depth of our Legacy Survey + Pan-STARRS cutouts.If this system turns out to be a dual quasar, the physical separation would be 9.3 kpc.We suggest this system is a dual quasar, but designate this system as an Unclassified Quasar Pair (UQP) until deeper high-resolution imaging can be obtained.
PS J0557−2959 − The NTT/EFOSC2 resolved spectrum shows two quasars at redshift z = 2.04 with identical continua and broad emission lines.Even though our shallower r-band 2p2 data is well fitted with two PSFs, the same fit of Legacy Survey data shows a clear sign of a lensing galaxy in the z-band data.This residual subsequently disappears when including a Sersic profile, which falls between the quasars.We also note that the PSF separation is 0.783 when fitting without a galaxy, and 0.808 when fitting with a galaxy, the latter of which agrees better with the Gaia separation of 0.810 , providing more evidence for the existence of the lensing galaxy.We therefore conclude that this system is a lensed quasar.
PS J0909−0749 − The NTT/EFOSC2 resolved spectrum shows two quasars at redshift z = 1.08 with distinct differences in the ratio of the continuum and broad lines.There is no evidence for a lensing galaxy in the WFI 2p2 imaging.The simultaneous fit of Legacy Survey iz and Pan-STARRS grizy data with two point sources is good, and does not significantly improve when adding a Sersic profile.The Gaia separation is also recovered without the need for a model containing a galaxy.Moreover, the 38 similar Table 1.Confirmed quasars pairs.Are given the separation between the two images from Gaia, the spectroscopic redshift and the Gaia G-magnitudes.
When applicable, we also indicate the separation recovered from the models in Figs. 4 and  Given the slight differences in spectra and lack of lensing galaxy, we suggest this system is a dual quasar with physical separation 6.6 kpc.However we conservatively designate this system as an UQP, needing deeper high-resolution imaging.
PS J1019−1322 − The NTT/EFOSC2 resolved spectrum shows two quasars at redshift z = 2.32 with very similar emission line profile shapes.The flux ratio across wavelength is ∼2.2, with some small variations around the emission lines.There is no evidence for a lensing galaxy in the WFI 2p2 imaging.Similarly, the simultaneous fit of the Pan-STARRS grizy and Legacy Survey griz data with two point sources does not leave significant residuals, and the Gaia separation is recovered as well.Similar mocks from the OM10 catalogue can have galaxy magnitudes as faint as i ∼ 22.6.We suggest this system is a lensed quasar given the similarity of the spectra, however a deep image detecting the lensing galaxy is still needed.We designate this system as an UQP.
J0125-1012 J0557-2959 J0909-0749 J1019-1322 J1037+0018 J1225+4831 zir zir y zi zig irg irg PS J1037+0018 − The NTT/EFOSC2 resolved spectrum shows two very similar quasars at redshift z = 2.46.The flux ratio is approximately 5.5 across all wavelengths, and both spectra show a damped Lyman-alpha system.The available HSC images are poorly fit with two point sources, but the subsequent introduction of a Sersic profile produces much cleaner residuals.Subtracting the two point sources then reveals the lensing galaxy close to the fainter point source.We designate this system as a lens.
PS J1225+4831 − This system has also been discovered by Yue et al. (2023), who obtained Large Binocular Telescope LUCI1 imaging with an effective PSF of 0.35 (FWHM).The images show no evidence for a lensing galaxy, and they further note a possible flux ratio difference between the Hβ and [Oiii] lines.We present resolved NTT/EFOSC2 spectra of the system, which show two quasars at redshift z = 3.09 with a reasonably constant flux ratio of 3 across wavelength.Similar systems in the OM10 mock catalogue have lensing galaxies as faint as i ∼ 23.7.One possible explanation of these observations is that the lensing galaxy falls very close to one of the images and is being fit by the PSF subtraction of Yue et al. (2023), and the fainter component's spectrum is too noisy to conclude incompatibility with the brighter object's spectrum.Additionally, jointly fitting the available Pan-STARRS and Legacy Survey data produces slightly prominent residuals with a 2-PSFs model, but the color of the said residuals matches that of the point sources.We thus attribute the poor fit to a defect in our PSF model, and thereby cannot find evidence for a lensing galaxy either.If this system were a dual quasar, its physical separation would be 6.87 kpc.We designate this system as an UQP, needing deeper imaging to confirm the presence or absence of the lensing galaxy.
PS J1800+0146 − The NOT/ALFOSC resolved spectrum shows two very similar quasars at redshift z = 1.53.The flux ratio is roughly one across wavelength.Here the 116 similar OM10 mocks suggest a lensing galaxy as faint as i ∼ 21.6, which is relatively bright and was blindly picked up by starlets in the 2p2 data.The Gaia separation of 1.12 is also recovered by the fit.We designate this system as a lens.
PS J2122−1621 − The NTT/EFOSC2 resolved spectrum shows two very similar quasars at redshift z = 1.47, with a flux ratio of ∼ 1.2 across wavelength.Just like in the case of 999000999, page 7 of 11 A&A proofs: manuscript no.diffimg PS J1800+0146, the lensing galaxy could be fitted blindly with STARRED in the 2p2 data.The Gaia separation of 1.20 was recovered as well.We designate this system as a lens.

Conclusions
Our search using difference imaging in Pan-STARRS has led to the identification of four new doubly lensed quasars and four Unclassified Quasar Pairs that are either lenses or dual quasars, and one projected quasar pair.Notably, all of the lenses have narrow separation between 0.81 and 1.24 , and three have source redshifts above z = 1.5.This distribution is thereby starting to fill the gap pointed out by Lemon et al. (2023), where many small separation high redshift doubles are missing.Furthermore, Lemon et al. (2020) note that there is a quad bias when variability is used to select lensed quasars, which can be understood by the fact that the more variable quasars are the least massive ones (MacLeod et al. 2010): these have fainter absolute magnitudes and would need to be highly magnified to reach a detectable threshold.Thereby searching for lensed quasars through their variability introduces a high-magnification bias, and is complementary to traditional search methods for the construction of a magnitude-limited lensed quasar catalogue.So, given the bias for strong magnification that our technique introduces, it may seem surprising that we do not find more quads.This is simply due to the fact that all quads were already found at the Gaia limit (Lemon et al. 2018;Delchambre et al. 2019;Lemon et al. 2023).As we extend the search to beyond the Gaia limit, we can expect that this time, quads will be detected at an increased rate.
In conclusion, we have presented a difference imaging technique and applied it to Pan-STARRS imaging with a Gaia+WISE pre-selection, focusing on detecting extended variability as an efficient method for the selection of lensed quasars.This work underlines the potential of extended variability and difference imaging as an effective method for lensed quasar selection, with great promise of uncovering highly magnified lenses in future surveys such as the forthcoming Vera Rubin Observatory's Legacy Survey of Space and Time (Ivezić et al. 2019, LSST).Indeed, considering the anticipated scarcity of spectroscopic resources in the future relative to the vast number of candidates expected to be uncovered by LSST, attaining a high level of purity in candidate selection through their variability will become increasingly important.Filho de Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cientifico e Tecnologico and the Ministerio da Ciencia, Tecnologia e Inovacao, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey.The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energeticas, Medioambientales y Tecnologicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciencies de l'Espai (IEEC/CSIC), the Institut de Fisica d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig Maximilians Universitat Munchen and the associated Excellence Cluster Universe, the University of Michigan, NSF's NOIRLab, the University of Nottingham, the Ohio State Univer-

Acknowledgments
extract the Dark Energy Survey multi-epoch photometry of lensed quasars, quasar pairs, and quasar-star projections, and show that a simple cut in 999000999, page 1 of 11 arXiv:2307.13729v1[astro-ph.GA] 25 Jul 2023

Fig. 1 .
Fig. 1.Illustration of the difference-imaging process using a 65 stamp of Pan-STARRS data.Panel (A) shows an image obtained at any given epoch (I(x) in the text or target image) while panel (B) shows the best-seeing image of the time series, i.e. the reference image.Panel (C) shows the pixels that are used in the least square procedure.Note that saturated pixels are not considered.The blurred reference image is show in panel (D) and panel (E) shows the residuals after subtracting the blurred reference from the target image.The kernel used to blur the reference image is shown in panel (F).In the present case it is simply a low-pass filter coupled to a small translation accounting for the shift between the target and reference images.

Fig. 2 .
Fig.2.Pan-STARRS gri views (top) and r+i variability maps (bottom) of the candidates selected for spectroscopic follow up.Extended variability is, albeit sometimes barely raising above the noise, still well visible for all the candidates.Those that are found to be false positives are displayed in gray, whereas the quasar pairs are displayed in white.Note that the variability map of PS J0125−1012 has an artifact.

Fig. 3 .
Fig.3.Long slit spectroscopy of double-type candidates performed using EFOSC2 on the ESO/NTT or ALFOSC on the Nordic Optical Telescope.The candidates are arranged in ascending order of right ascension, with lenses or quasar pairs written in bold and stars or projected quasars in gray.Shaded regions indicate wavelength ranges affected by atmospheric absorption.The spectra were reliably de-blended using the procedure outlined in Section 4.1.

Fig. 4 .
Fig. 4. First row: MPIA 2.2-meter telescope R c -band data (4×320 seconds).Second row: modeling of the systems using the STARRED algorithm.Third row: residuals after subtractions of the point-source component of the STARRED model.Last row: residuals after subtraction of the complete STARRED model (PSFs + background).Each row shows a stack, but the data was jointly deconvolved.The red crosses indicate the positions of the PSFs, the separations of which are close to the Gaia separations down to 0.01 (Table1).At this depth and filter, the lensing galaxy is detectable only in the two brightest objects.

Fig. 5 .
Fig. 5. First row: Legacy Survey, Pan-STARRS (italicized) or HSC (bolded) data.Second row: Models with 2 PSFs, and also an additional Sersic profile for J0557-2959 and J1037+0018.Third row: Residuals after subtraction of the models.Fourth row: When the model has a Sersic profile, residuals after subtraction of the PSFs only.The bands are arranged by wavelength (redder first), highlighting the redder hue of the lensing galaxies seen in the fourth row compared to the quasar pairs.
This program is supported by the Swiss National Science Foundation (SNSF) and by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (COSMICLENS: grant agreement No 787886).TA acknowledges support from the Millennium Science Initiative ICN12_009 and the ANID BASAL project FB210003 2p2: The 2p2 imaging was made possible by an agreement between MPIA and EPFL, the ESO program ID is 0108.A-9005(A).