© The Authors 2023

1 Introduction

Measurements of the Hubble-Lemaitre constant (H₀) from the local and the early Universe have shown substantial differences (Riess et al. 2019; Verde et al. 2019), for instance, among the measurements of the cosmic microwave background (CMB; see Bennett et al. 2013; Planck Collaboration VI 2020) and baryon acoustic oscillations (BAO; see Addison et al. 2018; Abbott et al. 2019) as well as those carried out in the more local Universe based on supernovae (SNe; see Dhawan et al. 2018; Macaulay et al. 2019; Dainotti et al. 2021), the tip of the red giant branch (TRGB; see Freedman et al. 2019; Yuan et al. 2019), and Cepheid variables (Pietrzyński et al. 2019; Riess et al. 2019). This so-called ‘Hubble tension’ is considered one of the most significant crises of modern cosmology. Significant efforts have been devoted to solving the problem (Vagnozzi 2020; Adhikari 2022; Goicoechea & Shalyapin 2022; Niedermann & Sloth 2022), but the answer remains inconclusive.

Independent of all of the aforementioned methods, strong lensing time delays provide valuable measurements of H₀ (Birrer et al. 2019; Liao et al. 2019; Wong et al. 2020; Shajib et al. 2020) that may assist in understanding these discrepancies once potential selection bias and unignorable statistical uncertainties in the technique are fully calibrated. Achieving this will require a much larger sample of time delay systems than what is currently available (Shajib et al. 2018) via the gravitationally lensed QSO database¹ (GLQ, Lemon et al. 2019).

In addition to constraining the cosmological model, strong lensing time delay systems, typically multiply-imaged lensed QSOs, provide valuable insights into astrophysical problems such as constraining the distributions of dark and luminous matter of the lenses (Oguri et al. 2014; Suyu et al. 2014; Sonnenfeld & Cautun 2021; Van de Vyvere et al. 2022) as well as uncovering the properties of distant active galactic nuclei (AGN) and their host galaxies to a level of detail that is not possible without lens magnification (e.g., McGreer et al. 2010; More et al. 2015; Fan et al. 2019; Yue et al. 2022). In the case of the latter, microlensing caused by small structures within the lens have revealed fine-level details of AGN morphology such as accretion disk characteristics (Anguita et al. 2008; Motta et al. 2012; Braibant et al. 2014; Fian et al. 2021).

With increasing depth and sky coverage of large-scale surveys, the sample size of multiply imaged QSOs is predicted to grow remarkably. For instance, within third-generation surveys such as Dark Energy Survey (DES, Akhazhanov et al. 2022), KiDS (Kilo-Degree Survey, Kuijken et al. 2019), and Gaia (Lemon et al. 2023), there is a general forecast of about 2000 multiply-imaged QSOs systems (Oguri & Marshall 2010). The current mainstream strategy that has proven effective in seeking these lensed QSOs comprises two steps: (1) finding candidates with high completeness and (2) confirming the candidates with spectra to improve the purity. Previous studies have seen construction of several candidate catalogues, such as those from Agnello et al. (2015), Krone-Martins et al. (2018), Spiniello et al. (2018, 2019), Wu et al. (2022), Akhazhanov et al. (2022), and subsequent spectroscopic follow-ups have provided a number of confirmed lensed QSO samples (see e.g. Lemon et al. 2018, 2019). Sample sizes are set to increase even further with upcoming fourth-generation surveys such as Euclid (Laureijs et al. 2011), the Large Survey of Space and Time (LSST Ivezić et al. 2019), Roman (Eifler et al. 2021), and the Chinese Space Station Telescope (CSST; Cao et al. 2018).

In this paper, we describe our catalogue-based algorithm for finding the candidates of multiply-imaged lensed QSOs from the QSO candidate catalogue of He & Li (2022) extracted from the ninth data release of the Dark Energy Spectroscopic Instrument Legacy Imaging Surveys (DESI-LS; Dey et al. 2019), which covers ~20 000 deg² of the extragalactic sky in the g, r, z-bands. Specifically, we compiled a catalogue containing 971 multiply imaged lensed systems, 620 of which are new candidates not contained in the known lensed QSO databases of Dawes et al. (2022, D22 hereafter) or Lemon et al. (2023, L22 hereafter). Catalogues and images are available at the CDS.

This paper is organised as follows. Section2 introduces the datasets used in this project. The methodology for identifying strongly lensed QSO candidates is detailed in Sect. 3. The results, including candidate catalogues, the images and lens models of selected candidates, and the comparison with known lensed QSOs and other works are presented in Sect. 4. Finally, Sect. 5 presents our discussion and conclusions. In this paper, a fiducial cosmological model Ω_m = 0.26, Ω_DE = 0.74, h = 0.72, w₀ = −1 and w_a = 0 is assumed, matching that adopted by Oguri & Marshall (2010, OM10 hereafter). Unless otherwise stated, all magnitudes quoted in this paper are in the AB system.

Fig. 1 Magnitude distributions of the QCC in g, r, z (left panel) and W1, W2 (right panel). The kernel density estimate curves (solid lines) are generated by kdeplot in seaborn package with binsize = 0.5 and the default Gaussian kernel.

2 Datasets

