© ESO 2022

1. Introduction

Gravitationally lensed quasars provide a powerful means to study both galaxy evolution and cosmology, and appear when a foreground lensing galaxy deflects the light from a distant quasar source. The configuration of lensed images enables us to study the mass structures and substructures of lensing galaxies and to further probe galaxy evolution in the Universe (Suyu et al. 2012; Dalal & Kochanek 2002; Vegetti et al. 2012; Nierenberg et al. 2017; Gilman et al. 2019). The time delays between image pairs can be used to infer the Hubble constant, H₀, a crucial cosmological parameter in determining the size, age, and critical density of the Universe (Refsdal 1964; Chen et al. 2019; Wong et al. 2020). Given the current tension between H₀ measurements from the early and late Universe, independent measurements from lensed quasars are vital (Verde et al. 2019).

There are currently around 200 known gravitationally lensed quasars. The Cosmic Lens All-Sky Survey (CLASS; Myers et al. 2003) provided the largest statistical sample of radio-loud lensed quasars by identifying the flat-spectrum radio sources with multiple lensed images. From optical data, the SDSS Quasar Lens Search discovered 62 lenses based on both morphological and colour selection of spectroscopically confirmed quasars (Oguri et al. 2006, 2008, 2012; Inada et al. 2008, 2010, 2012). Jackson et al. (2012) extended the SQLS sample by using the better image quality from the UKIRT Infrared Deep Sky Survey (UKIDSS). Given the increasing depth and image quality of ongoing surveys, it is now possible to identify lensed quasars from imaging alone, with no spectroscopic pre-selection.

Chan et al. (2015) demonstrated such a search using CHITAH – an automated pixel- and mass-modelling code for finding lensed quasars in the Hyper-Suprime Cam (HSC) Survey (Aihara et al. 2018). Other image-based lens-finding strategies have relied on machine learning, ring-finders, and citizen science (Sonnenfeld et al. 2018, 2019, 2020; Wong et al. 2018; Chan et al. 2020; Jaelani et al. 2020, 2021), however these have mainly had success for lensed galaxies where the extended lensed arcs are distinguishable from other astrophysical objects. In the case of lensed quasars, often just two point sources are present, which outshine the light from the lensing galaxy and are outnumbered by visually similar contaminants like binary stars or quasar–star projections. Lemon et al. (2018, 2019) used simultaneous optical–infrared extracted colours combined with astrometric measurements from Gaia to avoid such contaminants. Thanks to its exceptional spatial resolution, the catalogued image positions in Gaia have been used to identify new lenses (Ducourant et al. 2018a,b; Krone-Martins et al. 2018).

Here, we present a search for lensed quasars in the Canada-France Imaging Survey (CFIS; Ibata et al. 2017). We analyse CFIS Data Release 2 (DR2) data, covering 2500deg². While the excellent image quality and depth of the r-band allow for detection of faint quasar images and lensing galaxies, the u-band provides an important colour indicator for removal of many common contaminants. Despite previous spectroscopic searches in this area of the sky, we expect many lenses to still be left undiscovered either due to lack of depth and/or stringent pre-selection of previous search techniques. The CFIS data, thanks to their depth and seeing, also consist in a major advantage to find these missing lenses.

This paper is organised as follows. The descriptions of the CFIS pixel data and the quasar catalogue selection are provided in Sect. 2. We describe the automated part of our lens search method in Sect. 3, and the visual inspection stage in Sect. 4. Prioritisation of high-grade systems for spectroscopy is outlined in Sect. 5, with the results of follow-up spectroscopy presented in Sect. 6. We present our results and discuss the individual systems in Sect. 7. In Sect. 8 we conclude our findings. All images are oriented with north up and east to the left. Optical magnitudes quoted in this paper are in the AB system, and infrared magnitudes in the Vega system. When required, a flat cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3 and Ω_Λ = 0.7 is used.

2. Data

We first describe the available CFIS imaging data (Sect. 2.1), and then the initial quasar catalogues in which we perform our selection (Sect. 2.2).

2.1. CFIS imaging data

The Canada-France Imaging Survey (CFIS) is a component of the Ultraviolet Near Infrared Optical Northern Survey (UNIONS) project, which is a collaboration of wide-field imaging surveys of the northern hemisphere. UNIONS consists of CFIS, conducted at the 3.6-m Canada-France-Hawaii Telescope (CFHT) on Maunakea in Hawaii, members of the Pan-STARRS team, and the Wide Imaging with Subaru Hyper Suprime-Cam of the Euclid Sky (WISHES) team. CFHT/CFIS is obtaining deep u- and r-bands; Pan-STARRS is obtaining deep i-band and moderate-deep z-band imaging; and Subaru/WISHES is obtaining deep z-band imaging. These independent efforts are, in part, working towards securing optical imaging to complement the Euclid space mission, although UNIONS is a separate collaboration aimed at maximizing the science return of these large and deep surveys of the northern sky. Ibata et al. (2017), Fantin et al. (2019), and Guinot et al. (2021) provided an additional illustration of CFIS imaging data. This CFIS component is mapping the northern sky in the u-band (CFIS-u, covering 8000 deg²) and r-band (CFIS-r, covering 4800 deg²). For the present work, we use all available data from CFIS DR2 with both u- and r-band coverage, i.e. ∼2500 deg² (see Fig. 1 in Savary et al. 2021). The full survey is expected to finish in 2023. CFIS-r provides exquisite image quality with median seeing of ∼0.6″ to a depth of 24.1 (10σ, 2″diameter aperture), and CFIS-u has median seeing of ∼0.8″ to a depth of 23.6.

