© K. Rojas et al. 2022

1. Introduction

Gravitational lensing is the phenomenon by which light rays are deflected by a gravitational field. In so-called “strong” lens systems, it is possible to observe multiple images, arcs, or rings of a distant source around a foreground galaxy, group, or cluster. In such cases, these systems can serve as important tools in the study of diverse and fundamental questions about the Universe. Some examples are the study of luminous and dark matter components of the deflector (Kochanek & Dalal 2001; Oguri et al. 2002; Davis et al. 2003; Jiménez-Vicente et al. 2015), measuring the Hubble constant H₀ using time delays (Falco et al. 1997; Vuissoz et al. 2007; Bonvin et al. 2017; Wong et al. 2020; Millon et al. 2020) and constraining the dark energy equation of state (Biesiada et al. 2010; Collett & Auger 2014; Cao et al. 2012, 2015). However, most of these applications are limited by the paucity of known systems, as only a few hundred such systems are confirmed. Therefore, an effort on the discovery and confirmation of more lenses is required.

Since the serendipitous discovery of the first lensed quasar (Walsh et al. 1979), new discovery methods have been developed on the basis of novel datasets and techniques. Recent searches have included: algorithms based on identifying lens features such as ArcFinder (Alard 2006), RingFinder, which searches for blue features blended with red light (Gavazzi et al. 2014), principal component analyses (PCA) of galaxies to search for lensed features in the residual images using machine learning (Joseph et al. 2014; Paraficz et al. 2016), and CHITAH, which evaluates point source configurations as possible lensed images using lens modeling (Chan et al. 2015). In recent years, the growing amount of available data has motivated the use of more automated techniques such as artificial neural networks (ANNs; Rosenblatt 1957) and, in particular, convolutional neural networks (CNNs; LeCun et al. 1989). These latter techniques are based on supervised machine learning algorithms capable of solving complex problems such as pattern recognition or image classification when a proper training set is provided.

The biggest challenge of using CNNs for lens finding is creating a robust training set that contains diverse lens systems for positive examples and non-lens galaxies, including some that can be mistaken as lenses like spirals, rings, and mergers, as negative examples. We currently lack sufficient numbers of known examples of both lens systems and common false positives. The only solution is to then simulate them as realistically as possible. This has already been addressed in several ways: fully simulating images using analytical profiles for both the lens and source (Jacobs et al. 2019a), using an analytical profile for the source but a real image of the lens Petrillo et al. (2019b), and using real data for both the deflector and background galaxy (Cañameras et al. 2020). The fully analytical approach has the advantage of having full control over all parameters to create a sample as varied as possible, but lacks the ability to mimic features of real images like artifacts, noise, and companion galaxies.

The main differences between previous searches depend on mock simulation methods, use of single or multi-band data, and the architecture design. Previous searches include: the Canada-France-Hawaii Telescope Legacy Survey (CFHTLS; Jacobs et al. 2017), the Kilo Degree Survey (KiDS; Petrillo et al. 2017, 2019a,b; He et al. 2020), the Dark Energy Survey (DES) Year 3 (Jacobs et al. 2019a,b), the Dark Energy Spectroscopic Instrument (DESI) Legacy Imaging Surveys (Huang et al. 2020, 2021), the Pan-STARRS 3π survey (Cañameras et al. 2020), and the VST Optical Imaging of the CDFS and ES1 fields (VOICE survey; Gentile et al. 2021). Overall, these studies have shown that CNNs are a promising tool for listing thousands of new lens candidates. However, these tools rely on a high degree of human visual inspection afterwards in order to compile the final candidate list. Improving the training process with realistic lenses and diverse types of galaxies is therefore of key importance when the next generation of surveys such as Euclid Space Telescope (Laureijs et al. 2011) and the Rubin Observatory Legacy Survey of Space and Time (LSST, LSST Science Collaboration 2009; Ivezić et al. 2019) start producing data. Current lens-finding efforts, including some performed in simulated data (Lanusse et al. 2018; Avestruz et al. 2019), serve as important preparation for their advent, as it is expected that over 100 000 new strong lensing systems will be discovered (Collett 2015) and a visual inspection process is not affordable.

In this work, our main aim is to find new strong-lensing systems. As a secondary goal, we want to characterize a subset of our false positives. We perform our search in the footprint of the Dark Energy Survey (DES; The Dark Energy Survey Collaboration 2005). The data description is presented in Sect. 2. Details of the simulation procedure of our training sets are in Sect. 3. We present the details of both the CNN training and validation in Sect. 4. We grade the parent sample using the CNN model and we select the ones with higher score to perform a dedicated visual inspection to identify the best lens candidates and subclassify false positives. In Sect. 5, we give our visual classification procedure for the best graded candidates from the CNN, presenting both the best lens candidates and a compilation of ring galaxy candidates that will help to improve future lens finding searches. Finally, in Sect. 6, we show the results of an automatic modeling tool on a sample of our best candidates.

2. Data selection

We used data from DES based on the Dark Energy Camera (DECam, Honscheid & DePoy 2008; Flaugher et al. 2015) on the Blanco 4-m telescope at Cerro Tololo Inter-American Observatory (CTIO), Chile. In brief, DECam is a 570 Megapixel camera with a field of view of 2.2 square degrees and a pixel size of 0.27″. Observations were performed in the optical grizY bands. The first DES Data release (DES-DR1, Abbott et al. 2018) contains images taken over the first three years of operation, covering an area of 5186 deg². The images from DR1 have been co-added and each filter re-scaled to have a fixed zero point of 30 mag.

We used the NOAO Data Lab (Fitzpatrick et al. 2016) service to build our sample from the des_dr1.galaxies catalog and selected a sample of luminous red galaxies (LRGs) in order to maximize the lensing cross-section (e.g., Turner et al. 1984). To do so, we applied the following cuts in color and magnitude: 1.8 < g − i < 5, 0.6 < g − r < 3, 18 < r < 22.5, g > 20, i > 18.2, where the magnitude used is the mag_auto column reported in the DES data release. Our color selection is summarized in Fig. 1 and is similar to the one adopted by Jacobs et al. (2019b). However, we slightly widened the g − i range and adopted a brighter magnitude limit in the r-band. This selection allows us to better account for the contamination of the lens light by the bluer color of any putative lensed source, but also increases the probability that other types of galaxies (e.g., mergers, spirals, or ring galaxies) can be selected in the sample. The result is a sample of 18 745 029 galaxies located in 10 388 coadded tiles from DES DR1, which we refer to as the parent sample in this work. We downloaded all the g, r, and i bands of these tiles, and generated cutouts around the galaxies of 50 × 50 pixels, corresponding to ∼13 × 13′′. When a point spread function (PSF) is required for the simulation process and modeling, PSFEx (Bertin 2011) is run on the relevant coadd tile, extracting a model of the PSF from the FITS image, allowing us to retrieve a PSF at any position on the tile.

Fig. 1. Color-color diagram of our parent sample of ∼19 million galaxies. The plot displays where most galaxies in our selection are located. The green shaded area shows the density of galaxies from the parent LRG sample. The green solid lines shows the 1, 2, 3, and 6σ contours. The lens candidates and contaminants from this work are shown in overlay: Sure Lens (purple stars), Maybe Lens (open red triangles), and ring galaxies (open cyan circles). See Sect. 5.3 for details.

3. Simulated galaxy-scale lenses

Since our search for galaxy-scale strong lenses is based on a CNN, we required a training set that mimics as closely as possible both the lenses we want to find and also the common non-lensed systems. When simulating lenses, we adopted an approach where the training set is data-driven, in the sense that both the images of the lenses and the background sources are obtained from real data. The outline of our procedure is as follows. First, we selected a sample of high-redshift background galaxies with available high-resolution imaging and accurately measured colors (Sect. 3.1). Then, we selected a sample of LRGs that will act as deflectors and matched them to the background sources to create pairs of lens-source suitable for simulations (Sect. 3.2). Finally, we lensed the source light using the lens equation, whose parameters (deflection angles) are defined by the lens, to produce a simulated image of a lens system that combines images from both samples (Sect. 3.3).

3.1. Background galaxies selection

Our aim is to obtain realistic source galaxies, imaged at HST-resolution, with color information and a high signal-to-noise ratio (S/N). Such a catalog of background sources has already been compiled by Cañameras et al. (2020). The galaxies in this catalog were selected from the Galaxy Zoo catalog (Willett et al. 2017) and are also included in the COSMOS2015 photometric catalog (Laigle et al. 2016). All the objects categorized as galaxies in these catalogs were picked, with no previous selection in color or magnitude, as the depth is limited by the Galaxy Zoo selection, namely, down to F814W ∼ 23.5. Stars and artifacts were manually removed from the sample, but also extended galaxies and galaxies with nearby companions were removed, leaving a final sample of 52 696 objects. Spectroscopic redshifts were obtained from several follow-up surveys (Lilly et al. 2007; Comparat et al. 2015; Silverman et al. 2015; Le Fèvre et al. 2015; Tasca et al. 2017; Hasinger et al. 2018), as well as photometric redshifts from Laigle et al. (2016) for the ones lacking spectra.

