Finding Multiply Lensed and Binary Quasars in the DESI Legacy Imaging Surveys

The time delay between multiple images of strongly lensed quasars is a powerful tool for measuring the Hubble constant (H 0). To achieve H 0 measurements with higher precision and accuracy using the time delay, it is crucial to expand the sample of lensed quasars. We conduct a search for strongly lensed quasars in the Dark Energy Spectroscopic Instrument (DESI) Legacy Imaging Surveys. The DESI Legacy Surveys comprise 19,000 deg2 of the extragalactic sky observed in three optical bands (g, r, and z), making it well suited for the discovery of new strongly lensed quasars. We apply an autocorrelation algorithm to ∼5 million objects classified as quasars in the DESI Quasar Sample. These systems are visually inspected and ranked. Here, we present 436 new multiply lensed and binary quasar candidates, 65 of which have redshifts from Sloan Digital Sky Survey Data Release 16. We provide redshifts for an additional 18 candidates from the SuperNova Integral Field Spectrograph.


Introduction
Multiply lensed quasars are powerful cosmological probes.In particular, quantifying the time delay between quasar images can deliver an independent measurement of the Hubble constant, H 0 (e.g., Refsdal 1964;Treu & Marshall 2016;Suyu et al. 2017).As tension persists between the value of H 0 inferred from the cosmic microwave background observations assuming the flat ΛCDM model (Planck Collaboration et al. 2020) and direct late-Universe measurements of H 0 (e.g., Riess et al. 2019Riess et al. , 2022)), the development of a precise, accurate, and independent measurement of H 0 becomes ever more important.Careful analysis of small samples of lensed quasars have already yielded measurements of H 0 with competitive precision (e.g., Wong et al. 2020).For lensed quasars to make significant progress toward resolving the current tension, statistical uncertainties need to be reduced to below 1% (Treu & Marshall 2016), while effectively addressing systematic uncertainties (Birrer et al. 2020;Shajib et al. 2020).A larger and more diverse sample of lensed quasars can improve measurements of H 0 to subpercent precision (Rathna Kumar et al. 2015;Treu & Marshall 2016;Birrer & Treu 2021).
Beyond the study of cosmological parameters, lensed quasars offer opportunities for advancement in black hole physics, especially in understanding the co-evolution of supermassive black holes and their host galaxies.Local Universe (z < 0.1) observations have found that supermassive black hole mass (M BH ) and various properties of the host galaxy (e.g., bulge luminosity, stellar velocity dispersion, and stellar mass) are tightly correlated, suggesting that central black holes and their host galaxies are physically coupled (Kormendy & Richstone 1995;Magorrian et al. 1998;Ferrarese & Merritt 2000;Gebhardt et al. 2000;Marconi & Hunt 2003;Häring & Rix 2004).To understand the coupling of central black holes and their host galaxies at earlier epochs, it is important to observe high-redshift galaxies beyond the local Universe.Strong gravitational lensing amplifies high-redshift quasars and magnifies their host galaxies.This, coupled with high-resolution imaging, allows for analysis of quasar and host galaxy properties through modeling and reconstruction (Ding et al. 2017(Ding et al. , 2021)).Despite this potential for lensed quasars in black hole physics, the relatively small number of known systems limits their impact, making the discovery of more lensed quasars important.
There are now ∼200 confirmed lensed quasar systems (e.g., Weymann et al. 1979;Inada et al. 2012;More et al. 2016;Agnello et al. 2018;Lemon et al. 2018Lemon et al. , 2019;;Jaelani et al. 2021) and ∼130 lensed quasar candidates (e.g., Sergeyev et al. 2016;Agnello et al. 2018;Chan et al. 2020), a development that can be mostly attributed to large-scale surveys such as the Sloan Digital Sky Survey (SDSS; e.g., Inada et al. 2003), the Kilo-Degree Survey (e.g., Spiniello et al. 2018), the Dark Energy Survey (DES; e.g., Anguita et al. 2018), PanSTARRS (e.g., Rusu et al. 2019), and Gaia (e.g., Lemon et al. 2018).As the search for lensed quasars pivots toward large-scale surveys, recent papers tend to include many candidates alongside new confirmed systems (e.g., Chan et al. 2020;Lemon et al. 2020).This trend can be attributed to the much increased size of the new large-scale surveys and high volume of new discoveries.Upon further examination of these candidates, many quasar pairs, in particular, turn out to be physical quasar-quasar binaries rather than the result of lensing (Kochanek et al. 1999;Mortlock et al. 1999).Quasar pairs that turn out to be physical binaries are also valuable, however, as physical binaries lend insight to how galaxy-galaxy interactions and mergers can enhance or even trigger quasar activity (Begelman et al. 1980;Djorgovski 1991;Di Matteo et al. 2005;Hopkins et al. 2008;Ellison et al. 2011;Liu et al. 2012;Bogdanović et al. 2022), or as the progenitor of binary black hole systems (e.g., Boroson & Lauer 2009).
We have searched the Dark Energy Spectroscopic Instrument (DESI) Legacy Imaging Surveys for multiply lensed quasars.Previous papers have searched the DESI Legacy Imaging Surveys for strongly lensed galaxies (Huang et al. 2020(Huang et al. , 2021;;Stein et al. 2022;Storfer et al. 2022), but this is the first search of the survey for lensed quasars in particular.We note that in this paper, "strongly lensed" and "multiply lensed" are used interchangeably.In Section 2, we give an overview of the DESI Legacy Surveys.We describe our methodology in Section 3. In Section 4, we present the results of our search, as well as the spectroscopic observations of a subset of our candidates.We discuss our results in Section 5 and conclude in Section 6.

Observations
The DESI Legacy Imaging Surveys (Dey et al. 2019) cover ∼19,000 deg 2 of the extragalactic sky visible from the northern hemisphere.The footprint is split into two contiguous parts by the Galactic plane, and each was observed with at least three passes in the grz bands (Figure 1).The Legacy Surveys are composed of a northern region and a southern region (see Figure 1).The 14,000 deg 2 southern region, with δ  + 32°, comprises the Dark Energy Camera Legacy Survey (DECaLS).For the 5000 deg 2 northern region, the Beijing-Arizona Sky Survey (BASS; Zou et al. 2017) carried out the gr-band observations and the Mayall z-band Legacy Survey (MzLS), the z-band observations.For DECaLS, the g, r, and z bands deliver an image quality with FWHM of approximately 1 29, 1 18, and 1 11, respectively, compared with 1 61, 1 47, and 1 01 in the MzLS/BASS subregion.Figure 1 shows the subregions of the Legacy Survey footprint and their depths.
The Legacy Surveys serve as the precursor to the DESI spectroscopic survey with the purpose of identifying targets for spectroscopic observations.The Legacy Surveys source catalogs are constructed using The Tractor (Lang et al. 2016), a forward-modeling algorithm that performs source extraction on pixel-level data.A small set of light profiles (point source, point-spread function, or PSF; round exponential galaxy with a variable radius; de Vaucouleurs profile, for elliptical galaxies; exponential profile, for spiral galaxies; and Sérsic profile) are fit to each source to determine the best-fit model. 7Important to our search are sources modeled as PSFs, which are likely quasars or stars.
The DESI Quasar Sample (Yèche et al. 2020) identifies the potential quasars among DESI targets.DESI quasar selection utilized the three optical bands (grz) of the Legacy Surveys as well as the W1 and W2 bands of the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010).The criteria for selection are restricted to objects with stellar morphology (PSFs), r-band magnitude <22.7 AB mag, and to targets not in areas with corrupted or masked pixels (e.g., targets in the vicinity of bright stars, globular clusters, or bright galaxies).Yèche et al. (2020) used two distinct methods for quasar target selection: color cuts and a machine-learning algorithm.The color cuts consist of W1 − W2 and r − W versus g − z, where W is the magnitude computed from the weighted average of W1 and W2 fluxes.The W1 and W2 bands are crucially important in distinguishing quasars from stars due to the infrared excess observed in quasars at all redshifts.In addition, a random forests (RF)-based algorithm has been trained on 230,000 known quasars within the Legacy Surveys footprint and 210,000 unclassified objects presumed to be stars.98% of quasar targets selected via the color-cut method are contained within the RF selection, which is more complete than the color-cut method at low and high redshifts.1) within the Legacy Surveys footprint are overlaid.
While both methods have been employed in selecting targets for the DESI main survey, the current DESI Quasar Sample corresponds to the RF selection.