For all candidates, we generate image cutouts of 8.2″ × 8.2″ (44 pixels on-a-side). We also produce models of the point spread function (PSF) and its spatial variations across co-added images which have been reduced, processed, and calibrated at CADC using an improved version of the MegaPipe pipeline (Gwyn 2008, 2019). The model PSF is exploited with PSFEx (Bertin 2013) from the stacked image, and each image of the local PSF is oversampled by a factor of two. We notice that the optimal PSF model should be obtained from the individual exposure images, although the imperfect PSFs only marginally affect the detection of lensed quasars. More accurate and precise PSF models are beyond the scope of this work. Weight images along with other data quality diagnostics are also produced for each candidate in each of the available bands. The pixel scale of co-added images is 0.186″.

2.2. Quasar catalogues

The goal of this work is to find convincing lensed quasar candidates within the CFIS footprint, while minimising the manual inspection time. We therefore choose to start our search from catalogues of systems likely to contain at least one quasar. In future searches, we expect to be able to use the CFIS data themselves as a way to select quasars, because the u and r-bands provide an efficient way to select quasars as objects with a u-band excess (see e.g. Nakoneczny et al. 2019).

The Million Quasars Catalog v6.4 (hereafter MILLIQUAS; Flesch 2019) provides a compilation of both confirmed and candidate quasars. It is composed of all spectroscopically reported quasars, radio and X-ray detections, and probable quasars based on AllWISE colours and photometry (Secrest et al. 2015). Restricting the original catalogue of 1 980 903 objects to those with CFIS-u and -r data leaves 253 299 systems.

While around 30% of the MILLIQUAS catalogue comprises quasar candidates selected from infrared colours (Secrest et al. 2015), the colour-selection targets isolated quasars. We expect that many lensed quasars will have bluer infrared colours either due to blending with the lensing galaxy, or due to having high-redshift (z ≳ 3) sources. We therefore choose to augment the MILLIQUAS catalogue with a selection based on a bluer WISE limit of W1–W2 > 0.5. To minimise contaminants, we cross-match these WISE detections (requiring W1 < 16.5) to Gaia detections, keeping only those systems with two Gaia detections within 3″ of the AllWISE catalogue detection (Secrest et al. 2015; Gaia Collaboration 2018). Finally, we impose a cut on the proper motion significance (PMSIG) of less than 15σ for each Gaia detection, as defined in Lemon et al. (2019). This leaves 9244 systems within the CFIS DR2 footprint, of which only 3015 currently have both CFIS-u and -r imaging. Including 253 299 systems in MILLIQUAS, we have 256 314 pre-selected quasar candidates in total to be classified as lensed quasar candidates.

3. Classification method

Our search algorithm relies upon identifying systems with multiple, similar-colour point sources, as is expected in lensed quasars. This straightforward approach can be effective when applied to systems known to contain at least one quasar point source, because in the case of a quasar+star projection, we expect different u − r colours (e.g. see Fig. 1 in Finlator et al. 2000). This separation of stars and quasars in colour-space is due to the u-band excess of quasars at redshifts z ≲ 2.7. At higher redshifts, the Lyman-alpha forest is redshifted into the u-band and we expect a higher rate of contamination from quasar+star projections. We note that the following algorithm is applicable to all imaging datasets with two or more bands; however we use CFIS as a testbed given its depth, excellent image quality, and the aforementioned use of the u-band.

The search algorithm is as follows:

1. Sharpening the point sources: Before comparing colours of multiple point sources, we require reliable identification of point sources in the ground-based imaging, where nearby point sources are often overlapping and blended. Conventional extraction algorithms are limited by this blending, and adjusting the deblending threshold parameters leads to multiple spurious detections from extended sources. However, we can identify positions of likely point sources in each of our candidate cutouts because we have the expected two-dimensional PSF profile at the location of the system. To achieve this, we measure the second derivative of the observed surface brightness and compare it with the second derivative of the PSF (see Fig. 7 in Cantale et al. 2016). In practice, we calculate the convolution of the Laplacian of the original image, D_ori, and the Laplacian of its corresponding PSF to obtain the sharpened image, D_shp, as:

$$\begin{array}{c}\hfill D_{\mathrm{shp}}=\mathrm{\Delta }{D}_{\mathrm{ori}}\ast \mathrm{\Delta }\mathrm{PSF},\end{array}$$(1)