To create high-resolution gri-images of our sources, we created cutouts combining the morphological information from HST/ACS F814W high-resolution images (Leauthaud et al. 2007; Scoville et al. 2007; Koekemoer et al. 2007) and the color information from Hyper Suprime Cam (HSC) ultra-deep stack images (Aihara et al. 2018). The detailed procedure to combine the information from these two surveys is described in Cañameras et al. (2020), who followed the steps described in Griffith et al. (2012). In summary, HST/ACS F814W images were aligned and rescaled as if they were observed in the HSC i-band. These HSC images were then resampled to the resolution of the HST/ACS F814W images and were multiplied by an illumination map obtained by dividing the HST/ACS F814W image by the HSC i-band image. Each galaxy stamp has a size of 10 × 10′′ and a pixel size of 0.03″, that is, the HST resolution and PSF but with the HSC observed colors. The morphology of the source is the same in each band. Since the HST PSF is much sharper than that of the ground-based DES images we do not deconvolve our stamps from the HST PSF, which would introduce noise and possible artifacts.

3.2. Lens-source association

Ideally, we would want spectroscopic redshift and velocity dispersions for each member of our LRG sample but the vast majority are lacking this information. To cope with this limitation, we performed a prediction of those parameters using a simple K-nearest-neighbors (KNN) algorithm, assuming that other galaxies with similar gri magnitudes will also have similar redshifts and velocity dispersions. For a reference data set where colors, redshifts and velocity dispersions are available, the KNN algorithm provides a match between the gri magnitudes and the redshift and velocity dispersion of the galaxies. It then predicts the redshift and velocity dispersions for new data based on the K-objects with the most similar colors. We trained the algorithm with 1 400 000 SDSS galaxies that match the color-magnitude cuts of the parent sample and have redshift and velocity dispersion measurements available. We tested the model on another set of 99 382 spectroscopically-confirmed SDSS galaxies, obtaining the predictions for the parameters from the ten nearest neighbors in the gri color space of the training set. We found that the rms scatter in the predictions was σ_z = 0.06 for the redshift, and σ_vel = 69 km s⁻¹ for the velocity dispersion. Finally, we used this model to predict the most likely redshift and velocity dispersion for each of our galaxies in the parent sample. The distributions of the predicted redshifts and velocity dispersions are shown in Fig. 2.

Fig. 2. Redshifts and velocity dispersion distributions used for simulations. Top panel: redshift distributions of the lenses (red) and sources (blue) in the simulated training set. Bottom panel: lens velocity-dispersion distribution, predicted from the K-nearest neighbor algorithm (red) and shifted to match the criteria of pairing lenses and sources (orange), as described in Sect. 3.2. We note that the actual lens velocity dispersions used in the simulations are shifted to higher values, so that lensing features can be seen even at the DES resolution.

We then paired the LRGs with source galaxies, requiring that our simulations have a uniform distribution in Einstein radii spanning 1.2′′ < θ_E < 3.0′′. We chose a conservative lower limit on θ_E because we noticed that, given the average seeing in the gri bands of 1.12″, 0.97″, and 0.88′′ respectively (Abbott et al. 2018), simulations with θ_E < 1.2 create lensing features that are too close or blended with the lens galaxy, and can easily be mistaken as non-lenses, that is, galaxies with extended disks. To evaluate θ_E we use the redshift of the source galaxy, the redshift of the lens, and the velocity dispersion of the lens. To match a source to each LRG, we first take a random lens galaxy from the parent sample and compute the Einstein radii for all source galaxies. We then formed lens-source pairs that produce Einstein radii falling only within our desired bounds and filling bins that produce an uniform distribution. In case where no lens-source pairs satisfied the θ_E conditions, we artificially increased the velocity dispersion of the lens galaxy up to 1.5 times its original value. If still no pair satisfied our criteria, we discarded the LRG. We note, as illustrated in Fig. 2, that this results in a high bias of the velocity dispersion distribution of the lenses. Although this procedure tends to produce lenses with dark matter halos larger than the predicted from the actual galaxy velocity dispersion, this ensures that the lensing features are clearly noticeable to the CNN.

Finally, we enforce the final Einstein radii distribution to be uniform. In other words, our training set is not representative of the true distribution of Einstein radii on the sky, but gives equal probability to all possible values, allowing for more discriminating power in our trained CNN.

3.3. Lensing simulation

We then combined these components to create realistic images of lenses. We adopted the singular isothermal ellipsoid (SIE) as our lensing mass model, which is defined by the Einstein radius (θ_E), the position angle (PA), and axis ratio converted into a complex ellipticity (e₁, e₂), and the central position (x₁, x₂). As mentioned before, the Einstein radius follows a uniform distribution, while the other parameters are acquired individually according to the light distribution of each lens galaxy. Our simulation toolbox uses the Python package Lenstronomy¹ (Birrer et al. 2015; Birrer & Amara 2018). The first step is to determine a simple but realistic representation for the mass of the lensing object. The Einstein radius is calculated using the lens and source redshifts as well as the lens velocity dispersion of the lens derived in Sect. 3.2. The ellipticity parameters and mass centroid are estimated by fitting an elliptical Sérsic profile to the DES r-band image of the LRG. We optimized the fitting procedure using 50 iterations of particle swarm optimization (PSO; Kennedy & Eberhart 1995), with 50 particles. This simple model provides us with parameters for a mass distribution that broadly follows the light distribution of the brightest object in the image. We note that limiting to 50 iterations results in ellipticity parameters that are not perfect and that naturally mimic the effect of deviations of the dark matter profile with respect to the light, without introducing extra complexity in the simulation pipeline. Even though it could produce a few lenses with exotic properties, for instance, very elliptical mass profiles or unusually large dark matter halos, it was found to be adequate for our goal of building realistic simulations in the vast majority of the cases.

The second step is to deflect the light rays from the source according to the lensing mass model. To ensure we can distinguish the final lensed source features against the lens galaxy light, we first increased the original source brightness by one magnitude. To decide where in the source plane our background galaxy is located, we selected a random position inside a square that encloses the caustic curves (curves that mark the location of the maximum magnification and delimit the region inside which a source will be multiply-imaged). Then, we performed a ray-tracing simulation to map the source image onto the image plane and we further convolved the resulting lensed source with the relevant stamp PSF. To convert this image into the DES characteristic pixel resolution, we down-sampled the pixels from 0.03′′ (HST) to 0.27′′ (DES), and re-scaled the flux to match the DES zero points in each filter. As a last step, we added the convolved, resized, and flux-normalized image of the lensed source to the original image of the LRG lens. The latter has, by construction, the right DES PSF and noise properties. Thus, our simulations preserve the characteristics of the original image, such as background noise, seeing, the presence of artifacts, and neighbouring galaxies or stars in the field of view.

To build the multi-band gri simulations, we used the same mass model for all bands, with its parameters derived only from the r-band, and lensed the source image in each band according to this model. We then added the lensed source in each band to the corresponding image of the lens taken from the DES images in the g, r, and i bands. Our final set of simulated galaxy-scale lenses consists of 100 000 systems with a uniformly distributed Einstein radius in the range 1.2′′ < θ_E < 3.0′′. Examples of these stamps are shown in Fig. 3, along with the non-lensed objects.

Fig. 3. Examples of cutouts used to train the neural network Top: examples of simulated lenses based on real DES images. Stamps are ordered by increasing Einstein radii (top-left to bottom-right). The top row corresponds to θE = 1.2–1.8″, the middle row to θE = 1.8–2.4″, and the bottom row to θE = 2.4–3.0″. Bottom: examples of LRGs used either as non-lenses during the training of the CNNs (see Sect. 4.1) or as objects onto which we inject a lensed source to build simulated lens systems, as shown in the top panel. All the cutouts are 50 pixels on-a-side, corresponding to 13′′..

4. Lens finding using CNN

Artificial neural networks (ANNs; Rosenblatt 1957) consist of an interconnected group of nodes which are typically organized into the so-called input, hidden, and output layers. In particular, CNNs (LeCun et al. 1989) – which are especially good at solving image classification problems (He et al. 2015) – have hidden layers that are of key importance, as they highlight the patterns in the data using a series of convolutional, pooling, normalization, and fully connected layers. The level of abstraction in the pattern features increases with the depth of the convolutional layers, assisting in the classification of objects into different classes. Here, we train a CNN to recognize strong lens systems against isolated red galaxies.

4.1. CNN training

