When Spectral Modeling Meets Convolutional Networks: A Method for Discovering Reionization-era Lensed Quasars in Multiband Imaging Data

In this work, we use the Dark Energy Survey (DES) Data Release 2 (Abbott et al. 2021) as the primary catalog. This survey is based on optical imaging collected using the Dark Energy Camera (DECam; Honscheid & DePoy 2008; Flaugher et al. 2015) on the Blanco 4 m telescope at the Cerro Tololo Inter-American Observatory, Chile. DES covers ∼5000 deg² of the southern Galactic cap in five bands, with a median point-spread function (PSF) FWHM of g, r, i, z, Y = 111, 095, 088, 083, and 090, respectively. The median coadded catalog depth is g_DES = 24.7, r_DES = 24.4, i_DES = 23.8, z_DES = 23.1, and Y_DES = 21.7 mag for a 195 diameter aperture and a signal-to-noise ratio (S/N) of 10. By default, we use the aperture-based magnitude (MAG_AUTO), where the aperture size varies depending on the extent of each source, as reported in the main table of DES. As an additional note, DES tile images have a pixel scale of 0263 and a fixed zeropoint of 30 mag.

We then complement our primary data with near-IR and mid-IR (MIR) photometry from the Visible and Infrared Survey Telescope for Astronomy Hemisphere Survey (VHS) Data Release 5 (McMahon et al. 2013, 2021) and the Wide-field Infrared Survey Explorer (unWISE version; Wright et al. 2010; Schlafly et al. 2019) catalogs, respectively, using a 2'' crossmatching radius. The VHS Petrosian J-magnitude (Jpmag) and unWISE photometric bands (W1 and W2) are extremely useful for distinguishing between high-z quasars, low-z galaxies, and MLT dwarfs (e.g., Andika et al. 2020, 2022). This complementary data is also an effective way of removing spurious sources, such as diffraction spikes or moving objects, which are often found in one survey, but not in the others. Finally, we correct all of the photometric measurements from Galactic reddening by making use of the dust map from Bayestar19 ¹² (Schlafly et al. 2019), following the Fitzpatrick (1999) extinction relation. As a side note, we crossmatch our main catalog with a list of MLT dwarfs, employing a 2'' crossmatching radius (Best et al. 2018; Carnero Rosell et al. 2019) and quasars (Flesch 2021) to keep track of the known objects.

We proceed to model the target sources by generating 900 mock quasars. Here, we distribute them in redshifts of 5.6 ≤ z ≤ 7.2 and rest-frame 1450 Å absolute magnitudes of −30 ≤ M₁₄₅₀ ≤ −20, following the quasar luminosity function from Matsuoka et al. (2018b), in the form of:

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Phi }}({M}_{1450},z)\\ \,=\,\displaystyle \frac{{10}^{k(z-6)}{{\rm{\Phi }}}_{* }}{{10}^{0.4({\alpha }_{\mathrm{qso}}+1)({M}_{1450}-{M}^{* })}+{10}^{0.4({\beta }_{\mathrm{qso}}+1)({M}_{1450}-{M}^{* })}}.\end{array}\end{eqnarray} \tag{ 1 }$$

We adopt faint- and bright-end slopes of α_qso = −1.23 and β_qso = −2.73, a redshift evolution term of k = −0.7, a break magnitude of M^∗ = −24.90, and normalization of Φ^∗ = 10.9 (Matsuoka et al. 2018b).

We then produce the aforementioned quasar spectra, following the Yang et al. (2016) prescription, implemented with the SIMQSO ¹³ simulation code (McGreer et al. 2013). SIMQSO is designed to produce synthetic quasar spectra that match ∼60,000 Sloan Digital Sky Survey (SDSS; York et al. 2000; Eisenstein et al. 2011; Dawson et al. 2013, 2016; Blanton et al. 2017) quasars at 2.2 < z < 3.5. The code also assumes that the quasar SEDs do not evolve significantly with redshift (Shen et al. 2019). However, we note that z ≳ 6 quasars frequently exhibit high-ionization broad lines with more dramatic velocity shifts than other z ≲ 5 quasars of similar luminosity (e.g., Mazzucchelli et al. 2017; Meyer et al. 2019; Schindler et al. 2020). Nonetheless, a large number of low-z SDSS quasars still provide a good reference for building a high-z quasar composite spectrum.