In this section, we describe the QSO candidate catalogue, which we mined for multiply-imaged QSO candidates (Sect. 2.1, and a reference sample of known lensed QSOs (Sect. 2.2), which we used for optimising the mining process.

2.1 The QSO candidate catalogue

The QSO candidate catalogue (QCC) of He & Li (2022) contains 24 440 816 objects in total. It was created from the point-like sources identified in DESI-LS DR9 using a random forest (RF; Breiman 2001) classification model. The catalogue includes the RA, Dec, the five band magnitudes (g, r, z, W1, and W2), and the probability of being a QSO given by the RF model. The magnitude distributions are shown in Fig. 1. The r-band magnitude of the QCC ranges from 18 to 26, with a mean of 22.44. Evaluated with the testing set that mimics the magnitude and colour distributions of the point-like sources of DESI-LS, the candidate catalogue’s completeness and purity are ~99% and ~25%, respectively.

DESI-LS covers an area of ~20 000 deg² in g, r, z and comprises three different sub-projects: The Dark Energy Camera Legacy Survey (DECaLS², Dey et al. 2019), Beijing-Arizona Sky Survey (BASS³, Zou et al. 2017, 2019), and The Mayall z-band Legacy Survey (MzLS⁴, Dey et al. 2019). The point-source sensitivities of these surveys are as follows: for DECaLS, the 5σ detection limits in AB mag for a point source in individual images are 23.95, 23.54, and 22.50 of g, r, z-bands; for BASS, they are 23.65 (g-band) and 23.08 (r-band); for MzLS, it is 22.60 (z-band). The DESI-LS DR9 catalogue also includes four mid-infrared bands at 3.4, 4.6, 12, and 22 μm (corresponding to W1, W2, W3, and W4, respectively) observed by the Wide-field Infrared Survey Explorer⁵ (WISE, Wright et al. 2010).

2.2 Reference sample

We adopted a reference sample containing 57 spectroscopically confirmed multiply imaged QSO systems to optimise the candidate selection process, such as fine-tuning selection thresholds and training the inspectors. The reference sample (referred to as RLQ hereafter) was acquired by cross-matching the GLQ and QCC. There are 111 objects (~50.5% of all objects in GLQ, composed of 93 pairs + 18 quads) having at least one matched member, and 57 of these 111 objects have at least two matched members in the QCC. The presence of multiple matched members is an indication that our group finder can successfully detect all of these multiply imaged QSO systems (see Sect. 3.1 for more details).

Our reference sample comprises only those objects with two or more matches. The remaining excluded 54 systems with only one match arising due to three possible situations: (1) in a given multiply-imaged QSO system, only one of the lensed images is labelled as the point-like source in the DESI-LS DR9 catalogue and the other lensed images are labelled as extended sources. Hence, the QCC did not include these multiple images in the first place (see the left panel of Fig. 2); (2) as shown in the middle panel of Fig. 2, not all of the multiple images of lensed QSOs are detected in all five-bands and, thus, they are not included in the QCC; (3) the multiple images are smeared by the imaging point spread function (PSF) and therefore have been labelled as a single point-like source in the DESI-LS DR9 catalogue and in the QCC as well (see the right panel of Fig. 2).

The catalogue of RLQ holds both the information provided in GLQ (including RA, Dec, Name, lens and source redshifts, and the number of images) and the labels and scores attributed in this work, such as the internal IDs given by the group finder, S_colour and Srf (defined in Sects. 3.2.1 and 3.2.2), and labels indicating whether the system exists in the sample for visual inspection (Sect. 3.4). The detailed RLQ catalogue and the corresponding description are published at the CDS.

Fig. 2 Demonstrating the three possible situations that cause a known multiply-imaged QSO to have only one match in the QCC. The colours of the cycles represent the morphology types in DESI-LS catalogues. Green circles represent point-like sources, whereas white circles represent extended sources. For WG0214-2105, three out of its four images are recorded in the catalogue, but only one is labelled as point-like. For Β1608+656, the image with the green pentagon marker lacks a g-band detection. In the case of SDSSJ0924+0219, three of the lensed images have been identified as a single point-like source and the fourth image has been identified as an extended object.

3 Methodology

Our procedure for mining multiply-imaged lensed QSOs in the QCC is broadly divided into three phases. Section 3.1 describes the first phase, where we use a grid-search algorithm on the angular positions of the members in the QCC to identify goups of candidate QSOs. The second phase cleans these groups depending on whether they contain two (Sect. 3.2) or more (Sect. 3.3) members. The third phase applies visual inspection to finalise the list of multiply imaged candidates (Sect. 3.4).

3.1 QSO group finder

Our QSO group finder selects QSO candidate groups within the QCC based on their angular positions. This comprises four steps:

Firstly, all sources in the QCC are put into grids generated by HEALPix (Gorski et al. 1999) with nside = 2¹⁵, leading to a grid bin size of 6.4 × 6.4 arcsec². The 24 440 816 sources in the QCC occupy 23 394 395 HEALPix grids. Secondly, the programme begins with grid bins that contain at least one QSO candidate. For each of these bins, the surrounding eight grid bins are checked for QSO candidates. Surrounding bins that contain QSOs are then connected to the central bin. After that, the surrounding bins of new groups are checked; if there are QSO candidates in the neighbouring pixels of all the members of the group, the new candidates are further added to the group. This process is then repeated until there are no more candidates found in the surrounding bins of the final group. Next, Once a group of bins is obtained in this way, the QSO candidates in the group are labelled as belonging to that QSO candidate group. Finally, for a QSO candidate group holding two members only, their angular distance must be greater than 0.5 but less than 10 arcsec; otherwise, the group is removed.

As a result, 562 206 QSO candidate groups were identified and placed into a candidate QSO group catalogue (QGC). As expected, the QGC includes all groups contained in the RLQ. The QGC provides the RA, Dec, the probability of being a QSO given by He & Li (2022), and the five-band magnitudes g, r, z, W1, W2 for each candidate of every group. We note that the QGC also gives the quantities S_colour and S_RF, which we define below.)

We divided the QGC into two parts, one part containing groups with two members and the other containing groups with more than two members. We labelled these PAIR and MUL, respectively. PAIR includes 462 608 systems, while MUL includes 99 538 systems. We also labelled those groups found in the RLQ as either RLQ-PAIR (47 systems) or RLQ-MUL (10 systems) depending on whether they belong to PAIR or MUL. The QGC and its corresponding description are published at the CDS.

3.2 QSO candidate groups in PAIR

PAIRs are approximately five times more prevalent in the QGC than MULs. We discarded all PAIRs having broad colour differences (see below) and/or low probabilities of their members being QSO candidates. As we discuss in Sect. 3.3, we implemented a different cleaning strategy for the MULs.

3.2.1 Selection with colour similarity

Since gravitational lensing is an achromatic process, differences in colour between multiple images can only arise through a range of scattering processes along the different paths between the observer and source. Unless a lens has strong dust gradients, the colours of lensed images are therefore generally very similar within a given multiply imaged system. As such, we can use colour similarity to eliminate PAIRs that are not likely double-image systems. Hence, we define the colour similarity of a group in PAIR, S_colour, as follows:

$$S_{\text{colour}}=\{\begin{array}{ll}1-\frac{1}{10}{\displaystyle {\sum }_{i=1}^{10}{\sigma }_i}\hfill & \text{if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\frac{1}{10}{\displaystyle {\sum }_{i=1}^{10}{\sigma }_i}<1\hfill \\ 0\hfill & \text{if\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\frac{1}{10}{\displaystyle {\sum }_{i=1}^{10}{\sigma }_i}\ge 1\hfill \end{array},$$(1)

where σ_i is the standard deviation of the ith colour (computed over both members of the PAIR) out of the ten unique colours provided by the g, r, z, W1, and W2 magnitudes in the QCC.

All PAIRs with a value of S_colour that was lower than a given threshold were rejected. We chose °.5 as the threshold in this work as a compromise between rejecting as many systems in the QGC as possible whilst minimising the rejection of confirmed RLQ-PAIR systems. This choice of threshold rejects ~36% of systems in the QGC while only losing ~4% of confirmed systems classified as RLQ-PAIR. We label the remaining 297 502 groups of QSO PAIR candidates as PAIR-CS. Two known lensed QSOs are missing: RXJ0911+0551 and SBS1520+530. These systems exhibit significantly different colours between their images, presumably due to very different levels of dust encountered along the path to each image. In particular, these systems have significant differences between images in the colours that include the W1 and W2 magnitudes. However, if we do not include these magnitudes, the rejection rate drops dramatically, which leads to a much higher workload for human inspection. We discuss this further in Sect. 5.

3.2.2 Selection with the probabilities of being QSOs

We used an additional metric, S_RF, that measures the overall probability of the candidates truly being QSOs, to remove the groups with possible false QSOs. For a given QSO candidate group, this score is defined as:

$$S_{\text{RF}}=\frac{1}{n}{\displaystyle \sum _{i=1}^{n}pro{b}_i}$$(2)

where n is the number of QSO candidates in the group. The quantity prob_i is provided in the QCC and gives the probability of a candidate being a true QSO. The score S_RF therefore represents the likelihood that a candidate group is a QSO group.

In this work, we choose a threshold of S_RF = 0.85 and reject any groups in PAIR-CS with a value of S_RF less than this. This removes two more confirmed multiply imaged QSOs (SBS0909+532 and DESJ0405-3308), leaving 43 out of the initial 47 confirmed systems in RLQ-PAIR. Overall, 102 468 (~22.2% PAIR, ~34.4% PAIR-CS) groups in PAIR are selected and labelled as PAIR-CS-RF.

Since we provide both S_RF and S_colour in our online catalogue, users of the data can apply different thresholds according to whether a higher recall rate or a better precision is required.

3.3 QSO candidate groups in MUL

Our procedure for cleaning PAIRs is not appropriate for candidate groups classified as MUL since it is overly sensitive to groups that contain additional contaminating members. As such, we define a different strategy for the MUL groups in the QGC.

First, for each group, members with S_RF < 0.85 were removed. Then, groups containing only one member after the above process are also discarded. The remaining 52 582 groups (~52.8% MUL) were labelled MUL-RF. Collectively, groups within MUL-RF contain a total of 134 976 QSO candidate members. Secondly, we explored the colour similarity of the members in each group to quantify the possibility of each group being a multiply imaged QSO system. The details of this step are as follows: (1) for a given group, S_colour, was calculated for all possible combinations of group members, from those containing two members to the largest combination that contains all members; (2) all combinations with S_colour < 0.5 were rejected; (3) finally, the combination containing the most members (and that with the highest S_colour if there is more than one of these) was retained and the others were rejected.

As a result, 45 905 groups (containing a total of 111 761 members) remain (~46.1% MUL, ~87.3% MUL-RF) and are labelled MUL-RF-CS. In MUL-RF-CS, all ten of the confirmed systems in RLQ-MUL are retained. Together, the MUL-RF-CS and PAIR-CS-RF contain 53 out of 57 RLQs; that is, the recovery rate is ~93%.

3.4 Visual inspection

The final phase of our mining procedure implements human visual inspection (VI) to improve the true positive rate and to grade the multiply imaged QSOs candidates. All inspectors were initially trained by visually analysing the images of the 53 confirmed lensed QSOs in PAIR-CS-RF and MUL-RF-CS. Then, all candidates in PAIR-CS-RF and MUL-RF-CS are inspected, yielding a total of 971 candidates. These 971 candidates are then graded into one of three grades, A, B, or C by two inspectors independently (the first and fourth author). An A grade is awarded to the most reliable systems with strong lensing features similar to the confirmed systems; B grades are awarded to systems with features that are less similar to the confirmed systems, such as a larger image separation. To help the inspectors improve their objectivity, the statistics, S_RF and S_colour, are referred to. Grade-C systems typically have a few lensing features; if they are lensed QSOs, their configurations are atypical.

The grading process is naturally subjective, so the grades of the two inspectors are listed simultaneously in our online catalogue to give an indication of human bias. For simplicity, the grading results presented in this paper are solely those given by the first author unless otherwise stated. The following features were considered during the grading: (1) the higher S_colour and S_RF, the better; 2) the existence of an apparent lens galaxy makes the system more plausible, especially if it has a red colour; 3) separations between images on opposite sides of the lens should usually be less than 3 arcsec and a [0.5, 1.5] arcsec separation makes a candidate more convincing; 4) for the pairs with apparent lens light, the angle between the two position vectors (measured from the lens centre to the image) is greater than 120 degrees (Chan et al. 2015). The candidates that do not satisfy this criterion are rejected or graded as C; 5) for the pairs without apparent lens light, the S_colour and S_RF influence the grading more strongly; 6) for quads, the configuration of the multiple images should be similar to those of the systems in the RLQ.

The resulting visually inspected catalogues are labelled PAIR-CS-RF-VI and MUL-RF-CS-VI. We combined the 971 candidates contained in total by these two catalogues into a single catalogue that we refer to as H22 hereafter. This catalogue is available at the CDS.

4 Results

In this section, we present the statistical properties of the candidates, comparing our final H22 candidate catalogue with existing, known multiply-imaged QSOs and two candidate catalogues determined by other studies.

4.1 Lensed QSO candidates in this work

Table 1 lists the properties given for every candidate lensed QSO system in our H22 catalogue. Included in this list of properties is an internal ID linking to the QGC, the RA and Dec of the candidate, the awarded grade, the evaluation metrics (S_RF, S_colour) and physical parameters such as image separation and number of images. Systems that have been previously modelled also contain the axis ratio, position angle and Einstein radius lens model parameters. We also include the labels ‘in_L22’ and ‘in_D22’ that indicate overlap with the candidate catalogues of Lemon et al. (2023) and Dawes et al. (2022) respectively. Where an overlap occurs, we give the grade/classification awarded by that catalogue in the properties ‘classification_by_L22’ and ‘grade_by_D22’.

We also provide an additional catalogue which we refer to as H22-details, which includes all images belonging to each candidate lensed QSO system in H22⁶). Table 2 lists the properties given for each image of each candidate lensed QSO system. These include the internal ID that links to the QGC (different images of one candidate share the same ID), RA, Dec, g, r, z, W1, W2, and the source redshift (if applicable) from the Sloan Digital Sky Survey (SDSS, Blanton et al. 2017) extended Baryon Oscillation Spectroscopic Survey (eBOSS, Dawson et al. 2016) DR16. The RA and Dec here are the coordinates of the images in candidate systems.

The distributions of S_colour and S_RF are plotted in Fig. 3 for H22, the QGC, PAIR-CS-RF, MUL-RF-CS, and the catalogue of confirmed systems, RLQ. A general observation is that the distributions of PAIR-CS-RF, MUL-RF-CS and H22 match those of the RLQ catalogue significantly better than the QGC distributions, although this is to be expected given the selection criteria we have applied based on S_colour and S_RF. In addition, the difference between the RLQ and QGC catalogues is much more pronounced for S_RF, indicating that this is a more efficient statistic to base the selection of lensed QSOs. It is also apparent from Fig. 3 that H22 more closely matches the distributions of RLQ than PAIR-CS-RF or MUL-RF-CS do, demonstrating that the final step of human visual inspection has made a significant improvement. However, as is mentioned in Sect. 3, our catalogue-based approach may miss the candidates that do not exist in the parent samples (the QCC), if the various selection criteria for compiling RLQs are considered, then the difference between H22 and RLQ is inevitable. This difference can be reduced when spectroscopic follow-ups remove the false positives in H22.

