© D. C.-Y. Chao et al. 2020

1. Introduction

Gravitationally lensed quasars are powerful tools to study astrophysics and cosmology. For example, they can be used to probe the 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 galaxies (e.g., Peng et al. 2006; Ding et al. 2017). With the time delays between multiple images, we can also measure the Hubble constant, H₀, which is crucial for testing cosmological models and studying dark energy (e.g., Refsdal 1964; Birrer et al. 2019; Wong et al. 2020; Chen et al. 2019).

Given the importance of lensed quasars, there have been several systematic searches of these objects. 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 a compact flat spectrum and then looking for the high-resolution, multiply-imaged components in the radio follow-up. In the optical, the Sloan Digital Sky Survey (SDSS; York et al. 2000) Quasar Lens Search (SQLS; Oguri et al. 2006; Inada et al. 2008, 2010, 2012; More et al. 2016a), which has found the largest sample of lensed quasars (∼60), started 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; next, these authors conducted follow-up observations to confirm these candidates after a visual inspection of the objects selected by the morphological or color selection. Both CLASS and SQLS start with spectroscopic searches, and further look into the morphology or images for lensed quasar candidates. 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, owing to the increased depth and much larger areal coverage, many ongoing and upcoming wide-field surveys provide a great pool of new lensed quasars and allow lens searches directly with images. With the aim of exploiting the unprecedented advantage of those surveys, many dedicated searches have been initiated with new lens search techniques. Some lens search techniques explore multiband catalogs and employ various cuts in 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. 2019). 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) shows that lensed quasars can also be found through citizen science.

With the exceptional resolution of Gaia (Gaia Collaboration 2016), it is possible to conduct quasar lens search by looking for multiple detections in Gaia or comparing the flux and position offsets from other surveys for objects in SDSS, Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chambers et al. 2016) (e.g., Lemon et al. 2017, 2018, 2019; Ostrovski et al. 2018), Dark Energy Survey (DES; Sánchez 2010) (e.g., Agnello & Spiniello 2019; Agnello et al. 2018a), or Kilo-Degree Survey (KiDS; de Jong et al. 2013) (Spiniello et al. 2018). Using the ability of Gaia to carry out multiple detections to search for lensed quasars could also be applied to the ATLAS¹ footprint (Agnello et al. 2018b). We can also combine Gaia detections with other lens searches to find lensed quasars. While the search for lensed quasars with small separations benefits most from Gaia’s exceptional resolution, Gaia usually can detect only one or two multiple image(s) of a lens system, which is limited by its shallow depth, i ≲ 20.7.

In the near future, thousands of lensed quasars are expected to be detected by the Rubin Observatory Legacy Survey of Space and Time (LSST; Ivezić et al. 2019) (Oguri & Marshall 2010, hereafter OM10). With further imaging and spectroscopic analysis, we are able to discover and confirm these LSST detections, enlarging the sample size of lensed quasars by at least an order of magnitude. Moreover, the LSST is a cadenced survey that will reveal the time variability of variable objects, such as lensed quasars with significant variability². 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 method based on time variability benefits lens searches given the enormous number of objects in the LSST.

Inspired by Kochanek et al. (2006), we developed a lens search algorithm for lensed quasars based on their time variability. The Hyper Suprime-Cam (HSC; Miyazaki et al. 2012; Aihara et al. 2018) transient survey (Yasuda et al. 2019) provides a great opportunity to simulate the difference images that could be used 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. 2018) with similar image quality expected for the 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 these objects as lensed 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). In spite of the higher number of doubles than quads, quads provide more constraints on the lens potential distribution and stellar mass fraction. Furthermore, the configuration of quads is so unique that it is unlikely to be mistakenly identified, while many objects such as quasar binaries are easily mistaken as doubles³.

The organization of the paper is as follows. In Sect. 2, we present a new simulation pipeline of time-varied lensed images, and apply the new simulation pipeline to the HSC transient survey in Sect. 3. The lens search algorithm is detailed in Sect. 4. The lens search performance is shown in Sect. 5 and we conclude in Sect. 6.

2. 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 a lens search algorithm based on time variability. We used the mock configuration from the catalog created by OM10 to generate the simulated lensed quasars. In order to have realistic simulations, we simulated 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. The flow chart of this simulation pipeline is shown in Fig. 1.

Fig. 1. Flow chart illustrating the simulation of time-varying lensed quasars images. For quasars, mock quasar light curves are generated based on their source redshifts zs, and being shifted in time and magnified the mock quasar light curves on the image plane (see Fig. 2). For quasar hosts, their Sérsic light profiles are created on the source plane and their lensed images are produced on the image plane assuming a SIE model. Lens galaxies are selected from a real survey based on their lens redshifts zd 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 all the components for one lens system (image of lensed quasar, image of lensed quasar host, and image of lens galaxy) are produced, these components are convolved with a realistic PSF and are added to obtain one complete image.

2.1. Mock lenses from OM10

Oguri & Marshall (2010; OM10) predicted the distribution of lensed quasars and produced a mock catalog of lensed quasars. For each lens system, OM10 provides the source position (η₁, η₂), source redshift z_s, unlensed i-band magnitude m_s of the quasar, lens redshift z_d, velocity dispersion σ, ellipticity e, and position angle (PA) θ_e of the lens galaxy. Also, OM10 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