The primary component of our quasar spectral model is continuum emission constructed with a power-law function with four breaks. In the rest wavelength range of 1100–5700 Å, the corresponding continuum slope α_ν follows a Gaussian distribution with a mean of μ(α_ν ) = −1.5. Meanwhile, for the wavelength ranges of 5700–10850 Å, 10850–22300 Å, and > 22300 Å, we change the slopes to μ(α_ν ) = −0.48, −1.74, and −1.17, respectively. We use a dispersion value of σ(α_ν ) = 0.3 for each of the aforementioned slopes, following Yang et al. (2016).

Next, we add a series of UV-to-optical Fe emissions to the spectral model, based on the Vestergaard & Wilkes (2001), Tsuzuki et al. (2006), and Boroson & Green (1992) templates for rest-frame wavelengths of < 2200 Å, 2200–3500 Å, and 3500–7500 Å, respectively. After that, the emission lines from the broad- and narrow-line regions are appended to the model, following the equivalent widths and FWHM line distributions of SDSS quasars. Furthermore, the IGM absorption by the Lyα forest is also implemented in the simulated spectra (Songaila & Cowie 2010; Worseck & Prochaska 2011). On top of this, we add the Lyα damping wing model, following the Miralda-Escudé (1998) formalism, assuming a fixed value of 3 Mpc for the proximity zone size and a randomly drawn neutral hydrogen fraction value in the range of 0%–10% (e.g., Euclid Collaboration et al. 2019; Eilers et al. 2020; Andika et al. 2022). Finally, internal reddening was applied to each mock quasar spectrum, using the Calzetti et al. (2000) dust model, with a randomly chosen E(B − V) value in the range of −0.02 to 0.14. We should note that the negative reddening values are to account for quasar models with bluer continua than encompassed by the original templates. After this, the photometry is calculated from the simulated spectra, and the corresponding uncertainties are derived from the observed magnitude and error relations taken from each survey (e.g., Yang et al. 2016).