Figure 4 displays some example grade-A systems. We show some examples of systems also found by D22, some confirmed systems in RLQ, and, finally, some examples of new systems not found in D22 or RLQ.

In Fig. 5, we compare the distributions of S_colour and S_RF in H22 with those of D22, L22 and the ‘rediscovered’ systems in RLQ that are found in H22. These are all statistically similar. The figure also shows the distributions of S_colour and S_RF split by grade in H22. As expected, the grade-A candidates match RLQ best, while grade-Bs and Cs have similar distributions and differ more from RLQ.

The redshift distribution and g-band magnitude distribution are shown in Figs. 6 and 7, respectively, and compared with OM10, L22, and D22. The source redshifts are acquired from the QSO catalogue of eBOSS DR16 (Lyke et al. 2020, DR16Q hereafter). Altogether, 195 out of 971 entries of H22 have at least one match in DR16Q and, thus, the corresponding source redshifts are available. In source redshift space, our samples are mostly distributed within [0,3.5]. Our method tends to select the samples at lower redshifts (<2.0), missing high-redshift lensed QSOs because of the depth limit of DESI-LS and WISE. The g-band magnitude distribution reveals that the luminosity distribution of the QSOs in H22 deviates from the power-law trend given by OM10 around g ≃ 20.5, indicating again that our method becomes insufficient beyond the depth limits of the observations (mainly limited by WISE data).

Moreover, among the 195 candidates with redshift detections, 158 candidates have one redshift detection, and 37 (23 As + 5 Bs + 9 Cs) systems have two. The redshift overlaps between DR16Q and H22 give us some clues about the false positive rates of H22. Among 23 grade-As, 21 have similar redshifts (difference <0.02), but for 5 grade-Bs, the result is 3 and for 9 grade-Cs, the result is 5. Assuming that the candidates with two similar redshift detections are more likely to be the true lensed QSOs, this suggests that the false positive rate is increasing from A to B to C, although the trend is of course subject to a relatively large Poisson error.

Distributions of image separations plotted in Fig. 8 show that there are significant differences between grades. The image separation distribution of grade-As peaks at ~0.8 arcsec, while those of grade B and C peak at ~2.2 arcsec. The dissimilarity is due to the criteria applied during human visual inspection (Sect. 3.4) which assigns a lower grade to larger image separations. The distribution of grade A image separations agrees with OM10, but those of grade B and C do not. Contamination by stars is more likely at larger image separations and therefore the false positive rate in the grade A lenses is expected to be lower.

Fig. 3 Distributions of Scolour (upper panel) and SRF (lower panel) for systems in QGC, PAIR-CS-RF, MUL-RF-CS, H22 (all candidates), and RLQ.

Fig. 4 Examples of grade A candidates in H22. The top row shows candidates also found in RLQ, the middle row shows candidates also found in the candidate catalogue of D22 and the bottom row shows new systems neither contained in RLQ nor D22. Each image is titled with the internal ID (linking the system to QGC) and gives the values of SRF, Scolour, and image separation.

4.2 New lensed QSO candidates

Our catalogue of candidate multiply imaged lensed QSOs contains 620 new systems not contained in the RLQ catalogue, L22 or D22. In this section, we discuss the key differences between our new catalogue, H22, and the RLQ catalogue and those of D22 and L22.

4.2.1 Compared to RLQ

Compared to RLQ, we have identified 918 extra candidates out of 971 and labelled them as H22-new-RLQ hereafter (see Fig. 4).

We note that some of these 918 candidates are also found in D22and L22; we discuss these overlaps in Sects. 4.2.2 and 4.2.3.

Comparisons of the distributions of S_colour and S_RF between H22 and RLQ are shown in Fig. 5. At the upper end of S_colour, H22 shows a more significant difference to RLQ, which indicates that there are items with low S_colour in H22. This is likely due to false positives which spectroscopic follow-up would reject. In contrast, the S_RF statistic of RLQ shows a much more similar distribution to our catalogue, implying that S_RF is a more distinguishing criterion in the identification process.

The redshift and magnitude distributions of strongly lensed images in the candidate systems are displayed in Figs. 6 and 7. We find that H22 covers broader redshift and magnitude ranges than RLQ, for example, H22’s redshift range beyond 1.8 and magnitude range beyond g ≃ 19. Comparing the max image separation distributions between H22 and RLQ (Fig. 8), we find that H22 has a wider distribution than RLQ. Notably, the lower end is filled with grade-A candidates (whose distribution matches the theoretical predictions of OM10), while grade-B and grade-C candidates dominate the upper end (whose distributions differ from OM10).

The majority of systems in H22-new-RLQ are pairs. A few possible quads (110720, 2484493, 11419327, and 11125158) exist. Thus, the quads-to-pair ratio (~1/30) is unexpectedly lower than the theoretical prediction of ~1/6 given by OM10. Multiple causes may lead to such an issue. One possibility is that the quadruple systems are more prone to the effects of PSF smearing than duals; when the max image separations are comparable to the PSF size, it is harder to distinguish quads than it is to distinguish duals, which likely causes quads to be under-respresented in the QCC (see Fig. 2 and the corresponding discussion). It is likely that adopting image-based deep learning approaches similar to those developed for searching galaxy-galaxy strong lensing systems (e.g., Petrillo et al. 2017; Lanusse et al. 2018; He et al. 2020; Li et al. 2021; Huang et al. 2021; Rojas et al. 2022) would help retrieve these missing quads. Another possibility is that a low quads-to-duals ratio could be caused by the dual candidates having a higher false positive rate than the quad candidates. In principle, the probability of the emergence of two-image systems of non-lenses is significantly higher than that of four-image systems, especially when a requirement of visual inspection is that the quads must have typical lensed-image configurations. Hence, our candidate catalogue is unsuitable for estimating the double-to-quad ratio but reveals the issues that need to be improved.

Fig. 5 Distributions of Scolour (top two panels) and SRF (bottom two panels) comparing our new candidate systems (H22) with those of D22, L22, and RLQ (left panels) and splitting our new candidates by awarded grade (right panels).

Fig. 6 Comparison of redshift distributions of H22, D22, L22, and OM10. We note that the redshifts of L22are from Lemon et al. (2023), while those of H22 and D22are from SDSS eBOSS DR16.

Fig. 7 Comparison of magnitude distributions of H22, D22, and OM10 shown in the left panel. Since L22did not provide g-band magnitudes, it is not plotted here. Distributions split by awarded grade in H22 are shown in the right panel.

Fig. 8 Comparison of max image separation of systems in H22, D22, L22, RLQ, and OM10. The panel on the right splits this distribution for H22 by awarded grade.

4.2.2 Compared to D22