The training set, consisting of 50 × 50 pixel cutouts in each of the gri-bands, is composed of two equal subsets: the first being 100 000 simulated lens cutouts from Sect. 3, and the other containing LRGs that were not used in the simulation process. We labeled our data with: 1 for lenses and 0 for non-lenses. We kept 20% of each sample as a validation set. Before training the CNN, we pre-processed our data by normalizing each image brightness to range between 0 and 1. We also augmented our sample by flipping each image horizontally and vertically. Data augmentation increases the probability that the network correctly classifies different orientations of the same image, but it does not transform the CNN into a rotationally invariant one. To achieve this, a different architecture must be used that is not explored in this work. The training process is performed using the Keras Deep Learning framework².

Our CNN uses a model from the EfficientNet family (Tan & Le 2020), which has been designed to achieve better performance than other CNNs. The network of this model uses a compound coefficient to scale the depth, width, and resolution, which are key parameters for obtaining better accuracy and efficiency. In particular, the EfficientNet implementation in Keras counts with 8 different variants B0-B7, whose depth, width, and resolution parameters have been carefully selected and tested to produce good results. The complexity and requirements of the models increase as we move from B0 to B7. As running a more complex model also implies the use of more computational resources, we decided to use an EfficientNet-B0, whose architecture is described in Tan & Le (2020), and is sufficient for our classification task and the characteristics of our data. After this CNN model, we added a sequence of fully connected hidden layers. The network has a total of 4 182 205 trainable parameters.

During training, the neural network learned how to grade images of galaxies and distinguish between lenses and non-lenses. At each iteration the network analyzed subsets of 32 images. When all iterations are completed through the entire training set, it is counted as one epoch. Within each epoch, the accuracy and loss of the model is monitored using the validation set. The next step is to minimize a binary cross-entropy loss function using a stochastic gradient descent optimizer (Adam) with a learning rate of 0.0001, and stop the training if either the loss value does not improve by more than 0.0001 over 10 epochs – or when 100 epochs are reached.

4.2. Evaluation of the CNN performance

The network provides a score, S_CNN, between 0 and 1, for each processed image. This means that those images classified as lenses obtain S_CNN ∼ 1 while non-lenses obtain S_CNN ∼ 0. The training process was performed in a single GPU Nvidia GTX 1080 Ti in about 8 hours. It converged, within our criteria above, after 57 epochs and achieved a 99.9 (99.8) percent accuracy in the training (validation) sets and a loss of 0.01 (0.02). This near-perfect accuracy achieved in the training set might be understood as overfitting, thus, to evaluate this possibility we compared the loss and accuracy learning curves for the training and validation sets (Fig. 4). For both loss and accuracy we see that after ten epochs the training set reaches a stable point with minimal changes, while the validation set follows the same trend with a small gap showing less accuracy and more loss than the training set, as is expected. The lack of overfitting signs, that is, the training loss continues decreasing or the validation loss starts increasing again after several epochs, leading us to the conclusion that our model is able to learn and generalize this classification problem.

Fig. 4. Learning curves for the accuracy (top) and loss (bottom) for the training (blue) and validation (red) sets, as a function of epoch. In the case of loss, we crop the y-axis, but the initial training set loss is 0.76.

In order to evaluate the performance of the CNN, we built two test sets. The first one contains 40 000 cutouts with the same characteristics as the training set, namely, half made up of simulated lenses and half of LRGs. The latter has 636 cutouts where half are known lenses or lens candidates (visually selected to have noticeable lensing features) and half are LRGs not seen by the CNN during training. The known lenses are taken from the Master Lens Database³ and the candidates from Jacobs et al. (2019a,b). The purpose of this second test set is to have a more realistic idea of the performance of the CNN in grading real strong lens systems instead of simulations. The distribution of S_CNN for both test sets (Fig. 5) shows that objects labeled as lenses are concentrated around S_CNN > 0.9 and non-lenses around S_CNN < 0.1, as expected.

Fig. 5. CNN score (SCNN) distribution for the different datasets. Top panel: SCNN distribution of both test sets: lens simulations (dashed line) and real confirmed and candidate lenses (solid line). Both datasets contain images labeled as lenses (purple lines) and real LRGs labeled as non-lenses (red lines). The two test-sets are normalized to their corresponding maximum value in the distribution. We shifted the x-axis of the second test set distribution by 0.02 for clarity. Bottom panel: SCNN distribution for the objects in the parent sample with scores above 0.5. We crop the x-axis for visualization as 99% of the sample is below 0.5.

To evaluate the number of lenses correctly identified, we used a receiver operating characteristic (ROC) curve, (Fig. 6) which shows the true positive rate (TPR) against the false positive rate (FPR), naturally both functions of the decision threshold applied to the score. It illustrates the performance of a binary classifier in discriminating between the two classes as the decision threshold is varied. The first test set shows a very good performance reaching an accuracy of 99.7% and a loss of 0.02. From the ROC curve we see that choosing S_CNN = 0.5(0.9) gives a TPR = 99.8% (99.4%) and a FPR = 0.21% (0.12%). On the other hand, in the second test set the performance of the network decreases obtaining an accuracy of 89.6%, and a loss of 0.44, with TPR = 76.1% (65.7%) and FPR= 0.31% (>0.01%) for S_CNN = 0.5(0.9). Thus, while the accuracy in the second test set is still high and the network did not grade any LRG above 0.9, the loss and TPR are significantly worse than for the dataset with similar characteristics to the training set. We think that this decrease in the performance of the CNN is because it was trained to recognize lens simulations and (despite having been created in a fully data-driven way) lacks the diversity and uniqueness of some strong lens systems (e.g., multiple deflectors, distortions produced by substructures or external sources, etc.). For example, most of the false negatives in this second test set are compact lens systems or have lensing features that are too faint to be properly recognized. Nevertheless, we found that our model is able to generalize and accomplish the goal of successfully classifying a high percentage of strong lens systems, although we are aware that in a realistic scenario we misclassify more objects compared to the simulations (as Figs. 5 and 6 show).

Fig. 6. ROC curve for the test set containing simulations (test set 1 in red) along with the confirmed and candidate lenses (test set 2 in blue); both datasets contain real LRGs as non-lensed examples. The FPR is plotted on a logarithmic scale to aid visualization. The TPR and FPR for SCNN = 0.5 (green) and SCNN = 0.9 (black) are also shown for each set.

When we applied the CNN to our parent sample, we found that 98.6% of cutouts obtained S_CNN ≤ 0.1, 133 322 obtained S_CNN ≥ 0.5 (Fig. 5, bottom panel), and 76 582 cutouts obtained S_CNN ≥ 0.9. The choice of S_CNN ≥ 0.9 is driven both by the ROC shown in Fig. 6, and also the resulting number of candidates being reasonable for human inspection.

5. Visual inspection

The 76 582 cutouts scored above S_CNN = 0.9 by the CNN were visually inspected by seven of the authors of this work (K.R., E.S., B.C., F.C., C.L., J.C., and G.V.).

5.1. Visualization tools and guidelines

We created two visualization tools⁴: one to quickly select lens candidates from many objects displayed simultaneously in a mosaic configuration and one to visually inspect each individual object in more detail and classify them into specific categories.

The mosaic tool simultaneously displays 100 color cutouts, each of which the user can mark for selection. The user can choose a random seed for displaying the images in a random fashion on the grid, to avoid all users seeing each object at the same location in the grid. This has the objective of preventing any possible bias from the position of the object on the mosaic coupled with the different level of concentration when looking at many mosaics in a row. This turned out to be very efficient, as illustrated by the “heat-map” of user grid selections displayed in Fig. 7, which are fairly flat, with a small bias towards selecting more objects from the top, bottom, and left row for this particular example. With the mosaic tool, we classified our sample into only two categories, namely, objects that we determined to be displaying potential strong-lensing features and the rest that we discarded from any following step of the visual inspection.

Fig. 7. Heat-map for the normalized mean number of times that each cell in the mosaic was clicked on, among all seven visual inspectors during phase 1 of the visual inspection. The values in each cell were obtained by calculating the mean of the total number of clicks per-cell among the seven users; we normalized these values by a factor of 24.78 that represents the mean number of clicks in a cell for this specific classification.

The second visualization tool allows us to inspect one by one all lens candidates selected with the mosaic tool. In doing so, we display the gri color stamps, which allow the user to change the display scale and color map. With this tool, we can classify each object into one of four categories: 1- sure lens, 2- maybe lens, 3- single arc, and 4- non-lens. In addition, we can define five sub-categories for objects classified as non-lens: 1- ring, 2- spiral, 3- elliptical, 4- disk, and 5- merger.