To obtain the deflectors for the lens modeling, we first search for a sample of spectroscopically confirmed galaxies using the SDSS Data Release 17 (Abdurro'uf et al. 2022) catalog, via the CasJobs ¹⁴ webpage. We retrieve all sources classified as "GALAXY" by the SDSS pipeline and limit our search to those having reliable velocity dispersion measurements σ_v with less than 30% error, which is one important parameter for calculating the lensing effect later. Moreover, we only select the galaxies with redshifts of 0.3 < z < 2.0, considering that the majority of the lensing optical depth for z ≳ 6 sources comes from z ≲ 1.5 early-type lens galaxies. (Wyithe et al. 2011; Mason et al. 2015; Pacucci & Loeb 2019; Yue et al. 2022a). After this, we crossmatch these sources to the DES, VHS, and unWISE catalogs, with a searching radius of 1'', to obtain their corresponding optical/IR magnitudes, when available. In the end, we obtain a sample of 109,808 galaxies and query their images accordingly via the DESaccess ¹⁵ webpage.

As the next step of our method, we need to obtain the spectra of the deflectors that will later be used for constructing the mock galaxy–quasar lens SEDs. Based on the multiband photometric information of the aforementioned galaxies, we proceed to model their spectra using the Bagpipes ¹⁶ code (Carnall et al. 2018). Several priors on the associated physical parameters need to be supplied to Bagpipes. In this case, we first parameterize the star formation histories (SFHs) of the galaxies using a double-power-law function:

$$\begin{eqnarray}&&\mathrm{SFR}({\rm{t}})\ \propto \ {\left[{\left(\displaystyle \frac{{\rm{t}}}{\tau }\right)}^{{\alpha }_{\mathrm{gal}}}+{\left(\displaystyle \frac{{\rm{t}}}{\tau }\right)}^{-{\beta }_{\mathrm{gal}}}\right]}^{-1},\end{eqnarray} \tag{ 2 }$$

where the falling (α_gal) and rising (β_gal) slopes range between 0.01 and 1000 in logarithmic space (Carnall et al. 2019). We then vary the peak times of star formation (τ) between 0 Gyr and the universe's age at the observed redshifts. Next, the prior on the SFH normalization is varied in the range of M_formed = 10–10¹⁵ M_*/M_⊙, where M_* is the stellar mass. Also, we vary the galaxy metallicity in the range of Z/Z_⊙ = 0.0–2.5, where Z_⊙ = 0.02 is solar metallicity.

Following Carnall et al. (2019), we add the Charlot & Fall (2000) dust attenuation model, in the form of A_λ ∝ λ⁻ⁿ . Here, we impose a Gaussian prior on the slope of the attenuation curve n, with a mean and standard deviation of 0.3 and 1.5, respectively. The attenuation in the V-band A_V is set to range from 0 to 8. In addition to this, we utilize a constant value for the maximum lifetime of stellar birth clouds, of t_BC = 0.01 Gyr, and for the ratio of attenuation between stellar birth clouds and the wider interstellar medium, of = 2. It should be noted that the priors on the galaxy photometric redshifts are set to match the ones derived from SDSS-based spectroscopy, with slight variation (i.e., Δz = ± 0.015). Finally, using the σ_v information from the SDSS catalog, we convolve the spectral model with a Gaussian kernel in velocity space. As a result, we produce a sample of SEDs of galaxies ranging from UV to IR wavelengths, and their associated photometry. An example of a lens galaxy fitted with the Bagpipes code is presented in Figure 1.

Zoom In Zoom Out Reset image size

Finding lensed quasar candidates using catalog-level photometry requires knowledge of the shape of their SEDs—i.e., the addition of foreground galaxy and background quasar emissions. Here, we combine the simulated quasars and galaxies described in the previous sections to produce the mock lensed quasar spectra.

To describe the mass profile of the lenses, we assume a singular isothermal sphere (SIS; Schneider 2015) model. The Einstein radius can be inferred from the 1D stellar velocity dispersion σ_v in the potential of the mass distribution, using the formula

$$\begin{eqnarray}&&{\theta }_{{\rm{E}}}=4\pi \displaystyle \frac{{\sigma }_{v}^{2}}{{c}^{2}}\displaystyle \frac{{D}_{\mathrm{ds}}}{{D}_{{\rm{s}}}},\end{eqnarray} \tag{ 3 }$$

where the speed of light is c, while D_ds and D_s are the angular diameter distances between the deflector and the source and the observer and the source, respectively.

It is important to note that galaxy mass distributions in nature are not perfectly symmetric. The symmetry-breaking is often caused by the ellipticity of the mass distribution or external shear forces—e.g., the tidal gravitational fields of adjacent galaxies. Thus, this will alter the SIS lens characteristics, and might produce more than two images. However, we refrain from using a more complex lens model, like the singular isothermal ellipsoid (SIE; Barkana 1998) model, which usually requires accurate measurements of deflector axis ratios and position angles. Nevertheless, SIS appears to reproduce lens systems quite well, and the typical image separation remains in the same order of magnitude as predicted by the SIE model (Schneider 2015). For details of the SIS lens calculation, we refer the reader to Appendix A.

We then produce the mock lensed quasar spectra by pairing each simulated galaxy with a randomly chosen quasar model created beforehand. Assuming that β is the true angular position of the source, the quasar is randomly placed in a specified region behind the lens, within 001 ≤ β ≤ 15. After this, the source image is projected onto the lens plane, where the deflection angle and magnification are calculated based on the corresponding lens configuration. We accept the mock lens if it contains a strong lensing effect, with a magnification factor of μ ≥ 2 within θ_E ≤ 1'', and the lensed quasar flux has >5σ detection in the DES Y band. Our choice of the β and θ_E ranges is motivated by Pacucci & Loeb (2019), who predict that a significant fraction of quasars at z ≳ 6 are strongly gravitationally lensed, with small image separations—i.e., Δθ = 2θ_E ≤ 1''.

We show in Figure 2 the distributions of the lens galaxy redshifts, velocity dispersions, Einstein radii, and i-band magnitudes used for the simulation. In terms of redshift, the number of galaxies peaks at z ≈ 0.5. Then, for the i-band magnitudes, we see an increase in the deflector numbers up to i_DES ≈ 19.5, before a rapid decrease toward the faint end. Hence, our training data set is skewed toward brighter, more massive lens galaxies. This phenomenon is primarily caused by how SDSS selects target galaxies for its spectroscopic surveys, following the criteria established by Dawson et al. (2013) and Prakash et al. (2016), to study the universe's large-scale structure. Most of the target galaxies are luminous early-type galaxies at z < 1, which are excellent probes for constraining the baryon acoustic oscillation signal and, subsequently, the universe's expansion rate (e.g., Woodfinden et al. 2022; Zhao et al. 2022). Note that the faint limits for galaxies chosen for spectroscopic observations are i = 19.9 for SDSS III, extending to i = 21.8 for SDSS IV (Prakash et al. 2016).

Zoom In Zoom Out Reset image size

Finally, of the initial 109,808 foreground galaxies and 900 background quasars, we obtain 20,823 surviving lens configurations that meet our criteria. Naturally, the lensed quasar spectra can be obtained simply by adding the fluxes of the foreground galaxy and the magnified background quasar, as presented for one example in Figure 3. However, in some cases, the lens galaxy could reside far away from the quasar image, resulting in two close-by objects, so that the photometry extracted from the quasar image only has partial coverage of the lens galaxy flux. To take this into account, we produce more SED templates, by arbitrarily scaling the lens galaxy flux contribution in the range of 10%–100%, before adding it to the magnified quasar flux.

Zoom In Zoom Out Reset image size

As the next step in our lensed quasar search methodology, we select sources detected in the DES Y band that have counterparts in the VHS J band and the unWISE W1 band, within a searching radius of 2''. Note that there are some cases where the magnitude values of the other bluer DES bands—i.e., g, r, i, and z—are not available, or the sources are simply not detected in those bands, i.e., with S/Ns < 3. In this case, we replace the catalog magnitudes with their associated 3σ upper limits, where we infer these limits based on the measured flux uncertainties reported in the DES table. We then apply the following flag and color cuts:

$$\begin{eqnarray}&&\mathrm{DES}\,{\mathtt{IMAFLAGS}}\_{\mathtt{ISO}}\ (g,r,i,z,Y)=0,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\mathrm{DES}\,{\mathtt{FLAGS}}\ (g,r,i,z,Y)\lt 4,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\rm{S}}/{\rm{N}}({Y}_{\mathrm{DES}})\geqslant 5,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{z}_{\mathrm{DES}}-{Y}_{\mathrm{DES}}\gt 0.5,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&-0.6\lt {Y}_{\mathrm{DES}}-{J}_{\mathrm{VHS}}\lt 1.0,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&-0.5\lt {Y}_{\mathrm{DES}}-{\rm{W}}1\lt 2.5.\end{eqnarray} \tag{ 9 }$$

