# Models of the strongly lensed quasar DES J0408-5354

A. Agnello, H. Lin, L. Buckley-Geer, T. Treu, V. Bonvin, F. Courbin, C.

Lemon, T. Morishita, A. Amara, M.W. Auger, et al.

► To cite this version:

A. Agnello, H. Lin, L. Buckley-Geer, T. Treu, V. Bonvin, et al., Models of the strongly lensed quasar DES J0408-5354. Mon.Not.Roy.Astron.Soc., 2017, 472 (4), pp.4038-4050. 10.1093/mnras/stx2242. hal-01645195

HAL Id: hal-01645195

https://hal.science/hal-01645195v1

Submitted on 24 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Models of the strongly lensed quasar DES J0408–5354

A. Agnello,¹★ H. Lin,² L. Buckley-Geer,² T. Treu,³† V. Bonvin,⁴ F. Courbin,⁴

C. Lemon,⁵ T. Morishita,³ A. Amara,⁶ M. W. Auger,⁵ S. Birrer,⁶,³ J. Chan,⁷,⁸

T. Collett,⁹ A. More,¹⁰ C. D. Fassnacht,¹¹ J. Frieman,² P. J. Marshall,¹²

R. G. McMahon,⁵,¹³ G. Meylan,⁴ S. H. Suyu,⁸ F. Castander,¹⁴ D. Finley,²

A. Howell,¹⁵ C. Kochanek,¹⁶ M. Makler,¹⁷ P. Martini,¹⁶ N. Morgan,¹⁸ B. Nord,²

F. Ostrovski,⁶,¹⁹ P. Schechter,²⁰ D. Tucker,² R. Wechsler,²¹ T. M. C. Abbott,²²

F. B. Abdalla,²³,²⁴ S. Allam,² A. Benoit-Lévy,²³,²⁵,²⁶ E. Bertin,²⁵,²⁶ D. Brooks,²³

D. L. Burke,¹²,²⁷ A. Carnero Rosell,²⁸,²⁹ M. Carrasco Kind,³⁰,³¹ J. Carretero,³²

M. Crocce,³³ C. E. Cunha,¹² C. B. D'Andrea,³⁴ L. N. da Costa,²⁸,²⁹ S. Desai,³⁵

J. P. Dietrich,³⁶,³⁷ T. F. Eifler,³⁸ B. Flaugher,² P. Fosalba,¹⁴ J. García-Bellido,³⁹

E. Gaztanaga,¹⁴ M. S. Gill,²⁷ D. A. Goldstein,⁴⁰,⁴¹ D. Gruen,¹²,²⁷ R. A. Gruendl,³⁰,³¹

J. Gschwend,²⁸,²⁹ G. Gutierrez,² K. Honscheid,¹⁶,⁴² D. J. James,²²,⁴³ K. Kuehn,⁴⁴

N. Kuropatkin,² T. S. Li,²,⁴⁵ M. Lima,²⁸,⁴⁶ M. A. G. Maia,²⁸,²⁹ M. March,³⁴

J. L. Marshall,⁴⁵ P. Melchior,⁴⁷ F. Menanteau,³⁰,³¹ R. Miquel,⁴⁸,⁴⁹

R. L. C. Ogando,²⁸,²⁹ A. A. Plazas,³⁸ A. K. Romer,⁵⁰ E. Sanchez,⁵¹ R. Schindler,²⁷

M. Schubnell,⁵² I. Sevilla-Noarbe,⁵¹ M. Smith,⁵³ R. C. Smith,²² F. Sobreira,²⁸,⁵⁴

E. Suchyta,⁵⁵ M. E. C. Swanson,³¹ G. Tarle,⁵² D. Thomas⁹ and A. R. Walker²²

Affiliations are listed at the end of the paper

Accepted 2017 August 30. Received 2017 August 29; in original form 2017 January 30

## ABSTRACT

We present detailed modelling of the recently discovered, quadruply lensed quasar J0408–5354, with the aim of interpreting its remarkable configuration: besides three quasar images (A,B,D) around the main deflector (G1), a fourth image (C) is significantly reddened and dimmed by a perturber (G2) which is not detected in the Dark Energy Survey imaging data. From lens models incorporating (dust-corrected) flux ratios, we find a perturber Einstein radius $0.04 \text{ arcsec} \lesssim R_{E, G2} \lesssim 0.2 \text{ arcsec}$ and enclosed mass $M_p(R_{E, G2}) \lesssim 1.0 \times 10^{10} \text{ M}_\odot$. The main deflector has stellar mass $\log_{10}(M_*/M_\odot) = 11.49_{-0.32}^{+0.46}$, a projected mass $M_p(R_{E, G1}) \approx 6 \times 10^{11} \text{ M}_\odot$ within its Einstein radius $R_{E, G1} = (1.85 \pm 0.15) \text{ arcsec}$ and predicted velocity dispersion 267–280 km s$^{-1}$. Follow-up images from a companion monitoring campaign show additional components, including a candidate second source at a redshift between the quasar and G1. Models with free perturbers, and dust-corrected and delay-corrected flux ratios, are also explored. The predicted time-delays ($\Delta t_{AB} = (135.0 \pm 12.6) \text{ d}$, $\Delta t_{BD} = (21.0 \pm 3.5) \text{ d}$) roughly agree with those measured, but better imaging is required for proper modelling and comparison. We also discuss some lessons learnt from J0408–5354 on lensed quasar finding strategies, due to its chromaticity and morphology.

**Key words:** gravitational lensing: strong–methods: statistical–techniques: image processing–catalogues.

<footnote>* E-mail: aagnello@eso.org</footnote>

<footnote>† Packard Fellow.</footnote>
# 1 INTRODUCTION

Strongly lensed quasars are interesting astrophysical objects for diverse purposes (Courbin, Saha & Schechter 2002). The morphology of the multiple images, accompanied by arcs or rings tracing the lensed host galaxy, enables the description of the mass profile of the lens/deflector galaxy, which typically sits at redshift $z_1 \approx 0.5–1$ (e.g. Oguri, Rusu & Falco 2014). Thanks to magnification, the source can be superresolved, well beyond what is possible for un lensed distant quasars. Astrometric and flux ratio ‘anomalies’ among the multiple images are signatures of luminous and/or dark substructure surrounding the lens (Dalal & Kochanek 2002; Nierenberg et al. 2014), as well as faint features such as extended discs or isophotal twist, boxiness or disciness (Möller, Hewett & Blain 2003; More et al. 2009; Vegetti et al. 2012; Gilman et al. 2016; Hsueh et al. 2016). When the source luminosity varies over time, the time delay between different images can be measured (e.g. Schechter et al. 1997; Tewes, Courbin & Meylan 2013; Bonvin et al. 2016) and used to measure cosmological distances (as originally envisioned by Refsdal 1964, for lensed supernovae) and hence the expansion rate of the Universe, yielding low-redshift ($z_1$) constraints on cosmological parameters that are independent of local distance-scale calibrations (cf. Treu & Marshall 2016; Suyu et al. 2016, and references therein).

Image configuration has a central role for these studies. Systems with four images of the source quasar (hereafter *quads*) provide more information on the mass profiles of the deflector. In contrast, systems with two well-separated images (or *doubles*) can generally be more easily monitored for time variability with ground-based long-cadence observations, since fewer point sources must be blended within the same region. Systems in a *fold* configuration, where two of the quasar images are close to one another, are an interesting transition case: while enabling a robust lens mass reconstruction (see Ding et al. 2017, for a discussion of specific systems), they are also easier to follow up for time delays than other configurations with comparable Einstein radii. In particular, in a *fold* configuration the source lies close to the caustic separating the double and quad regimes, with a *merging pair* of two of the images, thereby giving a highly stretched view of the quasar host near its centre (More et al. 2009; Rusu et al. 2014; Agnello et al. 2016).

Wide-field surveys offer a significant opportunity to discover new systems with suitable configuration, to be followed up for ancillary data. In particular, the Dark Energy Survey (hereafter DES: Sánchez & DES Collaboration 2010) has opened a new window for lens searches in the Southern hemisphere, thanks to a combination of large footprint, depth and good image quality of the Dark Energy Camera (Flaugher et al. 2015; Dark Energy Survey Collaboration et al. 2016).

Here, we detail the first models of a new quasar lens, J0408−5354 (RA = 62.091333, DEC = −53.900266). This lens was discovered by Lin et al. (2017) in the Y1A1 release of DES (Diehl et al. 2014, Drlica-Wagner et al. 2017, in preparation), through a visual inspection of blue objects near red galaxies. Its multiband images show four compact sources, compatible with being point-like, given the DES point-spread-function (PSF), around a luminous red galaxy, as shown in Fig. 1. A spectroscopic confirmation campaign (Lin et al. 2017) shows that the three bright, blue point sources are images of the same source quasar at redshift $z_s = 2.375$, with absorption features at $z_1 = 0.597$ that can be attributed to the lens galaxy. The fourth compact source to the South-West (Fig. 1) is redder than the other confirmed quasar images. Detailed modelling is required to determine whether the anomalous colour is given by

dust extinction, microlensing or an additional red galaxy along the line of sight.