In order to achieve a more consistent classification among users, we all agreed to follow the same guidelines for the four main categories. “Sure lens” is selected when the cutout shows a clear strong lensing configuration without the help of a higher resolution image. This means that several clear multiple images can be identified or that there are signs of a counter-image. “Maybe lens” is chosen if the object shows a promising lensing-like configuration but a clear identification of multiple images is not possible visually. This category also includes cases where several objects or a single arc-like object lie on one side of the central galaxy but no clear counter-image can be distinguished on the other side. In this case, high-resolution imaging or spectroscopy will be required to decide whether it is a false positive or a genuine lens. When there is a single image object or a single arc far away from the central galaxy with signs of tangential distortion, the cutout goes to the category of “single arc”. Finally, everything else that does not fit these categories is classified as “non-lens”.

5.2. Visual selection procedure

We used both tools in four different phases to ensure that we have a clean sample of not only potential lens candidates but also a sub-sample of contaminants, such as ring galaxies, which are a source of confusion for CNNs and a matter of debate among visual inspectors. The different steps carried out to perform the visual inspection are described below.

First, for the lens and ring galaxy selection we used the mosaic tool. We selected from among the 76 582 cutouts all objects that presented signs of lensing features or looked like ring galaxies in one category and we discarded the rest. An average of 2478 cutouts were selected per visual inspector, the normalized mean distribution of clicks per cell are shown in Fig. 7. A total of 9210 objects was selected by at least one user, while 89 of them were selected by all the users unanimously (see Fig. 8, top panel).

Fig. 8. Number of objects classified by a certain number of visual inspectors in phase 1 (top panel: to select lenses and ring galaxies) and phase 2 (bottom panel: to select only the ring galaxies from phase 1). Objects selected by the seven visual inspectors represent a 100% in agreement among users, while the sum along the different bins give us the union, namely, the number of objects selected by at least one user. The exact number in each bin is shown at the top of each bar.

Second, for the ring galaxies selection we used the mosaic tool to select only ring galaxies from the 9210 objects. A mean of 230 cutouts were selected per visual inspector, but only 71 were classified by all seven, while a total of 1445 were selected by at least one (see Fig. 8, bottom panel).

Next for the lens systems classification we visually inspected all the 9210 objects selected in phase 1 using the one by one visualization tool, looking specifically for lens systems. Here, we showed again the classified ring galaxies from the previous step as a consistency check (users should re-classify them as rings, or at least not classify them as lenses). We classified each object into: “sure lens”, “maybe lens”, “single arc”, and “non-lens”. Optionally, if a non-lens was clearly identified by the user as a spiral, merger or ring galaxy, the object was sub-classified into the corresponding category. From this visual inspection, we obtained a total of 275 “sure lens”, 2666 “maybe lens”, 2602 “single arc”, and 9125 “non-lens”, classified at least by one visual inspector. In a unanimous agreement among the seven visual inspectors, we counted only 6 “sure lenses”, 1 “maybe lens”, 1 “single arc”, and 4716 “non-lens”. On the other hand, K.R., E.S., B.C., F.C., and J.C. sub-classified 359 ring galaxies, 22 mergers, and 49 spirals with an agreement of 50% among the visual inspectors. In Table 1, we summarize the individual classification by category and subcategory of each user and in Fig. 9 we show examples of objects classified in each category and sub-category, with the exception of “sure lens”.

Fig. 9. Example of objects classified as “maybe lens” (ML) or “single arc” (SA) and sub-classified as part of the “ring galaxy”, “spiral”, or “merger”, or “sure lens” category is shown in Fig. B.1. At the top of each image is the name of each system, while at the bottom, we show the CNN score (SCNN) and visual inspection score (VIS) obtained in the corresponding category.

Finally we implemented a group visual inspection. This last step was included due to the lack of agreement among users for the main classification categories (sure lens and maybe lens), as Fig. 10 shows. In this figure, we can see, for example, that user 2 classified 1396 objects in both categories; this user stands as the one that classified more cutouts as candidates, but the overlap with other users is not higher than 321 (along with user 4 and user 5). This step was performed by K.R., E.S., B.C., F.C., J.C, and G.V. altogether. Using the one-by-one visualization tool, we revised 2690 objects selected as “sure” or “maybe” lens (or both) by at least one visual inspector. The aim was to obtain a final selection of potential candidates that can be suitable for follow-up high-resolution imaging and spectroscopic confirmation and to avoid spending telescope time on false positives. We classified them into the two main categories as follows: 81 sure lenses and 296 maybe lenses. This represents 0.5% of the sample with S_CNN ≥ 0.9 and 0.002% of the initial LRG selection sample.

Fig. 10. Correlation among different classifiers for the categories “sure lens” and “maybe lens” after phase 3. The values in the diagonal represent the total number of objects classified into both categories for each user, while in the adjacent rows and columns the number of objects that both users classified into the same category.

An extra visual inspection step was performed for the cutouts classified by the CNN with scores between 0.8 < S_CNN < 0.9 (hereafter referred to as bin80) by K.R., E.S, and B.C. with the purpose of quantifying how many objects we could have missed by selecting only those with S_CNN > 0.9. Similarly to the previous analysis, we first used the mosaic visualization tool to inspect the 17 779 cutouts, selecting 190 potential lens candidates, which were then inspected one by one. A set of 24 objects were classified as lenses by at least one visual inspector and 115 as maybe lenses. Finally, K.R., E.S., B.C, F.C, J.C, and G.V. conducted a group visual inspection of the 190 initially selected candidates to compile a final sample with 9 sure lenses and 19 maybe lenses that were then added to our candidate list. In total, only 0.2% of the data visually inspected in the bin80 was considered as a lens candidate, while for all the objects with S_CNN ≥ 0.9 we selected 0.5% of the cutouts in the categories of “sure” or “maybe” lenses. Furthermore, taking into account that the amount of cutouts classified in the bin80 is about four times smaller than those with S_CNN ≥ 0.9, we concluded that the number of expected candidates with 0.8 < S_CNN < 0.9 was very low, showing that we reached a point of diminishing returns which would make the additional human visual inspection of images classified with S_CNN < 0.8 ineffective.

5.3. Final catalogs

As a final product, we present two main catalogs: one containing lens candidates and one containing ring galaxy candidates. We assigned a visual inspection score (VIS) to each candidate, computed using the percentage of visual inspectors that classified it into a certain category. In the case of lens candidates, we used the percentage of users that classified a system as either a “sure” or a “maybe” lens. We summed up these percentages to obtain a “strong lensing percentage”. Then, we considered this percentage as the final visual inspection score for lens systems (VIS_L). In the case of ring galaxies, we had candidates from step 2 (using the mosaic tool) or step 3 (or both), based on the one-by-one method). We averaged the percentage of users who classified each object as a ring in each step and we present this as the final visual inspection score for Rings (VIS_R). If the candidate was selected only by one of the tools, the final score obtained is the one corresponding to that classification (i.e. not the average).

The ultimate catalog of lens systems can be split in two categories, “sure lens” (SL) with 90 systems (Fig. B.1) that display prominent lensing features and counterpart images and “maybe lens” (ML) with 315 systems that show promising lensing features but for which more evidence, such as higher resolution imaging and spectra, is needed (see Fig. 9 for examples). From the figures, we can conclude that a large portion of our candidates are group- and cluster-scale lenses. This is mostly because to be able to identify them in ground-based data, most of them should have an Einstein radius above 1′′.

The CNN and visual inspection scores of both the SL and ML candidates is shown in Fig. 11. Here, we clearly see that most SL systems are clustered towards the upper right corner, indicating that in general they obtained a high score from both methods, while very few of them had either CNN scores below 0.95 or visual inspection scores below 0.5. On the other hand, a large majority of ML objects did not receive a high visual inspection score, including two that originally were rejected by visual inspectors, but upgraded after the group visual inspection. Several of the ML objects still got very high scores from the CNN, indicating that the visual inspection step is needed to refine the final catalog.

Fig. 11. CNN score (SCNN) against the visual inspection score (VISL), with their respective distributions for the final catalog of lens candidates containing 405 systems. The distribution of the category “sure lens” is in purple and “maybe lens” is in red, while in the scatter plot “sure lens” systems are represented by stars and “maybe lens” systems by triangles.

In order to identify lens candidates that were not previously published, we cross-matched our final catalog with available astronomical databases such as Vizier (Ochsenbein et al. 2000), Simbad (Wenger et al. 2000), the Master Lens database, and other lens-finding works, including Wong et al. (2018), Jacobs et al. (2019a,b), Petrillo et al. (2019b), Cañameras et al. (2020), Jaelani et al. (2020), Huang et al. (2020, 2021). As a result, we found that our catalog contains 219 previously identified candidates (74 SL, 145 ML), including at least 5 spectroscopically-confirmed systems and 186 new candidates (16 SL, 170 ML). The detailed information for these systems can be found in Table 2, available at the CDS.