Equations (4) and 5 are the primary flags for selecting well-behaved objects—i.e., there are no problems in the source extraction process, with missing pixels, etc.—from the DES photometry catalog (Abbott et al. 2021). These criteria also ensure that our targets will have catalog entries in all DES bands. On the other hand, the criteria in Equations (6)–(9) are determined empirically, based on the photometric features of the mock lenses created in Section 2.4.

As seen in Figure 4, there is a substantial overlap between the lenses and contaminants in the color space where z_DES − Y_DES ≤ 0.5. Thus, we choose to limit our color cut to z_DES − Y_DES > 0.5, which consequently focuses our search on the lens systems where the background quasars are located, at z = 6.6–7.2. More details of the quasar selection function in our lensed quasar search criteria can be found in Appendix B. To quantify our selection completeness, we start by applying a z_DES − Y_DES > 0.5 color cut to the remaining 20,823 simulated lenses produced in the previous section. As a result, only 3657 of them pass this criterion. We then use this number of remaining mock lenses as the basis for estimating our selection completeness in recovering lens systems having quasars at z = 6.6–7.2. Further cuts are made by applying the criteria in Equations (6)–(9), which results in 3650 surviving mock lenses, equivalent to around 99% completeness. This choice also rejects a substantial fraction (≳70%) of the contaminants—i.e., low-z galaxies and MLT dwarfs. Therefore, these are the candidates on which we will conduct SED modeling for further selection.