$$\begin{array}{c}\hfill \kappa ({\theta }_1,{\theta }_2)=\frac{{\theta }_{\mathrm{Ein}}\sqrt{1-e}}{2}\frac{\lambda (e)}{\sqrt{{\theta }_1^2+{(1-e)}^2{\theta }_2^2}},\end{array}$$(1)

where (θ₁, θ₂) 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

$$\begin{array}{c}\hfill \varphi ({\theta }_1,{\theta }_2)=\frac{\gamma }{2}({\theta }_1^2-{\theta }_2^2)cos2{\theta }_{\gamma }+\gamma {\theta }_1{\theta }_{2}sin2{\theta }_{\gamma },\end{array}$$(2)

where γ is the magnitude of the external shear and θ_γ is the orientation of the external shear.

2.2. Quasar

For the point sources, that is, the lensed quasars without the (lensed) host galaxies, we used the source redshift z_s, unlensed i-band magnitude m_s, lensed image positions, magnifications, and time delays provided by OM10 to generate the simulated image for each lensed quasar.

We generated the source light curve of quasars via CARMA_PACK (Kelly et al. 2009, 2014) based on z_s and m_s from OM10. The CARMA_PACK package was developed to quantify stochastic variability, especially for quasar variability, and the quasar light curves are generated with a damped-random-walk model based on magnitude m and the fitting relation on redshifts z in Kelly et al. (2009). We then shifted (in time) and magnified the generated source light curve according to the time delay and magnification of each lensed image⁴.

Figure 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). 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 generated 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.

Fig. 2. Illustration of generating light curves for a lensed quasar without its (lensed) host galaxy. Left: positions of the four multiple images (A, B, C, and D) and the lens galaxy (cross). Top right: mock quasar light curve. Bottom right: the lensed light curves of the four multiple images in the left panel, for the period marked by the blue band in the top right panel.

2.3. Quasar host galaxy

Since OM10 only provides the configuration of point sources and there is no definite morphological relation between quasar and the host galaxy with little scatter, we randomly drew the morphological properties based on the Sérsic profile for the host galaxy of quasar and used GLAFIC (Oguri 2010) to generate the lensed image of host galaxy. We conservatively adopted broad ranges for the parameters in the Sérsic profile of the host galaxy (based on, e.g., Blanton & Moustakas 2009; Peng et al. 2006; 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 adopted 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 are consistent with the lensed quasar image positions.

2.4. Lens galaxy

For the lens galaxies, we assigned the light profile to the lens galaxy for each lens system with fitted profile of galaxy from a real survey because only the mass model of lens galaxy is given by OM10. We denoted such a real survey “S_imag”, where S_imag can be any survey devoted to finding lens systems, such as the SDSS or the HSC survey.

We began with a catalog 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, we do not have to cross-match S_imag with another survey. From the cross-matched catalog, we selected “matched galaxies” for each lens galaxy with the closest values to the lens redshift z_d, OM10 and the velocity dispersion σ_OM10 from the mass model provided by OM10. We first selected the galaxies with

$$\begin{array}{c}\hfill |d_z|<0.01\end{array}$$(3)

and

$$\begin{array}{c}\hfill |d_{\sigma }|<\mathrm{\Delta }\sigma ,\end{array}$$(4)

where d_z = z_spec − z_d, OM10, d_σ = σ_spec − σ_OM10, and Δσ is the error in σ_spec from S_spec. We then got the matched galaxies from the selected galaxies with the smallest values of d_z, σ, where d_z, σ is defined as

$$\begin{array}{c}\hfill d_{z,\sigma }=\sqrt{d_z^2+d_{\sigma }^2}.\end{array}$$(5)

Ideally, for one lens galaxy in OM10, we could have ten galaxies with the smallest d_z, σ from the cross-matched catalog as matched galaxies 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 matched galaxies. For such cases, we took all the available galaxies passing the selection criteria (Eqs. (3) and (4)) as matched galaxies. In this work, we used the data from the SDSS Data Release 14 (DR14; Abolfathi et al. 2018) as S_spec.

The velocity dispersion σ_SDSS is limited in the cross-matched catalog when z_SDSS >  1.04, so we applied only Eq. (3) to select matched galaxies for the lens galaxy with 1.04 <  z_d, OM10 <  1.57; we selected up to ten galaxies as the matched galaxies with the smallest values of d_z. Owing to the lack of objects with z_SDSS >  1.57 in the cross-matched catalog, we did not simulate the lens systems with z_d, OM10 >  1.57.

We then obtained the image for each matched galaxy from S_imag, and fit the light distribution of the matched galaxy with the Sérsic profile. For the fitting, we define

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _i^{N_{\mathrm{p}}}\frac{{(I_{i,\mathrm{fit}}-I_{i,\mathrm{data}})}^2}{{\sigma }_{i,\mathrm{data}}^2},}\end{array}$$(6)

and the reduced χ² as

$$\begin{array}{c}\hfill {\chi }_{\mathrm{reduced}}^2=\frac{{\chi }^2}{N_{\mathrm{p}}},\end{array}$$(7)

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 ranked the matched galaxies by ${\chi }_{\text{reduced}}^2$. The smaller ${\chi }_{\text{reduced}}^2$, the higher is the rank.