Our second catalog is composed of ring galaxy candidates classified during two different steps of the visual inspection process. We identified 1445 ring galaxy candidates during the second visual inspection step using the mosaic tool, while 985 galaxies were classified into this category by at least one visual inspector using the one by one tool in the third step. A cross-match between these two selections gave an intersection of 854 galaxies for a total of 1 576 ring galaxies selected by at least one user using either of the two methods. The final catalog was built by selecting the objects picked by at least 50% of the visual inspectors using the mosaic or the one by one tool, resulting in 539 ring galaxy candidates. In Fig. 9 we present the six top-graded candidates and in Table 3, available at the CDS, we detail the information for the full sample.

Finally a catalog using the classification “single arc” (SA) will be created after a more detailed analysis of these objects, but this is beyond the scope of the current work. We expect that this catalog could serve as a probe for such works as Birrer (2021).

5.4. Lens-source decomposition

For ground-based observations of most strong lensing systems, the light from the source and lens galaxies are blended. In order to better visualize our 90 lens candidates in the SL catalog, we designed a prototypical automated procedure for deblending the light from lens and source galaxies based primarily on their colors.

Due to the complexity of the light profile of the lensed sources we chose to represent them in a non-parametric way using undecimated isotropic wavelets (starlets, Starck et al. 2007), as implemented in the Multi-band morpho-Spectral Component Analysis Deblending Tool MuSCADeT⁵ (Joseph et al. 2016). Starlets are a family of functions that allow free-form modelling of images at various spatial scales and present advantages for modelling smooth galaxy profiles as discussed in the MuSCADeT paper. Estimating both color and morphology of sources requires a large number of parameters, larger than the number of pixels in the starlet-decomposed image, making it a degenerate problem. To overcome this we use a combination of the scarlet⁶ (Melchior et al. 2018) and MuSCADeT algorithms. In both methods, multi-band images are modeled as sums of factorised components, where each object, i, in an image has a 2D surface brightness, S_i, with as many pixels as there are in the image bands, and a spectrum, A_i, with as many entries as there are bands (see details in Melchior et al. 2018; Joseph et al. 2016), such that:

$$\begin{array}{c}\hfill Y={\displaystyle \sum _{i<o}}{A}_i{S}_i+N,\end{array}$$(1)

where Y is a multi-band cube of images, o is the number of objects in the scene, and N is the noise map.

The strategy implemented in MuSCADeT only allows for crude estimates of source colors, based on principal component analysis of pixel fluxes. Instead, the scarlet software is able to estimate the colors of each source in the field provided that the morphology is constrained to be a monotonic profile. Monotonicity of galaxy profiles from the center out does not suit the description of complex lensed sources, hence the need for MuSCADeT to model strongly lensed galaxies in a non-parametric way coupled with sparse regularization. In short, MuSCADeT is used to model the 2D profile of galaxy images, including a complex lensed source, while scarlet recovers the colour of the objects.

Scarlet requires detection of the brightest pixel of each source to model, which is a challenging and ill-defined problem in the case of strongly lensed galaxies, where lensed features are often multi-modal and strongly blended with the deflector’s light. In order to circumvent this issue and make sure we capture (lensed) sources with a bluer spectrum than the central LRGs, we allow scarlet to model one source with Starlets, initialized with a “blue” spectrum. This allows scarlet to capture blue features with complex morphologies that might not have been detected due to blending, while limiting degeneracies with other sources. The blue normalised spectra used for initialisation are empirically set to [0.4, 0.4, 0.2], where each of the three values reflect the relative contributions to g, r, and i bands respectively. Other non-lensed sources are detected using the sep⁷ package (Bertin & Arnouts 1996). We thus run sep on a filtered version of the images. The filtering is done by computing the starlet decomposition of an image and setting to zero the coefficients that contribute to low frequencies before reconstructing the image. This amounts to a high pass filtering that favours peak detection. The position of the brightest pixels of objects detected by sep are fed as entries to scarlet. For each object detected with sep, then scarlet estimates a spectrum (flux in each band). From these spectra obtained with scarlet, we select the bluest and reddest spectra by finding those that maximize the scalar product between the normalized spectra [0.667, 0.333, 0] (for blue) and [0, 0.333, 0.667] (for red). This ensures that two components with different colors are extracted, with the expectation that the red component features the morphology of the LRG and its neighbours, while the blue component extracts the morphology of the lensed star-forming background galaxies.

The summary of the procedure for deblending strong gravitational lens candidates is as follows:

The detection of sources in the image is done using the source extraction package, sep (Bertin & Arnouts 1996) on a starlet-filtered version of the image where only the first two levels of the starlet decomposition are used. We set the detection threshold to 1 noise standard deviation of the noise upon running sep. This may seem aggressive and as if it would potentially lead to a shredding of the objects upon deblending, but the smoothness of the images imposed by the PSF prevents such an effect. Furthermore, shredding is not an issue as we only intend to capture the spectra of the reddest and bluest objects.

The initialization of scarlet sources takes one extended source per detected object plus one starlet component with blue spectra. Scarlet uses a target PSF of Gaussian profile with a standard deviation of 0.5 pixel. This is the target resolution to which all bands are uniformly deconvolved to.

The next step is to run Scarlet for up to 200 iterations and extract spectra for each source in the field of view by simply measuring the flux in each source. Then we identify the bluest and reddest sources through scalar product with predefined red and blue spectra. We run MuSCADeT with the red and blue spectra for 200 iterations. The threshold for starlet reconstruction is set to 5 sigma of the standard deviation of the noise. This means that the model reconstructs features that are 5 sigma above noise levels.

Finally, we extract red and blue components by computing the difference between the multi-band images and the model for each MuSCADeT component. The details of the processing can be found in the notebook that was used to generate these images: https://github.com/herjy/Lensdeblend/releases/tag/prototype.

The results for our best lens candidates are shown in Fig. B.1 and display for each system the red residuals, namely, the data from which the blue model has been subtracted, and the blue residuals, namely, the data from which the red model has been subtracted. In the following, we refer to these as R_r and R_b, respectively, defined as:

$$\begin{array}{cc}& R_{b,j\in \{g,r,i\}}=Y_j-A_{r,j}{S}_{r,j},\hfill \end{array}$$(2)

$$\begin{array}{cc}& R_{r,j\in \{g,r,i\}}=Y-A_{b,j}{S}_{b,j},\hfill \end{array}$$(3)

where A_r, jS_r, j and A_b, jS_b, j are the models for the red and blue components in each band j.

The results in Fig. B.1 show that the lens and source light can be deblended efficiently without fitting any analytical profile. The effectiveness of the method to deblend the profiles comes mostly from the spectral decomposition of the objects and on their representation on an array of pixels to which we apply sparse regularization with wavelets (starlets). This procedure is well suited to automated use in a pipeline but assumes that lensed sources are significantly bluer than the lens light. This is the case for most of our lenses as by construction our lens finding method is based on a preselection of objects that favors such a configuration. Still, we do have objects where the lens-source color contrast deviates significantly from our assumption. In this case the deblending works less efficiently and we see leakage of flux between the lens and source. Another case of leakage, leading to sub-optimal deblending can be observed in systems where the image contains sources with colors different from those of the lenses or sources. In this case, since the whole image is modeled as two fields of light, the spectra of the color components tends to offset towards an average spectrum that better matches all the colors in the patch. This can be observed in systems DES J013522-423223, DES J024911+004848, and DES J010826-262019, where the blue components contain light from the lens galaxies and contain objects with colors different from that of the main deflector. These shortcomings are motivation enough for further refinement of our deblending, in particular with focus on using scarlet to better model individual, non-lensed sources, which is beyond the scope of this paper. Finally, it is important to emphasize that our light deblending confirms our visual grading and does not discard any of our best candidates.

6. Model

We developed an automated modeling pipeline in order to further explore the highest rated lens candidates obtained from the visual inspection. Our candidate sample is very heterogeneous, containing galaxy, group, and cluster scale systems. Thus, in order to perform this automatic modeling we split the sample and selected through visual inspection only the images in which there appeared to be a single lens galaxy as a deflector. The 52 images selected for modeling are labeled with an “M” in the mosaics of Fig. B.1. This pipeline allows us to efficiently model large samples of lens candidates acquired in current and future lens finding efforts and to explore the model parameter distributions in search of meaningful trends.

6.1. Automated modeling pipeline

We modeled the images using single elliptical Sérsic profiles for the light distributions of both the deflector and the source. For the mass distribution of the deflector, we used a singular isothermal ellipsoid profile (SIE), along with an additional external shear component (γ_ext). The simplicity of these profiles allows us to model many lens candidates efficiently, while still fitting most images well enough for us to observe large-scale trends in the properties of the sample. The pipeline supports multi-band fitting, so we fit the DES lens candidates using images in the g, r, and i bands. We used a separate elliptical Sérsic profile for each of the three photometric bands when fitting the deflector and source light components, although we fixed the center positions between bands and we added priors to bound the semi-minor and semi-major axes. The deflector mass profile is shared across all bands.