where Δ is the Laplacian operator and is the convolution operator. The effect of this convolution is to reduce the size of the PSF FWHM by around 25%, allowing more efficient detection of small separation point sources. However, a by-product of this convolution on any PSF is a negative annulus around the sharpened PSF. We replace any negative values with zeros in D_shp, providing the final image, denoted as $D_{\text{shp}}^{\text{mask}}$, for the subsequent point source detection. In Fig. 1 we illustrate this sharpening on a known lens. As we have imaging data in two filters well-separated in colour, we can separate the lens and the lensed images based on the colour measurements, as described below.

Fig. 1. Point source detection with and without the convolution of Laplacian. (a) Measured PSF, (b) r-band data, (c) the sharpened image of PSF, and (d) the sharpened image of r-band data. The sharpened images are obtained by convolving the Laplacian of the image with the Laplacian of the corresponding PSF (see Eq. (1)). We replace the negative values with zeros in (c) and (d) due to the byproduct of this convolution. The red crosses in (b) and (d) indicate the point source detection using DAOStarFinder with 2D Gaussian kernels fitted from (a) and (c), respectively. The sharpened image (d) with improved dynamical range allows us to detect blended images better as shown in the red central cross.

2. Identifying the point sources: After obtaining the sharpened images, including their sharpened PSFs, we use the python module DAOStarFinder¹ for locating the possible point sources; see Figs. 1 and 2. DAOFIND, used in DAOStarFinder, searches images for local maxima with a peak amplitude greater than a given threshold and with a size and shape similar to the defined 2D Gaussian kernel (Stetson 1987). During this process, we supply the 2D Gaussian fit of the local PSF to DAOStarFinder². The result of this process is to return positions of identified point sources. In this step, we do not consider the flux measurements even though they are part of the output, as the sharpened images and PSFs create a bias in the flux estimation. We illustrate this in Figs. 1b and d.

Fig. 2. Output of the classification method on the known lens SDSSJ1537+3014. (a) CFIS composite image made from the u and r bands, (b) SED of the detected point sources, (c) u-band data, (d) u-band residuals upon subtraction of blue group point sources, (e) r-band data, and (f) r-band residuals upon subtraction of the blue point sources. Circles and diamonds in (a) show the positions of the blue and red point sources, respectively, with single-band detections overlaid in panels (c) and (e).

3. Grouping the point sources by position: The identification of point sources is run separately on each band, and can lead to different detections per system. This can be due to low signal-to-noise ratio in one band, differing PSF widths, or a very red or blue spectral energy distribution (SED). We collect all detections across both bands, and group any double detections into one detection when their distances are less than 1.5 pixels in the oversampled image (i.e. 0.14″), as shown in Fig. 2. The final position is taken as the mean position of each group.

4. Removing systems with one or no identified point sources: We remove any system with one or no point sources detected in the preceding step, as it is unlikely to contain the multiple point sources seen in a lensed quasar.

5. Measuring the fluxes of point sources: To extract the colours of each component, we first estimate the flux in each band by performing a least-squares fit of PSFs placed at their derived positions from step 3. The relevant χ² is:

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{i,j}\frac{{[\sum _{k=1}^n{P}_k(i,j)-D_{\mathrm{ori}}(i,j)]}^2}{\mathrm{var}(i,j)},}\end{array}$$(2)

where i = 1…N_x and j = 1……N_y are the pixel indices of the image cutout of dimensions N_x × N_y, P_k(i, j) is the k-th scaled PSF placed at the measured position, and var(i, j) is the pixel variance of D_ori(i, j). As we oversample our cutouts by a factor of two, N_x = N_y = 88. We note that, if detected, the fluxes of lensing galaxies are likely underestimated, but this is partly cancelled out when deriving the colour of the extended component.

6. Removing the faint point sources: To avoid the false detection on the noise peaks, we abandon the point sources beyond the observation depth of each tile in both bands. In cases where we miss faint lensed images, we extend the observation limit in CFIS to about 0.5 mag larger than the depth of 10σ point source.

7. Grouping the blue point sources: A main contaminant of our sample is single quasars with a nearby or coincident star or galaxy. Therefore we aim at separating the quasars from other objects, expecting quasars to be the bluest objects in any cutout. We first identify the bluest point source, i.e. the smallest u − r value, denoted as u_b − r_b, and collect any other detections with a u − r value within 1 mag of this bluest detection, i.e. |(u_k − r_k)−(u_b − r_b)| < 1 mag. We assign these detections, along with the original, as the blue group and any remaining point sources are assembled as the red group. As shown in Figs. 2a and b as an example, circles 1 and 2 (in blue and orange respectively) are classified in the blue group and diamond 3 (in green) is classified in the red group. We note that sometimes an object can be composed of only one group, particularly when no galaxies are present. Also the blue group may be arbitrarily red, for example when a quasar is dusty or at high redshift.

8. Selecting systems with multiple blue sources: As typical lensed quasar systems contain two or four images, we classify systems with 2 ≤ n_img ≤ 5 as lens candidates, where n_img is the number of detections in the blue group. The upper limit of n_img = 5 is empirically chosen to reduce contamination from dense star-forming galaxies or high-density stellar fields.