After we had the fitted light profile for each matched galaxy, we compared 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 matched galaxy, if |e_fit − e_OM10|< 0.2, we picked the matched galaxy as the lens galaxy; otherwise, we checked e_fit of the matched galaxy in the next lower rank. For the lens galaxy that has no matched galaxy with the satisfied value of e_fit, we picked the matched galaxy with the smallest value of ${\chi }_{\text{reduced}}^2$. In total, we simulated 2033 quads with lens redshift z_d, OM10 <  1.57 from OM10⁵.

2.5. Simulated image with all components

Once we had the individual simulated parts of a whole lens system (the lensed quasar, lensed host galaxy, and lens galaxy), we convolved these individual parts with the point spread functions (PSFs) in S_imag for each epoch and added the parts to become one complete image. The complete image in an epoch t is expressed as

$$\begin{array}{c}\hfill I_t(x,y)=(I_{\mathrm{q},t}(x,y)+I_{\mathrm{h}}(x,y)+I_{\mathrm{d}}(x,y))\otimes {\mathrm{PSF}}_t(x,y),\end{array}$$(8)

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 generated complete images for all the epochs in S_imag, we obtained the images of the mock lensed quasar in the time series of S_imag.

3. 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.

3.1. HSC transient survey

The HSC transient survey was done in the COSMOS (Scoville et al. 2007) field as part of the HSC-SSP (Subaru Strategic Program; Aihara et al. 2018; Miyazaki et al. 2018; Komiyama et al. 2018; Kawanomoto et al. 2018; Furusawa et al. 2018) observation from November 2016 to April 2017 with a pixel size of 0.168 arcsec, covering 1.77 deg² in the Ultra-Deep layer to a median depth per epoch of 25.9 mag in i band. Since OM10 has only i-band magnitude, we simulated lenses only in i band. The HSC transient survey has 13 epochs (nights of observations) over six months for i band (Table 3).

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 co-adding multiple exposures from multiple epochs during March 2014 to April 2016 (Table 4). Only the exposures with seeing better than 0.7 arcsec were used to create 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 co-adding several warped images, and each warped image corresponds to a distortion-corrected image of the sky from a single exposure. For each epoch, the difference imaging was performed on every warped image by subtracting the deep reference image to produce the warped difference images, and all the warped difference images were then co-added to create the deep difference image. The HSC transient survey convolves the exposures used in creating the deep reference image (Table 4) with a kernel to match the PSF of each warped image when performing the difference imaging process.

To generate realistic deep difference image of the mock lenses, we selected the locations that have no object within 5 arcsec, which we call “empty regions”, and we injected 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 was done, we used the difference imaging method mentioned above to create the deep difference image of the mock lenses. Practically, we selected 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 were used to create the deep reference image (Table 4), and 175 warped images were used to create the deep difference images for all the 13 epochs in i band.

3.2. Convolution with the HSC PSFs and lens injection

In order to obtain realistic deep difference image, we injected 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 image 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 co-added, 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 simulated the images for the whole lens system as mentioned in Sect. 2, accounting for the quasar variability in the epoch, and then we created 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 injected the complete images into the corresponding warped images.

For the lens injection into the deep reference image, we first weighted 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 convolved the weighted image with the PSF from the deep reference image and further injected 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), can be described as

$$\begin{array}{c}\hfill I_{\mathrm{ref}}(x,y)=({\displaystyle \sum _t\frac{I_{\mathrm{q},t}(x,y)·N_t}{N_{\mathrm{total}}}+I_{\mathrm{h}}(x,y)+I_{\mathrm{d}}(x,y))\otimes {\mathrm{PSF}}_{\mathrm{ref}}(x,y),}\end{array}$$(9)

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 simulated 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 did 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 Sect. 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 applied the deep difference imaging method to obtain the deep difference image for the mock quads (bottom row in Fig. 3). As mentioned in Sect. 3.1, we first get the warped difference images of injected lenses as intermediate products, and the deep difference images are created by co-adding 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 can see the light from both the lensed quasar and lens galaxy in each panel, which appear as “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 can 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 because the pipeline might not perform perfect subtraction for all the non-variable components.

Fig. 3. Examples of lenses injected into the HSC transient survey. Top: single-epoch images co-added by the pipeline of the HSC transient survey. Bottom: corresponding difference images of the top row produced by the image subtraction in the pipeline of the HSC transient survey. The size of each image cutout is 10″ × 10″.

4. Search method: Spatial extent in difference image

With the multiple-image feature and the variable brightness, strongly lensed quasars are expected to exhibit multiple point-like image residuals in the difference image shown in Fig. 3. However, lensed quasars are often not deblended by the survey pipelines, and most of these show extended morphology. Since most astrophysical variable sources are isolated and point-like (e.g., variable stars, unlensed quasars, unlensed supernovae), targeting spatially extended or multiple point-like variable sources in the difference image is an effective approach to find lensed quasars, as previously noted by Kochanek et al. (2006). We developed a lens search algorithm for lensed quasar via a method that quantifies the extendedness of an object in the difference image and selects objects that are large in spatial extent as lens candidates. We present the method that is based on the extendedness in Sect. 4.1, and an enhancement of the method with a secondary criterion in Sect. 4.2.

4.1. 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 collected its deep difference images from all the 13 epochs (Table 3) and created 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