The modeling pipeline is entirely written in Python and it makes use of the Lenstronomy lens modeling package (Birrer et al. 2015; Birrer & Amara 2018). For the parameter optimization, we used a particle swarm optimization (PSO; Kennedy & Eberhart 1995), and to estimate the variances in the sample we used a Markov chain Monte Carlo (MCMC) sampler. For each image, the pipeline first performed a chain of pre-sampling PSOs before running the sampling with the MCMC. The MCMC is performed using an affine-invariant MCMC ensemble sampler (Goodman & Weare 2010; Foreman-Mackey et al. 2013), which was implemented using the emcee⁸ Python package.

In order to obtain realistic results on the parameters, we introduced priors that punitively discourage extremes in some of the model parameter values. While we cannot make any definite assumptions about the position angles of the Sérsic or SIE profiles, we can use Gaussian priors on the ratio between the semi-minor and semi-major axes, q. The Gaussian prior was centered on a value of q̄ = 0.8, in accordance with the distributions of 138 269 galaxies from the Galaxy And Mass Assembly (GAMA) database that were modeled in Kelvin et al. (2012). This prior is applied first to the r-band, and then to allow only small variations between bands in the light components of the model, the Gaussian prior of the other two bands is centered on the result obtained in the r-band. We also used a similar Gaussian prior method to constrain the deflector mass eccentricity and position angle to values close to those of the deflector light. Lastly, we also applied a prior distribution for the effective (half-light) radius, R_eff, and Sérsic index n_s parameters of the source light. The source priors we used came from the Sérsic parameter distributions of 56 062 galaxies from the COSMOS survey. This data was used as a training set in the development of the GalSim⁹ software. (Rowe et al. 2015). We show these distributions in Fig. 12.

Fig. 12. Prior distributions for the effective radius (Reff) and the Sérsic index (ns) used for constraining source light parameter values.

When modeling the lens candidates, it is common for image cutouts to contain neighboring objects in the field of view that are unrelated to the lens system. Light contamination from these “satellites” can be mistaken as having originated from the lensed source if not masked properly. In the literature, this problem is handled differently by various authors. Shajib et al. (2020) excluded systems contaminating satellites in their sample, and modeled only isolated lenses from the SLACS survey (Auger et al. 2009). Nightingale et al. (2018) did not pre-select isolated lenses but, instead, these authors masked all the pixels outside of a circular region with a fixed radius of 3.9′′.

In our case, we designed the pipeline to be flexible in handling a large variety of lens system configurations and sizes. The steps of our masking procedure are illustrated in Fig. 13, and begins with applying filters in order to identify the brightest regions in the image as well as their centroid locations. We first applied a Laplacian of Gaussian (LoG) filter to detect areas with rapid changes in flux. Next, we took all remaining pixels with flux less than a threshold of six times the rms background and set them to zero. This results in a final filtered image with only the areas of the image containing the most light. We find the centroid locations of these areas by finding the local maxima, or peaks, in the final filtered image. For our masking algorithm we made use of both these peak locations as well as the pixel values – these are labeled with black and red markings in the bottom left panel of Fig. 13, respectively. The peak locations are used first for determining the center of the lens system, that is, the position of the deflector galaxy, and assume that this is the peak detected object nearest to the image center. Because the deflector is assumed to be an LRG, we use the reddest available band (i-band) for this step. Next, we take the detected peaks in the bluest band (i.e., the g-band) to estimate the lens system size. This is because the source galaxies in lens systems are usually younger and more active galaxies, meaning that the lensed source light will be more prominent in the bluest image band. We assume that the second closest detected peak to the center is the first of the lensed images of the source. We also assume that the furthest lensed image from the deflector is not more than 1.5′′ further out than the nearest one. Therefore, our estimated lens system radius is the distance from the deflector to the closest lensed source object plus 1.5″. We show this estimation as a black circle enclosing the lens system in the bottom middle panel of Fig. 13.

Fig. 13. Illustration of the automated masking procedure using an example DES image. In the upper row, we show (from left to right) the original image, the image after applying an LoG filter, and the result of setting pixels with flux below the threshold to zero. In the bottom row we show in the leftmost frame the original image annotated with the remaining pixels from the filtering step (red “+” marks), along with the detected peaks (black “x” marks). In the middle we show the estimated size of the lens system with a black circle, as well as the detected bright pixels that are considered contaminant light. These pixels are then used for the mask, and we show the areas covered in the rightmost panel by setting the corresponding pixels to a large constant value.

Using the estimated size of the lens system from the g-band image and the location of the deflector obtained from the i-band, we create a circular mask for each band that is centered on the deflector location and only covers detected bright pixels outside of the circular region with our determined size. The mask itself is a boolean array with the same shape as the original data, and has the value of zero at any pixel that is to be ignored in Lenstronomy computations and ones everywhere else. In the bottom right panel we illustrate the coverage of the mask by setting all of the “ignored” pixels of the original image to a large constant value.

On average, our pipeline took 4.3 hours to model a gri DES system. This includes reading data, the masking process, and performing the modeling sequence to find the best parameters that describe the lens candidate.

6.2. Modeling Results

Using our automated pipeline, we modeled 52 of the systems in the SL catalog that appeared to have only a single galaxy as a deflector. We show in Fig. B.2 a sequence of images for visualizing the modeling results in the r-band, including the corresponding image, a reconstructed image, normalized residuals, convergence map, and the reconstructed source light. In Table B.2, we present the best model parameters obtained for each system, and we show the obtained distributions for the Einstein Radii, the external shear, and the effective radius and Sérsic index for the lens and source light, in the histograms in Fig. 14.

Fig. 14. Model best-fit parameter distributions for the lens mass, lens light, and source light profiles. Results for which reduced χ2 ≤ 1.5 are shown in green and those with reduced χ2 > 1.5 are displayed in red. Top left: Einstein radii in arcseconds. Top right: External shear strength, γext. Middle left: Sérsic half-light radii of the lens light. Middle right: Sérsic indices of the lens light. Bottom left: Sérsic half-light radii of the source light. Bottom right: Sérsic indices of the source light.

We obtained acceptable fits for 41 systems, which represent 79% of the sample, and we observed 11 failures in the fitting, which we define as fits with mean reduced χ² per pixel above χ² = 1.5. In the lens mass components, we observed Einstein Radii, R_E, distributed between ∼1′′ and ∼3.5′′. For the external shear strengths, we observed γ_ext ≲ 0.47 for all lenses, except for one in which the fit failed. The distribution of the values have a peak at 0.14 , these are typical shear values for strong lens systems (Keeton et al. 1997). For the effective radii, R_eff and Sérsic indices, n_s, of the lens light profiles in the r-band, we observed peaks at R_eff ∼ 2′′ and n_s ∼ 5, respectively. Because the CNN searches selected lens systems from a catalogue of LRGs, we expect to obtain deflector light parameters that are typical for LRGs, and that is indeed what we recover. For the parameter distributions of the source light, we observe the effective radii and Sérsic indices peaking at R_eff ∼ 0.2′′ and n_s ∼ 1, respectively. This is also the expected behavior for smaller, low-mass galaxies that are usually dominating the lensed galaxy source population.

When modeling these lenses, the primary source of failures lies in the masking procedure. For example, the estimated size for the lens system is either slightly too small or too large, resulting in parts of the lensed source light being masked or neighboring contaminants not being masked and instead treated as lens features. This has happened for four systems with failed models (DES J060653−585843, DES J015216−583842, DES J032216−523440, and DES J051047−263222) and for two considered to be characterized by acceptable fits (DES J034713−453506 and DES J040822−532714). Since it is common for images to contain companion objects very close to the lens systems, there is a small margin for error in determining the lens system size. For the system, DES J012042−514353, the contaminant is actually residing among the lensed images of the source, a situation which cannot simply be handled with a more precise measurement of the lens size. A method would be needed for better untangling the contaminant light from the lens features. Finally, there are two lens systems (DES J010553−053419, DES J041809−545735) in which the contaminant light distributions were spread out enough that the mask failed to adequately cover all this light. In general, we need to improve our masking procedure to avoid these problems during an automatic fitting of the lens. In the meantime, for all of the systems for which the masking algorithm did not perform well, we recreated masks by hand and performed the modeling a second time. These results are shown in the rows directly below the original results for the specific system and both sets of results are enclosed in a red dashed box in Fig. B.2. Each time, we see a significant improvement in the residuals after using the better mask. On the other hand, the rest of the models that are considered as failures (DES J010659−443201, DES J021159−595624, DES J024803−061606, and DES J202855−523118) do not show an obvious reason for it, but are likely due to the compactness of the system, faint lens features, or the complexity in the shape of the source. For these cases we need further investigation to find a general solution to improve their models in the automatic pipeline.