In this follow-up paper, we aim to shed light on the lensing nature of J0408−5354, expanding upon the discovery paper. First, we model the DES images to obtain object positions and spectral-energy distributions (SEDs). The multiband SEDs of the point-sources can be used to quantify chromatic effects (such as microlensing or dust extinction), while the SED of the lens galaxy is used to estimate its stellar mass. The positions of deflector and images are used as inputs to gravitational lens models, whose results are then used to estimate the dark matter content of the lens and verify whether an additional galaxy, lying very close to the reddened compact source along the line of sight, is needed to reproduce the observed flux ratios. We will show that based on the data available so far, the most plausible interpretation of the system consists of a main deflector galaxy and a satellite producing four images of a background lens quasar. The satellite deflector is very well aligned with one of the images, suppressing its flux and contaminating its colours.

This paper is structured as follows. In Section 2, we detail the multiband model results of the DES *grizY* images. A comparison of different lens models is given in Section 3. We conclude in Section 4, including a discussion of the significance of J0408−5354 for different quasar lens searches, and briefly summarize in Section 5. Whenever needed, a standard flat $\Lambda$CDM cosmology is adopted with $\Omega_\Lambda = 0.7$ and $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$.

# 2 SYSTEM CONFIGURATION

J0408−5354 consists of point-like and extended objects (Fig. 1, 0.263 arcsec $px^{-1}$), which are blended in the DES segmentation maps. In order to obtain robust SED measurements, in this section we forward-model the *grizY* image cutouts as a superposition of objects, to recover robust magnitudes and relative positions with realistic uncertainties. The effective exposure time, as calculated by the reduction pipeline of single-epoch images is 540 s (resp. 225 s) for *griz* (resp. *Y*) cutouts. The on-target<fnref n="1" /> PSF full width at half-maximum (FWHM) ranges from 0.56 arcsec in *Y* band to 0.59 arcsec in *g* band. We deployed two independent versions of our own codes to model the cutouts; likelihood is explored via MCMC, whence the maximum-a-posteriori values and their uncertainties are obtained.

Follow-up imaging observations have been conducted with the Wide-Field Imager (WFI) on the 2.2-m telescope in La Silla, to measure the time-delays between the light curves of different images (Courbin et al. 2017). A co-add and optimal deconvolution (following Magain et al. 1998) of the best-seeing images obtained so far, shown in Fig. 2, reveals a more complex structure: besides G1, A, B, G2/C and D, at least three additional ‘blobs’ are visible (G3,G4,G5), as well as a nearly complete Einstein ring with radius $\approx 1.6$ arcsec. Better data are needed to ascertain the nature of this ring and whether G3, G4, G5 are physically connected to it. The $R_c$-band image in Fig. 2 has pixels of 0.12 arcsec per side and point-sources with a FWHM=0.2 arcsec, allowing to locate the position angle (p.a.) of G1 to $\approx 30$ deg E of N. We will discuss these aspects further in the following sections.

<footnote n="1">¹ The DES data processing and calibration system is described by Bertin (2011), Mohr et al. (2012) and references therein; its implementation on the Y1A1 release is described by Drlica-Wagner et al. (2017, in preparation).</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Multiband images of J0408−5354 in *grizY*, from DES single-epoch data with best image quality, plus colour-composites (*gri*, *riz*, *izY*) in the last three columns. The data are shown in the first line with overlaid best-fitting positions, the best-fitting model (as detailed in Section 2) in the second line and the residuals in the third line. An extra source between A and D is visible in the residuals, indicated as ‘G3’ in Fig. 2. Most of the residuals, besides G3, are due to PSF mismatch (around image A) and by blending of B and G2/C. North is up and East is left.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** WFI *R*<sub>c</sub>-band image of J0408−5354 after co-add and deconvolution of the best-seeing images. The central lens galaxy is *G1*. Image A, farther from *G1*, is the first to reach the observer, followed by image B, C and D; lens models will need a perturber *G2* near the location of image C (see Section 3); three additional blobs are visible, marked by blue circles and denoted as G3,G4,G5. *G3* sits on a nearly complete ring of radius ≈1.6 arcsec, indicated by the blue arrow.</figcaption></figure>

## 2.1 Image models

As illustrated in Figs 1 and 2, the system consists of a red galaxy (*G1*) surrounded by three blue point-like objects (A, B, D) and a red-der and compact object (G2/C). As will be quantified in Section 3, if the system is a genuine quad, then G2/C would be a saddle-point image ‘C’, merging with B in a fold-like configuration.

Follow-up Gemini South spectra (programme ID GS-2017A-LP-5, PI E. Buckley-Geer) were taken in 2017 March and April, using the same setup as in Lin et al. (2017), except that B and G2/C were

targeted on the same slit. This slit setup permits spatial decomposition of the spectra and resolves the previous issue of contamination of the G2/C spectrum by light from the much brighter image B. The results confirm that G2/C has the same quasar emission lines as B and the other blue images, albeit with strong differential reddening (Fig. 3).

Given the ordering of stationary points in the Fermat potential of a fold configuration (Saha & Williams 2003), the shortest arrival time corresponds to image A (minimum), followed by B (minimum) and C (first saddle-point) and then D (second saddle-point). For this reason, ‘C’ will be alternatively denoted as the *first saddle-point image* hereafter. Throughout this paper, we will treat this fourth image as an independent object, i.e. will not use its properties directly in constraining the lens models.

The DES cutouts are modelled as the superposition of a galaxy with a Sérsic (1968) profile for G1, and four point sources for A,B,D and G2/C. Different choices for the PSF are available, as it can be adopted from the DES PSF reconstruction or explicitly modelled as a superposition of analytic profiles. Each of these leads to a slight PSF mismatch on pixel-scales, but does not change the results appreciably. In order to test the robustness of the results, we opted for: (i) a model with the DES-reconstructed PSF; and (ii) a model with a Moffat profile (Moffat 1969) fit to a nearby star to determine a parametric PSF. In the model, we impose that the relative displacements of all components (with respect to image A) are the same in every band. The model then comprises: the position angle $\phi_1$, Sérsic index $n_s$ and half-light radius $R_{\text{eff}}$ of G1; the *grizY* positions of A; the relative displacements of G1, B, G2/C and D; and the *grizY* magnitudes of all objects. The Moffat PSF model (ii) includes G3.

The inferred parameters with their uncertainties are listed in Table 1. The differences in estimated magnitudes between the two choices of PSF are an outcome of PSF mismatch (and secondarily pixel scale), and the uncertainties are systematics-dominated. Unfortunately, the depth and image quality of the survey cutouts are not sufficient to constrain the shape parameters of G. Nevertheless, the multiband magnitudes of G1 (marginalized over everything else) are still well constrained. The (broad-band) SEDs of G1 and the four images are shown in Fig. 4. The colours of image G2/C can be obtained by adding a standard reddening law (Cardelli, Clayton &
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Follow-up spectra of images B and G2/C, obtained with a similar Gemini South setup as the discovery spectra of Lin et al. (2017). Both images have the same emission lines, but G2/C is considerably dimmed and reddened with respect to image B.</figcaption></figure>

<table><caption>**Table 1.** Positions (relative to image A) and magnitudes of the objects in J0408-5354, from a joint model of the DES grizY single-epoch images with best image quality, adopting the DES-reconstructed PSF (upper sub-table) or a parametric fit to a nearby star (lower sub-table). Image A is at (RA, DEC) = (62.091323, -53.900289). All the positions have an uncertainty of 0.25 × 10⁻⁴ deg = 0.09 arcsec, smaller than half the DES pixel size (0.27 arcsec), with zero covariance between δRA and δDEC. The naming scheme is illustrated in Fig. 2. With the current depth and image quality, there are degeneracies in the fitted parameters of G3 and those of other components, primarily G1. The zY magnitudes of the ‘blue plume’ G3 are quoted as upper limits.</caption><thead><tr><th>obj.</th><th>δRA(arcsec)</th><th>δDEC(arcsec)</th><th>g</th><th>r</th><th>i</th><th>z</th><th>Y</th></tr></thead><tbody><tr><td>A</td><td>0.0</td><td>0.0</td><td>20.07 ± 0.07</td><td>20.16 ± 0.07</td><td>20.16 ± 0.07</td><td>19.96 ± 0.10</td><td>20.04 ± 0.10</td></tr><tr><td>B</td><td>-6.34</td><td>1.85</td><td>19.98 ± 0.07</td><td>19.95 ± 0.07</td><td>19.74 ± 0.10</td><td>19.28 ± 0.08</td><td>19.34 ± 0.10</td></tr><tr><td>G2/C</td><td>-6.43</td><td>0.75</td><td>22.68 ± 0.20</td><td>21.98 ± 0.15</td><td>21.46 ± 0.15</td><td>20.91 ± 0.12</td><td>20.56 ± 0.16</td></tr><tr><td>D</td><td>-3.12</td><td>2.91</td><td>20.90 ± 0.07</td><td>20.94 ± 0.10</td><td>20.73 ± 0.12</td><td>20.42 ± 0.10</td><td>20.77 ± 0.13</td></tr><tr><td>G1</td><td>-3.31</td><td>1.48</td><td>22.18 ± 0.20</td><td>20.65 ± 0.03</td><td>19.77 ± 0.04</td><td>19.31 ± 0.03</td><td>19.12 ± 0.05</td></tr><tr><td>A</td><td>0.00</td><td>0.00</td><td>20.08 ± 0.01</td><td>20.15 ± 0.01</td><td>20.15 ± 0.02</td><td>19.90 ± 0.07</td><td>19.95 ± 0.14</td></tr><tr><td>B</td><td>-6.35</td><td>1.86</td><td>19.86 ± 0.01</td><td>19.79 ± 0.01</td><td>19.66 ± 0.02</td><td>19.29 ± 0.07</td><td>19.25 ± 0.15</td></tr><tr><td>G2/C</td><td>-6.42</td><td>0.69</td><td>23.16 ± 0.11</td><td>21.61 ± 0.05</td><td>20.92 ± 0.06</td><td>20.82 ± 0.09</td><td>20.45 ± 0.10</td></tr><tr><td>D</td><td>-3.13</td><td>2.96</td><td>20.86 ± 0.02</td><td>20.98 ± 0.02</td><td>20.90 ± 0.03</td><td>20.34 ± 0.07</td><td>20.51 ± 0.15</td></tr><tr><td>G1</td><td>-3.31</td><td>1.58</td><td>22.61 ± 0.16</td><td>20.52 ± 0.06</td><td>19.51 ± 0.06</td><td>19.34 ± 0.07</td><td>19.12 ± 0.08</td></tr><tr><td>G3</td><td>-1.10</td><td>1.63</td><td>22.09 ± 0.16</td><td>21.80 ± 0.17</td><td>21.50 ± 0.21</td><td>>21.20</td><td>>20.85</td></tr></tbody></table>

Mathis 1989, using $R_V = 3.1$ and $E(B-V) = 0.3$) to the SED of image B, but the overall magnitudes need an additional ‘grey’ dimming of 0.8 mags; we also sum the small contribution of a putative galaxy G2 3.5 mags fainter than G1, in order to better reproduce the zY-band fluxes. We will return to these points in Section 3. These results do not depend on the chosen PSF.