Zoom In Zoom Out Reset image size

The next part of our preselection method involves SED modeling, using four different spectral templates: (i) the lensed quasars derived from the previously created mock lenses; (ii) the elliptical galaxies at z ≤ 3; (iii) the MLT dwarfs; and (iv) the z ≥ 4 unlensed quasars. It is crucial to note that templates (ii)–(iv) are empirically derived from observations and described in detail in Appendix C. The main goal of this SED modeling is to calculate the probability of each target being a lens or a contaminant, as well as its associated photometric redshift. The SED modeling is implemented using the EAZY photometric redshift code (Brammer et al. 2008). Here, EAZY will step through a grid of spectral templates, fit them to the photometry of the candidates, and try to find the best SED templates. We choose as solutions the best models with the smallest reduced chi-square (${\chi }_{\mathrm{red}}^{2}$), which can be calculated for each template i as

$$\begin{eqnarray}&&{\chi }_{\mathrm{red},i}^{2}=\displaystyle \sum _{n=1}^{N}{\left(\displaystyle \frac{{\mathrm{data}}_{n}-{f}_{n}({\mathrm{model}}_{i})}{\sigma ({\mathrm{data}}_{n})}\right)}^{2}/(N-1),\end{eqnarray} \tag{ 10 }$$

where the number of photometric data points and degrees of freedom are denoted with N and (N − 1), respectively.

The sources with a high probability of being lensed quasars are chosen based on the calculated ${\chi }_{\mathrm{red}}^{2}$ of the lens (${\chi }_{\mathrm{red},{\rm{L}}}^{2}$), galaxy (${\chi }_{\mathrm{red},{\rm{G}}}^{2}$), MLT dwarf (${\chi }_{\mathrm{red},{\rm{D}}}^{2}$), and quasar (${\chi }_{\mathrm{red},{\rm{Q}}}^{2}$) templates, along with their associated ratios. In addition to this, we also use the estimated quasar photometric redshift z_qso to search specifically for z ≳ 6 unlensed quasars. Hence, the following criteria are employed to find high-z lensed and unlensed quasar candidates:

$$\begin{eqnarray}&&\displaystyle \frac{{\chi }_{\mathrm{red},{\rm{L}}}^{2}}{{\chi }_{\mathrm{red},{\rm{D}}}^{2}}\lt 0.2\quad \mathrm{and}\quad \displaystyle \frac{{\chi }_{\mathrm{red},{\rm{L}}}^{2}}{{\chi }_{\mathrm{red},{\rm{G}}}^{2}}\lt 0.2,\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{or}\\ \ \ {z}_{\mathrm{qso}}\gt 6.0\quad \mathrm{and}\quad \displaystyle \frac{{\chi }_{\mathrm{red},{\rm{Q}}}^{2}}{{\chi }_{\mathrm{red},{\rm{D}}}^{2}}\lt 0.3.\end{array}\end{eqnarray} \tag{ 12 }$$