7. Conclusions

We used DES-DR1 to search for galaxy-scale strong lensing systems using a CNN that carries out a binary classification of optical images in the g, r, and i bands. In doing so, we targeted massive galaxies, that is, LRGs, which were selected using a wide color-magnitude cut accounting for realistic color contamination by the putative background star-forming blue galaxies.

The design of our training set was data-driven in the sense that real DES images of LRGs were used to mock the light distribution of the lens plane. Real images of galaxies from the COSMOS HST were used to mock lensed sources. This helps ensure diversity in colors and morphologies for the sources and lenses, but also preserves the sky background characteristics, galaxies, or stars acting as companions, as well as any artifacts in the images.

We used these data-driven simulations as positive examples to train a CNN while using a portion of the LRG sample as negative examples, despite some previous searches that included as negative examples other types of galaxies. After analyzing the results from the visual inspection, we determined that the lack of representation of other types of galaxies was not important, as these were not the most relevant source of false positives. In fact, LRGs with bluish satellites near the line-of-sight are our most important contaminants. The CNN was trained and validated using a total of 200 000 images, half of them being mocked lens systems labeled as 1 and the other half being LRGs labeled as 0. Evaluating our model on a test set built from images with the same characteristics of the training set gave us an accuracy of 99.7%. On the other hand, a small test set built with 300 lens candidates and the same proportion of LRGs gives us a more realistic evaluation reaching an accuracy of 89.6%.

Applied to the 18 745 029 LRGs drawn from our color-magnitude selection, we obtained 76 582 images with CNN scores above or equal to 0.9, that several authors visually inspected. To do so, we created guidelines to separate them into different categories: “sure lens” (SL), “maybe lens” (ML), “single arc” (SA), “non-lens” (NL), and subcategories: “ring galaxy,” “spiral galaxy,” and “merger” for objects falling in the “non-lens” category. To perform the classification, we used a mosaic visualization tool displaying 100 images at once, as well as a one-by-one visualization tool that displayed the color composite image and each band for one object at a time. We classified 0.5% of the 76 582 images as lens candidates, with 81 falling in the SL category and 296 in the ML category. Additionally, we inspected the 17 779 cutouts with a CNN score in the range 0.8 < S_CNN < 0.9, with only 0.2% of the images classified as lens candidates, that is, 9 SL objects and 19 ML. The visual inspection of these low-score lenses allowed us to conclude that the reward for inspecting images with scores below 0.9 was very poor compared with the amount of work. We therefore did not consider systems with even lower scores at all.

From our visual inspection, we created two main catalogs: a lens candidates catalog and a ring galaxy candidates catalog, the latter being our main source of contaminants. The first catalog contains a total of 405 lens candidate systems: 90 SL and 315 ML. Out of these, 186 were totally new systems and 219 were identified (but not necessarily confirmed) in previous searches. We deblended the lens and source light for our 90 SL systems using the MuSCADeT software, which does not involve any profile fitting, but uses the color contrast between the lens and source together with sparse regularization. This was successful in deblending most of the cases, where there were clear differences in the colors of the lens and source. The second catalog contains 539 ring galaxy candidates. We expect to use this ring catalog in the future to improve the training of machine learning algorithms in the recognition between lenses and ring galaxies. Still, 539 objects is not a sufficient sample for training CNNs and, thus, furthers work with, for instance, generative adversarial networks, is likely to be needed.

Finally, we selected, from the SL category, the 52 systems that apparently had one well-defined galaxy as a deflector to test an automated modeling pipeline. The relatively simple SIE + γ_ext and elliptical Sérsic profiles used in the modeling appear to be sufficient in describing these lens systems, and additional complexity is not necessary for the purposes of this automated modeling pipeline, at least with the image quality of DES-DR1. We successfully modeled 41 of these systems, while the other 11 failed mainly due to problems in the masking algorithm, especially in the estimation of the lens system size. To address these failures, we plan improvements in a future version of the pipeline including the use of the decomposed images from MuSCADeT to initialize the code and find the correct position and size of the system.

The outcomes of this lens finding work in DES-DR1 include a catalog with 405 meticulously selected lens candidates that can serve as a start for spectroscopic confirmation. In our selection, we did our best to privilege quality of the candidates over their quantity. The methods and tools studied, developed, and presented here do have room for improvement, but they serve as a preview of what can soon be achieved for the future generation of surveys, such as LSST, Roman telescope, and the Euclid mission.

Acknowledgments

This work is supported by the Swiss National Science Foundation (SNSF) and by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (COSMICLENS: grant agreement No 787886). R.C. thanks the Max Planck Society for support through the Max Planck Research Group of S. H. Suyu. This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (LENSNOVA: grant agreement No 771776). G.V. has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodovska-Curie grant agreement No. 897124. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI: https://doi.org/10.26093/cds/vizier). The original description of the VizieR service was published in Ochsenbein et al. (2000). This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France (Wenger et al. 2000). This material is based upon work supported by the National Science Foundation under Cooperative Agreement 1258333 managed by the Association of Universities for Research in Astronomy (AURA), and the Department of Energy under Contract No. DE-AC02-76SF00515 with the SLAC National Accelerator Laboratory. Additional funding for Rubin Observatory comes from private donations, grants to universities, and in-kind support from LSSTC Institutional Members. This work rely on the following packages: numpy (van der Walt et al. 2011), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), astropy (Astropy Collaboration 2018) and https://github.com/herjy/scarletextensions/releases/tag/paper and csv. A notebook detailing the procedure of deblending can be found in the following repository: https://github.com/herjy/Lensdeblend/releases/tag/prototype.

References

Appendix A: Comparison with Jacobs et al. (2019b,a)

We attempted to compare our results with (Jacobs et al. 2019b,a, hereafter J19AB). Using DES, these two searches and obtained a total of 1256 candidates falling in three categories defined by the authors: “definitely,” “probably,” and “possibly” a lens. Making a comparison with our work is more than challenging as our procedures are very different in color and magnitude selections, cutout size, simulated training sets, CNN architectures, and visual inspection. Being aware of these methodological differences, we compared the results for the 693 candidates in J19AB that are in our parent sample. Our CNN gave a score above 0.8 to 262 of these objects. Among these, 39 are classified as "definitely," 98 are "probably," and 125 are "possibly" a lens, according to the J19AB classification. After our visual inspection we found that 50 are on our list of SL (split up into 28 "definitely," 21 "probably," and 1 "possibly" in the J19AB classification) and 83 in our ML classification (8 "definitely," 51 "probably," 24 "possibly"). This means that 129 candidates in the J19AB list that were also selected by our CNN did not pass our final visual inspection criteria. These objects are shown in Fig. A.1. Aside from a few objects in the first three rows on the top of the figure, we still find that no other object shows sufficient evidence of strong lensing. And indeed, J19AB graded them all as "possibly." The main difference between our work and J19AB for these 129 objects is that we discard them while J19AB still include them in their list, hence leading to a list of candidates that is more extensive than ours.

In addition, 505 candidates from J19AB obtained a score below 0.8 from our CNN and were not visually inspected. These are split into 6 in the "definitely" category of J19AB, 113 as "probably," and 386 as "possibly." We display 56 examples of these objects in Fig. A.2, noting that if we were to visually inspect them now, very few objects would pass our criteria, for the following reasons: 1) the bluish features are too close to the central galaxy that could be mistaken as a star forming galaxy or 2) there is no evidence for lensing at all because the lensing features are outside our stamp size or because of too low S/N. Determining more specific reasons why our CNN discarded these candidates would require a more detailed study that remains well beyond the scope of our work.

Fig. A.1. Images for the 129 candidates in J19AB that we discarded after visual inspection. Our CNN score is indicated below each object as well as the category allocated by J19AB, in which A stands for “definitely,” B for “probably,” and C for “possibly” a lens.

Fig. A.2. Example of cutouts of 56 candidates in J19AB that our CNN graded below 0.8. Our CNN score is indicated below each object as well as the category allocated by J19AB, in which A stands for “definitely,” B for “Probably,” and C for “Possibly” a lens.

Appendix B: Additional material

Fig. B.1. Images for the 90 lens candidates in the category Sure lens and their corresponding decomposition performed with MuSCADeT. In the first and fourth columns we have the gri-composite image of the system, the name is on the top, while the CNN score and the visual inspection score (VISL) are displayed at the bottom of each image. Additionally, we marked with an "M" those that we modeled in Sect. 6. Columns 2 and 5 show the subtraction of the blue model from the respective data. Columns 3 and 6 show the subtraction of the red model from the respective data, isolating the lensing features.

Fig. B.1. continued.

Fig. B.1. continued.

Fig. B.1. continued.

Fig. B.1. continued.

Fig. B.1. continued.