This method was developed alongside simultaneous tests on CFIS data of known lenses and contaminant systems. In particular, we used a set of 18 optically bright, known lensed quasars with CFIS u and r band imaging (from the lens compilation of Lemon et al. 2019)³. We recover 15 of these 18 lensed quasars, i.e. a recovery rate of 83%, and summarise these tests and the three failures in Appendix A. We now apply this automated search algorithm to the 256 314 pre-selected quasar candidates in Sect. 2, retaining 10 914 systems, i.e. reducing the initial catalogue to 4.3% of its original size. We compare the colours and flux ratios of this sample to the 15 recovered lenses in the CFIS footprint (see Fig. 3). The colour difference is measured using the brightest two potential lensed images (i.e. from the bluest group following step 7). We note that there is a limit on the absolute colour difference, as imposed during our search algorithm (i.e. |(u_k − r_k)−(u_b − r_b)| < 1 mag).

Fig. 3. Flux-ratio against colour difference for systems passing the algorithm detailed in Sect. 3, with known lenses overlaid as yellow stars. In the event of more than two detections, fluxes are based on the brightest detections within the blue group. The horizontal dashed line shows the flux ratio criterion we apply to further reduce the number of candidates.

As shown in Fig. 3, there are many candidate systems with large flux ratios that are not compatible with lens models and not seen in known lenses. Inspecting several of these high-flux-ratio systems indeed shows that they are mostly associated to companions near bright stars, or star-forming galaxies. We therefore apply a cut in the flux ratio, at an absolute magnitude difference of 2.5 mag in r-band, reducing the number of candidates by 20%. We retain 7815 candidates for the next step of visual inspection.

4. Visual inspection

We propose two phases for the visual inspection, namely ‘rapid’ and ‘refined’, performed by two authors (JC and CL). During the phase of rapid inspection, we adopt a simple binary classification: ‘possible lens’ and ‘non-lens’. For each candidate, we inspect an image of the system with the results of the analysis pipeline of Sect. 3. An example of such an image is shown as Fig. 2. We include the residual images after subtracting the point sources associated to the blue group in order to highlight any potentially faint compact lensing galaxies. The majority of the candidates are classed as non-lens for one or several of the following reasons:

1. Bad pixels in the cutout which lead to spurious point source detections from our pipeline;

2. Obvious extended galaxies or non-point sources, such as the host galaxy of a low-redshift AGN or a galaxy merger;

3. Noisy u-band data;

4. Wide-separation point sources with no apparent lensing galaxy.

Each case is illustrated in Fig. 4. In hindsight, the third case could have been avoided by increasing the threshold during the classification method; however, these are easily removed during our rapid inspection phase. Of the 7815 inspected systems, CL/JC classifies 514/640 as ‘possible lens’. At this stage, we check that both inspectors recovered all of the 15 known optically bright lensed quasars in the sample (see Appendix A for the classification of known lenses by the selection method of Sect. 3). Given the sufficient number for inspecting the candidates in detail, i.e the phase of refined inspection, we take the intersection of the two classifications as 239.

Fig. 4. Examples of objects classified as non-lens during the rapid visual inspection because of bad pixels (case 1), mergers/low-redshift AGNs (case 2), noise peaks (case 3), or wide binaries (case 4).

For each of these of 239 systems, we carefully inspect any existing data, including SDSS images and spectra, DECaLS images⁴ (Dey et al. 2019), unWISE images (Meisner et al. 2018), and Gaia DR2 proper motions (converted to the PMSIG parameter defined in Lemon et al. 2019), and astrometric excess noise (AEN; Gaia Collaboration 2018). A large PMSIG (≳10σ) is a strong indicator that the system contains a star, and a high AEN (≳10 mas) suggests that a candidate is possibly a star-forming galaxy (see Figs. 1 and 2 in Lemon et al. 2019). We grade the candidates with values of between 0 and 10, where:

• 10–9:

  very likely to be a lens, requiring resolved spectroscopy for confirmation,

• 8–7:

  blue point sources at modest separation with similar quasar-like colours and without an obvious lensing galaxy,

• 6–5:

  same as the grades 8–7 but with slightly high PMSIG, AEN, and/or more different colours,

• 4–0:

  a system that likely contains a star, star-forming galaxy, or a very low-redshift quasar.

Our final grade for each system (G_av) is the average of the two grades resulting from the rapid and refined inspections outlined above. There are 26 candidates graded G_av ≥ 8, detailed in Table 1 with their CFIS cutouts and colours in Fig. 5. Before spectroscopic follow-up (Sect. 6), we prioritise the candidates with detailed light modelling (Sect. 5).