D22⁷ is another catalogue of lensed QSO candidates extracted from DESI Legacy Imaging Surveys, created by Dawes et al. (2022). This study uses an independent approach to ours and contains 436 candidates. To understand the differences between the selection effects of H22 and D22, we compare distributions of S_colour, S_RF (Fig. 5), source redshift (Fig. 6), g-band magnitude (Fig. 7), and max separation (Fig. 8).

There are 292 matching systems between H22 and D22(which can be extracted by requiring ‘in_D22=True’ in H22). According to our grading, 133 are grade-A, 62 are grade-B, and 97 are grade-C. However, according to D22’s grading, 76, 83, and 133 out of 292 matches are grade-Α, B, and C respectively. We include the grade labels of D22in H22. Given that the grades are assigned by human inspectors independently, the above difference are a reflection of the subjective nature of this process.

The 679 systems in H22 not found in D22are labelled as H22-new-D22. These contain 151 grade-As, 221 grade-Bs, and 307 grade-Cs. However, there are 144 candidates proposed in D22but missed by H22. Specifically, 43 systems are rejected by our selection by S_colour and S_RF; S_colour causes 38 rejections while S_RF causes 5. The other 100 systems are not present in the QGC from the start. To summarise, 100 out of 144 mismatches are caused by different parent samples, 43 are caused by different selection methodologies and only one is caused by human inspection.

The comparisons of the distributions of S_colour and S_RF between H22 and D22in are shown in Fig. 5. We find that the distributions of D22are closer to RLQ than H22. If we select only the grade-Α systems from H22, the distributions become similar to D22, suggesting that D22has a slightly higher purity than H22. Figures 6 and 7 demonstrate that H22 is deeper than D22in both magnitude and redshift, mostly due to the different parent samples adopted in the two works. D22used the ‘DESI QSO Sample (Yèche et al. 2020)’ which has a hard magnitude cut of r = 22.7 for the QSO detections. In contrast, QCC does not have any cuts in magnitude. To provide some quantification of the effect of this, there are ~52% candidates in QCC that are fainter than r = 22.7. Regarding the distributions of image separation (Fig. 8), D22’s separation is smaller than H22’s but larger than the grade-A systems in H22. This indicates that H22-new-D22(especially grade-B and grade-C systems in H22) covers multiply imaged lensed QSO systems with large image separations missed by D22and thus potentially improves the completeness of the sample of multiple image QSOs systems from DESI Legacy Imaging Surveys.

4.2.3 Compared to L22

L22 is a catalogue that contains the spectroscopic follow-up results of 175 systems selected from a multiply imaged QSO candidate catalogue given by Lemon et al. (2019) based on Gaia DR2, which is an important update to GLQ. It takes advantage of the high astrometric precision of Gaia to discover many new lensed QSOs, including confirmed lensed QSOs, nearly identical QSOs (NIQ), projected QSOs, and QSO pairs. Following the definitions in Lemon et al. (2023), the two QSOs of an NIQ system have similar spectra but cannot be confirmed as a strong lensing system because of the absence of a lens galaxy image. Nevertheless, the NIQs should be considered very promising candidates. ‘Projected QSO’s’ occur when the redshifts of two QSOs are different, while ‘QSO pairs’ are cases where the QSOs have similar redshifts but different spectra. For lensed QSO searching, QSO pairs and projected QSOs are contaminations.

There are 44 matches between L22and H22. The 948 systems in H22 not found in L22are labelled as H22-new-L22. In the 44 matches, 17 systems are labelled as lenses in L22, 20 are NIQs, 2 are projected QSOs, and 3 are QSO pairs. In other words, 5 systems are false positives, and 39 are confirmed lenses or promising candidates. The comparisons of the distributions of S_colour and S_RF between H22 and L22are shown in Fig. 5. The figure shows that compared to H22 and D22, L22has the most similar distributions to RLQ. This is unsurprising because most of the systems in L22are spectroscopically confirmed lenses and NIQ.

Figure 7 shows that L22is shallower than H22 and D22(i.e., lower mean redshift). This is due to their selection from shallower Gaia data (Gaia Collaboration 2018). The image separations of L22are similar to those of RLQ, which again is a reflection of the fact that L22primarily consists of spectroscopically confirmed lensed QSO systems (Fig. 8).

4.2.4 New candidates compared to RLQ, D22, and L22

We have identified 620 new candidates not found in any of RLQ, D22, or L22. We name the catalogue of these new candidates H22-new, which can be achieved by requiring ‘in_RLQ’=False, ‘in_D22’=False, and ‘in_L22’=False simultaneously in the online catalogue. In Fig. 9, we compare the distributions of S_colour, S_RF, g-band magnitude, and image separation between H22-new and the new systems found in H22 but not RLQ (H22-new-RLQ), D22(H22-new-D22), or L22(H22-new-L22).

Considering the distribution of S_colour, H22-new candidates have lower scores than the other samples, because the high-score systems are more likely to be rediscovered in other datasets, that is, RLQ, D22, and L22. Similar trends are also visible in the distribution of S_RF. In terms of the g-band magnitude, H22-new candidates occupy a fainter region than the other samples; there are more matches with the other datasets at brighter fluxes due to higher signal-to-noise ratios. Regarding image separations, H22-new candidates have larger separations on average compared to the other samples, meaning that the systems with small image separations are more likely to be rediscovered in the other datasets.

In summary, H22-new is a sample of multiply imaged QSO candidates not found in any of RLQ, D22, or L22. Although some candidates with high confidence (included in D22and L22) are excluded in H22-new, the 101 grade-As in H22-new are still valuable candidates worth spectroscopic follow-up. H22-new also includes 214 grade-Bs and 305 grade-Cs. These samples contain more false positives than in grade-As, but mining multiply imaged QSOs systems from them can still enhance the sample size of multiply imaged QSOs. Hence, spectroscopic confirmations are also worth pursuing for grade-Bs and Cs, albeit at a lower priority.

4.3 Basic lens modelling of selected candidates

Lens modelling provides a useful tool for enhancing the grading of candidate lensed QSO systems. A statistical goodness-of-fit of a lens model can be used to improve the confidence of a candidate being a confirmed lensed system. Thus, we performed a basic lens modelling of the candidates with visible lens light in H22-new-RLQ to validate our identification algorithm. As such, we selected 52 candidates for modelling. To simplify the modelling, we choose the candidates with visible light centres and assume the mass centres align with the light centres. This reduces the number of model parameters. For example, doubly imaged QSOs have six independent data points, including the positions and fluxes of the two images, however, if we do not fix the centre of mass, there are eight free model parameters for the singular isothermal ellipsoidal (SIE, Kormann et al. 1994) model, including five lens parameters and three source parameters.

