Lensed quasar search via time variability with the HSC transient survey

Gravitatonally lensed quasars are useful for studying astrophysics and cosmology, and enlarging the sample size of lensed quasars is important for multiple studies. In this work, we develop a lens search algorithm for four-image (quad) lensed quasars based on their time variability. In the development of the lens search algorithm, we construct a pipeline simulating multi-epoch images of lensed quasars in cadenced surveys, accounting for quasar variabilities, quasar hosts, lens galaxies, and the PSF variation. Applying the simulation pipeline to the Hyper-Suprime Cam (HSC) transient survey, we generate HSC-like difference images of the mock lensed quasars from Oguri&Marshall's lens catalog. We further develop a lens search algorithm that picks out variable objects as lensed quasar candidates based on their spatial extent in the difference images. We test our lens search algorithm with the mock lensed quasars and variable objects from the HSC transient survey. Using difference images from multiple epochs, our lens search algorithm achieves a high true-positive rate (TPR) of 90.1% and a low false-positive rate (FPR) of 2.3% for the bright quads with wide separation. With a preselection on the number of blobs in the difference image, we obtain a TPR of 97.6% and a FPR of 2.6% for the bright quads with wide separation. Even when difference images are only available in one single epoch, our lens search algorithm could still detect the bright quads with wide separation at high TPR of 97.6% and low FPR of 2.4% in the optimal seeing scenario, and at TPR of $\sim94%$ and FPR of $\sim5%$ in typical scenarios. Therefore, our lens search algorithm is promising and could find new lensed quasars in ongoing and upcoming cadenced surveys, such as the HSC transient survey and the Large Synoptic Survey Telescope. [abridged]


Introduction
Gravitationally lensed quasars are powerful tools to study astrophysics and cosmology. For example, they can be used to probe substructure of dark matter (e.g., Mao & Schneider 1998;Metcalf & Madau 2001;Dalal & Kochanek 2002;Nierenberg et al. 2014;Gilman et al. 2019), to examine the formation and evolution of supermassive black holes (Fan et al. 2019), and to study the properties of quasar host galaxy (e.g., Peng et al. 2006a;Ding et al. 2017). With the time delays between multiple images, we could also measure the Hubble constant, H 0 , which is crucial for testing cosmological models and studying dark energy (e.g.; Refsdal 1964;Schechter et al. 1997;Suyu et al. 2014;Bonvin et al. 2017;Birrer et al. 2019;Wong et al. 2019;Chen et al. 2019;Rusu et al. 2019).
Given the importance of lensed quasars, there have been several systematic searches of lensed quasars. The Cosmic Lens All Sky Survey (CLASS; Myers et al. 2003;Browne et al. 2003) discovered a large sample of lensed quasars in radio wavelengths by first selecting radio sources that have compact flat-spectrum and further looking for the high-resolution multiply imaged components in the radio follow-up. In the optical, the SDSS 1 Quasar Lens Search (SQLS; Oguri et al. 2006;Inada et al. 2008Inada et al. , 2010Inada et al. , 2012More et al. 2016a), which have found the largest sample of lensed quasars (∼60), started 1 SDSS, Sloan Digital Sky Survey (York et al. 2000) Article number, page 1 of 15 arXiv:1910.01140v1 [astro-ph.GA] 2 Oct 2019 A&A proofs: manuscript no. lens_search_time_variability with spectroscopic confirmed SDSS quasars. To these, they applied a morphological selection for potential lens candidates with narrow separation and a color selection for potential lens candidates that are deblended in the SDSS, and then conducted follow-up observations to confirm after a visual inspection of the objects selected by the morphological or color selection. In the present stage, further astrophysical or cosmological studies with lensed quasars are still limited by the small sample size of lensed quasars, given the fact that lensed quasars are rare.
Fortunately, many ongoing and upcoming wide-field surveys provide a great pool of new lensed quasars, owing to the increased depth and much larger areal coverage. Moreover, with the aim to exploit the unprecedented advantage of those surveys, many dedicated searches were initiated with new lens search techniques. Some lens search techniques explore multi-band catalogues and employ various cuts in the magnitude and color spaces to find lensed quasars (e.g., Agnello et al. 2015;Agnello 2017;Ostrovski et al. 2017;Williams et al. 2017;Rusu et al. 2018). Some other techniques look for lensed quasars by examining their configuration or recognizing their pattern in images (e.g., Agnello et al. 2015;Chan et al. 2015). Although not specific to lensed quasars, Space Warps (Marshall et al. 2016;More et al. 2016b) show that lensed quasars could also be found through citizen science.
With the exceptional resolution of Gaia 2 , one could conduct quasar lens search by looking for multiple detection in Gaia or comparing the flux and position offsets from other surveys for objects in SDSS, Pan-STARRS 3 (e.g., Lemon et al. 2017Lemon et al. , 2018Lemon et al. , 2019Ostrovski et al. 2018), DES 4 (e.g., Agnello et al. 2018a), or KiDS 5 ). Using Gaia multiple detection to search for lensed quasars could also be applied to the ATLAS 6 footprint (Agnello et al. 2018b). One could also combine Gaia detections with other lens searches to find lensed quasars.
In the near future, thousands of lensed quasars are expected to be found in LSST 7 (Oguri & Marshall 2010, hereafter OM10), enlarging the sample size of lensed quasars by at least an order of magnitude. Moreover, LSST is a cadenced survey that will reveal the time variability of variable objects, such as lensed quasars. In a cadenced survey, difference imaging is helpful in eliminating non-variable objects. As first suggested by Kochanek et al. (2006), looking for extended objects in the difference images is an effective means to find lensed quasars, since almost all the other variable objects (non-lensed variable objects) are point sources. Therefore, a time-variability based method will benefit lens searches given the enormous number of objects in the LSST.
Inspired by Kochanek et al. (2006), we develop a lens search algorithm for lensed quasars based on their time variability. The HSC 8 transient survey (Yasuda et al. 2019) provides a great opportunity to simulate the difference images that could be used 2 Gaia (Gaia Collaboration et al. 2016) 3 Pan-STARRS, Panoramic Survey Telescope and Rapid Response System (Chambers et al. 2016) 4 DES, Dark Energy Survey (Sánchez & Des Collaboration 2010) 5 KiDS, Kilo-Degree Survey (de Jong et al. 2013) 6 ATLAS, Asteroid Terrestrial-impact Last Alert System (Tonry et al. 2018) 7 LSST, Large Synoptic Survey Telescope (Ivezic et al. 2008) 8 HSC, Hyper-Suprime Cam (Miyazaki et al. 2012;Aihara et al. 2018) for the development of the lens search algorithm and to test the lens search performance. The HSC transient survey is an ongoing cadenced survey of the Subaru Telescope (Miyazaki et al. 2018a) with similar image quality expected for LSST.
In this work, we build a simulation pipeline for producing time-varied images of lensed quasars, and apply the simulation pipeline to the HSC transient survey to generate HSC-like difference images of mock lensed quasars. With the HSC-like difference images of lensed quasars, we develop a lens search algorithm that picks out variable objects with large spatial extent in the difference images and classifies them as lens quasar candidates. We further test the performance of our lens search algorithm in the HSC transient survey. Although the simulation and the search algorithm could also be applied to lenses with two-image configurations (double), we focus on lenses with four-image configuration (quad). Quads provide more constraints than doubles for further studies, despite that there are more doubles than quads.
The organization of the paper is as follows. In Sec. 2, we present a new simulation pipeline of time-varied lensed images, and apply the new simulation pipeline to the HSC transient survey in Sec. 3. The lens search algorithm is detailed in Sec. 4. The lens search performance is shown in Sec. 5 before we conclude in Sec. 6.