## 2.2 Lens stellar mass

The grizY SED inferred for the main deflector galaxy G1 can be used to estimate its stellar mass. We used the public<fnref n="2" /> version of FAST (Kriek et al. 2009). Motivated by Treu et al. (2010), we adopt a Salpeter stellar IMF, which is expected for massive early-type galaxies. A direct measurement of the lens velocity dispersion would enable an IMF-independent determination of the stellar mass (Auger et al. 2009). The best-fitting model is shown in Fig. 4. With the uncertainties from the SED modelling, we obtain $log_{10}(M_*/M_\odot) = 11.49^{+0.46}_{-0.32}$. We will compare this to the results of lens models in the next section. Once again, while we used the more uncertain grizY magnitudes (to remain on the cautionary side),

the inferred stellar mass does not depend appreciably on the PSF chosen to model the cutout.

## 3 LENS MODELS

The three images A, B, D have compatible SEDs, as is also confirmed by their long-slit spectra by Lin et al. (2017). We can then use their positions relative to G1 to model this system as a gravitational lens, obtaining estimates of the total mass (within the Einstein radius) and predicted time-delays between different images. Since G2/C is substantially redder than the other components, we do not include it in the lens model, but rather compare its properties with those predicted by the lens model fit to the other components. The technicalities of the lens model are described in Appendix A.

Conservatively, we adopt 0.2 arcsec positional uncertainties on A, B, D and 0.3 arcsec on G1, G2/C, about twice as large as those from the cutout modelling of Section 2 (relying solely on the DES cutouts). This allows us to explore a wide family of lens models and draw some general conclusions, in particular on the flux ratios allowed by different models. In one case, we also allow the positional uncertainties to be those given directly by the cutout modelling (last line of Table 2). The inferred lens model parameters for all models are given in Table 2. We stress that we are not using the smaller uncertainties from the WFI deconvolution, in order to highlight the

<footnote n="2">² Available at: http://w.astro.berkeley.edu/~mariska/FAST.html</footnote>
<table><caption>**Table 2.** Inferred lens model parameters in the case of a Singular Isothermal Ellipsoid (SIE, first row), a Singular Isothermal Sphere plus external shear (SIS+XS, second row), or the same plus a small perturber G2 where G2/C is observed in the DES cutouts, adopting 0.2 arcsec for the positional uncertainties of A, B, D and 0.3 arcsec for those of G1 and G2. For reference, the notation is fully explained in the main text. The number of degrees of freedom of each model corresponds simply to the number of parameters that are not kept fixed in this table. The lens p.a. of G1 (which may be different from that of its starlight) is quoted in 'mathematical notation' N of W, corresponding to ≈30 deg E of N. The perturber Einstein radius $\theta_{E,p}$ is not an additional parameter, being inferred directly from $b_p$ and $s_p$. Models with a sub-critical G2 ($s_p > b_p/2$) are not excluded. (a) This line shows the average parameters and standard deviations obtained when all uncertainties on positions are set to 0.1 arcsec. (b) The two bottom lines list the lens parameters form a model with G1 plus one (upper) or two (lower) perturbers with free positions, and flux ratios as additional constraints (Section 3.3.2). Letters in square brackets indicate the image next to which the perturbers lie.</caption><thead><tr><th></th><th>θ<sub>E,1</sub></th><th>q</th><th>φ<sub>1</sub><br>(N of W)</th><th>γ<sub>s</sub></th><th>φ<sub>s</sub><br>(N of W)</th><th>b<sub>p</sub></th><th>s<sub>p</sub>/b<sub>p</sub></th><th>θ<sub>E,p</sub></th></tr></thead><tbody><tr><td>SIE</td><td>(1.98 ± 0.08) arcsec</td><td>0.63 ± 0.06</td><td>-60.0 ± 2.0</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SIS+XS</td><td>(1.87 ± 0.08) arcsec</td><td>[1.00]</td><td>-</td><td>0.13 ± 0.3</td><td>29.8 ± 3.3</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SIS+XS+G2</td><td>(1.73 ± 0.15) arcsec</td><td>[1.00]</td><td>-</td><td>0.11 ± 0.03</td><td>18.4 ± 10.1</td><td>(0.33 ± 0.23) arcsec</td><td>0.24 ± 0.21</td><td>(0.26 ± 0.13) arcsec</td></tr><tr><td>SIS+XS+G2<sup>(a)</sup></td><td>(1.72 ± 0.10) arcsec</td><td>[1.00]</td><td>-</td><td>0.10 ± 0.02</td><td>16.5 ± 7.2</td><td>(0.35 ± 0.19) arcsec</td><td>0.30 ± 0.20</td><td>(0.22 ± 0.08) arcsec</td></tr><tr><td>SIE+XS+1p<sup>(b)</sup></td><td>(1.90 ± 0.10) arcsec</td><td>0.67 ± 0.04</td><td>-58.3 ± 1.4</td><td>0.016 ± 0.012</td><td>20.0 ± 3.0</td><td>(0.051 ± 0.008) arcsec</td><td>0.24 ± 0.17</td><td>(0.033 ± 0.013) arcsec [C]</td></tr><tr><td>SIE+XS+2p<sup>(b)</sup></td><td>(1.50 ± 0.10) arcsec</td><td>0.93 ± 0.02</td><td>-59.4 ± 2.3</td><td>0.004 ± 0.032</td><td>19.2 ± 2.2</td><td>(0.40 ± 0.07) arcsec<br>(0.44 ± 0.06) arcsec</td><td>0.14 ± 0.07<br>0.05 ± 0.04</td><td>(0.34 ± 0.07) arcsec [C]<br>(0.42 ± 0.06) arcsec [D]</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Top: grizY magnitudes of the multiple components; red (resp. blue) symbols indicate the galaxy G1 (resp. other compact images A,B,D), while the fainter SED with purple symbols corresponds to G2/C. The SED of image B, once reddened, needs an additional dimming of ≈0.8 mag in all bands to coincide with that of G2/C (black line), to which we also sum the contribution of a galaxy 3.5 mag fainter than G1 as discussed in Section 4. Bottom: Spectrum of the main deflector galaxy G1 from the best-fitting FAST model, yielding $\log_{10}(M_{\star}/M_{\odot}) = 11.49_{-0.32}^{+0.46}$. The observed photometry is given by the dark-green symbols.</figcaption></figure>

robustness of some conclusions that held already with DES-quality data. However, when ellipticity is included in the lens model (defined as 'SIE' below), its p.a. agrees well with that from the WFI images shown in Fig. 2.

The images A, B, D are mapped to the source plane according to the lens equation

$$ \theta_s = \theta_{im} - \alpha - \Gamma\theta_{im}, \quad (1) $$

where $\theta = (\delta x, \delta y)$ is the angular displacement relative to the best-fitting G1 centre from Section 2, the external shear matrix $\Gamma$ is defined as

$$ \Gamma = \gamma_s \begin{pmatrix} \cos(2\varphi_s) & \sin(2\varphi_s) \\ \sin(2\varphi_s) & -\cos(2\varphi_s) \end{pmatrix} \quad (2) $$

and $\alpha$ depends on how we describe the deflections by lensing galaxies. When describing lens galaxies, we use parametric models for their convergence profiles $\kappa = \Sigma/\Sigma_{cr}$, where $\Sigma_{cr} = c^2 D_s/(4\pi G D_t D_{ls})$ accounts for the dimensional dependence on angular-diameter distances. In particular, we use a Pseudo-Isothermal Ellipsoidal Mass Profile (PIEMD; Kassiola & Kovner 1993). This model provides a good representation of the gravitational potential of lens galaxies (e.g. Treu 2010) and the deflection angles $\alpha$ in coordinates $(X, Y)$ aligned with the principal axes of the iso-density ellipsoids