$$\begin{array}{c}\hfill I_{\mathrm{mask}}(i,j)=\{\begin{array}{cc}1,\hfill & \mathrm{if}|I(i,j)|>3\sigma (i,j)\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\end{array}$$(10)

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. Figures 4a and b 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

$$I_{\text{eff}}(i,j)=\{\begin{array}{ll}1,\hfill & \text{if}{\displaystyle \sum _{i^{\prime }=i-1}^{i+1}{\displaystyle \sum _{j^{\prime }=j-1}^{j+1}{I}_{\text{mask}}}}(i^{\prime },j^{\prime })>2\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}$$(11)

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 denoted by the white pixels show the spatial extent of the objects, after removing noise peaks in the 3σ masks. In the local extremum detection, pixels with values larger than all the neighboring pixels or smaller than all the neighboring pixels in the difference images within the effective regions are selected; the local extrema, indicated in cyan, correspond to the “blobs” (see Sect. 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 as a 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″.

where I_eff(i, j) is the pixel value of (i, j) in the effective region. Figure 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,

$$\begin{array}{c}\hfill A_{\mathrm{eff}}={\displaystyle \sum _{i=1}^{N_x}\sum _{j=1}^{N_y}{I}_{\mathrm{eff}}(i,j),}\end{array}$$(12)

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 can discard the effect from the noise peaks in the outskirts and quantify the extendedness of a 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_{\text{eff},t}^m$ 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):

$$\begin{array}{c}\hfill p\%\mathrm{of}M\mathrm{HSC}\mathrm{variables}\mathrm{with}{A}_{\mathrm{eff},t}^m<A_{\mathrm{eff},t}(p).\end{array}$$(13)

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_{\text{epoch}}^m(p)$. For example, if a HSC variable has $A_{\mathrm{eff},t}^{m\prime }>A_{\mathrm{eff},t}(p_{\mathrm{thrs}})$ at a threshold of percentile p_thrs%, for t = t₁, t₂,  and t₃, we say $N_{\mathrm{epoch}}^{m\prime }(p_{\mathrm{thrs}})=3$. Given a threshold of percentile p_thrs%, we set another threshold, N_thrs, and a HSC variable is classified as a candidate for lensed quasar if $N_{\mathrm{epoch}}^{m\prime }(p_{\mathrm{thrs}})>N_{\mathrm{thrs}}$.

We also perform the calculation in Eqs. (10)–(12) for each injected mock lensed quasar. Figure 4e shows an example of injected mock lensed quasar in the same epoch as Fig. 4a. Figures 4f and g are respectively the 3σ mask and the effective region of the mock lensed quasar in Fig. 4e. As shown in Figs. 4c and g, the injected mock lensed quasar has a much larger area of the effective region.

4.2. 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 either have 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”, as shown in Figs. 4a and e. In the difference image, these blobs are like local extrema in the zero-background. Owing 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 Sect. 4.1 with a preselection of HSC variables based on a criterion on the number of local extrema.

We first define a local extremum for the difference image. A pixel (i, j) is a local extremum, if

$$\begin{array}{c}\hfill I_{\mathrm{mask}}(i,j)=1\end{array}$$(14)

and

$$\begin{array}{c}\hfill I_{\mathrm{eff}}(i,j)=1,\end{array}$$(15)

with its neighboring pixels (i′,j′) satisfying

$$\begin{array}{c}\hfill I_{\mathrm{eff}}(i^{\prime },j^{\prime })=1\end{array}$$(16)

and

$$\begin{array}{c}\hfill I(i,j)>I(i^{\prime },j^{\prime })\mathrm{or}I(i,j)<I(i^{\prime },j^{\prime }),\end{array}$$(17)

for i′=i − 1, i, i + 1 and j′=j − 1, j, j + 1, except (i′,j′) = (i, j). The local extrema are indicated in cyan in Figs. 4d and h. Figures 4d and h only keep the pixel values 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 indicated in magenta in Fig. 4h. Even though 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 local extrema, $N_{\text{extrm}}^m$, 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 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 because this is the regime in which the difference imaging pipeline runs well. Given a criterion, N_crit, a HSC variable is discarded from the lens classification if $N_{\mathrm{extrm}}^{m\prime }\le N_{\mathrm{crit}}$, and a HSC variable is kept in the lens classification if $N_{\mathrm{extrm}}^{m\prime }>N_{\mathrm{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_{\mathrm{eff},t}^{m\prime }$, m′ = 1, …, M′), and further calculate the percentile p% for the area of the effective region from these M′ HSC variables, $A_{\text{eff},t}^{\prime }(p)$, in each epoch t, which is similar to Eq. (13). We notice that, at a same percentile p%, $A_{\text{eff},t}^{\prime }(p)$ is generally larger than A_eff, t(p), that is, the percentile p% from all the HSC variables before the preselection. This is because the HSC variables with fewer local extrema are also the HSC variables with fewer spatially extended blobs counting for the extendedness. Thus when we reject the HSC variables with a 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, although we would lose a few candidates, particularly narrowly separated lensed quasars (Sect. 5.2). We illustrate our lensed quasar candidate selection procedure (Sects. 4.1 and 4.2) in Fig. 5.

Fig. 5. Flow chart illustrating the lensed quasar candidate selection procedure. The effective region (Eqs. (10)–(12)) and number of blobs (Eqs. (14)–(17)) are calculated for each variable object from the HSC transient survey (HSC variable). The HSC variables with the number of blobs larger than Ncrit are preselected. The percentile of the effective region for each epoch from all the HSC variables that remain after the preselection on the number of blobs is calculated. Finally, two thresholds, pthrs% and Nthrs, are set; the HSC variables are selected that have an effective region larger than pthrs% percentile for more than Nthrs epochs, as lensed quasar candidates.

5. 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. Based on the estimation of the lensed quasar population in OM10, we expect that at most 1–2 quad(s) lie in the HSC transient survey according to the survey depth and the covering area. The 12910 HSC variables we use to test our lens search algorithm come from a ∼0.64 deg² sub-sky area in the HSC transient survey (1.77 deg²), therefore all of these HSC variables are likely to be false positives. Like the 2033 injected mock lensed quasars, we take the deep difference images from all the 13 epochs for each HSC variable.

5.1. Classification based only on spatial extent

The lens search algorithm performance is measured by the true-positive rate (TPR) and the false-positive rate (FPR), which are defined as

$$\begin{array}{c}\hfill \mathrm{TPR}=\frac{N_{\mathrm{TP}}}{N_{\mathrm{P}}}\end{array}$$(18)

and

$$\begin{array}{c}\hfill \mathrm{FPR}=\frac{N_{\mathrm{FP}}}{N_{\mathrm{N}}},\end{array}$$(19)

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. In this case, 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. The quantities N_N and N_FP only comprise the HSC variables and they do not correspond to any specific subgroup.

Figure 6 shows the receiver operating characteristic (ROC) curves, TPR against FPR, for the four subgroups. The ROC curves quantitatively show the lens search performance, and the best lens search performance would give points in the top left corner, indicating high TPRs and low FPRs. We can further quantify the lens search performance using the minimum distance between the top left corner ((TPR, FPR) = (100%,0%)) and the ROC curve, d. The smaller the value of d, the better the performance of the lens search algorithm is. The result in Fig. 6 is from the method in Sect. 4.1, by varying p_thrs% at N_thrs = 9⁶. 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 diamonds, triangles, and circles in Fig. 6 indicate TPRs and FPRs at p_thrs% = 90%,95%, and 97.5%, respectively, for each subgroup of the injected mock lensed quasars. We list the values of TPRs, FPRs, and d in Table 5, and the numbers of pixels for A_eff(p_thrs) at p_thrs% = 90%,95%, and 97.5% in Table 6.

Fig. 6. Lens search method ROC curves based only on spatial extent. A lens is considered to be bright if its third brightest image has m3rd <  22.0 mag, and faint if 22.0 mag ≤ m3rd <  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 pthrs% from 5% to 99.5% at Nthrs = 9 (see Sect. 5.1 for details). The diamonds, triangles, and circles indicate the points on the ROC curves when pthrs% = 90%,95%, and 97.5%, respectively. The lens search algorithm in this work could detect the bright lensed quasars with wide separation 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.

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 these objects 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 can 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 Sect. 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%. The false positives that our lens search method can possibly detect are very bright objects, such as very bright galaxies or very bright variable stars. Those objects are so bright that the transient data pipeline might not be able to perform the image subtraction properly, and their difference images might not represent the correct brightness change.

We examine our lens search algorithm’s sensitivity to m_3rd with p_thrs% = 95% and N_thrs = 9 in Fig. 7. The top panels in Fig. 7 show the TPRs, and the bottom panels show the number of mock lenses in each m_3rd bin. Our lens search algorithm can 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 to 90%. We show again that our lens search algorithm with the method from Sect. 4.1 has similar sensitivity to the bright lensed quasars with both wide and narrow separation; this is because 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%. Therefore, a bright lensed quasar typically has an area of the effective region larger than the percentile of p_thrs% and can be identified as a lens candidate regardless of its separation. The TPRs gradually drop when m_3rd ≥ 21.0 mag. Our lens search algorithm has poorer performance on the faint lensed quasars because the difference images of the faint lensed quasars have fewer pixels with values larger than 3σ (Eq. 10). Because the brightness change of quasar is generally small, the area of the effective region in each epoch, $A_{\text{eff},t}^m$ (Eqs. (12) and (13), see Sect. 4.1 for more detail), is also small even for the lenses with wide separation, which makes the faint lensed quasars unable to pass the given thresholds, p_thrs% and N_thrs.

Fig. 7. Sensitivity to the quasar image brightness m3rd 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 m3rd bin, when pthrs% = 95% and Nthrs = 9. The bottom panels show the number of mock lenses in each m3rd bin. Our lens search algorithm can identify more than 90% of the lensed quasars with m3rd <  21.0 mag.

5.2. 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 Sect. 4.2 for details). In Fig. 8, we plot the ROC curves by varying p_thrs% for N_crit = 2. As shown in Fig. 8, after the preselection is applied (N_crit = 2), our lens search algorithm can capture the bright lensed quasars with wide separation at (TPR, FPR) = (97.6%,2.6%) with thresholds, p_thrs% = 55% and N_thrs = 4.⁷ 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) gives a similar performance with much looser constraints on p_thrs% and N_thrs. The crosses, triangles, and circles in Fig. 8 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). The crosses in Fig. 8 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 list the values of TPRs, FPRs, and d (the minimum distance between the ROC curves and the top left corner) from the lens search performance with the preselection by the number of blobs in Table 7. The numbers of pixels for A_eff(p_thrs) at p_thrs% = 55% and 75% after the preselection are listed in Table 8.