Simulation of time-varying lensed quasars images
We present a simulation pipeline for creating realistic images of mock lensed quasars in a series of time. The mock images produced by this pipeline are useful for developing the time-variability-based lens search algorithm. We use the mock configuration from the catalog created by OM10 to generate the simulated lensed quasars. In order to have realistic simulations, we simulate not only the light from the variable object -the lensed quasar in this work -but also the light from the (lensed) host galaxy and lens galaxy. This also takes into account possible residuals on the difference image due to potential imperfection of the difference imaging method.

Mock lenses from OM10
OM10 have predicted the distribution of lensed quasars, and produced a mock catalog of lensed quasars. For each lens system, OM10 provides the source position (η 1 , η 2 ), the source redshift z s , the unlensed i-band magnitude m s of the quasar, the lens redshift z d , the velocity dispersion σ, the ellipticity e, and the position angle (PA) θ e of the lens galaxy. OM10 also provides the lensed image positions, the magnifications, and the time delays calculated by glafic (Oguri 2010), with the assumption of a singular isothermal ellipsoid (SIE) model (Kormann et al. 1994) for the lens and an external shear accounting for the effect of the lens environment (e.g. Kochanek 1991;Keeton et al. 1997;Witt & Mao 1997). The convergence κ from the SIE model of the lens is given by where (θ 1 , θ 2 ) is the coordinate of the lens, θ Ein is the Einstein radius in arcsec, e is the ellipticity, and λ(e) is the dynamical normalization defined in Oguri et al. (2012). The lens potential of the external shear φ is given by where γ is the magnitude of the external shear and θ γ is the orientation of the external shear.