Methodology
We assemble a catalog of known lensed quasars within the Legacy Surveys footprint, composed primarily of three separate, earlier catalogs: the Master Lens Database (Moustakas et al. 2012), the Gravitationally Lensed Quasar Database,8 and the VizieR Online Data Catalog: Gaia GraL.II.Known Multiply Imaged Quasars (Delchambre et al. 2019).We have since added dozens of confirmed lenses and candidates from more recent publications (Sonnenfeld et al. 2018(Sonnenfeld et al. , 2020;;Chan et al. 2020;Jaelani et al. 2021;Stern et al. 2021).Several of these papers present both candidates and confirmed systems (e.g., Agnello et al. 2018).In total, the catalog consists of ∼300 known lensed quasars and candidates (see Table 1).The compilation of this known lens catalog is used to inform decisions regarding our autocorrelation algorithm and subsequent human inspection.
To search for lensed quasars in the DESI Legacy Surveys, we perform autocorrelation on the aforementioned DESI Quasar Sample, which contains over 5 million targets.About 30% (94) of the ∼300 known lensed quasars within the Legacy Surveys footprint have two or more objects in the DESI Quasar Sample and therefore are "discoverable" by our algorithm (henceforth, "discoverable known" systems; see Figure 1, Table 1).Figure 2 shows the distribution of image separations among "discoverable known" systems.For the remaining known systems (∼70%), a subset are the candidates from the SuGOHI project (Sonnenfeld et al. 2018(Sonnenfeld et al. , 2020;;Chan et al. 2020;Jaelani et al. 2021), which, though promising lensed candidates, appear qualitatively different from the other known systems in the Legacy Surveys.9Excluding the SuGOHI candidates, the distribution of image separations for the rest10 is similar to the 94 "discoverable known." To err on the side of completeness, we identify and group multiple quasar targets within a 10″ radius of each other.We later reduce our separation cut by half, because more than 95% of the "discoverable known" systems have quasar targets separated by less than 5″.We perform autocorrelation within each of the Legacy Survey "bricks.""Bricks" are approximately 0.26 × 0.26 deg 2 (3600 pixels × 3600 pixels) in size and typically contain 15-20 quasar targets each.Performing autocorrelation within each brick risks missing groupings of quasars that straddle the boundaries of two or more bricks.With a 5″ image separation cut, we would miss some systems whose center happens to lie within 2 5 of a brick boundary.However, only some such systems are lost, depending on their position angle.Accounting for position angle, 0.3% is a generous upper limit to the fraction of lensed quasar candidates we would have missed due to not searching for a nearby target in a neighboring brick.
We developed our own algorithm in Python.The algorithm proceeds brick by brick (we made use of the mpi4py package to run in parallel across bricks).In each brick, the algorithm first groups pairs of quasars within 10″ of each other and then  recursively connects multiple pairs of quasar targets within the same candidate system (such as would occur in "triples" and "quads").11We test our algorithm on the 94 "discoverable known" systems and simulations, achieving a 100% success rate.We then deploy our autocorrelation algorithm on the ∼5 million quasar targets across the Legacy Surveys footprint.

Results
In this section, we present the results from the autocorrelation and subsequent human inspection.The 10″ radius search found over 27,000 recommendations (quasar targets, each of which is within 10″ of at least one other).Due to the abundance of recommendations and because more than 95% of the aforementioned "discoverable known" lenses (see Table 1) have quasar targets that are separated by less than 5″, we choose to impose a cut on the recommendations to those with separation <5″.This results in almost 6000 recommendations.Based on this and the DESI Quasar Sample surface density, we estimate that 2.4 recommendations will be chance superpositions.These will likely make it into our sample.
These recommendations are subjected to human inspection and graded as A, B, C, or nonlens.Of the 94 "discoverable known" systems, 72 are confirmed (henceforth, "discoverable confirmed" systems; see Table 1).The grading criteria are informed by the appearance of these "discoverable confirmed" systems in Legacy Surveys images.Figure 3 shows an assortment of these confirmed lenses.The following are common features of known lensed quasars that guided the visual inspection and grading: 1. quasar targets with similar color that are <2 5 apart, but sometimes slightly farther; 2. the presence or hint of a redder (putative lensing) galaxy between the quasar images; 3. in the case of unclear lensing galaxy light, one of the quasar images appearing dimmer and/or redder, which may be due to its closer proximity to the lensing galaxy; Because lensed quasars typically display a combination of these features, there is great diversity in their appearance.Furthermore, for some of the candidate systems, The Tractor has identified additional point sources, providing evidence that some systems we initially thought were doubles are actually quads (for 17 of the 24 quads identified by our algorithm, they each have only two objects from the DESI Quasar Sample, but three or four point sources from The Tractor).
To begin, C.D. and C.S. made "first pass" selections according to the above criteria.As a "second pass," C.D., C. S., and X.H. together assigned an integer score between 1 and 4 to the "first pass" selections, while coauthor A.D. did the same independently.These two scores were then averaged.This grading scheme is similar to that of Huang et al. (2021).The final grades for the candidates can be broken down as follows: 1. 3.5: Grade A. We are highly confident that these candidates are lensed quasars.Many of them have angular separations larger than the seeing and the putative lens is visible (or there is a hint of lens light) and the quasar images have similar colors.Others are close pairs and quads.2. = 3.0: Grade B. These have characteristics that are similar to those of the Grade As, but weaker.The lens light is often not obvious in close pairs.For some systems, the putative counterimages are somewhat redder, possibly due to contamination from the lens light, or just fainter.3. = 2.5 or 2.0: Grade C.These systems generally have large image separations and are fainter than the As and Bs.Many of them do not have a discernible lensing galaxy.For the few cases with a possible lens or lens light, they have an atypical lensing configuration.
In total, we have identified 530 candidate lensed quasar systems.As stated before, all 94 "discoverable known" systems have been rediscovered by our algorithm as recommendations.We therefore have found 436 new candidate systems.Table 1 gives a breakdown of the relevant types of systems discussed in this paper.
Figure 4 shows the magnitude distribution of our newly discovered candidates compared to that of "discoverable confirmed" systems.Our candidates include many more faint objects than the confirmed lensed quasars found in previous surveys.Figure 5 shows our new candidates plotted over a depth map of the Legacy Surveys.Among the new candidates, there are 102 As, 118 Bs, and 216 Cs. Figure 6 shows 10 notable systems, with nine doubles and one potential quad.
Figure 7 and Table 2 show the first 80 of our new candidates by ascending R.A., grouped by grade.All 436 new candidates can be found on our project website: https://sites.google.com/usfca.edu/neuralens.

Gaia Proper Motions and Parallax
Gaia Early Data Release 3 is used to provide further checks on the quality of our candidates; specifically, the low significance of both proper motion (PM) and parallax is an important indicator that a candidate system is composed of quasars and not Milky Way stars.380 (∼87%) of our candidates have Gaia PM and parallax information for at least one quasar image, which we obtained from the Gaia Archive. 12e follow a different definition of PM significance (PMSIG) compared to Lemon et al. (2019).By our definition, Additionally, we take into consideration parallax significance (PXSIG), where PXSIG PX .

PX s =
Figure 8 shows the distributions of PMSIG and PXSIG for the 380 new candidates and 72 "discoverable confirmed" systems with Gaia information.Informed by the PMSIG and PXSIG of the "discoverable confirmed" systems, we decide that PMSIG < 8 and PXSIG < 3.5 are acceptable for our candidates.

Redshifts for Our Candidate Systems
For the vast majority of our systems, the putative lensing galaxy is much fainter than the lensed images.As a result, we are only able to find photometric redshifts for 27 lenses (Section 4.2.1).For the lensed sources (quasars), we identified 65 systems with spectroscopic redshifts from SDSS and we have additionally obtained spectroscopy using the SuperNova Integral Field Spectrograph (SNIFS; Aldering et al. 2002) for a subset of our systems (Section 4.2.2).Six of our candidates have photometric redshifts for the foreground galaxy and spectroscopic redshifts for the background quasar.These are presented in Section 4.2.3.

Photometric Redshifts of Lensing Galaxies
For redshift information on our lensing galaxies, we include photometric redshifts for 27 of our candidate systems from Zhou et al. (2021).In cases where the centroid of The Tractorextracted source is ambiguous (i.e., straddling the putative lens and a lensed image), we did not include the photometric redshift.In cases where the putative lens and at least one lensed image are very close, the photometric redshift may be less reliable.

Spectroscopic Redshifts from SNIFS and SDSS
We have observed a subset of our high-quality candidates with r-band mag 20.0 on the SNIFS instrument on the University of Hawaii 88 inch (2.2 m) telescope (UH88), located on Maunakea.The spectrograph has two channels, blue (320 nm to 560 nm) and red (520 nm to 1μ m), with resolutions of R ∼1000 and R ∼1300, respectively.Our candidates were observed on the nights of 2020 September 21, October 12, October 19, and November 9 as well as 2021 June 11, June 14, July 12, and November 8 UTC.SNIFS has a 6.4 × 6 4 field of view.Given this field of view, we do not need to impose an image separation cut for our targets.SNIFS splits its field of  view via a 15 × 15 microlens array into 225 samples, each 0.43 × 0 43 (Lantz et al. 2004).The spectra were extracted from each spatiospectral data cube using a circular aperture with a radius of 2σ and a surrounding annual sky region spanning 5-7σ, where σ is the second moment of the PSF of the lensed quasar candidate system.For this pilot project, we did not attempt to disentangle the light from the individual sources.
We fit for the SNIFS spectra using the publicly available SDSS Data Release (DR) 5 spectral template for quasars. 13Of the 29 candidates observed, 10 had spectra with an insufficient signal-to-noise ratio (S/N; due to poor observing conditions), one was a Milky Way star,14 and 18 were confirmed as quasars.For these 18, the best-fit redshifts range from z = 0.75 to 2.60 (Figure 9).We later found that DESI-237.7387+ 02.3629 was observed by SDSS with a redshift consistent with the value from SNIFS (both are included in the online table on our project website).
For additional redshift information, we also include spectroscopic redshifts from SDSS DR16 for 65 quasars in our candidate systems. 15The redshifts of these quasars range from 0.36 to 3.55, and 61 of them have z > 0.5.Figure 11 shows redshift distributions for our candidates.The high-redshift distribution of the quasars is consistent with them being lensed.
We note that DESI-037.5191-07.0795appears to be a quasar with broad absorption-line (BAL) activity. 16We identify three potential rare quasar-quasar lensing systems:   image separation 3 69.17Last, we have found a potential double-source lensing system, DESI-181.2111+ 44.4764.The two quasar images have SDSS redshifts of 1.1438 and 1.8095, respectively, with an image separation of 3 99.In between these two images, there is clearly a galaxy foreground to both quasar images (based on color), which we consider to be the putative lens.For 11 of the 18 targets shown in Figure 9, the S/N is sufficiently high to see plausible evidence of multiple sources in wavelength slices from SNIFS data cubes that correspond to emission features.In Figure 10 doublet) can be observed for both quasars and galaxies, the first five are all features much more likely to be present for a quasar than a galaxy. 18Ten of the 11 systems have three to five of the five features predominantly seen in quasars, as mentioned above.For the last system, it is true that the emission features of Hβ and the [O III] doublet could be due to the quasar and its host galaxy, but it is still more plausible that the features are from multiple quasar images, because the contours are drawn around objects from a quasar catalog.Each slice is 6 4 on the side, corresponding to the SNIFS field of view.The shifting of the emission features from the longest to the shortest wavelength is due to atmospheric differential refraction (ADR).The ADR effect is greater for images further south, as expected, given the higher airmass from Maunakea.These wavelength slices show that in Notes.The first 80 new candidate systems arranged in ascending R.A., all of which are doubles.Thus, each system has two rows, corresponding to the two quasar images.To avoid confusion, the two rows for each system have the same system name.In the table on our project website, we provide a cutout centered on the individual quasar image in each system.Therefore, there will be no ambiguity.The columns are arranged as follows.Column (1): name of the system (R.A. and decl. in decimal format).Column (2): human inspection grade.Column (3): spectroscopic redshift of quasars from SDSS DR16 or SNIFS (with an asterisk).Column (4): photometric redshift of the putative lenses from Zhou et al. (2020).Columns (5), (6), and (7): the g-, r-, and z-band magnitudes, respectively.Column (8): parallax significance.Column (9): PM significance.Column (10) average separation of images in arcseconds.This portion of the table only shows the first 80 new candidates, while the full online version has the complete set of 436 candidates.
(This table is available in its entirety in machine-readable form.)these 11 systems, multiple quasar images are at the same redshift, suggesting that they are either lensing systems or binary quasars.
The SNIFS pilot study shows that to fully confirm a high fraction of these candidate systems, deeper observations with higher resolution are needed.Adaptive-optics-assisted integral field spectrographs (e.g., Very Large Telescope MUSE and Keck OSIRIS) would be particularly useful to obtain redshifts of the lens and multiple lensed images.Such observations would be also needed for lens modeling and constraining H 0 .

Candidates with Both Foreground and Background Redshifts
The following six candidate systems have photometric redshifts for the foreground galaxy and spectroscopic redshifts for the background quasar, with three Grade As, two Bs, and one C.We note that the lensing galaxies are ellipticals, which typically have reliable photometric redshifts (Zhou et al. 2020). 19The redshifts and image separations for these systems are consistent with lensing.However, for full confirmation, spectroscopic and high-resolution imaging and/or modeling are needed.We provide a brief description for each of them below.
DESI-028.0797-24.8105:this Grade A candidate is one of the visually impressive systems shown in Figure 6.It is also one of the first 80 systems arranged by R.A. (Table 2).It was observed by SNIFS (Figure 9) and we determined the quasar redshift to be 2.42.The photometric redshift for the foreground galaxy, the putative lens, is 0.286 ± 0.110.One of the two quasar images has the second highest PMSIG among the quasar images for these six systems, at 2.57, but it has a low PXSIG of 0.69.The quasar image separation is 1 93.DESI-060.4504-25.2439:this Grade A system is also one of the visually impressive systems shown in Figure 6.Its redshift was determined to be 1.33 from SNIFS observation (Figure 9).The putative lensing galaxy has a photometric redshift of 0.465 ± 0.159.All images have relatively low PXSIG and PMSIG, consistent with zero.The quasar image separation is 1 45.DESI-251.5695+ 44.0197: this is the third of the three Grade A systems.SDSS provides a quasar redshift of 2.0757, and the photometric redshift for the foreground galaxy is 0.620 ± 0.048.We did not find Gaia parallax or PM information for the images in this system.The image separation is 3 43.DESI-141.7488+ 06.3905: this is a Grade B system.The SDSS redshift for one quasar image is 0.6671 and the foreground galaxy has z phot = 0.620 ± 0.048.One of the images has the highest PMSIG among the quasar images for these six systems, at 3.87, but the PXSIG is only 0.04.This candidate has the largest image separation among the six systems, at 4 96.this is also a Grade B system.One of the quasar images and the foreground galaxy have a redshift of 1.6920 from SDSS and a photometric redshift of Figure 8.The distributions of significance of both PM and parallax in the images of our newly discovered candidates and of "discoverable confirmed" systems (see Table 1).
Figure 9. Legacy Surveys image cutouts alongside spectra from SNIFS on UH88 for 18 observed candidates with sufficient S/N ratio.We perform a correlation between the spectra and the SDSS quasar spectral template.Here we report the best-matched redshift.For each system, the spectrum is shown in orange, with the error spectrum in magenta shading.Redshifted templates are shown in blue, with the black lines representing emission features.DESI-237.7387+ 02.3629 was previously observed by SDSS with consistent redshift.Note that, for some systems, only the blue channel of SNIFS was available at the time of observation.
Figure 10.Each column shows an observed wavelength slice from SNIFS, corresponding to emission features in ascending wavelength, from Lyα to Mg II (see the text).We impose the contour from the Legacy Survey images (set to 80% of the brightest pixel for each image) on each of the slices.For each system, the extent, shape, and orientation of the contour matches the brightest pixels in the SNIFS slices, suggesting the presence of two or more source images.For each of these 11 systems (one per row), these wavelength slices show that there are multiple quasar images at the same redshift.The emission features of the last row do not align with the rest, given its lower redshift.Note the effects of the ADR (see the text).

Discussion
In this section, we assess the fraction of our Grade A candidates that are confirmed as quasars (Section 5.1) and compare the double-to-quad ratio for our candidates with forecast results (Section 5.2).

Contamination Estimation
For SNIFS targets, we selected visually convincing systems with r-band mag 20.0, which approximately corresponds to the brighter half of the distribution (Figure 4).Not surprisingly, a strong majority of these targets are Grade A systems: 15 of the 19 targets with sufficient S/N.One of them turned out to be a star (as mentioned in Section 4.2.2, this was subsequently removed from our candidate list), with the rest confirmed to be quasars.Among the 102 Grade A candidates, as with the rest of our candidates, approximately half of them (55) have at least one image with r-band mag <20.0.Based on the SNIFS targets, we expect 93% of them, or 51, to be actual quasars.

Double-to-quad Ratio
For the 436 new lensed quasar candidates discovered in the Legacy Surveys using the DESI Quasar Sample plus the 94 "discoverable known" systems from the same sample, 506 are doubles and 24 are quads.This amounts to a double/quad ratio of ∼21 to 1, much higher than expected.Oguri & Marshall (2010;hereafter, OM10), while estimating the double/quad ratio of detectable lensed quasars (those with a minimum image separation of (2/3) * θ PSF ) in various surveys according to their depth, predict that the DES will contain a double/quad ratio of ∼6 to 1 (they specifically cite the percentage of quads).DES is part of the DECaLS subregion of the Legacy Surveys.By filtering our candidates appropriately, we may more fairly compare our double/quad ratio to that estimated in DES by OM10.
In order to compare with the OM10 prediction, we restrict our sample as follows: 1.The MzLS/BASS subregion has worse seeing in the gr bands than DECaLS.Given our search algorithm, this will have a significant effect on the double/quad ratio in MzLS/BASS compared to DECaLS.Because most quads are not perfect Einstein crosses and are instead asymmetrical, often with three close-by (even blended) images and one fainter counterimage, this would imply only two objects would be identified by the DESI Quasar Sample.Our algorithm would then identify them as a double, resulting in the loss of a disproportionate number of quads.Thus, for comparison with the OM10 ratio, we consider the double/quad ratio for our candidates in DECaLS only.2. OM10 only considered systems with image separation >(2/3) * θ PSF , so using θ PSF = 1 18 in the r band in DECaLS, we further restrict our candidates in DECaLS to those with image separation >0 79. 3. We also restrict our candidates to those graded as As.
This decision is motivated primarily by our higher confidence in As compared to Bs and Cs.While we are confident in the quality of our B-and C-grade candidates, we are less confident in their typing as doubles or quads.
With the above restrictions-that is, A-grade candidates with image separation >0 79 within the DECaLS footprint-our double/quad ratio is ∼11.4 to 1 (80 to 7), which still differs from that predicted for DES (∼6 to 1).Poisson noise alone does not seem to fully account for the discrepancy between our ratio and OM10ʼs.One possible explanation for the discrepancy is that some of the doubles are physical quasar-quasar binaries, rather than doubly lensed quasars (e.g., Lemon et al. 2020).We note in Section 1 that physical binaries are also important for quasar physics.Follow-up spectroscopic and high-resolution observations are needed to determine the fraction of our candidates that are physical binaries.Traditionally, quads have been more highly valued because they more tightly constrain lens modeling (Suyu et al. 2017).Doubles, though, generally have much more precise time-delay measurements, because their time delays tend to be longer than those of quads (OM10).Current measurements of H 0 using doubles are becoming more competitive.For example, using 27 doubly lensed quasars, taking into account various systematic effects, Harvey (2020) reported a 4% measurement of H 0 .This is important, as doubles will dominate over quads among the lensed quasars yet to be discovered.Magnification bias strongly favors quads (Wallington & Narayan 1993), so deeper searches for lensed quasars are expected to yield more doubles than quads (e.g., Agnello et al. 2018;Lemon et al. 2019Lemon et al. , 2020;;Chan et al. 2020).With future surveys slated to produce ever higher double/quad ratios-the Vera C. Rubin Legacy Survey of Space and Time is expected to deliver a double/quad ratio of ∼7 to 1 (OM10)-doubly lensed quasars can play a significant role toward 1% precision H 0 measurements.

Conclusions
We have carried out a search for multiply lensed quasar systems in the DESI Legacy Imaging Surveys.We first apply an autocorrelation algorithm to the 5 million targets in the DESI Quasar Sample (Yèche et al. 2020).Then, guided by the visual appearance of the confirmed systems in the Legacy Surveys footprint, we inspect ∼6000 systems with image separations <5″ found by our algorithm.Among these, we discover and report 436 new lensed quasar candidates, with 102 Grade As, 118 Bs, and 216 Cs.We have found quasar redshifts for 65 systems from SDSS DR16, which range from 0.36 to 3.55.Of these, 61 have z > 0.5.We have obtained spectra using SNIFS on the University of Hawaii 2.2 m telescope for 18 additional quasars (one of which was previously observed by SDSS with consistent redshift).The best-fit redshifts for these 18 range from z = 0.75 to 2.60.The high-redshift distribution of the quasars is consistent with them being lensed.Based on the SNIFS observations, we estimate ∼93% of the brighter half of our Grade A candidates are actual quasars.Since our candidates are discovered from the DESI Quasar Sample, all of them will be spectroscopically observed by DESI.Among our candidates, we have identified the following: one system with BAL activity; three potential rare quasar-quasar lensing systems; a possible double-source lensing system; and six candidates with both foreground and background redshifts that are consistent with lensing.
These are very promising findings and represent a significant addition to the existing sample of lensed quasars.Oguri & Marshall (2010) predicted a double/quad ratio in DES of 6 to 1.After appropriate restrictions are applied, the corresponding ratio for our search results is ∼11.4 to 1.Even after accounting for Poisson noise, it is possible that some of the doubles among our candidates are physical binaries, which are important for gaining a deeper understanding of quasar physics.These discoveries will contribute a large number of systems for the time-delay measurement of H 0 with high accuracy and precision.

Figure 1 .
Figure 1.The DESI Legacy Imaging Surveys footprint in an equal-area Aitoff projection in equatorial coordinates.The dotted line separates the north (MzLS/BASS) and south (DECaLS) subregions.Slightly above δ = 32°, there is a small amount of overlap between MzLS/BASS and DECaLS.Patches with different shades of blue indicate the depth in the z band (5σ detection significance)."Discoverable" lensed quasars and other known lensed quasars (see Section 3,Table 1) within the Legacy Surveys footprint are overlaid.

Figure 2 .
Figure 2. Histogram of separations between lensed quasar images for the 94 "discoverable known" systems (see the text).

Figure 3 .
Figure 3.In this figure, we show representative configurations of the "discoverable confirmed" lensed quasar systems.The left panel consists of doubles: the first column shows systems that require deblending; the second column shows small-image-separation systems (<2 5); and the third column shows large-imageseparation systems (2 5).The right panel consists of quads: the first column shows systems that require deblending and the second column shows systems that do not.

Figure 4 .
Figure4.Distributions of g-, r-, and z-band magnitudes of the images of new candidates and "discoverable confirmed" systems.Note that the distributions of our candidates are fainter than the "discoverable confirmed" systems.

Figure 5 .
Figure 5. Grade A, B, and C candidates overlaid on the DESI Legacy Surveys (see Figure 1) are shown as red, orange, and yellow dots, respectively (those without Gaia data as triangles; see Section 4.1).Note that the 94 "discoverable known" systems are not included.Candidates with spectroscopic redshifts from SDSS DR16 (see Section 4.2.2) are outlined by a black diamond.

Figure 6 .
Figure 6.Ten visually impressive candidates.The naming convention is R.A. and decl. in decimal format.The top right corner of each image indicates its grade.North is up, and east is to the left.DESI-118.5478+ 10.4849 is a potential quad and the other nine are doubles.Note that observation by the blue channel of SNIFS did not detect quasar features in DESI-090.5672-43.5945,possibly due to the low S/N of the spectrum.

Figure 7 .
Figure7.The first 80 of our newly discovered candidates by ascending R.A., grouped by grade.For the naming convention, see the Figure6caption.The top right corner of each image indicates its grade.North is up, and east to the left.Additional information about each system can be found in Table2.

Figure 11 .
Figure 11.The redshift distributions for the lenses and quasars of our candidate systems.The redshifts of our quasars are from SDSS DR16 and, for 17 additional systems, from the SNIFS instrument on the University of Hawaii 88 inch (2.2 m) telescope.Photometric redshifts for the lenses are from Zhou et al. (2020).

Table 1
Lensed Quasar System Types

Table 2
Eighty of the New Candidates Arranged in Ascending R.A. (Also See Figure7)