The mass model adopts the SIE profile, which has been widely applied to describe the mass distribution of galaxies (e.g., Bolton et al. 2012; Sonnenfeld et al. 2013). The SIE profile requires the five parameters (x_l, y_l, q, phi, rein), which are the position coordinates, the axis ratio, the position angle, and the Einstein radius, respectively. We fix the mass centre, x_l, y_l, of each lens according to the peak of its light distribution as mentioned above. Three parameters are used to model the quasar source: two for the unlensed position (r_s) and one for the flux (f_s) on the source plane. We denote the observed position and flux of the i–th image of a lensed QSO by r_i and f_i, respectively. In our modelling, we fit to both the image fluxes and the image positions by varying the lens and source parameters. Quantitatively, we vary the model positions ($r_i^P$) and fluxes ($M_i^P{f}_s$) of the i-th image until the best fit with the corresponding observed parameters is found, that is, when $r_i^P\approx r_i$ and $M_i^P{f}_{\text{s}}\approx f_i$. Here, $M_i^P$ is the model magnification factor at the i–th lensed image. In practice, the above objective is achieved by using a simulated annealing algorithm that minimises the following penalty function,

$$\begin{array}{l}{\chi }^2={\chi }_{\text{position}}^2+{\chi }_{\text{flux}}^2\hfill \\ ={\displaystyle \sum _i\frac{{\left|r_i-r_i^P\right|}^2}{{\sigma }_i^2}+{\displaystyle \sum _i\frac{{\left(f_i-M_i^P{f}_s\right)}^2}{{\sigma }_{f,i}^2}}}\hfill \end{array},$$(3)

where σ_i and σ_f,i are the position and flux measurement uncertainties at the ith image, respectively. The calculations of $r_i^P$ and $M_i^P$ are performed with the open-source software lenstronomy (Birrer et al. 2015, 2021; Birrer & Amara 2018), and the python code used for the lensing analysis in this work is publicly available⁸. To evaluate the goodness of fit, we adopt the Bayesian information criterion (BIC, Liddle 2007), given by

$$\text{BIC=ln}\left(n\right)k-2\text{\hspace{0.17em}ln}\left(\widehat{L}\right),$$(4)

where $\widehat{L}=\exp \left(-{\chi }_{\text{minimised}}^2\right)$; n is the number of input data points and k is the number of free parameters of the model. For pairs, n = k = 6, while for the two quads (with internal IDs of 110720 and 2484493), n = 9, k = 6, since our approach missed one of the quadruply lensed images as shown in the first two panels of Fig. 4. A lower BIC indicates a better agreement between the model and the observations.

Consequently, we determine lens parameters (q,phi,rein) of 52 systems, including 50 pairs and 2 quads. The distribution of BICs of the above models is shown in Fig. 10. Examples of three different systems with different BICs are also plotted in Fig. 10 to demonstrate the agreement between the best-fit models and observations. There is no significant disagreement between models and observations, even for the two quad systems with a larger BIC. The mean BIC of grade-A, B, and C two-image samples are 14.94, 18.47, and 135.12, respectively, which suggests that the outcomes of our human grading procedure are in good agreement with the modelling.

We also present three sample systems in the bottom panels in Fig. 11, including one double-image system and two quadruple-image systems. Notably, the system with internal ID 110720 (also reported by Huang et al. 2021, as DESI-055.7976-28.4777), which has the largest BIC, seems to be a triple-image system, but the best-fit model presents a quad system that leads to the “worst” fitting case in this work. The primary reason is that the parent catalogue probably misses the fourth image due to the brightness limit or morphology classification. Hence, the quick lens modelling process can decrease the candidate sample’s false positive rates and improve the integrity of individual systems.

Fig. 9 Comparisons between H22-new-RLQ, H22-new-D22, H22-new-L22, and H22-new of distributions of Scolour (upper left), SRF (upper right), g-band magnitude (lower left), and image separation (lower right).

5 Discussion and conclusions

In this work, we have designed a catalogue-based approach for mining multiply imaged lensed QSO candidates from large catalogue datasets. The approach adopts three stages: (1) identifying groups of candidate QSOs; (2) refining the groups using a statistical measure of the similarity of group member colours and the likelihood of members being QSOs, according to He & Li (2022), which includes both the removal of entire groups and removal of individual members for groups with more than two members; (3) employing human inspection to grade candidate systems and further reject unlikely systems.

We applied the approach to the catalogue of He & Li (2022) which comprises 24 440 816 QSO candidates drawn from the DESI-LS DR9 source catalogue. Initially, 562 206 QSO candidate groups were identified, further refined by colour similarity and QSO likelihood into a set of 102 468 candidate groups containing two members and a set of 45 905 candidate groups containing more than two members. Subsequent human inspection resulted in a final total of 971 candidate multiply imaged QSOs. Of these, 620 candidates are new in the sense that they have not been identified in existing studies. Relevant images and catalogues of our candidate systems are available at the CDS.

The redshift range of the 971 candidate lensed QSOs is (0, 3.5], peaking at ~2. Their g-band apparent magnitudes span the range 17 ≤ g ≤ 24, peaking at g ~ 20.5. We estimated the recovery rate of our mining from the 57 known lensed QSOs contained in the catalogue of He & Li (2022); our catalogue of 971 contains 53 of these, indicating a recovery rate of approximately 93%. We note that the use of these 57 known systems in optimising our mining strategy likely results in selection bias towards discovering similar systems. Therefore, this recovery rate only approximately represents the true recovery rate, which will be lower in practice.

Our candidate catalogue includes grades awarded by two independent human inspectors, according to how likely it is that the candidates are actually lensed QSO systems. This grading is based on the visual inspection of images, relying mainly on image configurations but also turning to colour similarity and QSO likelihood in unclear cases. The grading process is subjective, and so the inclusion of both awarded grades in our catalogue gives an indication of subjectivity. There are 284 grade-As, 283 grade-Bs, and 404 grade-Cs. Grade-As show high similarity to the distribution of separations of the multiply imaged QSOs systems in OM10, while grade-Bs and grade-Cs generally have larger separations. We also implement lens modelling on selected candidates to validate the candidates, and the results show reasonable outcomes, especially for grade-A systems. Hence, we consider the grade-A systems as prioritised targets for follow-up spectroscopic campaigns.

Compared to the studies of D22and L22, ours adopts a catalogue-based approach and starts with a larger parent sample. In particular, our parent sample is deeper than those of D22and L22. D22utilises a relatively shallow parent sample (r < 11.7) and about 1/3 of D22candidates have been identified in our catalogue of 971 candidates. The remaining ~1/3 is caused by different parent samples (contributing about 68% difference) and selection methodologies (contributing about a 31% difference). Our method also applies colour similarity thresholding, which improves purity but sacrifices completeness; for example, two out of the 57 known systems are rejected after our selection based on colour similarity. There are 44 systems in common between H22 and L22. Among them, 39 candidates are confirmed lenses or promising candidates and 5 are false positives.