$$ \alpha_X = - \frac{b}{\sqrt{1 - q^2}} \arctan \left( \frac{X \sqrt{1 - q^2}}{s + \sqrt{q^2(s^2 + X^2) + Y^2}} \right), \quad (3) $$

$$ \alpha_Y = - \frac{b}{\sqrt{1 - q^2}} \operatorname{arctanh} \left( \frac{Y \sqrt{1 - q^2}}{q^2 s + \sqrt{q^2(s^2 + X^2) + Y^2}} \right), \quad (4) $$

are fully analytic, together with the convergence and the Fermat potential. The expression in coordinates $(x, y)$ in West-North orientation requires just rotations in the coordinates and deflections, for which we choose the lens long-axis p.a. $\phi_l$ as positive N of W. The spherical ($q \to 1$) and core-less ($s/b = 0$) limit reduces to the Singular Isothermal Sphere (SIS), for which $b$ is also the Einstein radius $R_E$ enclosing a mean convergence of 1. In the Singular Isothermal Ellipsoid (SIE) case ($q < 1$, $s/b = 0$), with the above notation we have $R_E = b/\sqrt{q}$ as the ellipsoidal coordinate of the contour enclosing $\langle \kappa \rangle = 1$. In the case where $q = 1$ but $s/b > 0$, the Einstein radius is $R_E = b\sqrt{1 - 2s/R_E}$, which means that the PIEMD can be sub-critical ($\kappa < 1$ everywhere) when $s > b/2$. The
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Output magnifications from a model with one SIE in the lens plane; the green contours represent the 68percent and 95percent quantiles of the marginalized posterior (no parameters held fixed). While the magnifications of B and D relative to A are in qualitative agreement with the SED fit results, the predicted image C should be almost as bright as image B and appreciably brighter than image A. This is not observed even after differential reddening is added to fit the colours (Section 2), and so it cannot be solely the result of dust extinction.</figcaption></figure>

Einstein radius can be used to estimate the lens velocity dispersion via<fnref n="3" />

$$ \sigma_{\text{sis}} = c \sqrt{\frac{R_E D_s}{4\pi D_l D_{ls}}} = 203 (\theta_E / 1 \text{ arcsec})^{1/2} \text{ km s}^{-1} \quad (5) $$

while the projected mass within $R_E$ is

$$ M_p(R_E) = \pi \Sigma_{cr} R_E^2 = 2.0 * 10^{11} (\theta_E / 1 \text{ arcsec})^2 M_{\odot}, \quad (6) $$

regardless of the lens model. Here and in what follows, $\theta_E = R_E/D_l$ is the Einstein radius in angular units, the same as for the lens strength parameter b.

## 3.1 Models with one deflector

For the first models, we describe the lensing mass distribution as given solely by G1. The first model (SIE) comprises simply a SIE representing G1. The second model (SIS+XS) adopts an SIS for G1, with the addition of external shear with non-null $\gamma_s$. The resulting parameters are listed in Table 2.

Both the SIE and SIS+XS models reproduce the positions of images A,B,D and predict a saddle-point image ‘C’ near the position G2/C found in Section 1 (Fig. 6), whose relative position can vary from model to model, still within one or two DES pixels. The inferred Einstein radius $\theta_{E,l}$ of G1 is slightly less than half the A-to-B image separation (≈2.2 arcsec), due to quadrupole contributions to the deflection either by ellipticity or by shear.<fnref n="4" /> The lens velocity dispersion and mass within $R_E$ can be estimated as (286 ± 6) km s⁻¹ and (7.9 ± 0.6)10¹¹ M☉ [resp. (280 ± 6) km s⁻¹ and (7.0 ± 0.6)10¹¹ M☉] for the SIE (resp. SIS+XS) model.