Here, the values in Equation (11) are derived by fitting the simulated photometric data of the mock lenses created in Section 2.4, which are empirically optimized to maximize the purity and completeness of our lensed quasar selection. On the other hand, the criteria in Equation (12) are inferred empirically, by modeling the SEDs of known MLT dwarfs and z ≳ 6 quasars in the DES catalog (Carnero Rosell et al. 2019; Flesch 2021). Finally, applying the SED fitting to the previously surviving mock lenses yields 3650 remaining objects, corresponding to a completeness of around 99% for selecting z = 6.6–7.2 lensed quasars. An example of a lens candidate selected using our SED modeling is presented in Figure 5, and a summary of our selection steps is reported in Table 1.

Zoom In Zoom Out Reset image size

Since only one galaxy–quasar lens system has been found at z ≳ 6 (Fan et al. 2019), we need to mock up additional lens images. To do this, we adopt a data-driven approach to construct the training data set, which is outlined as follows. First, we create images of mock lenses based on the deflector galaxies, simulated quasars, and lens configurations used in Section 2.4. Then, through the lens equation whose parameters have been defined, we paint lensed point-source lights over the actual galaxy images, to construct realistic galaxy–quasar lens models. In this case, the lensing effect and the ray-tracing simulation are generated according to the PyAutoLens ¹⁷ lensing code (Nightingale et al. 2021). The resulting deflected quasar lights are then convolved with a Gaussian PSF model at the location of the lens and coadded with the original DES galaxy images. As a reminder, we accept the mock image if it contains a strong lensing case with a configuration of 001 ≤ β ≤ 15 and μ ≥ 2, within θ_E ≤ 1''. To better illustrate this, we present examples of mock color images that have been created by combining the izY-band data in Figure 6, and we show the employed simulation workflow in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We categorize the input data for training the CNN classifier into four classes: (i) the galaxy–quasar lenses created in the previous section; (ii) the DES galaxies that are not selected in the lensing simulation; (iii) a sample of MLT dwarfs, randomly drawn from the Carnero Rosell et al. (2019) catalog; and (iv) a sample of spectroscopically confirmed quasars at z ∼ 1.5–7.0, from the Flesch (2021) database. Classes (i) and (ii) each contain 10,000 sources, while classes (iii) and (iv) each consist of 8000 objects. Hence, 36,000 sources in total are used for the CNN inputs.

The images for the training inputs are based on grizY bands of DES cutouts, having a size of 48 × 48 pixels—i.e., equivalent to angular sizes of 126 × 126. These images are then min–max normalized, so the fluxes range between 0 and 1. Note that the relative fluxes between bandpasses are conserved, so the spectral colors of the corresponding sources are preserved. After this, we apply data augmentation to the images through random ± π/2 rotations, 4 pixel translations, and horizontal or vertical flips. This technique subsequently increases the amount of training data, while enhancing the likelihood that the network will correctly classify multiple orientations of the same image.

Inspired by classical CNN architectures (e.g., LeCun et al. 1989; Sultana et al. 2019), we start building the network with three convolutional layers, followed by one fully connected layer of 128 neurons (see Figure 8). The convolutional layers each have a stride of 1 × 1, the "same" padding, and a kernel size of 3 × 3 × C, with C = 32, 64, and 64 for the first, second, and third layers, respectively. For each convolutional layer, we also add a max pooling layer with of size 2 × 2, a stride of 2 × 2, and the same padding. The dropout regularizations are utilized on both the convolutional (drop rate = 0.2) and fully connected layers (drop rate = 0.5). We also add L2-norm regularization, with a weight decay of 10⁻⁴. The learning rate is set to 10⁻³ initially, while the weight and bias of each neuron are generated randomly and then updated throughout the training. The rectified linear unit functions are employed everywhere, except for the final layer, which utilizes the softmax activation. After passing through all of the convolutional layers, the data cube is flattened and processed by the fully connected layer to obtain the four outputs—i.e., the probabilities of a candidate being a lensed quasar, a galaxy, an MLT dwarf, and a normal quasar. All the training processes and the CNN modeling are implemented using the TensorFlow ¹⁸ deep learning framework (Abadi et al. 2016; Developers 2022).