Fig. 8. ROC curves of the lens search method using both spatial extent and number of blobs. The subgroups of lensed quasars are the same as the subgroups in Fig. 6. The ROC curves are plotted by varying pthrs% from 0% to 99.5% at Nthrs = 4 with the preselection on the number of blobs Ncrit = 2. The crosses, triangles, and circles indicate the points on the ROC curves when pthrs% = 0%,55%, and 75%, respectively. With the preselection Ncrit = 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 Ncrit = 2. The ROC curves are zoomed-in with the FPR (x-axis) spanning from 0 to 0.05 in the small panel.

We further look into each subgroup of lensed quasars to compare 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 discards 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. 6 in which 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). This 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. Figure 9 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 this is because, 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. 10, zoomed-in with the x-axis spanning from 0 to 0.05. The crosses and triangles in Fig. 10 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 can 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 starts to lead to more decrease in TPR compared to FPR, and is not as advantageous.

Fig. 9. Sensitivity to the quasar image separation θLP of the lens search method using both spatial extent and the number of blobs. Bright: the group of lensed quasars with the third brightest image brightness m3rd <  22.0 mag. Faint: the group of lensed quasars with 22.0 mag ≤ m3rd <  24.0 mag. The top panels show the TPRs of each θLP bin, when pthrs% = 0% and Nthrs = 4, with the preselection by the number of blobs Ncrit = 2. The bottom panels show the number of mock lenses in each θLP bin.