Fig. 5. Lens candidates with Gav ≥ 8 or with follow-up spectroscopy with ALFOSC on the NOT. Each candidate contains ‘SED’, ‘colour’, and ‘Data−PS’ from light modelling. The labels in ‘SED’ and ‘colour’ are the same as Figs. 2a and b. The number in square brackets indicates the grade of visual inspection [Gav]. ‘Data−PS’ illustrates the residual of point source subtraction. The numbers of point sources and galaxies for the fit are indicated at the top of the panels as NPS and NGAL, respectively. The results of light modelling for each component are listed in Table B.1. We label the confirmed lenses with green stars, the probable lenses with yellow stars, the NIQs (nearly identical quasars) with blue diamonds, the quasar pairs with yellow diamonds, the contaminants with red crosses, and the inconclusive systems with yellow squares. The systems with NOT spectra are labelled with cyan triangles.

5. Light modelling

A lensed quasar system should contain a lensing galaxy surrounded by multiple lensed images. In most cases, a lensing galaxy is fainter than a quasar and is barely to be seen in ground-based imaging. Although our classification method is able to highlight lensing galaxies by removing lensed quasar images (see Fig. 2 as an example), imperfect PSFs from CFIS data sometimes blend lens light and residuals, or residuals of double point sources can appear to mimic a lensing galaxy. We therefore attempt to optimally subtract the obvious point sources (or single point source and galaxy) in search of a lensing galaxy (or counter image) in the deep, sharp CFIS imaging for the most promising systems resulting from visual inspection.

A Moffat profile (Moffat 1969) is used for the PSF (and hence also quasar images), and a Sérsic profile (Sérsic 1963) for lensing galaxies. We manually decide the numbers of point sources and galaxies for each system, and fit the u and r bands simultaneously. The PSF parameters are fitted directly from the point sources in our candidates. Results from modelling the best candidates are shown in Fig. 5 and Table B.1.

6. Spectroscopic follow-up

Long-slit follow-up spectra were taken with the Alhambra Faint Object Spectrograph and Camera (ALFOSC) mounted on the 2.56-m Nordic Optical Telescope (NOT) on the nights of 16 and 17 April 2021 with seeings of ∼1.2″ and ∼0.75″ respectively. To reliably identify a priori unknown quasar redshifts, we require multiple broad lines. We therefore used grism #4 with a wide wavelength coverage of 3200 Å–9600 Å and dispersion of 3.3 Å pixel⁻¹. Slit position angles were calculated from the best-fit image positions from the CFIS light modelling (Sect. 5). The slit width was 1″ and no atmospheric dispersion corrector was used to optimise throughput.

The 2D spectra were bias-subtracted and cosmic rays masked, and a model for the sky background based on the 5″ either side of the object was subtracted. We extract separate spectra using a forward-modelling process, as described fully in Lemon et al. (2022). Briefly, we forward model the quasar images as Moffat profiles onto a pixelised grid, and derive the 1D spatial spread function after applying a slit of 1″ in width. The Moffat parameters, image separation, and spatial variation along the slit are found in each wavelength bin using MCMC sampling. A second-order polynomial fit to these parameters then allows a model at each pixel in wavelength bins, followed by a least-squares fit to determine the flux. The residuals are inspected for goodness-of-fit and the spectra are binned into equal noise bins to facilitate comparison of spectra.

We prioritise our systems according to the visual inspection and the light modelling. In total we, obtained spectra for 18 systems, four of which have grades G_av < 8 due to our observing constraints. The spectra are shown in Fig. 6, and determined redshifts and classification are listed in Table 1.

Fig. 6. Extracted 1D spectra of candidates observed with ALFOSC on the NOT. Redshifts are measured from the average of the peaks of broad emission lines. The CFIS colour image is shown to the right of each spectrum, with the slit position overlaid as a green box.

7. Result

Combining our spectroscopy and light modelling, the results for our best 30 lens candidates are summarised in Table 1. We note that redshifts are measured from either NOT or other existing spectrographs.

In this work, we confidently confirm five new lensed quasars: J0943+3532, J1111+3804, J1146+3254, J1201+4013, and J1555+3118. One further system, J1559+3839, is a highly probable lensed quasar; however, we did not obtain spatially resolved spectroscopy for this system. There are six systems that we designate as potential lensed quasars, as they show nearly identical quasar (NIQ) spectra in two components. Two systems show double quasar spectra at the same redshift, but different profiles in emission lines. We therefore categorise these latter two as quasar pairs. We also identify two contaminants in our best list, because we can only fit the light distribution with a single point source and a galaxy. Our results for the remainder of the sample are inconclusive.