Models with just one central deflector predict that image ‘C’ should be about as bright as image ‘B’, even with relatively large adopted uncertainties on the image positions (0.2" instead of 0.09 arcsec). This is summarized in Figs 5 and 7, and in Table 3.

## 3.2 Models with perturbers

The first saddle-point predicted by models with one deflector would fall near the position of G2, which however is appreciably redder

than the other images and significantly fainter than predicted even in the reddest bands. Extinction as measured in other lensed quasars (Dai et al. 2006; Mediavilla et al. 2005) does not differ substantially to that measured in the Milky Way and Magellanic Clouds (e.g. Cardelli et al. 1989). However, while the simple addition of a standard reddening law<fnref n="5" /> to the SED of image B can reproduce the colours of image G2/C, it still requires a ‘grey’ dimming of ≈0.8 mag in each band to match its overall magnitudes as in Fig. 5.

### 3.2.1 Fixed perturber at G2/C

Since G2/C lies close to image B, the differential reddening should be produced by a local overdensity, such as a small galaxy, whose lensing effect can also alter the magnification of image C. In general, saddle-points of the Fermat potential are *suppressed*, i.e. dimmed, by the presence of nearby perturbers, whereas minima fluctuate less (Schechter & Wambsganss 2002; Keeton 2003).

For this reason, we add a galaxy at the location of G2/C, which we describe as a PIEMD with $q = 1$. The addition of a perturber at a fixed position increases the number of parameters by two (core size and Einstein radius), making the model underconstrained. However, we can still rely on the priors on positions given by Section 2, and examine the range of parameter configurations that are compatible with the observed image configuration.

For simplicity, and due to the lack of an independent redshift measurement, we place the perturber in the same plane of the main lens G1. In general, models of lenses with four images have degeneracies among the monopole and quadrupole parameters (Kochanek, Schneider & Wambsganss 2006). As verified above, the SIS+XS and SIE models do not differ appreciably in the output image positions and magnifications (Tables 2, 3 and Fig. 7).

The inferred lens parameters of the new model (with a main lens G1 and a perturber G2), given in Table 2, suggest a fairly small (≈0.2 arcsec) Einstein radius and do not rule out a sub-critical perturber. Even when $s_p < b_p/2$, the formation of multiple images and their magnifications depend on the distance between the perturber and image C; even though this does not occur in this configuration, the constraint that multiple images are not produced by the perturber (or not observed) could be used if sharper and deeper data are available. Similarly to the findings of Nierenberg et al. (2014) on a different lens, limits on how massive the perturber can be are

<footnote n="3">³ The numerical prefactors in the second equalities are specific to the redshifts $z_s$, $z_l$ of source and deflector in this particular case.</footnote>

<footnote n="4">⁴ We should caution the reader that shear angles obtained in these models depend on the convention chosen for the orientation and on how to write the shear contribution in equation (1). With this convention, the shear direction corresponds to mass overdensities perpendicular to it. Different choices can yield to shear angles perpendicular to these, i.e. aligned with the p.a. of G1.</footnote>

<footnote n="5">⁵ With $R_V = 3.1$ and $E(B-V) = 0.3$, blueshifting the DES wavebands to the lens rest frame.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Time-delay contours for the case with one SIS plus external shear (left panel) or with the addition of a perturber at G2 (right). Models generally predict the fourth image position within one pixel-length in each direction from G2. This has a magnification comparable to that of image B if no perturber is present nearby. Being a saddle-point image, its magnification is easily suppressed by the presence of a small perturber at G2.</figcaption></figure>

<table><caption>**Table 3.** Inferred logarithmic magnifications for the three models fit to positions of A,B,D and G1, with one SIE (top) or SIS+XS (middle) in the lens plane, or with the addition of a perturber at G2/C (bottom). The first line of each block is $log_{10}(\mu)$ from the best-fitting model, while the second line shows the mean and standard deviation from the MCMC posterior. The last column lists the predicted displacement of image C, in terms of West-ward and North-ward displacements from the best-fitting position of G1 from Section 2 (identified with $\delta x = \delta y = 0$). The positional uncertainties are systematics-dominated, as the predicted position (especially $y_C$) can change appreciably across models.</caption><thead><tr><th>Model</th><th>log<sub>10</sub>μ(A)</th><th>log<sub>10</sub>μ(B)</th><th>log<sub>10</sub>μ(D)</th><th>log<sub>10</sub>μ(C)</th><th>x<sub>C</sub> - x<sub>G1</sub> (")</th><th>y<sub>C</sub> - y<sub>G1</sub> (")</th></tr></thead><tbody><tr><td>SIE</td><td>0.47<br>0.45 ± 0.09</td><td>0.89<br>0.93 ± 0.17</td><td>0.52<br>0.45 ± 0.18</td><td>0.76<br>0.86 ± 0.20</td><td>1.65 ± 0.05</td><td>-0.89 ± 0.03</td></tr><tr><td>SIS+XS</td><td>0.52<br>0.64 ± 0.13</td><td>1.05<br>1.09 ± 0.15</td><td>0.58<br>0.71 ± 0.17</td><td>1.14<br>0.93 ± 0.19</td><td>1.60 ± 0.05</td><td>-0.70 ± 0.05</td></tr><tr><td>SIS+XS+G2</td><td>0.64<br>0.77 ± 0.19</td><td>0.82<br>0.97 ± 0.14</td><td>0.51<br>0.62 ± 0.16</td><td>0.49<br>0.70 ± 0.25</td><td>1.58 ± 0.03</td><td>-0.47 ± 0.07</td></tr></tbody></table>

given simply by the requirement that the other images (A,B,D) are not shifted by the perturber beyond their measured uncertainties.

With the addition of G2 in the lens model, the predicted image 'C' (Fig. 6) is suppressed by the presence of the small perturber, making it slightly fainter than image A (Fig. 8). With a small perturber at ≈0.2 arcsec from image C West-ward and South-ward, its SED can be easily reddened even though it lies very close to B. The small separation between C and G2 makes them hardly distinguishable even in the Gemini acquisition image of Lin et al. (2017), whose PSF has a quoted FWHM≈0.5 arcsec. Within this model, the lens velocity dispersion of G1 is (267 ± 12) km s⁻¹, and its projected mass within $R_E$ is (6.0 ± 1.0)10¹¹ M☉. Even though G2 is not excluded to be sub-critical, we can still estimate its velocity dispersion and enclosed projected mass as (95 ± 17) km s⁻¹ and ≲ 1.0 × 10¹⁰ M☉, respectively.

### 3.2.2 Perturbers with free positions

A perturber (G2) at the location of G2/C, with constraints set only by the relative image positions, is even too effective at saddle-point suppression. In fact, the magnifications predicted by that model (Table 3) are to be meant before dust-extinction, which in this case

is substantial, as mentioned in Section 2. Then, in order to properly assess the presence of a perturber, one must correct the magnitudes from Section 2 for dust-extinction, microlensing and delays in the light curves, and use those as constraints to the lens models. The long-slit spectra (Lin et al. 2017) show some microlensing, in that flux ratios among emission lines are different from flux ratios among the continua, but the amount is small (≈10percent) and, to our aim, negligible in the wavelength range 6000 < $\lambda$ < 10 000 Å. Differential reddening can be inferred quite robustly from the differences in colours among different images, as mentioned above. In order to correct for variability, one needs shifted light curves from a monitoring campaign. In what follows, we use shifted magnitudes from the 2.2m-WFI campaign (Courbin et al. 2017), correct them for differential extinction and neglect microlensing effects. The grey dimming on image C is robust since its delay from image B is small.

Here, we analyse two models with perturbers with free positions, Einstein radii and core radii. In the first case, one perturber near the fold pair (B,C) is added to the primary deflector (G1), and a magnitude difference $m_C - m_B = 2.5log_{10}(\mu_B/\mu_C) = (0.85 \pm 0.1)$ is added as constraint. In the second case, a second perturber is added near image D and the dust-corrected, delay corrected magnitude differences $m_D - m_A = (0.45 \pm 0.05)$, $m_A - m_B = (1.95 \pm 0.05)$
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Output magnifications from a model SIE (top), SIS+XS (middle), and SIS+XS+G2 (bottom, see the text), in logarithmic units. The SIE and SIS+XS models produce similar results, particularly for the predicted ordering of magnifications. Adding a perturber (G2) at the observed location of G2/C dims image C considerably.</figcaption></figure>

<table><caption>**Table 4.** Time-delays predicted by different lens models, adopting flat-ΔCDM cosmology with $H_0 = 70\,\text{km s}^{-1}\,\text{Mpc}^{-1}$. Models using 0.2 arcsec positional errors have large time-delay uncertainties and posteriors that are very skewed from those of the best-fitting solution. The delays from models with free perturbers are in coarse agreement with those subsequently measured by Courbin et al. (2017), but sharper and deeper imaging data are needed for proper modelling and comparison.</caption><thead><tr><th>Model</th><th>Δt<sub>AB</sub><br>(d)</th><th>Δt<sub>BD</sub><br>(d)</th><th>Δt<sub>BC</sub><br>(d)</th></tr></thead><tbody><tr><td>SIE</td><td>125</td><td>20</td><td>8</td></tr><tr><td>SIS+XS</td><td>100</td><td>19</td><td>2.5</td></tr><tr><td>SIS+XS+G2</td><td>87</td><td>29</td><td>6</td></tr><tr><td>SIE+XS+1p</td><td>135.0 ± 12.5</td><td>21.0 ± 3.5</td><td>1.8 ± 0.3</td></tr><tr><td>SIE+XS+2p</td><td>140.0 ± 12.0</td><td>14.0 ± 2.4</td><td>0.0 ± 2.0</td></tr></tbody></table>

(Courbin et al. 2017) are additionally used as constraints. The model with two perturbers is explored because a model with just one perturber predicts higher magnifications for image D than are observed. As a model for G1, we adopt an SIE+XS with the following priors: uniform in $1.3 \, \text{arcsec} < b_l < 2.2 \, \text{arcsec}$, $0.0 < \gamma_s < 0.14$, $0.5 < q < 0.99$; Gaussian in $\phi_l = (-60.0 \pm 5.0)$ deg, $\phi_s = (20.0 \pm 5.0)$ deg. Uniform priors are used for the parameters of the perturbers. The inferred parameters are appended in the two bottom lines of Table 2, and positions sampled from the posteriors are displayed in Fig. 9. The ‘butterfly’ pattern in the allowed perturber positions is common to other lenses (e.g. Nierenberg et al. 2014).

### 3.3 Time-delays

From the lens models, we can also give some forecasts on the expected delays between the arrival times

$$t_i = \frac{(1 + z_l)D_l D_s}{c D_{ls}} \left[ \frac{1}{2} |\theta_{im,i} - \theta_s|^2 - \Phi \right], \quad (7)$$

where the projected potential $\Phi$ is analytic in all models chosen.

Their values are approximately the same across different models, albeit with ≈20 d differences, and the ordering is always the same: the first image is A, followed by B, C shortly after, and finally D. The ordering is general and does not depend on whether a perturber is included in the model, being determined by the configuration of critical points (e.g. Saha & Williams 2003). Within the SIS+XS+G2 model, we have $\Delta t(AB) = 85$ d, $\Delta t(BC) = 6$ d, and $\Delta t(BD) = 29$ d,

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Output magnifications from a model with SIS+XS plus a small perturber (G2) at the observed location of G2/C; again, the green contours represent the 68 percent and 95 percent quantiles of the marginalized posterior. Within this class of models, both images C and D are slightly fainter than image A and significantly fainter than image B, even before dust extinction is accounted for. This saddle-point suppression is excessive on image C, due to the proximity of G2. For this reason, we explore models where one or two perturbers have free position parameters, and flux ratios (corrected for dust and light curve delays) are used as constraints (see the text, Fig. 9).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Positions of G1 (red), perturbers (grey) and model-predicted quasar images (blue) as sampled from the posteriors with one (top) or two (bottom) perturbers with free positions. The black star-symbols mark the image positions (relative to G1) from Section 2.</figcaption></figure>

where $\Delta t(i, j) = t_j - t_i$ is positive when the arrival-time of image i is shorter than that of image j. The quoted values have large uncertainties, due to the wide degeneracies in the lens models, and their marginalized posterior is offset from the values from the best-fitting model. Models with free perturbers have a more symmetric posterior, and smaller uncertainties in the predicted time-delays. Delays like these are ideal for ground based monitoring, because they are long enough to yield 1–3 per cent precision with daily cadence, yet short enough that one or two observing seasons are sufficient. The delays predicted by models with free perturbers are in coarse agreement with those accurately measured by an ongoing, high-cadence 2.2m-WFI monitoring campaign (Courbin et al. 2017). However, the current data and models are not sufficient to allow for a proper comparison, primarily because of long-term microlensing effects and uncertainties in the lens potential.

# 4 DISCUSSION

We have modelled J0408−5354 to obtain the photometry of its individual components, the stellar mass of the main lens galaxy and lens parameters for a choice of plausible models. The predicted time-delays and image configuration make this system amenable to follow-up for time-delay cosmography, as well as for studies of the quasar host near the central engine and substructure near the quasar images.

With the current data quality, there are vast degeneracies in the lens model parameters, which however can be easily relieved with high-resolution imaging data. This will also help locate the perturber G2 responsible for the reddening *and* dimming of image C. The occurrence of both cases would not be uncommon, as seen e.g. for the lens B1608+656 (Myers et al. 1995; Fassnacht et al. 1996; Suyu et al. 2009). A direct measurement of the lens velocity dispersion, together with the recently measured time-delays (Courbin et al. 2017), would yield a direct measurement of the angular-diameter distance to the lens via $D_l \propto c^3 \Delta t/\sigma^2$ (see Paraficz & Hjorth 2009; Jee et al. 2016, for a general discussion).

## 4.1 System configuration

J0408−5354 consists of three blue point-like images of the source quasar at $z_s = 2.375$, and two redder components of which G1, the main lens, is at $z_l = 0.597$, whereas the nature of G2/C is unclear, given its colours and the significant contamination from spectra of other components (Lin et al. 2017). We have modelled this system as a superposition of an extended galaxy (G1) plus four compact sources A, B, D, G2/C and obtained deconvolved SEDs. In particular, the SED of G1 suggests a stellar mass $M_* \approx 3.2 \times 10^{11} M_{\odot}$ for the lens, within $\approx 0.4$ dex uncertainty. There is a degree of systematic uncertainty in the positions and fluxes of B and G2/C, given their proximity and the fact that B is more than a magnitude brighter than G2/C in bluer bands. Some faint residuals are given by PSF mismatch on pixel scales, regardless of the choice of PSF. The WFI images with best seeing, once deconvolved, show what could be an additional ring with radius $R \approx 1.6$ arcsec, which is appreciably smaller than that inferred from lens models based on the positions of images A,B,D (as summarized below). If it corresponds to a second source, it can map to a unique location only if it lies at a redshift between those of quasar and deflector.

## 4.2 Lens model properties

The positions of images A, B, D relative to G1 have been used to explore lens models of J0408−5354. Models with one main lens, adopted as SIE or SIS plus external shear, predict an Einstein radius $\approx 1.9$ arcsec and a lens p.a. $\approx 30$ deg East of North, or equivalent quadrupole contributions from $\gamma_s = 0.1$ shear. Both models, while successful at reproducing the positions of images A, B, D, would predict a saddle-point image where G2/C lies, and about as bright as image B, which is not observed even in band Y or after differential extinction corrections. Models with a small perturber at the location of G2 reproduce the same image positions, but are able to suppress image C by about a magnitude even before dust is accounted for. The dust-corrected and delay-corrected flux ratios, as well as the images measured from the DES cutouts, are well reproduced by two small perturbers near the saddle-points.

The projected mass within the Einstein radius is $M_p(R_E) \approx 6 \times 10^{11} M_{\odot}$, about twice the stellar mass estimated from the SED of G1. A proper evaluation of the dark matter fraction in the lens,
however, would require a measurement of the effective radius of G1. When the perturber has non-null Einstein radius, its enclosed mass is $M_p(G2) \approx 1.3 \times 10^{10} M_*$.
The contribution of a small galaxy with magnitudes $m(G2) = m(G1) - 2.5\log_{10}(M_p(G2)/M_p(G1))$ is barely noticeable in *gri* bands. It can be well reproduced by reddening and offsetting the SED of image B by $\approx 0.8$ mag, and makes the *zY*-band magnitudes of image G2/C in complete agreement with the values measured from Section 2 (black line in Fig. 4).

The estimated time-delays between images A,B,D (from lens models with free perturbers) are somewhat in agreement with what is measured by Courbin et al. (2017) through a 2.2m-WFI dedicated monitoring campaign. Still, given the uncertainties on image positions and few constraints, the derived uncertainties are appreciable and higher-resolution imaging data will be required to tighten the model-predicted uncertainties on the delays.

If indeed two sources are present at different redshift, J0408−5354 can also be used to measure Dark Energy cosmological parameters via the ratio of distance ratios $D_s/D_{ls}$ to the different sources (e.g. Paczynski & Gorski 1981; Soucail, Kneib & Golse 2004; Collett et al. 2012), besides time-delay cosmography to measure $H_0$. The only other system with time delays and multiple source-planes that is known and studied to date is the galaxy cluster MACSJ1149.5+2223 (Treu et al. 2016).

## 4.3 Relevance of J0408−5354 for lens searches

The photometry and configuration of J0408−5354 make it an interesting testbed for different techniques of lensed quasar candidate selection. These, in turn, have implications for substructure studies, as the composition of lens-selected or source-selected samples affects the sensitivity to substructure, especially for lens searches that are tailored on simple lenses or on systems dominated by ‘isolated’ quasar SEDs.

Like the serendipitous quad of More et al. (2016), J0408−5354 was originally found by visual inspection of objects selected solely on *gri* survey properties, instead of relying on hybrid infrared ‘excess’ colours (Warren, Hewett & Foltz 2000) that have been used to target quasars (Maddox & Hewett 2006; Maddox et al. 2012; Peters et al. 2015) or lenses (Ofek et al. 2007) and applied in other lens searches in DES (Agnello et al. 2015a; Ostrovski et al. 2017). After the initial discovery via the blue-near-red search of Lin et al. (2017), different teams have examined their own search methods. Here we provide a summary of the different findings.

### 4.3.1 Cutout classification: CHITAH

CHITAH (Chan et al. 2015) examines the image cutouts of objects to detect at least two blue compact sources and a red galaxy, evaluating how plausible the configuration is as a strong lens via the corresponding source-plane $\chi^2$. This approach relies on the requirement that the blue images have very similar SEDs, distinct from the lens SED.

When applied to the *grizY* cutouts of J0408−5354, it did not flag this system as a possible quad since the fourth image G2/C is significantly redder than the others. However, based on A, B and G1, it did classify this system as a possible double. These findings suggest that pixel-based automatic recognition, such as CHITAH or LENSTRACTOR<fnref n="6" /> could be made more flexible by accounting for possible SED variations of the predicted images.

### 4.3.2 Target selection: data mining

The first technique used to select lensed quasars in the DES relied upon Artificial Neural Networks (ANNs) trained on SDSS *griz* and WISE (Wright et al. 2010) W1, W2 bands of four main classes of objects (Agnello et al. 2015a). Despite the success of the first discovery results (Agnello et al. 2015b), further improvements could be made for wider application to DES, as discussed in Appendix B. With these new ANNs, J0408−5354 was automatically flagged as an extended quasar with $z_s > 1.75$, one of the two classes (besides ‘lens’) to be retained for visual inspection.<fnref n="7" /> Despite the improvement in the ANNs and the *blind* re-discovery of J0408−5354, there is considerable scatter in the SDSS-DES translated magnitudes, which can cause some interesting objects to slip out of the selection boundaries (and false positives to leak in). The outlier selection method (Agnello 2017), in which J0408−5354 is rediscovered as a $>3\sigma$ outlier among quasars and with low probability to be a galaxy, is somewhat immune from this issue, as are Population Mixture classifications (Williams et al. 2016; Ostrovski et al. 2017).

## 5 SUMMARY

J0408−5354 has an interesting fold-like image configuration, with three well-separated images (A,B,D) and a fourth one (C) in a merging pair with the brightest image (B). Besides the three, clearly identifiable blue images of the source quasar, a fourth component G2/C is fainter than simple lens-model predictions and appreciably red. While image B is already redder than the farthest image A, with $\Delta(Y - r) \approx 0.65$ compatible with a simple (Cardelli et al. 1989) reddening law with $E(B - V) = 0.3$, image G2/C is further reddened (additional $E(B - V) = 0.3$) and also requires a grey dimming of $\approx 0.8$ mag in every band.

A small perturber ($R_{E,p} \approx 0.2$ arcsec, $M_p \approx 1.0 \times 10^{10} M_{\odot}$) near the location of G2/C explains both the needed reddening and dimming over the whole *grizY* range. The data on positions and (dust-corrected, delay-corrected) flux ratios allow for two perturbers near the saddle-point images C and D, albeit with very uncertain masses (Table 2) due to the current scarcity of constraints.

The image separation makes this system particularly apt to time-delay measurements. The B-A delay $\approx 120-140$ d predicted by the full models (SIE+XS+pert.) is in good agreement with the results of a dedicated 2.2m-WFI monitoring campaign (Bonvin et al. 2017, in preparation). The lens mass within the Einstein radius $R_E = 1.73$ arcsec is $M_p \approx (6.0 \pm 1.0) \times 10^{11} M_{\odot}$, about twice the stellar mass $M_* \approx 3 \times 10^{11} M_{\odot}$ of the main galaxy G1.

The chromaticity and morphology of J0408−5354 mean that different search techniques, while successfully flagging it as a lens candidate, are triggered by different features. Also, the peculiar colours and configuration of the quasar images are a powerful reminder that automated search techniques should be flexible enough to encompass these systems, in view of homogeneous lens-selected or source-selected samples for follow-up science. Oguri & Marshall (2010) estimated 1146 quasar lenses within a depth of $i = 23.6$ in the 5000 deg$^2$ final DES footprint, of which 14 percent are quads. Past and ongoing lens searches show that a suite of complementary techniques are needed to maximize the number of detected lenses, especially at magnitudes fainter than $i \approx 19$.

The composition of J0408−5354, with a primary (massive) lens and a small perturber and a merging image-pair, make it both an

<footnote n="6">⁶ Available at https://github.com/davidwhogg/LensTractor</footnote>

<footnote n="7">⁷ In particular, the blend D+G1 with catalogue ID = 3070264166, RA = 62.0904688061, DEC = −53.8996413857</footnote>
interesting system for follow-up and a rather peculiar system to model. Spectroscopic and high-resolution imaging observations would enable more accurate models, both for cosmography and for substructure studies, and a highly magnified view of the source quasar and its host.

## ACKNOWLEDGEMENTS

AA and TT acknowledge support by the Packard Foundations through a Packard Research Fellowship and by the National Science Foundation through grant AST-1450141.

This paper was written as part of the STRong lensing Insights into the Dark Energy Survey (STRIDES) collaboration, a broad external collaboration of the Dark Energy Survey, http://strides.astro.ucla.edu

We are grateful to the anonymous referee for a 'fresh' and thorough reading of the paper, which helped clarify several points.

Funding for the DES Projects has been provided by the DOE and NSF(USA), MISE(Spain), STFC(UK), HEFCE(UK). NCSA(UIUC), KICP(U. Chicago), CCAPP(Ohio State), MIFPA(Texas A&M), CNPQ, FAPERJ, FINEP (Brazil), MINECO(Spain), DFG(Germany) and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne Lab, UC Santa Cruz, University of Cambridge, CIEMAT-Madrid, University of Chicago, University College London, DES-Brazil Consortium, University of Edinburgh, ETH Zürich, Fermilab, University of Illinois, ICE (IEEC-CSIC), IFAE Barcelona, Lawrence Berkeley Lab, LMU München and the associated Excellence Cluster Universe, University of Michigan, NOAO, University of Nottingham, Ohio State University, University of Pennsylvania, University of Portsmouth, SLAC National Lab, Stanford University, University of Sussex, and Texas A&M University. The DES Data Management System is supported by the NSF under Grant Number AST-1138766. The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2012-39559, ESP2013-48274, FPA2013-47986 and Centro de Excelencia Severo Ochoa SEV-2012-0234. Research leading to these results has received funding from the ERC under the EU's 7th Framework Programme including grants ERC 240672, 291329 and 306478.

## REFERENCES

* Agnello A., 2017, MNRAS, 471, 2013

* Agnello A., Kelly B. C., Treu T., Marshall P. J., 2015, MNRAS, 448, 1446

* Agnello A. et al., 2015, MNRAS, 454, 1260

* Agnello A., Sonnenfeld A., Suyu S. H., Treu T., Fassnacht C. D., Mason C., Bradač M., Auger M. W., 2016, MNRAS, 458, 3830

* Auger M. W., Treu T., Bolton A. S., Gavazzi R., Koopmans L. V. E., Marshall P. J., Bundy K., Moustakas L. A., 2009, ApJ, 705, 1099

* Bertin E., 2011, in Evans I. N., Accomazzi A., Mink D. J., Rots A. H., eds, ASP Conf. Ser. Vol. 422, Astronomical Data Analysis Software and Systems XX. Astron. Soc. Pac., San Francisco, p. 435

* Bonvin V., Tewes M., Courbin F., Kuntzer T., Sluse D., Meylan G., 2016, A&A, 585, A88

* Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245

* Chan J. H. H., Suyu S. H., Chiueh T., More A., Marshall P. J., Coupon J., Oguri M., Price P., 2015, ApJ, 807, 138

* Collett T. E., Auger M. W., Belokurov V., Marshall P. J., Hall A. C., 2012, MNRAS, 424, 2864

* Courbin F., Saha P., Schechter P. L., 2002, in Courbin F., D. Minniti D., eds, Lecture Notes in Physics, Vol. 608, Gravitational Lensing: An Astrophysical Tool. Springer-Verlag, Berlin, p. 1

* Courbin F. et al., 2017, A&A, preprint (arXiv:1706.09424)

* Dai X., Kochanek C. S., Chartas G., Mathur S., 2006, ApJ, 637, 53

* Dalal N., Kochanek C. S., 2002, ApJ, 572, 25

* Dark Energy Survey Collaboration et al., 2016, MNRAS, 460, 1270

* Diehl H. T. et al., 2014, in Peck A. B., Benn C. R., Robert L. Seaman R. L., eds, Proc. SPIE Conf. Ser. Vol. 9149, Observatory Operations: Strategies, Processes, and Systems V. SPIE, Bellingham, p. 91490V

* Ding X. et al., 2017, MNRAS, 465, 4634

* Fassnacht C. D., Womble D. S., Neugebauer G., Browne I. W. A., Readhead A. C. S., Matthews K., Pearson T. J., 1996, ApJ, 460, L103

* Flaugher B. et al., 2015, AJ, 150, 150

* Gilman D., Agnello A., Treu T., Keeton C. R., Nierenberg A. M., 2017, MNRAS, 467, 3970

* Hsueh J.-W., Fassnacht C. D., Vegetti S., McKean J. P., Spingola C., Auger M. W., Koopmans L. V. E., Lagattuta D. J., 2016, MNRAS, 463, L51

* Jee I., Komatsu E., Suyu S. H., Huterer D., 2016, JCAP, 4, 031

* Kassola A., Kovner I., 1993, in Surdej J., Fraipont-Caro D., Gosset E., Refsdal S., Remy M., eds, Proc. 31st Liege International Astrophysical Colloquia, Gravitational Lenses in the Universe. Institut d'Astrophysique, Universite de Liege, Liege, p. 571

* Keeton C. R., 2003, ApJ, 584, 664

* Kochanek C. S., Schneider P., Wambsganss J., 2004, in Meylan G., Jetzer P., North P., eds, Proc. 33rd Saas-Fee Advanced Course, Part 2 of Gravitational Lensing: Strong, Weak & Micro. Springer-Verlag, Berlin

* Kriek M., van Dokkum P. G., Labbé I., Franx M., Illingworth G. D., Marchesini D., Quadri R. F., 2009, ApJ, 700, 221

* Lin H. et al., 2017, ApJL, 838, L15

* Maddox N., Hewett P. C., 2006, MNRAS, 367, 717

* Maddox N., Hewett P. C., Péroux C., Nestor D. B., Wisotzki L., 2012, MNRAS, 424, 2876

* Magain P., Courbin F., Sohy S., 1998, ApJ, 494, 472

* Mediavilla E., Muñoz J. A., Kochanek C. S., Falco E. E., Arribas S., Motta V., 2005, ApJ, 619, 749

* Moffat A. F. J., 1969, A&A, 3, 455

* Mohr J. J. et al., 2012, in Radziwill N. M., Chiozzi G., eds, Proc. SPIE Conf. Ser. Vol. 8451, Software and Cyberinfrastructure for Astronomy II. SPIE, Bellingham, p. 84510D

* Möller O., Hewett P., Blain A. W., 2003, MNRAS, 345, 1

* More A., McKean J. P., More S., Porcas R. W., Koopmans L. V. E., Garrett M. A., 2009, MNRAS, 394, 174

* More A. et al., 2016, MNRAS, 465, 2411

* Myers S. T. et al., 1995, ApJ, 447, L5

* Nierenberg A. M., Treu T., Wright S. A., Fassnacht C. D., Auger M. W., 2014, MNRAS, 442, 2434

* Ofek E. O., Oguri M., Jackson N., Inada N., Kayo I., 2007, MNRAS, 382, 412

* Oguri M., 2010, PASJ, 62, 1017

* Oguri M., Marshall P. J., 2010, MNRAS, 405, 2579

* Oguri M., Rusu C. E., Falco E. E., 2014, MNRAS, 439, 2494

* Ostrovski F. et al., 2017, MNRAS, 465, 4325

* Paczynski B., Gorski K., 1981, ApJ, 248, L101

* Paraficz D., Hjorth J., 2009, A&A, 507, L49

* Peters C. M. et al., 2015, ApJ, 811, 95

* Refsdal S., 1964, MNRAS, 128, 307

* Rusu C. E., Oguri M., Minowa Y. Iye M., More A., Inada N., Oya S., 2014, MNRAS, 444, 2561

* Saha P., Williams L. L. R., 2003, AJ, 125, 2769

* Sánchez E., DES Collaboration, 2010, J. Phys.: Conf. Ser., 259, 012080

* Schechter P. L., Wambsganss J., 2002, ApJ, 580, 685

* Schechter P. L. et al., 1997, ApJ, 475, L85

* Sérsic J. L., 1968, Cordoba, Argentina: Observatorio Astronomico

* Soucail G., Kneib J.-P., Golse G., 2004, A&A, 417, L33

* Suyu S. H., Marshall P. J., Blandford R. D., Fassnacht C. D., Koopmans L. V. E., McKean J. P., Treu T., 2009, ApJ, 691, 277

* Suyu S. H. et al., 2016, MNRAS, 468, 2590

* Tewes M., Courbin F., Meylan G., 2013, A&A, 553, A120

* Treu T., 2010, ARA&A, 48, 87

* Treu T., Marshall P. J., 2016, A&AR, 24, 11
Treu T., Auger M. W., Koopmans L. V. E., Gavazzi R., Marshall P. J., Bolton A. S., 2010, ApJ, 709, 1195

Treu T. et al., 2016, ApJ, 817, 60

Vegetti S., Lagattuta D. J., McKean J. P., Auger M. W., Fassnacht C. D., Koopmans L. V. E., 2012, Nature, 481, 341

Warren S. J., Hewett P. C., Foltz C. B., 2000, MNRAS, 312, 827

Williams P., Agnello A., Treu T., 2017, MNRAS, 466, 3088

Wright E. L. et al., 2010, AJ, 140, 1868

## APPENDIX A: LENS MODELLING SPECIFICS

Regardless of the model specifics, all images must map to the same source-position. For each choice of the lens model parameters, a source at $\theta_s$ in the source plane corresponds to images $\theta_i$ in the image plane, and the goodness of fit can be described by the image-plane $\chi^2$

$$ \chi_{ip}^2 = \sum_{i=1}^{3} \frac{|\theta_i - \theta_{im,i}|^2}{\delta_i^2} = \sum_{i=1}^{3} \frac{|\mathbf{A}_i(\theta_s - \theta_{s,i})|^2}{\delta_i^2}, \quad (A1) $$

where $\theta_{im,i}$ and $\theta_{s,i}$ are the measured image-positions and their model-predicted source-plane positions for images A,B,D, $A_i = \partial\theta_{im,i}/\partial\theta_{s,i}$ is the magnification tensor around each image and $\delta_i$ is the positional uncertainty on image $i$. The second equality relies on the fact that, near a reasonable lens solution, we can linearize the lens equation around the measured image positions. Then, within this approach, the best-fitting source-position and corresponding image-plane $\chi^2$ and image-positions can be found by straightforward linear operations. Its validity has been tested extensively by Oguri (2010). Based on those functional tests and ours, the lens models are accurate provided one iterates $\approx 30$ times between the linearized best-fitting $\theta_s$ and $\theta_{im,i}$, ensuring that all model image-positions are re-mapped to the same source-point.

Writing the $\chi^2$ as above relies on a Gaussian distribution of the measured image positions, with isotropic positional uncertainties, and is equivalent to drawing image positions with infinite precision from Gaussians $G(\theta_{im,i}, \delta_i)$, considering (for each choice) a highly-penalized image-plane $\chi^2 = p\chi_{ip,1}^2$ in the lens model with

$$ \chi_{ip,1}^2 = \sum_{i=1}^{3} |\theta_i - \theta_{im,i}|^2 \quad (A2) $$

and $p \gg \delta_i^{-2}$. This allows us to generalize the lens model likelihood to image configurations that do not have isotropic and Gaussian uncertainties. In particular, we can draw the relative displacements of G1, B and D with respect to image A as given by the likelihood explored in Section 2, which we call $\mathcal{L}_{SED}$. At very high values of $p$, the only parameter combinations that are explored are those that correspond to all image positions mapping back to the same source position, because other configurations are heavily penalized.

Another hypothesis underlying this approach is that the measured image position uncertainties are simply given by the extraction of Section 2, so that each image carries a weight proportional to its (squared) magnification in the $\chi^2$. This does not account for systematic uncertainties in the image positions given by the proximity of different objects and PSF mismatch. This problem is evident for the brightest image B, which would instead carry the highest weight in $\chi_{im}^2$. We then opt for a penalized source-plane $\chi^2$ of the form

$$ \chi_{sp}^2 = p \sum_{j=1}^{3} |\theta_{s,j} - \langle \theta_s \rangle|^2, \quad (A3) $$

where $\langle \theta_s \rangle = (\theta_{s,A} + \theta_{s,B} + \theta_{s,C})/3$ for each choice of the model parameters, and consider the lens-model likelihood as

$$ \mathcal{L} \propto \mathcal{L}_{SED}(\theta) \times e^{-\chi_{sp}^2/2}. \quad (A4) $$

The penalty parameter $p$ is gradually increased, until all possible models are effectively producing images originating from the same source-position, within milli-arcsecond tolerance, and the model uncertainties are driven by $\mathcal{L}_{SED}$.

## APPENDIX B: MINING ACROSS SURVEYS

The original implementation of ANNs was based upon SDSS data and four main classes of objects. In order to be more widely applicable to DES, it was improved in three ways: (i) more object classes, including multiple redshift intervals for the ‘quasar’ class to distinguish low-redshift contaminants from higher-redshift objects; (ii) less restrictive colour-cuts, that would otherwise exclude known lenses with higher $g-i$ or lower $W1-W2$; and (iii) accounting for the differences in photometry between SDSS and DES via a cross-calibration valid for blue extended objects.<fnref n="8" /> The best-fitting regressions have

$$ \begin{aligned} g_{\text{des}} &= g_{\text{sdss}} + 0.05, & r_{\text{des}} &= r_{\text{sdss}} + 0.088, \\ i_{\text{des}} &= i_{\text{sdss}} + 0.112, & z_{\text{des}} &= z_{\text{sdss}} + 0.159, \end{aligned} \quad (B1) $$

for the $p \ll f$ magnitudes, and

$$ \begin{aligned} g_{\text{des}} &= g_{\text{sdss}} + 0.165 - 0.092(g_{\text{des}} - r_{\text{des}} - 0.4) \\ r_{\text{des}} &= r_{\text{sdss}} + 0.118 - 0.215(g_{\text{des}} - r_{\text{des}} - 0.4) \\ i_{\text{des}} &= i_{\text{sdss}} + 0.04 - 0.2(i_{\text{des}} - z_{\text{des}}) \\ z_{\text{des}} &= z_{\text{sdss}} + 0.078 - 0.044(z_{\text{des}} - Y_{\text{des}} - 0.17) \end{aligned} \quad (B2) $$

for the model magnitudes. There is considerable scatter (0.11–0.18 mag) in the translated magnitudes, given by the extendedness of the objects and different depth and image quality between SDSS and DES. This means that interesting candidates (resp. contaminants) can leak out of (resp. within) the hyperplanes defining class boundaries as identified by the ANN classification.

<footnote n="1">¹European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching bei München, Germany</footnote>

<footnote n="2">²Fermi National Accelerator Laboratory, Batavia, IL 60510, USA</footnote>

<footnote n="3">³Department of Physics and Astronomy, PAB, 430 Portola Plaza, Box 951547, Los Angeles, CA 90095-1547, USA</footnote>

<footnote n="4">⁴Laboratoire d'Astrophysique, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, CH-1290 Versoix, Switzerland</footnote>

<footnote n="5">⁵Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="6">⁶Institute for Astronomy, Department of Physics, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093, Zurich, Switzerland</footnote>

<footnote n="7">⁷Institute of Astronomy and Astrophysics, Academia Sinica, PO Box 23-141, Taipei 10617, Taiwan</footnote>

<footnote n="8">⁸Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85741 Garching, Germany</footnote>

<footnote n="9">⁹Institute of Cosmology & Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK</footnote>

<footnote n="10">¹⁰Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan</footnote>

<footnote n="11">¹¹Department of Physics, University of California Davis, 1 Shields Avenue, Davis, CA 95616, USA</footnote>

<footnote n="12">¹²Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA 94305, USA</footnote>

<footnote n="8">⁸ We refer to Agnello et al. (2015b) for the definition of ‘blue and extended’ in this case.</footnote>
<footnote n="13">¹³Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="14">¹⁴Institut de Ciències de l'Espai, IEEC-CSIC, Campus UAB, Carrer de Can Magrans, s/n, 08193 Bellaterra, Barcelona, Spain</footnote>

<footnote n="15">¹⁵Department of Physics, University of California, Santa Barbara, CA 93106, USA</footnote>

<footnote n="16">¹⁶Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="17">¹⁷ICRA, Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, CEP 22290-180, Rio de Janeiro, RJ, Brazil</footnote>

<footnote n="18">¹⁸Staples High School, 70 North Ave, Westport, CT 06880, USA</footnote>

<footnote n="19">¹⁹CAPES Foundation, Ministry of Education of Brazil, Brasília - DF 70040-020, Brazil</footnote>

<footnote n="20">²⁰MIT Kavli Institute for Astrophysics and Space Research, 37-664G, 77 Massachusetts Avenue, Cambridge, MA 02139</footnote>

<footnote n="21">²¹Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA</footnote>

<footnote n="22">²²Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, Casilla 603, La Serena, Chile</footnote>

<footnote n="23">²³Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="24">²⁴Department of Physics and Electronics, Rhodes University, PO Box 94, Grahamstown 6140, South Africa</footnote>

<footnote n="25">²⁵CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014 Paris, France</footnote>

<footnote n="26">²⁶Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014 Paris, France</footnote>

<footnote n="27">²⁷SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA</footnote>

<footnote n="28">²⁸Laboratório Interinstitucional de e-Astronomia - LineA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil</footnote>

<footnote n="29">²⁹Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil</footnote>

<footnote n="30">³⁰Department of Astronomy, University of Illinois, 1002 W. Green Street, Urbana, IL 61801, USA</footnote>

<footnote n="31">³¹National Center for Supercomputing Applications, 1205 West Clark St., Urbana, IL 61801, USA</footnote>

<footnote n="32">³²Institut de Física d'Altes Energies, Universitat Autònoma de Barcelona, E-08193 Bellaterra, Barcelona, Spain</footnote>

<footnote n="33">³³Institut de Ciències de l'Espai, IEEC-CSIC, Campus UAB, Carrer de Can Magrans, s/n, 08193 Bellaterra, Barcelona, Spain</footnote>

<footnote n="34">³⁴Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA</footnote>

<footnote n="35">³⁵Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India</footnote>

<footnote n="36">³⁶Excellence Cluster Universe, Boltzmannstr. 2, 85748 Garching, Germany</footnote>

<footnote n="37">³⁷Faculty of Physics, Ludwig-Maximilians University, Scheinerstr. 1, 81679 Munich, Germany</footnote>

<footnote n="38">³⁸Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA</footnote>

<footnote n="39">³⁹Instituto de Física Teorica UAM/CSIC, Universidad Autonoma de Madrid, 28049 Madrid, Spain</footnote>

<footnote n="40">⁴⁰Department of Astronomy, University of California, Berkeley, 501 Campbell Hall, Berkeley, CA 94720, USA</footnote>

<footnote n="41">⁴¹Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA</footnote>

<footnote n="42">⁴²Department of Physics, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="43">⁴³Astronomy Department, University of Washington, Box 351580, Seattle, WA 98195, USA</footnote>

<footnote n="44">⁴⁴Australian Astronomical Observatory, North Ryde, NSW 2113, Australia</footnote>

<footnote n="45">⁴⁵George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA</footnote>

<footnote n="46">⁴⁶Departamento de Física Matemática, Instituto de Física, Universidade de São Paulo, CP 66318, CEP 05314-970, São Paulo, SP, Brazil</footnote>

<footnote n="47">⁴⁷Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA</footnote>

<footnote n="48">⁴⁸Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain</footnote>

<footnote n="49">⁴⁹Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra, Barcelona, Spain</footnote>

<footnote n="50">⁵⁰Department of Physics and Astronomy, Pevensey Building, University of Sussex, Brighton BN1 9QH, UK</footnote>

<footnote n="51">⁵¹Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, Spain</footnote>

<footnote n="52">⁵²Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="53">⁵³School of Physics and Astronomy, University of Southampton, Southampton SO17 1BJ, UK</footnote>

<footnote n="54">⁵⁴Universidade Federal do ABC, Centro de Ciências Naturais e Humanas, Av. dos Estados, 5001, Santo André, SP 09210-580, Brazil</footnote>

<footnote n="55">⁵⁵Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA</footnote>

<footnote>This paper has been typeset from a TeX/LATeX file prepared by the author.</footnote>