Fig. B.2. Modeling results for the 52 lens candidates that appear to have only a single lens galaxy acting as a deflector. 1st column: Observed DES image of the lens system in the r-band. 2nd column: Reconstructed image using best-fit model parameters. The black regions are "masked" pixels that are ignored in the modeling as they contain light from contaminant objects in the image. The red curves are the critical lines of the lens model. 3rd column: Normalized residual map showing the difference between the best-fit model and the original data. 4th column: Convergence map of the lens model. 5th column: Reconstructed source light profile (un-lensed). The caustic curves are shown in yellow.

Fig. B.2. continued.

Fig. B.2. continued.

Fig. B.2. continued.

Fig. B.2. continued.

Fig. B.2. continued.

Fig. B.2. continued.

Fig. B.2. continued.

Fig. B.2. continued.

Fig. B.2. continued.

Fig. B.2. continued.

Fig. B.2. continued.

Fig. B.2. continued.

All Tables

All Figures

Fig. 1. Color-color diagram of our parent sample of ∼19 million galaxies. The plot displays where most galaxies in our selection are located. The green shaded area shows the density of galaxies from the parent LRG sample. The green solid lines shows the 1, 2, 3, and 6σ contours. The lens candidates and contaminants from this work are shown in overlay: Sure Lens (purple stars), Maybe Lens (open red triangles), and ring galaxies (open cyan circles). See Sect. 5.3 for details.

In the text

Fig. 2. Redshifts and velocity dispersion distributions used for simulations. Top panel: redshift distributions of the lenses (red) and sources (blue) in the simulated training set. Bottom panel: lens velocity-dispersion distribution, predicted from the K-nearest neighbor algorithm (red) and shifted to match the criteria of pairing lenses and sources (orange), as described in Sect. 3.2. We note that the actual lens velocity dispersions used in the simulations are shifted to higher values, so that lensing features can be seen even at the DES resolution.

In the text

Fig. 3. Examples of cutouts used to train the neural network Top: examples of simulated lenses based on real DES images. Stamps are ordered by increasing Einstein radii (top-left to bottom-right). The top row corresponds to θE = 1.2–1.8″, the middle row to θE = 1.8–2.4″, and the bottom row to θE = 2.4–3.0″. Bottom: examples of LRGs used either as non-lenses during the training of the CNNs (see Sect. 4.1) or as objects onto which we inject a lensed source to build simulated lens systems, as shown in the top panel. All the cutouts are 50 pixels on-a-side, corresponding to 13′′..

In the text

	Fig. 4. Learning curves for the accuracy (top) and loss (bottom) for the training (blue) and validation (red) sets, as a function of epoch. In the case of loss, we crop the y-axis, but the initial training set loss is 0.76.
In the text

Fig. 5. CNN score (SCNN) distribution for the different datasets. Top panel: SCNN distribution of both test sets: lens simulations (dashed line) and real confirmed and candidate lenses (solid line). Both datasets contain images labeled as lenses (purple lines) and real LRGs labeled as non-lenses (red lines). The two test-sets are normalized to their corresponding maximum value in the distribution. We shifted the x-axis of the second test set distribution by 0.02 for clarity. Bottom panel: SCNN distribution for the objects in the parent sample with scores above 0.5. We crop the x-axis for visualization as 99% of the sample is below 0.5.

In the text

	Fig. 6. ROC curve for the test set containing simulations (test set 1 in red) along with the confirmed and candidate lenses (test set 2 in blue); both datasets contain real LRGs as non-lensed examples. The FPR is plotted on a logarithmic scale to aid visualization. The TPR and FPR for SCNN = 0.5 (green) and SCNN = 0.9 (black) are also shown for each set.
In the text

Fig. 7. Heat-map for the normalized mean number of times that each cell in the mosaic was clicked on, among all seven visual inspectors during phase 1 of the visual inspection. The values in each cell were obtained by calculating the mean of the total number of clicks per-cell among the seven users; we normalized these values by a factor of 24.78 that represents the mean number of clicks in a cell for this specific classification.

In the text

Fig. 8. Number of objects classified by a certain number of visual inspectors in phase 1 (top panel: to select lenses and ring galaxies) and phase 2 (bottom panel: to select only the ring galaxies from phase 1). Objects selected by the seven visual inspectors represent a 100% in agreement among users, while the sum along the different bins give us the union, namely, the number of objects selected by at least one user. The exact number in each bin is shown at the top of each bar.

In the text

	Fig. 9. Example of objects classified as “maybe lens” (ML) or “single arc” (SA) and sub-classified as part of the “ring galaxy”, “spiral”, or “merger”, or “sure lens” category is shown in Fig. B.1. At the top of each image is the name of each system, while at the bottom, we show the CNN score (SCNN) and visual inspection score (VIS) obtained in the corresponding category.
In the text

	Fig. 10. Correlation among different classifiers for the categories “sure lens” and “maybe lens” after phase 3. The values in the diagonal represent the total number of objects classified into both categories for each user, while in the adjacent rows and columns the number of objects that both users classified into the same category.
In the text

	Fig. 11. CNN score (SCNN) against the visual inspection score (VISL), with their respective distributions for the final catalog of lens candidates containing 405 systems. The distribution of the category “sure lens” is in purple and “maybe lens” is in red, while in the scatter plot “sure lens” systems are represented by stars and “maybe lens” systems by triangles.
In the text

	Fig. 12. Prior distributions for the effective radius (Reff) and the Sérsic index (ns) used for constraining source light parameter values.
In the text

Fig. 13. Illustration of the automated masking procedure using an example DES image. In the upper row, we show (from left to right) the original image, the image after applying an LoG filter, and the result of setting pixels with flux below the threshold to zero. In the bottom row we show in the leftmost frame the original image annotated with the remaining pixels from the filtering step (red “+” marks), along with the detected peaks (black “x” marks). In the middle we show the estimated size of the lens system with a black circle, as well as the detected bright pixels that are considered contaminant light. These pixels are then used for the mask, and we show the areas covered in the rightmost panel by setting the corresponding pixels to a large constant value.

In the text

Fig. 14. Model best-fit parameter distributions for the lens mass, lens light, and source light profiles. Results for which reduced χ2 ≤ 1.5 are shown in green and those with reduced χ2 > 1.5 are displayed in red. Top left: Einstein radii in arcseconds. Top right: External shear strength, γext. Middle left: Sérsic half-light radii of the lens light. Middle right: Sérsic indices of the lens light. Bottom left: Sérsic half-light radii of the source light. Bottom right: Sérsic indices of the source light.

In the text

	Fig. A.1. Images for the 129 candidates in J19AB that we discarded after visual inspection. Our CNN score is indicated below each object as well as the category allocated by J19AB, in which A stands for “definitely,” B for “probably,” and C for “possibly” a lens.
In the text

	Fig. A.2. Example of cutouts of 56 candidates in J19AB that our CNN graded below 0.8. Our CNN score is indicated below each object as well as the category allocated by J19AB, in which A stands for “definitely,” B for “Probably,” and C for “Possibly” a lens.
In the text

Fig. B.1. Images for the 90 lens candidates in the category Sure lens and their corresponding decomposition performed with MuSCADeT. In the first and fourth columns we have the gri-composite image of the system, the name is on the top, while the CNN score and the visual inspection score (VISL) are displayed at the bottom of each image. Additionally, we marked with an "M" those that we modeled in Sect. 6. Columns 2 and 5 show the subtraction of the blue model from the respective data. Columns 3 and 6 show the subtraction of the red model from the respective data, isolating the lensing features.

In the text

	Fig. B.1. continued.
In the text

	Fig. B.1. continued.
In the text

	Fig. B.1. continued.
In the text

	Fig. B.1. continued.
In the text

	Fig. B.1. continued.
In the text

Fig. B.2. Modeling results for the 52 lens candidates that appear to have only a single lens galaxy acting as a deflector. 1st column: Observed DES image of the lens system in the r-band. 2nd column: Reconstructed image using best-fit model parameters. The black regions are "masked" pixels that are ignored in the modeling as they contain light from contaminant objects in the image. The red curves are the critical lines of the lens model. 3rd column: Normalized residual map showing the difference between the best-fit model and the original data. 4th column: Convergence map of the lens model. 5th column: Reconstructed source light profile (un-lensed). The caustic curves are shown in yellow.

In the text

	Fig. B.2. continued.
In the text

	Fig. B.2. continued.
In the text

	Fig. B.2. continued.
In the text

	Fig. B.2. continued.
In the text

	Fig. B.2. continued.
In the text

	Fig. B.2. continued.
In the text

	Fig. B.2. continued.
In the text

	Fig. B.2. continued.
In the text

	Fig. B.2. continued.
In the text

	Fig. B.2. continued.
In the text

	Fig. B.2. continued.
In the text

	Fig. B.2. continued.
In the text