We note that there are three systems with X-ray detections, namely J1325+3629, J1555+3813, and J1559+3839 (a probable lens), which can be used to study the inner structure of quasars (Chartas 2000). One of the NIQs and a possible lens, J1111+3804, which can be a possible lens, shows a radio detection. If this system is confirmed as a lens, we can exploit the radio flux ratio anomalies induced by substructures in lensing galaxies. Estimating the number of substructures allows us to discern diverse dark matter models (Harvey et al. 2020). The redshift range in our sample of lens candidates spans from 0.24 to 3.36 (from 1.21 to 3.36 for confirmed lenses), showing that our classification method is still able to detect lensed quasars at high redshift (z ≳ 2.7), if there is no blue object located in the image cutout. Most of our candidates have small image separations, ranging from 0.6″ to 4.2″ (from 0.8″ to 2.5″ for confirmed lenses). The smallest separation is approximately equal to the seeing in CFIS-r, which can help to extend the sample size for lensing statistics. The lenses confirmed in this work are doubly lensed quasars, which can potentially serve as tools for the measurement of H₀ (Birrer et al. 2019). We do not find quadruply lensed quasars because they are missing in our quasar catalogues. The potential binary quasars found in this work, especially those with close separations (∼kpc), can be used for the study of low-frequency gravitational waves from dual supermassive black holes (Goulding et al. 2019; Casey-Clyde et al. 2022), and those at high redshifts (z ≳ 2) can even provide extra science cases with which to study the distant double quasars, such as dynamical evolution and successive mergers in the early Universe (Shen et al. 2021). We detailed a few interesting systems in the individual notes as follows.

7.1. J0112+3307 (inconclusive)

This system was given a high average grade of 9 (see Sect. 4), though was unfortunately not visible during our spectroscopic follow-up. We measure the redshift at z = 0.243 using the LAMOST spectra (Luo et al. 2015), as shown in Fig. C.1. Given the image separation and source redshift, we estimate the velocity dispersion σ_v = 250 km s⁻¹ (assuming lens redshift at 0.1). Although a galaxy component can be detected from the light modelling (see Fig. 5), it seems too faint to achieve such high velocity dispersion (Bolton et al. 2008). We are therefore unable to classify this system as a lens. High-resolution spectroscopic or imaging data are needed to confirm its nature.

7.2. J0943+3532 (lens)

Given the small separation of 0.831″ between the two components, the spectrum shows a single blended component, with emission lines of a quasar at z = 1.21. Gaia detects two point-source components (with small excess noise and proper motion) and the CFIS image modelling requires an additional extended component between these two point sources (see Fig. 5). We therefore classify this system as a lens.

7.3. J1006+3040 (NIQ)

The two components have the same CIII and MgII lines. We noticed that H and K of CaII absorption and OII emission lines appear in the yellow spectrum at z ∼ 0.3. However, we are not able to see an extended component after fitting two point sources (see Fig. 5). This system requires higher resolution observations to confirm its lens feature.

7.4. J1039+5126 (NIQ)

SiIV+OIV, CIV, and CIII emission lines fit the same redshift at z = 2.05. The slight deviation in CIV profiles could be attributed to microlensing. We therefore classify J1039+5126 as a NIQ.

7.5. J1111+3804 (lens)

Quasar source at high z = 3.02 and lensing galaxy can be seen in light modelling. An existing SDSS spectrum centred on the lower quasar confirms the source redshift at z = 3.018. The top component is redder in the spectrum (blue in Fig. 6), which suggests some light leaking from a galaxy. We can see the 4000 Å break at H and K lines of CaII and G at z ∼ 0.6. We also note that there is a strong radio detection for this system from VLASS.

7.6. J1146+3254 (lens)

This new lens is an example of the power provided by the improved depth and seeing of CFIS in finding lensed quasars. A single point source and galaxy subtracted from Pan-STARRS or DECaLS⁵ images does not show any counter image, whereas in CFIS it is immediately apparent. Our resolved spectra detect the CIV line of the counter image blended with the lensing galaxy, which we measure to be at z ∼ 0.3.

7.7. J1201+4013 (lens)

This system was originally considered by the SDSS quasar lens search as a potential lens, SDSS J120118.92+401318.1, and it was selected as a morphological candidate with an SDSS QSO spectrum at z = 1.936. Inada et al. (2010) ruled this out as a possible lens based on I-band imaging from Tek2k CCD at UH88 which was fit as a single point source+galaxy model. Thanks to the deeper and higher resolution spectra provided by Gaia, we selected this system as having two Gaia detections, and thus highly likely to have two point sources present in the system. Modelling of the CFIS data requires a counter image and a detection of the lensing galaxy as well (see Fig. 5). Resolved NOT spectra clearly show CIV and CIII emission from the fainter component, with a redder continuum, which is likely due to the lensing galaxy. We subtract the spectrum of A from B, and rescaled by the median flux ratio below 4800 Å. Excluding microlensing and differential reddening effects, this leaves the spectrum of the lensing galaxy, from which we can see absorption lines and a 4000 Å break associated to z ∼ 0.4.

7.8. J1219+3453 (NIQ)

Although two quasar images have the same redshift of z = 1.15 and similar spectra, we are not able to confirm the nature of the central object (see the red line in Fig. 6); it is well fit by a PSF and is therefore likely a star. High-resolution imaging is required to determine the presence or lack of a lensing galaxy. We therefore classify this system as a NIQ instead of a lens.

7.9. J1325+3629 (inconclusive)

This system shows the X-ray detection in Chandra imaging. The faint component is likely to be a contaminant, because its spectrum does not have apparent emission lines. However, the very small separation makes it difficult to separate the 2D spectrum. We need high-resolution imaging or better resolved spectroscopy.