Since our method directly applies to catalogue data, it is heavily dependent on the source extraction algorithm used which may not be optimised for finding strongly lensed QSO systems. Hence, we are generally biased against finding small image separation systems and quadruply imaged systems. This is true of the DESI-LS catalogue mined in this work. The PSF of DESI-LS results in about 64.0% of quads (or 41.5% of pairs) being missed when creating the initial list of QSO candidates because the quads are more likely to be smeared to one point source. Although we ensured a high recovery rate of known lenses, the selection by colour and QSO likelihood results in some lensed QSO being rejected, as shown in comparison with D22. Due to this process, 26 grade-A candidates (according to D22’s grading) were missed. Besides, mis-labelled objects in the DESI-LS catalogue lead to false positives, for instance, the candidate with internal ID 110720 in H22 is actually a galaxy-galaxy strong lensing system (Lemon et al. 2023; Schmidt et al. 2023). As such, in future work, we plan to develop a methodology that combines catalogue-based and image-based approaches to take advantage of both whilst avoiding their disadvantages.

To summarise, our work provides the largest catalogue of multiply imaged lensed QSO candidates to date, comprising 620 new lensed QSO candidates, of which over 100 are high-grade. With future large-scale spectroscopic follow-up of these from, for example, DESI and the 4-metre Multi-Object Spectrograph Telescope (de Jong et al. 2019), existing samples of lensed QSOs used for cosmological and astrophysical studies could therefore be greatly increased for much-improved statistical power.

Fig. 10 Distribution of BICs of the 52 lensing models. BICs of two-image systems are represented by histograms, while vertical orange lines indicate the BICs of the two quad systems. In the square panels, blue circles show the observed image positions, orange circles show the modelled image positions, and green diamonds show the source positions given by the lens models. The green and red curves indicate the caustic and critical curves of the corresponding lensing models.

Fig. 11 Modelled caustic curves for the three candidates in the third line of Fig. 4. The plots are centred at the lens light centre, which is determined by peak of the lens light. Symbols are the same as in Fig. 10.

Acknowledgements

We thank China-VO for providing DESI-LS data. We thank Huanyuan Shan, Xinzhong Er, Dongxu Zhang, Yunfei Chen for the extensive discussions. We thank astropy, healpix, pandas, seaborn, lenstronomy for providing convenient and reliable python packages. We acknowledge the support from the Ministry of Science and Technology of China (Nos. 2020SKA0110100); the science research grants from the China Manned Space Project (no. CMS-CSST-2021-A01); CAS Project for Young Scientists in Basic Research (no. YSBR-062). H.Z. acknowledges the supports by the National Natural Science Foundation of China (NSFC, Grant no. 12120101003), Beijing Municipal Natural Science Foundation under grant 1222028 and the science research grants from the China Manned Space Project with no. CMS-CSST-2021-A02. This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the DES collaboration. Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cientifico e Tecnologico and the Ministerio da Ciencia, Tecnologia e Inovacao, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energeticas, Medioambientales y Tecnologicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciencies de l’Espai (IEEC/CSIC), the Institut de Fisica d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig Maximilians Universitat Munchen and the associated Excellence Cluster Universe, the University of Michigan, NSF’s NOIRLab, the University of Nottingham, the Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University.

References

• Abbott, T. M. C., Abdalla, F. B., Alarcon, A., et al. 2019, MNRAS, 483, 4866
• Addison, G. E., Watts, D. J., Bennett, C. L., et al. 2018, ApJ, 853, 119
• Adhikari, S. 2022, Phys. Dark Universe, 36, 101005
• Agnello, A., Kelly, B. C., Treu, T., & Marshall, P. J. 2015, MNRAS, 448, 1446
• Akhazhanov, A., More, A., Amini, A., et al. 2022, MNRAS, 513, 2407
• Anguita, T., Schmidt, R. W., Turner, E. L., et al. 2008, A&A, 480, 327
• Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20
• Birrer, S., & Amara, A. 2018, Phys. Dark Universe, 22, 189
• Birrer, S., Amara, A., & Refregier, A. 2015, ApJ, 813, 102
• Birrer, S., Treu, T., Rusu, C. E., et al. 2019, MNRAS, 484, 4726
• Birrer, S., Shajib, A., Gilman, D., et al. 2021, J. Open Source Softw., 6, 3283
• Blanton, M. R., Bershady, M. A., Abolfathi, B., et al. 2017, AJ, 154, 28
• Bolton, A. S., Brownstein, J. R., Kochanek, C. S., et al. 2012, ApJ, 757, 82
• Braibant, L., Hutsemékers, D., Sluse, D., Anguita, T., & García-Vergara, C. J. 2014, A&A, 565, A11
• Breiman, L. 2001, Mach. Learn., 45, 5
• Cao, Y., Gong, Y., Meng, X.-M., et al. 2018, MNRAS, 480, 2178
• Chan, J. H. H., Suyu, S. H., Chiueh, T., et al. 2015, ApJ, 807, 138
• Dainotti, M. G., De Simone, B., Schiavone, T., et al. 2021, ApJ, 912, 150
• Dawes, C., Storfer, C., Huang, X., et al. 2022, ApJ, submitted [arXiv:2208.06356]
• Dawson, K. S., Kneib, J.-P., Percival, W. J., et al. 2016, AJ, 151, 44
• de Jong, R. S., Agertz, O., Berbel, A. A., et al. 2019, The Messenger, 175, 3
• Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168
• Dhawan, S., Jha, S. W., & Leibundgut, B. 2018, A&A, 609, A72
• Eifler, T., Miyatake, H., Krause, E., et al. 2021, MNRAS, 507, 1746
• Fan, X., Wang, F., Yang, J., et al. 2019, ApJ, 870, L11
• Fian, C., Mediavilla, E., Motta, V., et al. 2021, A&A, 653, A109
• Freedman, W. L., Madore, B. F., Hatt, D., et al. 2019, ApJ, 882, 34
• Gaia Collaboration (Brown, A.G.A., et al.) 2018, A&A, 616, A1
• Goicoechea, L. J., & Shalyapin, V. N. 2022 ApJ, submitted [arXiv:2201.01271v1]
• Gorski, K. M., Wandelt, B. D., Hansen, F. K., Hivon, E., & Banday, A. J. 1999, ArXiv e-prints [arXiv:astro-ph/9905275]
• He, Z., & Li, N. 2022, RAA, 22, 095021
• He, Z., Er, X., Long, Q., et al. 2020, MNRAS, 497, 556
• Huang, X., Storfer, C., Gu, A., et al. 2021, ApJ, 909, 27
• Ivezić C.Z., Kahn, S.M., Tyson, J.A., et al. 2019, ApJ, 873, 111
• Kormann, R., Schneider, P., & Bartelmann, M. 1994, A&A, 284, 285
• Krone-Martins, A., Delchambre, L., Wertz, O., et al. 2018, A&A, 616, A11
• Kuijken, K., Heymans, C., Dvornik, A., et al. 2019, A&A, 625, A2
• Lanusse, F., Ma, Q., Li, N., et al. 2018, MNRAS, 473, 3895
• Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv e-prints [arXiv:1110.3193]
• Lemon, C. A., Auger, M. W., McMahon, R. G., & Ostrovski, F. 2018, MNRAS, 479, 5060
• Lemon, C.A., Auger, M.W., & McMahon, R.G. 2019, MNRAS, 483, 4242
• Lemon, C., Anguita, T., Auger-Williams, M. W., et al. 2023, MNRAS, 520, 3305
• Li, R., Napolitano, N. R., Spiniello, C., et al. 2021, ApJ, 923, 16
• Liao, K., Shafieloo, A., Keeley, R. E., & Linder, E. V. 2019, ApJ, 886, L23
• Liddle, A. R. 2007, MNRAS, 377, L74
• Lyke, B. W., Higley, A. N., McLane, J. N., et al. 2020, ApJS, 250, 8
• Macaulay, E., Nichol, R. C., Bacon, D., et al. 2019, MNRAS, 486, 2184
• McGreer, I. D., Hall, P. B., Fan, X., et al. 2010, ApJ, 140, 370
• More, A., Oguri, M., Kayo, I., et al. 2015, MNRAS, 456, 1595
• Motta, V., Mediavilla, E., Falco, E., & Muñoz, J. A. 2012, ApJ, 755, 82
• Niedermann, F., & Sloth, M. S. 2022, Phys. Rev. D, 105, 063509
• Oguri, M., & Marshall, P. J. 2010, MNRAS, 405, 2579
• Oguri, M., Rusu, C. E., & Falco, E. E. 2014, MNRAS, 439, 2494
• Petrillo, C. E., Tortora, C., Chatterjee, S., et al. 2017, MNRAS, 472, 1129
• Pietrzyński, G., Graczyk, D., Gallenne, A., et al. 2019, Nature, 567, 200
• Planck Collaboration VI. 2020, A&A, 641, A6
• Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, ApJ, 876, 85
• Rojas, K., Savary, E., Clément, B., et al. 2022, A&A, 668, A73
• Schmidt, T., Treu, T., Birrer, S., et al. 2023, MNRAS, 518, 1260
• Shajib, A. J., Treu, T., & Agnello, A. 2018, MNRAS, 473, 210
• Shajib, A. J., Birrer, S., Treu, T., et al. 2020, MNRAS, 494, 6072
• Sonnenfeld, A., & Cautun, M. 2021, A&A, 651, A18
• Sonnenfeld, A., Treu, T., Gavazzi, R., et al. 2013, ApJ, 777, 98
• Spiniello, C., Agnello, A., Napolitano, N. R., et al. 2018, MNRAS, 480, 1163
• Spiniello, C., Agnello, A., Sergeyev, A. V., et al. 2019, MNRAS, 483, 3888
• Suyu, S. H., Treu, T., Hilbert, S., et al. 2014, ApJ, 788, L35
• Vagnozzi, S. 2020, Phys. Rev. D, 102, 023518
• Van de Vyvere, L., Gomer, M.R., Sluse, D., et al. 2022, A&A, 659, A127
• Verde, L., Treu, T., & Riess, A. G. 2019, Nat. Astron., 3, 891
• Wong, K. C., Suyu, S. H., Chen, G. C. F., et al. 2020, MNRAS, 498, 1420
• Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
• Wu, Q.-Q., Liao, S.-L., Ji, X., et al. 2022, Front. Astron. Space Sci., 9, 822768
• Yèche, C., Palanque-Delabrouille, N., Claveau, C.-A., et al. 2020, Res. Notes AAS, 4, 179
• Yuan, W., Riess, A. G., Macri, L. M., Casertano, S., & Scolnic, D. M. 2019, ApJ, 886, 61
• Yue, M., Fan, X., Yang, J., & Wang, F. 2022, ApJ, 925, 169
• Zou, H., Zhou, X., Fan, X., et al. 2017, PASP, 129, 064101
• Zou, H., Zhou, X., Fan, X., et al. 2019, ApJS, 245, 4