Quasar
For the point sources, i.e. the lensed quasars without the (lensed) host galaxies, we use the source redshift z s , the unlensed i-band magnitude m s , the lensed image positions, the magnifications, and the time delays provided by OM10, to generate the simulated image for each lensed quasar.
We generate the source light curve of quasars using carma_pack (Kelly et al. 2009(Kelly et al. , 2014) based on z s and m s from OM10. carma_pack is developed to quantify stochastic variability, especially for quasar variability, and the quasar light curves are generated with the fitting relation on redshifts z in Kelly et al. (2009). We then shift (in time) and magnify the generated source light curve according to the time delay and the magnification of each lensed image. 9 Fig. 2 is an illustration of the procedure to generate light curves for lensed quasars without the (lensed) host galaxies. In the left panel, we show the positions of the four multiple images of a symmetry quad lensed by a galaxy at z d = 0.23 from OM10 with source redshift z s = 3.26. The four multiple image positions, A, B, C, and D in the left panel, are relative to the lens galaxy, which is at the center for each lens system in OM10, (x, y) = (0, 0) (black cross in the left panel). Their magnifications, |µ A |, |µ B |, |µ C |, and |µ D |, and the time delays relative to Image B, ∆t AB , ∆t CB , and ∆t DB , are listed in Table 1 (where ∆t XB = t X − t B ). The source light curve that we generate with carma_pack based on the source redshift z s = 3.26 is shown in the top-right panel of Fig. 2. The bottom-right panel of Fig. 2 shows the light curves of the four multiple images (A, B, C, and D in the left panel), which are produced by shifting and magnifying part of the source light curve in the top-right panel (blue band) with the corresponding time delays and magnifications in Table 1.

Quasar host galaxy
Since OM10 provides only the configuration of point sources and there is no definite morphological relation between quasar 9 We do not explicitly include microlensing variability. The lensed quasar variability is changed but not canceled out by microlensing variability, so lensed quasars will still appear as multiple-point like variable objects in the difference images even when there is microlensing effect.

Property
Min. Max. log L host (L ) 9 13 e 0 0.8 θ e (deg) 0 180 r e (kpc) 1 10 n 1 5 Table 2: Ranges of the quasar host galaxy properties -luminosity L host , ellipticity e, PA θ e , effective radius r e , and Sérsic index n. and the host galaxy with little scatter, we randomly draw the morphological properties based on the Sérsic profile for the host galaxy of quasar and use glafic (Oguri 2010) to generate the lensed image of host galaxy. We conservatively adopt broad ranges for the parameters in the Sérsic profile of the host galaxy (based on, e.g., Blanton & Moustakas 2009;Peng et al. 2006b;Park et al. 2015;Bennert et al. 2010). The thresholds of the host galaxy properties used in generating the lensed images are listed in Table 2.
We put the host galaxy with its random Sérsic profile at the same source position as the (unlensed) quasar provided by OM10, and we adopt the same SIE model, external shear, and cosmology from OM10 for glafic to generate the lensed image of the host galaxy. Therefore, the lensed host galaxy image positions will be consistent with the lensed quasar image positions.

Lens galaxy
For the lens galaxies, we assign the light profile to the lens galaxy for each lens system with fitted profile of galaxy from a real survey, as only the mass model of lens galaxy is given by OM10. We denote such a real survey "S imag ", and S imag could be any survey where one is interested in finding lens systems, such as the Sloan Digital Sky Survey (SDSS) or the HSC survey.
We begin with a catalogue created by cross-matching the galaxies from S imag with another survey "S spec ", which provides spectroscopic redshift z spec and velocity dispersion σ spec . We note that if S imag itself has information about spectroscopic redshift and velocity dispersion, one does not have to cross-match S imag with another survey. From the cross-matched catalogue, we select "candidates" for each lens galaxy with the closet values to the lens redshift z d,OM10 and the velocity dispersion σ OM10 from the mass model provided by OM10. We first select the galaxies with and where d z = z spec − z d,OM10 , d σ = σ spec − σ OM10 , and ∆σ is the error in σ spec from S spec . We then get the candidates from the selected galaxies with the smallest values of d z,σ , where d z,σ is defined as Fig. 1: Flowchart illustrating the simulation of time-varying lensed quasars images. For quasars, we generate mock quasar light curves based on their source redshifts z s , shifting in time and magnifying the mock quasar light curves on the image plane (see Fig. 2). For quasar hosts, we create their Sérsic light profiles on the source plane, and produce their lensed images on the image plane assuming a SIE model. For lens galaxies, we pick galaxies from a real survey based on their lens redshifts z d and velocity dispersions σ, fitting the selected galaxies with Sérsic light profiles, and using the fitted Sérsic light profiles to produce their images on the image plane. Once we have all the components for one lens system (image of lensed quasar, image of lensed quasar host, and image of lens galaxy), we convolve them with a realistic PSF and add them to obtain one complete image.
as candidates for the lens light. However, for some lens systems in OM10, the number of selected galaxies from Eqs. 3 and 4 is not enough to have ten candidates. For such cases, we take all the available galaxies passing the selection criteria (Eqs. 3 and 4) as candidates. In this work, we use the data from the Sloan Digital Sky Survey Data Release 14 (SDSS DR14, Abolfathi et al. (2018)) as S spec .
The velocity dispersion σ SDSS is limited in the crossmatched catalogue when z SDSS > 1.04, so we apply only Eq. 3 to select candidates for the lens galaxy with 1.04 < z d,OM10 < 1.57; we select up to ten galaxies as the candidates with the smallest values of d z . Due to the lack of object with z SDSS > 1.57 in the cross-matched catalogue, we do not simulate lens system with z d,OM10 > 1.57. We then obtain the image for each candidate galaxy from S imag , and fit the light distribution of the candidate galaxy with the Sérsic profile. For the fitting, we define and the reduced χ 2 as where N p is the number of pixels used in fitting, I i,fit is the fitted intensity of light, I i,data is the observed intensity of light from S imag , and σ i,data is the noise from S imag . We rank the candidates by χ 2 reduced . The smaller χ 2 reduced is, the higher the rank is.
After we have the fitted light profile for each candidate galaxy, we compare the ellipticity from the fitted light profile, e fit , and the ellipticity from the mass model given by OM10, e OM10 , according to the ranking order. For one candidate, if |e fit − e OM10 | < 0.2, we will pick the candidate as the lens galaxy (Newton et al. 2011); otherwise, we check e fit of the candidate in the next lower rank. For the lens galaxy that has no candidate with the satisfied value of e fit , we pick the candidate with the smallest value of χ 2 reduced . In total, we simulate 2033 quads with lens redshift z d,OM10 < 1.57 from OM10.

Simulated image with all components
Once we have the individual simulated parts of a whole lens system (the lensed quasar, the lensed host galaxy, and the lens galaxy), we will convolve them with the point-spread-functions (PSFs) in S imag for each epoch and add the parts to become one complete image. The complete image in an epoch t is expressed as where I t (x, y) is the light distribution at position (x, y) of the complete image, I q,t (x, y) is the light distribution of lensed quasar, I h (x, y) is the light distribution of (lensed) host galaxy, I d (x, y) is the light distribution of lens galaxy, PSF t (x, y) is the PSF from S imag , and ⊗ represents the convolution. After we generate complete images for all the epochs in S imag , we will have the images of the mock lensed quasar in the time series of S imag .

Application to the Hyper Suprime-Cam Survey
In this section, we demonstrate our simulation in the HSC transient survey, and describe how we generate realistic difference images that will be used to develop the search algorithm through the pipeline of the HSC transient survey.
To obtain the difference image for each epoch, we need the reference image to subtract from the single-epoch image. In the HSC transient survey, the reference images are the deep reference images produced by coadding multiple exposures from multiple epochs during March 2014 to April 2016 (Table 4). Only the exposures with seeing better than 0.7 arcsec have been used in creating the deep reference images.
The HSC transient survey uses the method in Alard & Lupton (1998) and Alard (2000) for the difference imaging. The single-epoch image is created by coadding several warped images, and each warped image corresponds to a distortioncorrected image of the sky from a single exposure. For each epoch, the difference imaging has been performed on every warped images by subtracting the deep reference image to produce the warped difference images, and all the warped difference images are then coadded to create the deep difference image.
To generate realistic deep difference image of the mock lenses, we select the locations that have no object within 5 arcsec, which we call "empty regions", and we inject our mock lenses into these empty regions with the HSC pipeline (Bosch et al. 2018) version 4.0.5, for both the warped images in each epoch and the deep reference image. Once the injection is done, we use the difference imaging method mentioned above to create the deep difference image of the mock lenses. Practically, we select 75 empty regions in the "patch" of sky identified with the HSC Tract = 9813, Patch = (3, 4), to inject the mock lenses. In this patch, 44 warped images from five epochs are used to create the deep reference image (  all the 13 epochs in i-band.

Convolution with the HSC PSFs and lens injection
In order to obtain realistic deep difference image, we inject the mock lens into both the warped image and the deep reference image. For the warped images, the difference imaging method is sensitive to every warped images used to create the deep difference image for each epoch, so accounting for the variations in the PSFs of the warped images is crucial. Meanwhile, the deep reference image used for the subtraction is coadded, so weighting the variability from the dates in Table 4 plays a more important role. We detail the injection methods for the warped image and the deep reference image below.
For the lens injection into the warped images of one epoch, we first simulate the images for the whole lens system as mentioned in Sec. 2 accounting for the quasar variability in the epoch, and then we create the complete image by convolving with the PSFs of the empty region where the lens system is injected into the warped images. Since one epoch is composed of several warped images, the PSF of the same empty region changes across multiple warped images. Therefore, for one lens system in one epoch, the complete images are generated by convolving with different PSFs for different warped images, and we further inject the complete images into the corresponding warped images.
For the lens injection into the deep reference image, we first weight the images of the mock lens on an epoch basis by the number of the warped images used in one epoch, and unlike the lens injection into the warped images, we convolve the weighted image with the PSF from the deep reference image and further inject the mock lens directly into the deep reference image. In the deep reference image, the total light distribution at position (x, y), I ref (x, y), could be described as where t = 2014-03-28, ..., 2016-03-04 (Table 4), I q,t (x, y) is the light distribution of lensed quasar in t, N t is the number of warped images used in t, N total is the total number of warped images (N total = 44 in our work), I h (x, y) is the light distribution of (lensed) quasar host galaxy, I d (x, y) is the light distribution of lens galaxy, PSF ref (x, y) is the PSF from the deep reference image, and ⊗ represents the convolution. Injecting mock lenses into the deep reference image via this method takes into account the variability of lensed quasar across the multiple epochs used in the deep reference image.
We simulate the complete image for 2033 quads (z d,OM10 < 1.57) from OM10, and the size of each complete image is 10 × 10 with the lens at the center (top row in Fig. 3). The quads are injected into the empty region randomly, and the empty regions are used repeatedly for different quads. For simplicity, we do not add Poisson noise associated with the mock lens systems when injecting the quads into the empty region, given the small impact of the Poisson noise: We find that the sky background noise generally dominates over Poisson noise in single-epoch images, although the addition of Poisson noise would slightly affect the pixel counts for some bright or wide-separation lens systems (see Sec.5.1). We defer the inclusion of Poisson noise, and the propagation of noise from the single-epoch images into the difference images to future work. After the injection of the complete images, we apply the deep difference imaging method to obtain the deep difference image for the mock quads (bottom row in Fig. 3). As mentioned in Sec. 3.1, we first get the warped difference images of injected lenses as intermediate products, and the deep difference images are created by coadding the warped difference images. The top row in Fig. 3 are examples of the single-epoch images we obtain after the lens injection processed by the pipeline of the HSC transient survey, and we could see the light from both the lensed quasar and the lens galaxy in each panel, which appear like "bright spots". The bottom row in Fig. 3 are the corresponding deep difference images of the top row, and the "dark spots" (bright spots) are the spots that become fainter (brighter) than they were in the reference image. In contrast to the single-epoch images where we could only see bright spots and thus brightness changes are not obvious, brightness changes are clearly visible in the deep difference images through the dark spots (decreasing brightness) and the bright spots (increasing brightness). There are some residuals at the locations of the lens galaxies in the bottom panels of Fig. 3, which confirm that, simulating the light from all the components (quasar, host galaxy, and lens galaxy) instead of from only variable object (quasar) is important for realistic mock images, as the pipeline might not perform perfect subtraction for all the non-variable components.

Search method: spatial extent in difference image
With the multiple-image feature and the variable brightness, strongly lensed quasars exhibit multiple point-like image residuals in the difference image shown in Fig. 3. Therefore, lensed quasars appear as spatially extended in the difference image. Since most astrophysical variable sources are point-like (e.g., variable stars, unlensed quasars, unlensed supernovae), targeting spatially extended variable sources is an effective approach to find lensed quasars, as previously noted by Kochanek et al. (2006). We develop a lens search algorithm for lensed quasar via a method that (1) quantifies the extendedness of an object in the difference image and (2) selects objects that are large in spatial extent as lens candidates. We present the method that is based on the extendedness in Sec. 4.1, and an enhancement of the method with a secondary criterion in Sec. 4.2.

Quantification of the extendedness
The lens search algorithm starts from "HSC variables" in the HSC transient survey. In this work, we define a "HSC variable" as an object properly detected on the deep difference images at least twice in the HSC transient survey -the detections could be from two different epochs or two different bands. For each HSC variable, we first collect its deep difference images from all the 13 epochs (Table 3), and create the "3σ-mask" for each epoch by picking out the pixels with value larger than 3σ or smaller than −3σ. The pixel values in the 3σ-mask, I mask (i, j), are defined as where i = 1, ..., N x and j = 1, ..., N y are the pixel indices in both the deep difference image and the 3σ-mask of dimensions N x × N y (N x = N y = 59 in this work), I(i, j) are the pixel values in the deep difference image, and σ(i, j) are the estimated 1-σ uncertainties in the difference image from the HSC transient survey. Figs. 4a and 4b are examples of the deep difference image and the 3σ-mask from a HSC variable in an epoch, respectively. As indicated in Fig. 4b, 3σ-mask has several noise peaks in the outskirts that are not related to the HSC variable. Those noise peaks should not be counted in the extendedness of the HSC variable, so we further define the "effective region" by mock lensed quasar in the same epoch as Fig. 4a. Fig. 4f and 4g are respectively the 3 -mask and the e↵ective region of the mock lensed quasar in Fig. 4e. As shown in Fig. 4c and Fig. 4g, the injected mock lensed quasar has a much larger area of the e↵ective region.

Preselection by "number of blobs"
In addition to the extendedness, we could also apply the "blob-feature" as a secondary criterion to improve the search method. The image residuals of variable sources in the difference image have either positive pixel values or negative pixel values, depending on the brightness change, and the pixel values elsewhere should be zero, up to the noise level. The residuals with positive pixel values are like "white blobs", and the residuals with negative pixel values are like "black blobs". In the di↵erence image, these blobs are like "local extrema" in the "zero-background". Due to the multiple point-like image residuals, lensed quasars tend to have a larger number of "local extrema", compared to unlensed variable sources. Therefore, we enhance the lens search method in Sec. 4.1 with a preselection of HSC variable based on a criterion for the number of "local extrema".
We first define "local extremum" for the di↵erence image. A pixel (i, j) is a "local extremum", if and with its neighbouring pixels (i 0 , j 0 ) satisfying and The local extrema are marked by cyan in Fig. 4d and 4h. Fig. 4d and 4h keep the pixel values only for the pixels satisfying Eq. 14, and filter out the other pixels. We note that, with the definition in Eqs. 13-16, the local extrema do not completely correspond to the blobs. As shown in Fig. 4d, we lose the white blob in the left as a local extremum because the pixel supposed to be the local extremum, indicated as "edged extremum" in red, is located at the edge of the region that is e↵ectively related to the HSC variable. Furthermore, the local extrema could also come from the noise peaks, as marked by magenta in Fig. 4h. Despite that the local extrema and the blob are not in one-to-one  where I eff (i, j) is the pixel value of (i, j) in the effective region. Fig. 4c shows the effective region of the HSC variable in Fig. 4a.
Once we have the effective region, we determine the area of the effective region, A eff , by the sum of the pixel values in the effective region where i runs from 1 to N x , and j runs from 1 to N y . This area is equivalent to the number of pixels in the effective region where I eff = 1. As a result, we evaluate the extendedness of a HSC variable in the deep difference image for a given epoch by the area of this effective region. Doing so, we could discard the effect from the noise peaks in the outskirts and quantify the extendedness of HSC variable in the difference image.
We search for lensed quasars with all the 13 epochs. Assuming we have M HSC variables, we denote A m eff,t as the area of the effective region for a HSC variable in epoch t (Table 3), where m = 1, ..., M. For each epoch t, we compute the percentile p% for the area of the effective region from all the HSC variables, A eff,t (p): p% of M HSC variables with A m eff,t < A eff,t (p).
We note that, the computation of A eff,t (p) is based only on the HSC variables, independent of the mock lensed quasars. For each HSC variable m, we then count the number of epochs with the area of the effective region larger than A eff,t (p), N m epoch (p). For example, if a HSC variable has A m eff,t > A eff,t (p thrs ) at a threshold of percentile p thrs %, for t = t 1 , t 2 , and t 3 , we will say N m epoch (p thrs ) = 3. Given a threshold of percentile p thrs %, we set another threshold, N thrs , and a HSC variable will be classified as a candidate for lensed quasar if N m epoch (p thrs ) > N thrs .
We also perform the calculation in Eqs. 10-12 for each injected mock lensed quasar. Fig. 4e shows an example of injected mock lensed quasar in the same epoch as Fig. 4a. Figs. 4f and 4g are respectively the 3σ-mask and the effective region of the mock lensed quasar in Fig. 4e. As shown in Figs. 4c and 4g, the injected mock lensed quasar has a much larger area of the effective region.

Preselection by "number of blobs"
In addition to the extendedness, we could also apply the "blob-feature" as a secondary criterion to improve the search method. and 4h keep the pixel values only for the pixels satisfying Eq. 14, and filter out the other pixels. We note that, with the definition in Eqs. 13-16, the local extrema do not completely correspond to the blobs. As shown in Fig. 4d, we lose the white blob in the left as a local extremum because the pixel that is supposed to be the local extremum, indicated as "edged extremum" in red, is located at the edge of the region that is e↵ectively related to the HSC variable and thus fails the condition in Eq. 15. Furthermore, the local extrema could also come from noise peaks, as marked by magenta in Fig. 4h. Despite that the local extrema and the blobs are not in one-to-one correspondence, their numbers are still highly correlated.
Before we compute the percentile A e↵,t (p) at p% in Eq. 12, we count the number of the local extrema, N m extrm , for each HSC variable m (where m = 1, ..., M) in an epoch of choice. Depending on the imaging survey, this epoch could be for example the best-seeing epoch, or median-seeing epoch. In the HSC transient survey, the best-seeing epoch with a remarkable value of 0.42 arcsec is so good that the sharp images actually lead to significant artifacts in the deep di↵erence image that a↵ects the number of extrema. We therefore use the median-seeing epoch (with seeing of 0.72 arcsec) in the HSC transient survey for detecting the local extrema, as this is the regime where the di↵erence imaging pipeline runs well. Given a criterion, N crit , a HSC variable will be discard from the lens classification if N m 0 extrm  N crit , and a HSC variable will be kept in the lens classification if N m 0 extrm > N crit .
Assuming we have M 0 HSC variables with their number of the local extrema larger than N crit after the preselection (where M 0  M), we compute their area of the e↵ective region for each epoch t (A m 0 e↵,t , m 0 = 1, ..., M 0 ), and further calculate the percentile p% for the area of the e↵ective region from these M 0 HSC variables, A 0 e↵,t (p), in each epoch t, which is similar to Eq. 12. We notice that, at a same percentile p%, A 0 e↵,t (p) is generally larger than A e↵,t (p), the percentile p% from all the HSC variables before the preselection. The reason is that the HSC variables with fewer local extrema are also the HSC variables with fewer spatially extended blobs counting for the extendedness, so when we reject the HSC variables with small number of the local extrema, we also reject the HSC variables that are small in the area of the e↵ective region. Raising A e↵,t (p), the percentile p% for the area of the e↵ective region, with the preselection by the number of blobs makes the further lens classification more e cient and helps us to avoid large number of false candidates for lensed quasar, though we would lose a few candidates, particularly narrowly separated lensed quasars (Sec. 5.2).

Perfomance test
Sherry: since you've used "HSC variables" terminology in the previous section, let's use consistent terminology. So I've changed the "variable objects" below to "HSC variables".
In this section, we examine the performance of our lens search algorithm. To test the performance of the lens search algorithm, we need not only the deep di↵erence images of lensed quasars, but also the deep di↵erence images of the non-lensed objects. We exploit 12910 HSC variables from the HSC transient survey as the non-lensed objects to test our lens search algorithm. Like the 2033 injected mock lensed quasars, we take the deep di↵erence images from all the 13 epochs for each HSC variable.

Classification based only on spatial extent
The lens search algorithm performance is measured by the truepositive rate (TPR) and the false-positive rate (FPR), which are defined as and Local extremum detection Fig. 4: Examples of a HSC variable (top row) and an injected lensed quasar (bottom row). From left to right: difference image, the corresponding 3σ-mask, the corresponding effective region, and the corresponding local extremum detection. In 3σ-masks, the white pixels represent the pixels with values larger than 3σ or smaller than −3σ in the difference images. The effective regions marked by the white pixels show the spatial extent of the objects, after removing noise peaks in the 3σ-masks. In the local extremum detection, we pick out the pixels with values larger than all the neighbouring pixels or smaller than all the neighbouring pixels in the difference images within the effective regions; the local extrema indicated in cyan correspond to the "blobs" (See Sec. 4.2 for details) of the HSC variable/mock lensed quasar, the local extrema indicated in magenta are the noise peaks, and the pixel indicated as "Edged extremum" in red is a blob that should be picked out as local extremum but is missed because it is located at the edge of the effective region. The size of each image cutout is 10 × 10 .
as shown in Figs. 4a and 4e. In the difference image, these blobs are like "local extrema" in the "zero-background". Due to the multiple point-like image residuals, lensed quasars tend to have a larger number of "local extrema", compared to unlensed variable sources. Therefore, we enhance the lens search method in Sec. 4.1 with a preselection of HSC variables based on a criterion on the number of "local extrema".
We first define "local extremum" for the difference image. A pixel (i, j) is a "local extremum", if and with its neighbouring pixels (i , j ) satisfying and for i = i − 1, i, i + 1 and j = j − 1, j, j + 1, except (i , j ) = (i, j). The local extrema are marked by cyan in Figs. 4d and 4h. Figs. 4d and 4h keep the pixel values only for the pixels satisfying Eq. 15, and filter out the other pixels. We note that, with the definition in Eqs. 14-17, the local extrema do not completely correspond to the blobs. As shown in Fig. 4d, we lose the white blob in the left as a local extremum because the pixel that is supposed to be the local extremum, indicated as "edged extremum" in red, is located at the edge of the region that is effectively related to the HSC variable and thus fails the condition in Eq. 15. Furthermore, the local extrema could also come from noise peaks, as marked by magenta in Fig. 4h. Despite that the local extrema and the blobs are not in one-to-one correspondence, their numbers are still highly correlated.
Before we compute the percentile A eff,t (p) at p% in Eq. 13, we count the number of the local extrema, N m extrm , for each HSC variable m (where m = 1, ..., M) in an epoch of choice. Depending on the imaging survey, this epoch could be for example the best-seeing epoch, or median-seeing epoch. In the HSC transient survey, the best-seeing epoch with a remarkable value of 0.42 arcsec is so good that the sharp images actually lead to significant artifacts in the deep difference image that affects the number of extrema. We therefore use the median-seeing epoch (with seeing of 0.72 arcsec) in the HSC transient survey for detecting the local extrema, as this is the regime where the difference imaging pipeline runs well. Given a criterion, N crit , a HSC variable will be discarded from the lens classification if N m extrm ≤ N crit , and a HSC variable will be kept in the lens Article number, page 9 of 15 A&A proofs: manuscript no. lens_search_time_variability classification if N m extrm > N crit .
Assuming we have M HSC variables with their number of the local extrema larger than N crit after the preselection (where M ≤ M), we compute their area of the effective region for each epoch t (A m eff,t , m = 1, ..., M ), and further calculate the percentile p% for the area of the effective region from these M HSC variables, A eff,t (p), in each epoch t, which is similar to Eq. 13. We notice that, at a same percentile p%, A eff,t (p) is generally larger than A eff,t (p), the percentile p% from all the HSC variables before the preselection. The reason is that the HSC variables with fewer local extrema are also the HSC variables with fewer spatially extended blobs counting for the extendedness, so when we reject the HSC variables with small number of the local extrema, we also reject the HSC variables that are small in the area of the effective region. Raising A eff,t (p), the percentile p% for the area of the effective region, with the preselection by the number of blobs makes the further lens classification more efficient and helps us to avoid large number of false candidates for lensed quasar, though we would lose a few candidates, particularly narrowly separated lensed quasars (Sec. 5.2).

Performance test
In this section, we examine the performance of our lens search algorithm. To test the performance of the lens search algorithm, we need not only the deep difference images of lensed quasars, but also the deep difference images of the non-lensed objects. We exploit 12910 HSC variables from the HSC transient survey as the non-lensed objects to test our lens search algorithm. Like the 2033 injected mock lensed quasars, we take the deep difference images from all the 13 epochs for each HSC variable.

Classification based only on spatial extent
The lens search algorithm performance is measured by the truepositive rate (TPR) and the false-positive rate (FPR), which are defined as and where N TP is the number of correctly identified positive cases, N P is the number of total positive cases, N FP is the number of falsely identified negative cases, and N N is the number of total negative cases. Here, the positive cases are the quasar lenses, and the negative cases are the non-lensed objects, which are the HSC variables from the HSC transient survey.
We split the injected mock lensed quasars into four subgroups according to their brightness and quasar-image separation, and test the performance of the lens search algorithm individually for each subgroup. In this work, a lensed quasar is bright if the magnitude of the third brightest image m 3rd < 22.0 mag, otherwise this lensed quasar is faint (22.0 mag ≤ m 3rd < 24.0 mag); a lensed quasar has wide separation if the largest separation among the pairs of the lensed images θ LP > 1.5 , otherwise this lensed quasar has narrow separation (0.5 < θ LP ≤ 1.5 ). Therefore, the four subgroups are bright lensed quasars with wide separation, bright lensed quasars with narrow separation, faint lensed quasars with wide separation, and faint lensed quasars with narrow separation.
For each subgroup, N P is the total number of lensed quasars in the subgroup, and N TP is the number of lensed quasars in the subgroup that are classified as lens candidates by the lens search algorithm. For example, in the "bright-wide" group, N P is the total number of bright lensed quasars with wide separation, and N TP is the number of the bright lensed quasars with wide separation that are classified as lens candidates by the lens search algorithm. For all the four subgroups, N N is the total number of the HSC variables that we use to test the lens search algorithm (N N = 12910 in this paper), and N FP is the number of the HSC variables that are falsely classified as the lens candidates by the lens search algorithm. N N and N FP only comprise of the HSC variables and they do not correspond to any specific subgroup. Fig. 5 shows the receiver operating characteristic (ROC) curves, TPR against FPR, for the four subgroups. The result in Fig. 5 is from the method in Sec. 4.1, by varying p thrs % at N thrs = 9. 10 Increasing the value of p thrs % from 5% to 99.5% gives points along the ROC curves from the top-middle to the bottom-left corner. The best performance would give points in the top-left corner, indicating high TPRs and low FPRs. The nearer the curve is to the top-left corner, the better the performance of the lens search algorithm is. Diamonds, triangles, and circles in Fig. 5 mark TPRs and FPRs respectively at p thrs % = 90%, 95%, and 97.5%, for each subgroup of the injected mock lensed quasars. We summarize the values of TPRs and FPRs in Table 5.
We further explore the performance of our lens search algorithm in each subgroup when N thrs = 9. For the bright lensed quasars with wide separation, our lens search algorithm could identify them with a TPR = 90.1% and a FPR = 2.3%, at p thrs % = 95%. Although the ROC curve for the bright lensed quasars with narrow separation is slightly lower, our lens search algorithm could still capture these lensed quasars with a TPR = 84.0% at the same FPR (= 2.3%), when p thrs % = 95%. This indicates that our lens search algorithm with the method in Sec. 4.1 is sensitive to the bright lensed quasars regardless of the separation. Given the much lower ROC curves, the faint lensed quasars are mainly more difficult to detect by our lens search algorithm. For the faint lensed quasars with wide separation, the TPR hugely drops to 56.3% at p thrs % = 95%, while the nearest point to the top-left corner on the ROC curve is (TPR, FPR) = (89.0%, 13.0%) at p thrs % = 82%.
We examine our lens search algorithm's sensitivity to m 3rd with p thrs % = 95% and N thrs = 9 in Fig. 6. The top panels in Fig. 6 show the TPRs, and the bottom panels show the number of mock lenses in each m 3rd bin. Our lens search algorithm could detect all the lensed quasars with the third brightest image brighter than 20.5 mag (m 3rd < 20.5 mag), and the TPRs for the lensed quasars with m 3rd < 21.0 mag are larger than or equal  Table 5: TPRs and FPRs of our lens search algorithm at N thrs = 9. Fig. 5: ROC curves of the lens search method based only on spatial extent. A lens is considered to be bright if its third brightest image has m 3rd < 22.0 mag, and faint if 22.0 mag ≤ m 3rd < 24.0 mag; a lens is considered to have wide separation if the largest separation among the pairs θ LP > 1.5 , and narrow separation if 0.5 < θ LP ≤ 1.5 . The ROC curves are plotted by varying p thrs % from 5% to 99.5% at N thrs = 9 (see Sec. 5.1 for details). The diamonds, triangles, and circles are the points on the ROC curves when p thrs % = 90%, 95%, and 97.5%, respectively. For the bright lensed quasars with wide separation, our lens search algorithm could detect them with (TPR, FPR) = (90.1%, 2.3%). The lens search performance for the bright lensed quasars with narrow separation is similar to the bright lensed quasars with wide separation.
to 90%. We show again that our lens search algorithm with the method from Sec. 4.1 has similar sensitivity to the bright lensed quasars with both wide and narrow separation; the reason is that, most of the HSC variables are unlensed sources and have small areas of the effective region, resulting in small values for the percentile of p thrs %, so a bright lensed quasar typically has an area of the effective region larger than the percentile of p thrs % and would be identified as a lens candidate regardless of its separation. The TPRs gradually drop when m 3rd ≥ 21.0 mag. The reason our lens search algorithm has poorer performance on the faint lensed quasars is that, the difference images of the faint lensed quasars have fewer pixels with values larger than 3σ (Eq. 10), as the brightness change of quasar is generally small, so the area of the effective region in each epoch, A m eff,t , is also small even for the lenses with wide separation, which makes the faint lensed quasars unable to pass the given thresholds, p thrs % Fig. 6: Sensitivity to the quasar image brightness m 3rd of the lens search method based only on spatial extent. Wide: the group of lensed quasars with the largest separation among the pairs θ LP > 1.5 . Narrow: the group of lensed quasars with 0.5 < θ LP ≤ 1.5 . The top panels show the TPRs of each m 3rd bin, when p thrs % = 95% and N thrs = 9. The bottom panels show the number of mock lenses in each m 3rd bin. Our lens search algorithm could identify more than 90% of the lensed quasars with m 3rd < 21.0 mag. and N thrs .

Classification using both spatial extent and number of blobs
Now we explore the lens search performance with the preselection based on the number of blobs. We test how the ROC curves change with the criterion number of local extrema, N crit (see Sec. 4.2 for details). In Fig. 7, we plot the ROC curves by varying p thrs % for N crit = 2. As shown in Fig. 7, after the preselection is applied (N crit = 2), our lens search algorithm could capture the bright lensed quasars with wide separation at (TPR, FPR) = (97.6%, 2.6%) with thresholds, p thrs % = 55% and N thrs = 4. 11 Comparing to the lens search performance without the preselection by the number of blobs (N crit = 0), the lens search algorithm applying the preselection (N crit > 0) could have similar performance with much looser constraints on p thrs % and N thrs . Crosses, triangles, and circles in Fig. 7 represent TPRs and FPRs at p thrs % = 0%, 55%, and 75%, respectively. With the preselection by number of blobs, the ROC curves start from a much lower FPR, compared to the ROC curves without the preselection (N crit = 0). Actually, the crosses in Fig. 7 indicate that when N crit = 2, we have a much lower FPR (∼5%) at N thrs = 4, even when there is no constraint from p thrs %. We summarize the values of TPRs and FPRs from the lens search performance with the preselection by the number of blobs in Table 6.
We further look into each subgroup of lensed quasars for the comparison of the performance between the two conditions, with and without the preselection by number of blobs (N crit > 0 and N crit = 0). At small values of p thrs %, the TPRs are ∼100% for all the four subgroups when N crit = 0, while the TPRs for the four subgroups all drop when N crit = 2. This indicates that, although the preselection based on number of blobs hugely decrease the FPRs, the preselection also discard a part of the candidates for lensed quasars. The TPR decreases by ∼2.2% for the bright lensed quasars with wide separation, by ∼25.1% for the bright lensed quasars with narrow separation, and by ∼13.9% for the faint lensed quasars with wide separation. We notice that, unlike Fig. 5, where the ROC curve of the bright lensed quasars with narrow separation is higher than the ROC curve of the faint lensed quasars with wide separation, when N crit = 2, the ROC curve of the faint lensed quasars with wide separation becomes higher than the ROC curve of the bright lensed quasars with narrow separation (at FPR 0.02), which implies the blob-based preselection has a more significant impact on the lensed quasars with narrow separation.
To investigate the influence of the preselection, we examine the sensitivity to θ LP of the lens search algorithm setting p thrs % = 0%, at N thrs = 4, with a preselection N crit = 2. By setting p thrs % = 0%, we could study the impact from the preselection more clearly. Fig. 8 shows that, when there is no constraint from p thrs , our lens search algorithm with the preselection by number of blobs could detect the lensed quasars with θ LP > 1.5 at relatively stable TPRs (> 80%), and there is a great drop in TPRs when θ LP becomes smaller than 1.0 . The lensed quasars with narrow separation become even harder to capture when the preselection is applied, and the reason is that, as lensed images with narrow separation are blended together, the number of blobs decreases, making the lensed quasars with narrow separation difficult to meet the given criterion, N crit . We show the ROC curves for the bright lensed quasars with wide separation at different values of N crit (when N thrs = 4) in Fig. 9, zoomed-in with the x-axis spanning from 0 to 0.05. Crosses and triangles in Fig. 9 represent p thrs % = 0% and 55% respectively. As N crit increases, the ROC curve drops and slightly shifts to the left. Both the TPRs and the FPRs decrease when we raise N crit . However, FPR decreases faster than TPR at low values of N crit . In particular, we could decrease the FPR by a factor ∼2 and have a more efficient lens search when we raise N crit from 2 to 3, at the expense of lowering TPR by only ∼5%. Raising N crit beyond 3 would start to lead to more decrease in TPR compared to FPR, and not as advantageous.

Classification based on a single epoch
It is ideal and the most efficient to detect lensed quasars through only one epoch (in addition to the reference image). If the lens Fig. 7: ROC curves of the lens search method using both spatial extent and the number of blobs. The subgroups of lensed quasars are the same as the subgroups in Fig. 5. The ROC curves are plotted by varying p thrs % from 0% to 99.5% at N thrs = 4 with the preselection on the number of blobs N crit = 2. The crosses, triangles, and circles are the points on the ROC curves when p thrs % = 0%, 55%, and 75%, respectively. With the preselection N crit = 2, our lens search algorithm could identify the bright lensed quasars with wide separation at (TPR, FPR) = (97.6%, 2.6%), which is similar to the lens search performance when no preselection is performed. The ROC curves also start from a much lower FPR (crosses) with the preselection N crit = 2. The ROC curves are zoomed-in with the FPR (x-axis) spanning from 0 to 0.05 in the small panel.
search algorithm is able to identify the candidates for lensed quasars in a single epoch of an ongoing cadenced survey, possible spectroscopic follow-up could be conducted immediately for the confirmation and the lens model, and a further monitoring observation for measuring the time delays could start right away. Here, we study the performance of our lens search algorithm to capture lensed quasars in one single epoch.
If only one epoch is available, our lens search algorithm will employ only p thrs . Since the seeing of a single epoch might not be suitable for counting the local extrema, we do not apply the preselection using number of blobs. For one single epoch t , our lens search algorithm will classify a HSC variable m as a candidate for lensed quasar if where A m eff,t is the area of the effective region in epoch t for the HSC variable m , and A eff,t (p thrs ) is the percentile p thrs % for the area of the effective region from all the HSC variables in epoch t . We test the lens search performance in each epoch from the HSC transient survey, and show the ROC curves with varied values of p thrs % from three epochs with different seeings in Fig. 10. The result in the left panel of Fig. 10 Table 6: The values of TPRs and FPRs at N thrs = 4, when the preselection based on the number of blobs, N crit = 2, is applied. epochs, where the nearest point to the top-left corner on the ROC curve for the bright lensed quasars with wide separation is (TPR, FPR) = (97.6%, 2.4%) at p thrs % = 97.5%, which is similar to the performance when all the 13 epochs is used. The middle panel of Fig. 10 shows the average case of the lens search performance we have in one single epoch, with the nearest point to the top-left corner on the ROC curves, (TPR, FPR) = (94.2%, 5.0%) at p thrs % = 95.0%. In the right panel of Fig. 10 is the worst performance we have among the 13 epochs.
We further investigate the relation between the lens search performance and the seeing in one single epoch. For each epoch t, we define the distance d t as the distance between (TPR, FPR) = (100%, 0%) and the nearest point on the ROC curve of the bright lensed quasars with wide separation. The smaller d t is, the better the lens search performance is. We plot d t against the seeing for each epoch t in Fig. 11. d t becomes small when seeing is around 0.7 arcsec, indicating that our lens search algorithm would have good performance when the single epoch has seeing around 0.7 arcsec. In Fig. 11, we also see that d t suddenly increases when seeing is better than 0.5 arcsec, which verifies that when the see- ing is too good, our lens search performance will be heavily affected by the significant artifacts from the sharp images.

Conclusion
In this work, we present a comprehensive simulation pipeline of time-varied lensed images and develop a new algorithm for searching lensed quasars through their time variability. The simulation pipeline here is useful not only for lensed quasar search, but also for other studies, such as studying lensed supernovae and testing lens mass modelling. Our lens search method builds upon the one first proposed by Kochanek et al. (2006), and provides a practical way of selecting lensed quasar candidates through their difference images. We summarize the main results as follows: -Our simulation pipeline generates images of lensed quasar accounting for the quasar variability, quasar host, lens galaxy, and the PSF variation. The application of the simulation pipeline to the HSC transient survey yields HSC-like difference images of the mock lensed quasars. -The lens search algorithm we develop in this work quantifies spatial extent of the variable objects in difference images, and further selects the variable objects with large spatial extent as lensed quasar candidates. The test on a sample containing mock lensed quasars and HSC variables shows that our lens search algorithm could identify the bright lensed quasars with wide separation (m 3rd < 22.0 mag and θ LP > 1.5 ) at a high TPR (90.1%) and a low FPR (2.3%). In each panel, the ROC curves are plotted by varying p thrs %, and the diamonds, triangles, and circles are the points on the ROC curves at p thrs % = 90%, 95%, and 97.5%. The seeing of each epoch is also indicated in each panel. Our lens search algorithm has the best performance in the single epoch with a seeing of 0.72 , capturing the bright lensed quasars with wide separation at (TPR, FPR) = (97.6%, 2.4%) when p thrs % = 97.5%. In a single epoch with average seeing (∼0.98 ), our lens search algorithm could detect the bright lensed quasars with wide separation at (TPR, FPR) = (94.2%, 5.0%) when p thrs % = 95.0%. The worst lens search performance happens in the epoch with an exceptional seeing value of 0.48 due to artifacts appearing in the difference image pipeline when the seeing is "too" good. Fig. 11: The relation between the lens search performance and the seeing in each single epoch from the HSC transient survey. d t is defined as the minimum distance between (TPR, FPR) = (100%, 0%) and the ROC curve of the bright lensed quasars (m 3rd < 22.0 mag) with wide separation (θ LP > 1.5 ) for each epoch t. The smaller d t is, the better the lens search performance is. The 13 star symbols indicate the d t values for the 13 epochs, plotted against their corresponding seeing. Our lens search algorithm performs well when the single epoch has seeing ∼0.7 . The low lens search performance in the epoch with extraordinary seeing (< 0.5 ) is due to the significant artifacts in the difference image from the sharp image.
-With a preselection on number of blobs, our lens search algorithm could achieve an even higher TPR (97.6%) without a significant change in FPR (2.6%) for the bright lensed quasars with wide separation. The preselection is more sensitive to the lensed quasars with wide separation.
-Although our lens search algorithm mainly uses difference images from multiple epochs, it also works with only one single epoch. The lens search performance in one single epoch depends on the seeing. When the seeing is around 0.7 arcsec, we have the best lens search performance, (TPR, FPR) = (97.6%, 2.4%) for the bright lensed quasars with wide separation. If the seeing is substantially better or worse than 0.7 arcsec, our lens search performance become poorer.
Our lens search algorithm will be even more powerful when combined with other lens search techniques. While our lens search algorithm achieves a high TPR and a low FPR, the number of the input "variables" (e.g. the HSC variables in this work) are typically huge, and thus the absolute number of false-positive lens candidates are also large even with FPR of 1 − 2%, making further confirmation inaccessible. Although the lens search techniques in previous works use static approaches without any information from quasar variability, their exploration in catalogues and image configuration of lensed quasar are helpful for making our lens search algorithm more efficient. For example, we could use the techniques in the cuts of magnitude and color (e.g. Agnello et al. 2015;Williams et al. 2017) or Gaia multiplets to first select objects that are possible to be lensed quasar, and then we apply our lens search algorithm to examine their variable nature. Alternatively, after selecting the lensed quasar candidates by our lens search algorithm, we could apply the technique in Chan et al. (2015) to filter the candidates for lensed quasar. The variability across multiple bands could also be very useful for enhancing our lens search algorithm, while so far our lens search algorithm has exploited only i-band. We will explore the multi-band variability in the future work.
Since the lens search algorithm we have now is built on the HSC transient survey, we could start finding lensed quasars in the HSC transient survey right away (Chao et al. in prep.). More-Article number, page 14 of 15 over, by adapting our simulation pipeline for the time-varied lensed images, we could also apply the lens search algorithm to other cadenced surveys, such as the upcoming LSST. LSST is expected to have similar image quality as the HSC transient survey, but with a much larger area coverage. Around 2000 lensed quasars with four multiple images are predicted to be found in the LSST (Oguri & Marshall 2010). We expect our lens search algorithm to skillfully capture the new lensed quasars in the ongoing and upcoming cadenced surveys.