7.10. J1454+4146 (quasar pair)

Although two components have the same redshift, CIII and CIV profiles are not the same (see Fig. 6). We classify this system as a quasar pair.

7.11. J1537+3223 (NIQ)

Two components have very similar spectra at z = 1.825. However, because of the small separation, we are not able to see a lensing galaxy. We classify this system as a NIQ.

7.12. J1555+3118 (lens)

This system is a u-band dropout recovered by our selection at z = 3.360 with a lensing galaxy detected. It was originally selected in the Gaia double-WISE cross-match catalogue and not in the MILLIQUAS sample because of its marginally bluer WISE colours of W1 − W2 = 0.67 and W2 − W3 = 3.32, typical of higher redshift quasars (Secrest et al. 2015).

7.13. J1555+3813 (quasar pair)

Different Lyman alpha and CIV line profiles suggest that the two point sources are not from the same quasar source (see Fig. 6). In addition, we do not see an extended component between these two point sources (see Fig. 5). This system has a small image separation 0.749″ (6.04 kpc at z = 2.525). We notice that this system has X-ray detection in ROSAT.

7.14. J1559+3839 (probable lens)

Although we do not obtain the NOT follow-up spectra for this system, there is an SDSS spectrum with a QSO pipeline redshift of z = 4.37. However, closer inspection of the spectrum reveals the main emission line is likely MgII with apparent CIII, NeV, and NeVI at z = 1.334 (see Fig. C.1). This is also the classified redshift of the improved classification of Hewett & Wild (2010). The modelling residuals require a galaxy between the two point sources. We therefore classify it as a probable lensed quasar. The strong reddening of the background quasar can either be intrinsic, as in W2M J1042+1641 (Glikman et al. 2018) and MG 1131+0456 (Stern & Walton 2020), or due to the lensing galaxy. We also notice that this system has an X-ray detection in ROSAT.

7.15. J1739+4355 (NIQ)

Two quasars have similar spectra with a small image separation 1.538″ (12.4 kpc at z = 2.525). We notice that there is a big difference in Lyα or in the NV line (1240.22 Å). The deviation of Lyα could be attributed to microlensing, so we still classify this system as a NIQ/possible lens.

8. Conclusion

In this work, we present new lens candidates in CFIS DR2 (a component of UNIONS), covering 2500deg² with depths of u = 23.6 mag and r = 24.1 mag. We confirm the lens features based on light modelling and spectroscopic follow-up. We draw the following conclusions:

• We begin from the catalogues with MILLIQUAS and Gaia pairs, and then extract the image cutouts in CFIS DR2 with u and r bands. The total number of objects is 256 314.

• We designed an efficient classification method: searching for a system as a lens candidate that contains multiple blue point sources, from 2 to 5 in number. This method is able to recover 15 of the 18 known lenses.

• We classified 10 914 candidates with the classification method, and further removed 20% of candidates with the cut of r-mag difference with |r₁ − r₂|< 2.5 mag.

• We visually inspected 7815 candidates, finding 239 promising ones which we then graded from 0 (non-lens) to 10 (lens) according to criteria we outline in Sect. 4.

• We carried out a follow-up analysis of 30 high-grade candidates with light modelling and 18 of them are observed with ALFOSC on the NOT.

• We confirm five new lenses (J0943+3532, J1111+3804, J1146+3254, J1201+4013, and J1555+3118) and one probable lens (J1559+3839). Six systems are nearly identical quasars, which are potential lensed quasars. Two systems are classified as quasar pairs, showing double quasar spectra but different emission line profiles. We also identified two contaminants during the light modelling.

CFIS imaging data allow us to see more lenses of small separation (e.g. the smallest one is 0.831″ in this work), thanks to its high resolution in r band. The full CFIS sample will continue to be analysed as part of a future exploration. Upcoming data releases will contain deep and high-resolution imaging in the u, g, r, i and z bands from all the components of the UNIONS collaboration. The g and z bands obtained with Subaru/HSC will also allow us to make the most of the benefits of high angular resolution. The success of this work demonstrates that the classification method can be applied to the higher resolution imaging data in future, such as Euclid and the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST), and will therefore be exceptionally useful for finding the lenses with the narrowest separation on the sky.

Acknowledgments

This work is supported by the Swiss National Science Foundation (SNSF) and by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (COSMICLENS: grant agreement No 787886). Part of this work benefited from the Canadian Advanced Network for Astronomical Research (CANFAR) and Compute Canada facilities. R. G. thanks IoA and the Churchill College in Cambridge for their hospitality and acknowledges local support from the French government. G. V. has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodovska-Curie grant agreement No. 897124. This work is based on data obtained as part of the Canada-France Imaging Survey, a CFHT large program of the National Research Council of Canada and the French Centre National de la Recherche Scientifique. Based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA Saclay, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l’Univers (INSU) of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. This research is based in part on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. We are honored and grateful for the opportunity of observing the Universe from Maunakea, which has the cultural, historical and natural significance in Hawaii. Pan-STARRS is a project of the Institute for Astronomy of the University of Hawaii, and is supported by the NASA SSO Near Earth Observation Program under grants 80NSSC18K0971, NNX14AM74G, NNX12AR65G, NNX13AQ47G, NNX08AR22G, YORPD20_2-0014 and by the State of Hawaii.