All Tables

All Figures

	Fig. 1 Magnitude distributions of the QCC in g, r, z (left panel) and W1, W2 (right panel). The kernel density estimate curves (solid lines) are generated by kdeplot in seaborn package with binsize = 0.5 and the default Gaussian kernel.
In the text

Fig. 2 Demonstrating the three possible situations that cause a known multiply-imaged QSO to have only one match in the QCC. The colours of the cycles represent the morphology types in DESI-LS catalogues. Green circles represent point-like sources, whereas white circles represent extended sources. For WG0214-2105, three out of its four images are recorded in the catalogue, but only one is labelled as point-like. For Β1608+656, the image with the green pentagon marker lacks a g-band detection. In the case of SDSSJ0924+0219, three of the lensed images have been identified as a single point-like source and the fourth image has been identified as an extended object.

In the text

	Fig. 3 Distributions of Scolour (upper panel) and SRF (lower panel) for systems in QGC, PAIR-CS-RF, MUL-RF-CS, H22 (all candidates), and RLQ.
In the text

	Fig. 4 Examples of grade A candidates in H22. The top row shows candidates also found in RLQ, the middle row shows candidates also found in the candidate catalogue of D22 and the bottom row shows new systems neither contained in RLQ nor D22. Each image is titled with the internal ID (linking the system to QGC) and gives the values of SRF, Scolour, and image separation.
In the text

	Fig. 5 Distributions of Scolour (top two panels) and SRF (bottom two panels) comparing our new candidate systems (H22) with those of D22, L22, and RLQ (left panels) and splitting our new candidates by awarded grade (right panels).
In the text

	Fig. 6 Comparison of redshift distributions of H22, D22, L22, and OM10. We note that the redshifts of L22are from Lemon et al. (2023), while those of H22 and D22are from SDSS eBOSS DR16.
In the text

	Fig. 7 Comparison of magnitude distributions of H22, D22, and OM10 shown in the left panel. Since L22did not provide g-band magnitudes, it is not plotted here. Distributions split by awarded grade in H22 are shown in the right panel.
In the text

	Fig. 8 Comparison of max image separation of systems in H22, D22, L22, RLQ, and OM10. The panel on the right splits this distribution for H22 by awarded grade.
In the text

	Fig. 9 Comparisons between H22-new-RLQ, H22-new-D22, H22-new-L22, and H22-new of distributions of Scolour (upper left), SRF (upper right), g-band magnitude (lower left), and image separation (lower right).
In the text

Fig. 10 Distribution of BICs of the 52 lensing models. BICs of two-image systems are represented by histograms, while vertical orange lines indicate the BICs of the two quad systems. In the square panels, blue circles show the observed image positions, orange circles show the modelled image positions, and green diamonds show the source positions given by the lens models. The green and red curves indicate the caustic and critical curves of the corresponding lensing models.

In the text

	Fig. 11 Modelled caustic curves for the three candidates in the third line of Fig. 4. The plots are centred at the lens light centre, which is determined by peak of the lens light. Symbols are the same as in Fig. 10.
In the text