Zoom In Zoom Out Reset image size

The specific CNN architecture mentioned above was chosen after we employed the Hyperband algorithm to optimize the model hyperparameters (Li et al. 2018). To rapidly search for a model with high performance, the Hyperband tuner applies early stopping and the adaptive resource allocation method, utilizing a tournament bracket approach, similar to what is used in a sports competition. In principle, the algorithm trains several CNN architectures with randomly sampled configurations over a few cycles. In other words, we explored various models with different numbers of convolutional layers, kernel sizes, and dense units. We also tested the effects of different learning rates, in the range of 10⁻²–10⁻⁴, dropout rates of 0.1–0.5, convolutional filter numbers of 16–256, and L2-norm weight decays of 10⁻¹–10⁻⁵ on the performance of our CNN models. Ultimately, only the top half of best-performing models is carried forward to the next round, until the best-performing model is obtained. When compared to Bayesian optimization, random search, or grid search, this strategy maximizes the number of evaluated CNN configurations, and results in more effective resource consumption.

To assess network performance, we divide the images into three data sets: training (60% of the data), validation (20%), and testing (20%). Those data sets are further divided into random batches, with sizes of 256. Here, we use sparse categorical crossentropy ¹⁹ as the loss function. For each iteration, the networks predict the outputs for one batch (forward propagation), then continue the loops through all the batches from the training data set, to complete one epoch. The information about the current batch loss is then propagated back, to update the weights and biases of the neurons, which is performed by following a stochastic gradient descent algorithm to minimize the loss (i.e., the Adam optimizer; Kingma & Ba 2014). This optimization is performed for each epoch, initially for all of the training data set batches. The average loss for the whole training data set is then obtained. After this, the procedures are repeated for the validation data set, but without parameter updates to the neurons. In this way, the average loss for the validation data set is retrieved. By comparing the training and validation losses, we can check whether the model optimization is improved or whether overfitting happens. Note that overfitting occurs when a network fails to capture the features in the training data set and generalizes them for unseen data. Therefore, we randomly reorder our training data after each epoch, to improve generalization and produce networks with optimum accuracy. In the end, we can determine the optimum number of training cycles by locating the epoch with the lowest average validation loss over numerous runs.

For each image tensor that is fed to our CNN model, the classifier will provide probability scores of it being a lensed quasar, a galaxy, an MLT dwarf, and a normal quasar—i.e., P_lens, P_galaxy, P_dwarf, and P_quasar. The predicted class is then assigned by choosing which of the classes has the highest probability score. A value of P_lens = 1 means that there is a high chance that the classified image contains a lensed quasar. In contrast, P_lens = 0 means that the image is not a lensed quasar and is more similar to other contaminants. Note that the softmax activation function that we use in the last layer of our CNN model would produce P_lens + P_galaxy + P_dwarf + P_quasar = 1, by construction. In the end, our CNN training process converged after 82 epochs, obtaining 96.1% (96.7%) accuracy for the evaluation using the training (validation) data set, along with a loss value of 0.145 (0.131).

From a classical standpoint, the high accuracy attained in the training data set might have been caused by overfitting. Thus, we compare the accuracy–loss learning curves that were calculated by evaluating the CNN predictions on the training and validation data sets (see Figure 9). After decreasing for several epochs, the training and validation loss values stabilize and follow the same pattern. The lack of any overfitting signals—i.e., the training loss continues to decrease, but, in contrast, the validation loss begins to increase, after many epochs—makes us confident that our CNN classifier is capable of learning and generalizing.

Zoom In Zoom Out Reset image size