Fig. 10. Bright lensed quasars (m3rd <  22.0 mag) ROC curves with wide separation (θLP >  1.5″) of the lens search method with the preselection on the number of blobs at different Ncrit values.

5.3. Classification based on one single epoch

It is ideal and most efficient to detect lensed quasars through only one epoch (in addition to the reference image). If the lens search algorithm is able to identify 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. In this section, 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 employs 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 classifies a HSC variable m′ as a candidate for lensed quasar if

$$\begin{array}{c}\hfill A_{\mathrm{eff},t\prime }^{m\prime }>A_{\mathrm{eff},t\prime }(p_{\mathrm{thrs}}),\end{array}$$(20)

where $A_{\mathrm{eff},t\prime }^{m\prime }$ 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. 11. The result in the left panel of Fig. 11 is from the epoch with the best lens search performance among the 13 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. 11 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%. The worst performance we have among the 13 epochs is presented in the right panel of Fig. 11.

Fig. 11. Lens search performance based on a single epoch. The lens search algorithm is tested using difference images from only one single epoch for all the 13 epochs in the HSC transient survey; the ROC curves of the epochs with the best performance (left), the average performance (middle), and the worst performance (right) are shown. For lens search algorithm based on a single epoch, only pthrs is used as constraint. In each panel, the ROC curves are plotted by varying pthrs%; the diamonds, triangles, and circles indicate the points on the ROC curves at pthrs% = 90%,95%, and 97.5%. The seeing of each epoch is also indicated in each panel. The best lens search performance happened 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 pthrs% = 97.5%. In a single epoch with average seeing (∼0.98″), the lens search algorithm can detect the bright lensed quasars with wide separation at (TPR, FPR) = (94.2%,5.0%) when pthrs% = 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.

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, the better the lens search performance is. We plot d_t against the seeing for each epoch t in Fig. 12. The quantity d_t becomes small when seeing is around 0.7 arcsec, indicating that our lens search algorithm has good performance when the single epoch has seeing around 0.7 arcsec. In Fig. 12, we also see that d_t suddenly increases when seeing is better than 0.5 arcsec, which verifies that when the seeing is too good, our lens search performance is heavily affected by the significant artifacts from the sharp images.

Fig. 12. Relation between the lens search performance and the seeing in each single epoch from the HSC transient survey. The quantity dt is defined as the minimum distance between (TPR, FPR) = (100%,0%) and the ROC curve of the bright lensed quasars (m3rd <  22.0 mag) with wide separation (θLP >  1.5″) for each epoch t. The 13 star symbols indicate the dt 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.

6. Conclusions and discussion

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 in this work is useful not only for lensed quasar search, but also for other studies, such as studying lensed supernovae and testing lens mass modeling. Our lens search method builds upon the method 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%).

• 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 catalogs 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 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. For example, we could apply the method based on the spatial extendedness in Sect. 4.1 to the other bands and select the variable objects appearing as spatially extended in the difference images across multiple bands. Doing so could decrease the number of false-positive lens candidates from one single band. The thresholds on the effective region, p_thrs%, and the number of epochs, N_thrs, for the other bands depend on the survey, and we need the mock lenses for the other bands to investigate them. We will explore the multiband 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. Given the depth and the covering area, we estimate that 1 ± 1 quad(s) lie(s) in the footprint of the HSC transient survey up to June 2017, based on the lensing rates in OM10 and the correction for the faint end (i-band total brightness m_total >  18.1 mag) in Agnello & Spiniello (2019). There is already one known quad lying in this footprint (Anguita et al. 2009). Therefore, we expect the recovery of this known lensed quasar and a possible discovery of a new quad when we apply our lens search algorithm on ∼141 000 variable objects in the HSC transient survey. Most of the contaminants in the variable objects are very bright objects, such as bright variable stars or bright galaxies. (Chao et al., in prep.)