References

Appendix A: Test for known lenses

There are 18 optically bright known lensed quasars with CFIS DR2 imaging. We recover 15 of them via the number of blue source detections 2 ≤ n_img ≤ 5; see Fig. A.1 and Table A.1. The three failures are listed below:

1. SDSSJ0909+4449: overly wide separation (around 18″, Shu et al. 2018);

2. SDSSJ1452+4224: bad pixels in CFIS-u;

3. J1612+3920: omission of one lensed image identification.

Fig. A.1. Known lens test for the classification method in Sect. 3. The labels are same as Fig. 2 (a) and (b). The lenses that are not able to be recovered are highlighted with the red frames. SDSSJ0909+4449 has an overly wide image separation ∼18″. SDSSJ1452+4224 has bad pixels in CFIS-u image 1. J1612+3920 has a faint counter image hidden to the bottom-left of the lensing galaxy.

Appendix B: Result from light modelling

We measure the CFIS astrometry and photometry from light modelling in Table B.1. See details in Sect. 5. These candidates are the same as in Fig. 5 and Table 1.

Appendix C: Two unresolved spectra

We measure the redshifts of two additional systems: J0112+3307 and J1559+3839, from the public spectra of LAMOST and SDSS, respectively; see Fig. C.1.

Fig. C.1. One-dimensional spectra of J0112+3307 from LAMOST and J1559+3839 from SDSS.

All Tables

All Figures

Fig. 1. Point source detection with and without the convolution of Laplacian. (a) Measured PSF, (b) r-band data, (c) the sharpened image of PSF, and (d) the sharpened image of r-band data. The sharpened images are obtained by convolving the Laplacian of the image with the Laplacian of the corresponding PSF (see Eq. (1)). We replace the negative values with zeros in (c) and (d) due to the byproduct of this convolution. The red crosses in (b) and (d) indicate the point source detection using DAOStarFinder with 2D Gaussian kernels fitted from (a) and (c), respectively. The sharpened image (d) with improved dynamical range allows us to detect blended images better as shown in the red central cross.

In the text

Fig. 2. Output of the classification method on the known lens SDSSJ1537+3014. (a) CFIS composite image made from the u and r bands, (b) SED of the detected point sources, (c) u-band data, (d) u-band residuals upon subtraction of blue group point sources, (e) r-band data, and (f) r-band residuals upon subtraction of the blue point sources. Circles and diamonds in (a) show the positions of the blue and red point sources, respectively, with single-band detections overlaid in panels (c) and (e).

In the text

	Fig. 3. Flux-ratio against colour difference for systems passing the algorithm detailed in Sect. 3, with known lenses overlaid as yellow stars. In the event of more than two detections, fluxes are based on the brightest detections within the blue group. The horizontal dashed line shows the flux ratio criterion we apply to further reduce the number of candidates.
In the text

	Fig. 4. Examples of objects classified as non-lens during the rapid visual inspection because of bad pixels (case 1), mergers/low-redshift AGNs (case 2), noise peaks (case 3), or wide binaries (case 4).
In the text

Fig. 5. Lens candidates with Gav ≥ 8 or with follow-up spectroscopy with ALFOSC on the NOT. Each candidate contains ‘SED’, ‘colour’, and ‘Data−PS’ from light modelling. The labels in ‘SED’ and ‘colour’ are the same as Figs. 2a and b. The number in square brackets indicates the grade of visual inspection [Gav]. ‘Data−PS’ illustrates the residual of point source subtraction. The numbers of point sources and galaxies for the fit are indicated at the top of the panels as NPS and NGAL, respectively. The results of light modelling for each component are listed in Table B.1. We label the confirmed lenses with green stars, the probable lenses with yellow stars, the NIQs (nearly identical quasars) with blue diamonds, the quasar pairs with yellow diamonds, the contaminants with red crosses, and the inconclusive systems with yellow squares. The systems with NOT spectra are labelled with cyan triangles.

In the text

	Fig. 6. Extracted 1D spectra of candidates observed with ALFOSC on the NOT. Redshifts are measured from the average of the peaks of broad emission lines. The CFIS colour image is shown to the right of each spectrum, with the slit position overlaid as a green box.
In the text

	Fig. A.1. Known lens test for the classification method in Sect. 3. The labels are same as Fig. 2 (a) and (b). The lenses that are not able to be recovered are highlighted with the red frames. SDSSJ0909+4449 has an overly wide image separation ∼18″. SDSSJ1452+4224 has bad pixels in CFIS-u image 1. J1612+3920 has a faint counter image hidden to the bottom-left of the lensing galaxy.
In the text

	Fig. C.1. One-dimensional spectra of J0112+3307 from LAMOST and J1559+3839 from SDSS.
In the text