Another method for assessing the overall performance of the trained model is to use the receiver operating characteristic (ROC) curve and to calculate its corresponding area under the curve (AUC). In principle, this shows how well a binary classifier differentiates between two classes when the decision threshold is changed. Thus, we first define the positives (P) as the lenses and the negatives (N) as the nonlenses or contaminants. True positives (TP) are cases where the model correctly predicts the lenses, while true negatives (TN) occur when the nonlenses are correctly identified. On the other hand, false positives (FP) are where the classifier incorrectly predicts contaminants as lenses. In addition, we also define false negatives (FN), as cases where the model incorrectly rejects lenses. The ROC curve presents the false-positive rate (FPR) against the true-positive rate (TPR) for the validation data set, where

$$\begin{eqnarray}&&\mathrm{TPR}=\displaystyle \frac{\mathrm{TP}}{{\rm{P}}}=\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}\end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray}&&\mathrm{FPR}=\displaystyle \frac{\mathrm{FP}}{{\rm{N}}}=\displaystyle \frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{TN}}.\end{eqnarray} \tag{ 14 }$$

We then construct the ROC curve by progressively raising the probability threshold from 0 to 1. A perfect classifier will result in AUC = 1, while a classifier that only forecasts randomly will have AUC = 0.5. It is important to note that because we have four classes for categorizing the candidates—i.e., a multilabel classification—we need to binarize our CNN output using the "one versus all" methodology. Accordingly, we produce four ROC curves and present them in Figure 10, consisting of: (i) lensed quasars against galaxies and other point-source contaminants, shown by the solid blue line; (ii) galaxies against lens systems and other contaminants, shown by the dashed magenta line; (iii) MLT dwarfs against other sources, shown by the dashed yellow line; and (iv) quasars against lenses, galaxies, and others, shown by the dashed cyan line. Note that we created these ROC curves based on the CNN classifier evaluation of the previously unseen test data set, which resulted in high AUC values, indicating excellent performance.

Zoom In Zoom Out Reset image size

Based on the ROC curves, we use the geometric mean or G-mean ²⁰ metric to seek a balance between the TPR and FPR. The highest G-mean score indicates the best threshold of P_lens that maximizes TPR while minimizing FPR. In this case, we obtain a recommended threshold of P_lens > 0.13, which produces FPR = 0.01 and TPR = 0.96. However, we then decide to relax this probability threshold, instead using P_lens > 0.1 to classify the candidates as lenses, which results in FPR = 0.02 and TPR = 0.96. This choice is made to increase the final number of lensed quasar candidates, by accommodating candidates with lower lens probabilities. Below the aforementioned P_lens threshold, the number of candidates will grow exponentially, while their quality deteriorates, making visual examination more time-consuming and less effective. Without a reference sample of z ≳ 6 lensed quasars, it is impossible to determine the perfect threshold. In terms of the trade-off between completeness and purity, we need to find a balance where the number of candidates is manageable for follow-up observations.

As supplementary information, we performed an additional evaluation using an independent data set. This test data set was assembled based on the list of known z ≲ 3 lensed quasars taken from the Gravitationally Lensed Quasar Database ²¹ (Inada et al. 2012; Agnello et al. 2015, 2018; Lemon et al. 2018, 2019, 2020; Spiniello et al. 2018; Jaelani et al. 2021). Of the 220 lenses in the database, 34 of the known lenses have DES images. Combined with the contaminants generated in Section 3.2, we tested our CNN classifier on this new data set and successfully recovered 27 lenses—i.e., equivalent to a completeness of 79% and a purity of 11%. We refer the reader to Appendix D for more details on the resulting distribution of probability scores. Compared to the evaluation against the mock lenses, the performance of our CNN model decreases when it is tested against unseen real lens systems. This decreased performance might have been caused by the uniqueness of some of the strong lenses that were not covered by our simulation. In particular, some of the FNs in the test data set might have contained lensed arcs that were too faint to be identified, compact lenses, multiple distortions caused by substructures, or other things. Nonetheless, our CNN classifier can generalize and achieve a sufficiently high accuracy for our purpose. We should emphasize that the point of our CNN-based classification is to conduct a preselection before we go into the visual inspection process, not to assemble a final sample that is quantitatively pure or complete.