Moreover, 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. The 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 detected in the LSST (OM10). Although the LSST pixel size (0.2″) is slightly larger than the HSC transient survey (0.17″), the effective region is expected to be still distinguishable between potential lensed quasar candidates and non-lensed variable sources, as long as we adopt a large value for p_thrs% (p_thrs >  92.5%). Therefore, our lens search algorithm is still applicable for the LSST, if p_thrs is sufficiently large. The application of single-epoch-based lens search in Sect. 5.3 in the LSST should mostly depend on the seeing of the epoch that is used for lens search. However, we still need the mock lenses simulated for the LSST to investigate in detail the impact on our lens search performance from their image quality and cadence strategy. We expect our lens search algorithm to skilfully capture the new lensed quasars in the ongoing and upcoming cadenced surveys.

Acknowledgments

We thank Brandon Kelly, Takashi Moriya, Hiroko Niikura, Masaomi Tanaka, and Atsunori Yonehara for useful discussions. DCYC and SHS thank the Max Planck Society for support through the Max Planck Research Group for SHS. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Research, Innovation and Science. JHHC acknowledges support from the Swiss National Science Foundation (SNSF). This work was supported by JSPS KAKENHI Grant Numbers JP15H05892 and JP18K03693. ATJ is supported in part by JSPS KAKENHI Grant Number JP17H02868. This research made use of Astropy, (http://www.astropy.org) a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University.

References

Appendix A: Sensible separation or sensible time delay

We examine whether our lens search sensitivity has a dependence on the lensed quasar time delay by checking the longest time delay, Δt_max, in each lensed quasar. Figure A.1 shows the distribution of Δt_max in each θ_LP bin. The top panel in Fig. A.1 shows the mean value of Δt_max with 1σ standard deviation as the error bar from the lensed quasars in each θ_LP bin. The top panel in Fig. A.1 shows that both the group of non-selected lensed quasars and the group of selected lensed quasars by our lens search method (p_thrs% = 95%, N_thrs = 9) have similar distribution of Δt_max as all the lensed quasars, indicating that the lens search sensitivity has no preference on the lensed quasar time delay. On the other hand, the middle panel in Fig. A.1 shows again that the number of selected lensed quasars decreases as the θ_LP decreases, while the number of non-selected lensed quasars increases as the θ_LP decreases. Therefore, our lens search algorithm is sensitive to the lensed quasar separation not to the lensed quasar time delay.

Fig. A.1. Distribution of the longest time delay Δtmax for the non-selected lensed quasars and the selected lensed quasars. Top panel: mean values of Δtmax with 1 standard deviation as the error bars from the lensed quasars in each θLP bin. Middle panel: number of non-selected lensed quasars, the number of selected lensed quasars, and the number of total mock lenses in each θLP bin. Bottom panel: TPR in each θLP bin. The sensitivity of the lens search algorithm mainly comes from the separation (θLP), not the time delay (Δtmax).

Appendix B: Sensitivity to total (unresolved) brightness

For a more realistic study, we investigate the lens search sensitivity to the lensed quasars total brightness, m_total (i-band). Figure B.1 shows the sensitivity to m_total, when applying the lens search method based only on spatial extent at (p_thrs%,N_thrs) = (95%,9) (Sects. 4.1 and 5.1). As shown in Fig. B.1, regardless of the separation, our lens search algorithm could capture all the lensed quasars with m_total <  19.0 mag. We further examine the sensitivity for the lensed quasars with m_total >  19.0 mag in Fig. B.2. Figure B.2 shows that our lens search algorithm can still identify more than 70% of the lensed quasars with 19.0 mag <  m_total <  20.5 mag when θ_LP >  1.5″ (the largest separation among the pairs in one lens system), while the lens search performance is relatively poor for the lensed quasars with m_total >  19.5 mag when θ_LP ≤ 1.5″.

Fig. B.1. Sensitivity to the lensed quasar total brightness mtotal of the lens search method based only on spatial extent, when pthrs% = 95% and Nthrs = 9. The lens search algorithm could capture all the lensed quasars with mtotal 19.0 mag.

Fig. B.2. Sensitivity to the lensed quasar total brightness mtotal for the lensed quasars fainter than 19.0 mag. 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″.

All Tables

All Figures

Fig. 1. Flow chart illustrating the simulation of time-varying lensed quasars images. For quasars, mock quasar light curves are generated based on their source redshifts zs, and being shifted in time and magnified the mock quasar light curves on the image plane (see Fig. 2). For quasar hosts, their Sérsic light profiles are created on the source plane and their lensed images are produced on the image plane assuming a SIE model. Lens galaxies are selected from a real survey based on their lens redshifts zd 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 all the components for one lens system (image of lensed quasar, image of lensed quasar host, and image of lens galaxy) are produced, these components are convolved with a realistic PSF and are added to obtain one complete image.

In the text

	Fig. 2. Illustration of generating light curves for a lensed quasar without its (lensed) host galaxy. Left: positions of the four multiple images (A, B, C, and D) and the lens galaxy (cross). Top right: mock quasar light curve. Bottom right: the lensed light curves of the four multiple images in the left panel, for the period marked by the blue band in the top right panel.
In the text

	Fig. 3. Examples of lenses injected into the HSC transient survey. Top: single-epoch images co-added by the pipeline of the HSC transient survey. Bottom: corresponding difference images of the top row produced by the image subtraction in the pipeline of the HSC transient survey. The size of each image cutout is 10″ × 10″.
In the text

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 denoted by the white pixels show the spatial extent of the objects, after removing noise peaks in the 3σ masks. In the local extremum detection, pixels with values larger than all the neighboring pixels or smaller than all the neighboring pixels in the difference images within the effective regions are selected; the local extrema, indicated in cyan, correspond to the “blobs” (see Sect. 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 as a 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″.

In the text

Fig. 5. Flow chart illustrating the lensed quasar candidate selection procedure. The effective region (Eqs. (10)–(12)) and number of blobs (Eqs. (14)–(17)) are calculated for each variable object from the HSC transient survey (HSC variable). The HSC variables with the number of blobs larger than Ncrit are preselected. The percentile of the effective region for each epoch from all the HSC variables that remain after the preselection on the number of blobs is calculated. Finally, two thresholds, pthrs% and Nthrs, are set; the HSC variables are selected that have an effective region larger than pthrs% percentile for more than Nthrs epochs, as lensed quasar candidates.

In the text

Fig. 6. Lens search method ROC curves based only on spatial extent. A lens is considered to be bright if its third brightest image has m3rd <  22.0 mag, and faint if 22.0 mag ≤ m3rd <  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 pthrs% from 5% to 99.5% at Nthrs = 9 (see Sect. 5.1 for details). The diamonds, triangles, and circles indicate the points on the ROC curves when pthrs% = 90%,95%, and 97.5%, respectively. The lens search algorithm in this work could detect the bright lensed quasars with wide separation 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.

In the text

Fig. 7. Sensitivity to the quasar image brightness m3rd 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 m3rd bin, when pthrs% = 95% and Nthrs = 9. The bottom panels show the number of mock lenses in each m3rd bin. Our lens search algorithm can identify more than 90% of the lensed quasars with m3rd <  21.0 mag.

In the text

Fig. 8. ROC curves of the lens search method using both spatial extent and number of blobs. The subgroups of lensed quasars are the same as the subgroups in Fig. 6. The ROC curves are plotted by varying pthrs% from 0% to 99.5% at Nthrs = 4 with the preselection on the number of blobs Ncrit = 2. The crosses, triangles, and circles indicate the points on the ROC curves when pthrs% = 0%,55%, and 75%, respectively. With the preselection Ncrit = 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 Ncrit = 2. The ROC curves are zoomed-in with the FPR (x-axis) spanning from 0 to 0.05 in the small panel.

In the text

Fig. 9. Sensitivity to the quasar image separation θLP of the lens search method using both spatial extent and the number of blobs. Bright: the group of lensed quasars with the third brightest image brightness m3rd <  22.0 mag. Faint: the group of lensed quasars with 22.0 mag ≤ m3rd <  24.0 mag. The top panels show the TPRs of each θLP bin, when pthrs% = 0% and Nthrs = 4, with the preselection by the number of blobs Ncrit = 2. The bottom panels show the number of mock lenses in each θLP bin.

In the text

	Fig. 10. Bright lensed quasars (m3rd <  22.0 mag) ROC curves with wide separation (θLP >  1.5″) of the lens search method with the preselection on the number of blobs at different Ncrit values.
In the text

Fig. 11. Lens search performance based on a single epoch. The lens search algorithm is tested using difference images from only one single epoch for all the 13 epochs in the HSC transient survey; the ROC curves of the epochs with the best performance (left), the average performance (middle), and the worst performance (right) are shown. For lens search algorithm based on a single epoch, only pthrs is used as constraint. In each panel, the ROC curves are plotted by varying pthrs%; the diamonds, triangles, and circles indicate the points on the ROC curves at pthrs% = 90%,95%, and 97.5%. The seeing of each epoch is also indicated in each panel. The best lens search performance happened 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 pthrs% = 97.5%. In a single epoch with average seeing (∼0.98″), the lens search algorithm can detect the bright lensed quasars with wide separation at (TPR, FPR) = (94.2%,5.0%) when pthrs% = 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.

In the text

Fig. 12. Relation between the lens search performance and the seeing in each single epoch from the HSC transient survey. The quantity dt is defined as the minimum distance between (TPR, FPR) = (100%,0%) and the ROC curve of the bright lensed quasars (m3rd <  22.0 mag) with wide separation (θLP >  1.5″) for each epoch t. The 13 star symbols indicate the dt 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.

In the text

Fig. A.1. Distribution of the longest time delay Δtmax for the non-selected lensed quasars and the selected lensed quasars. Top panel: mean values of Δtmax with 1 standard deviation as the error bars from the lensed quasars in each θLP bin. Middle panel: number of non-selected lensed quasars, the number of selected lensed quasars, and the number of total mock lenses in each θLP bin. Bottom panel: TPR in each θLP bin. The sensitivity of the lens search algorithm mainly comes from the separation (θLP), not the time delay (Δtmax).

In the text

	Fig. B.1. Sensitivity to the lensed quasar total brightness mtotal of the lens search method based only on spatial extent, when pthrs% = 95% and Nthrs = 9. The lens search algorithm could capture all the lensed quasars with mtotal 19.0 mag.
In the text

	Fig. B.2. Sensitivity to the lensed quasar total brightness mtotal for the lensed quasars fainter than 19.0 mag. 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″.
In the text