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

HAL
open science

# A gravitationally lensed supernova with an observable two-decade time delay

Steven A. Rodney, Gabriel B. Brammer, Justin D. R. Pierel, Johan Richard, Sune Toft, Kyle F. O'Connor, Mohammad Akhshik, Katherine E. Whitaker

► To cite this version:

Steven A. Rodney, Gabriel B. Brammer, Justin D. R. Pierel, Johan Richard, Sune Toft, et al. A gravitationally lensed supernova with an observable two-decade time delay. Nature Astronomy, 2021, 5, pp.1118-1125. 10.1038/s41550-021-01450-9. insu-03711494

HAL Id: insu-03711494

https://insu.hal.science/insu-03711494v1

Submitted on 7 Oct 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_2.png" /></figure>

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

# A gravitationally lensed supernova with an observable two-decade time delay

Steven A. Rodney ID 1✉, Gabriel B. Brammer ID 2,3✉, Justin D. R. Pierel<fnref n="1" />, Johan Richard⁴, Sune Toft ID 2,3, Kyle F. O'Connor<fnref n="1" />, Mohammad Akhshik⁵ and Katherine E. Whitaker ID 2,6

When the light from a distant object passes very near to a foreground galaxy or cluster, gravitational lensing can cause it to appear as multiple images on the sky<fnref n="1" />. If the source is variable, it can be used to constrain the cosmic expansion rate² and dark energy models³. Achieving these cosmological goals requires many lensed transients with precise time-delay measurements⁴. Lensed supernovae are attractive for this purpose because they have relatively simple photometric behaviour, with well-understood light curve shapes and colours—in contrast to the stochastic variation of quasars. Here we report the discovery of a multiply imaged supernova, AT 2016jka ('SN Requiem'). It appeared in an evolved galaxy at redshift 1.95, gravitationally lensed by a foreground galaxy cluster⁵. It is probably a type Ia supernova—the explosion of a low-mass stellar remnant, whose light curve can be used to measure cosmic distances. In archival Hubble Space Telescope imaging, three lensed images of the supernova are detected with relative time delays of <200 d. We predict that a fourth image will appear close to the cluster core in the year 2037 ± 2. Observation of the fourth image could provide a time-delay precision of ~7 d, <1% of the extraordinary 20 yr baseline. The supernova classification and the predicted reappearance time could be improved with further lens modelling and a comprehensive analysis of systematic uncertainties.

We discovered AT 2016jka using data from the Hubble Space Telescope (HST) programme Resolved Quiescent Magnified Galaxies (REQUIEM, HST-GO-15663, principal investigator M.A.)⁶ (see a summary of observations in Supplementary Table 1). The REQUIEM project targets massive galaxies with low specific star formation rates that have been magnified by strong gravitational lensing. The brightest and most spectacular galaxy targeted by REQUIEM is MRG-M0138, a massive red galaxy (MRG) at redshift z=1.95 (ref. ⁷) behind the galaxy cluster MACS J0138.0-2155 (ref. ⁸). MRG-M0138 is quadruply lensed by a foreground galaxy cluster at z=0.338. During analysis of observations obtained on 13–14 July 2019 we discovered three point sources that were present in archival HST images from 18–19 July 2016, part of the programme (HST-GO-14496, principal investigator Newman) that first confirmed the MRG-M0138 galaxy as a strongly lensed object (Observations). Each point source is within 5 arcsec of one of the four MRG-M0138 images. None of the three point sources are present in the REQUIEM HST data in 2019 (Fig. 1). We infer that these are multiple images of a single astrophysical transient in MRG-M0138, most probably a supernova (SN).

To construct a lens model for the MACS J0138.0-2155 cluster we use the LENSTOOL software⁹,¹⁰ (Lens modelling). To avoid unintended bias, we kept the lens model development completely separate from the analysis of the SN. Only upon completion of both were the results combined for the analysis described here. The input model constraints are the positions and redshifts of the MRG-M0138 galaxy at z=1.95 (both the galaxy's centroid position and the SN location in each image) as well as a multiply imaged background galaxy at z=0.766, both having secure spectroscopic redshifts (Supplementary Table 3). We model the mass distribution in the cluster core as the combination of cluster-scale and galaxy-scale potentials (Supplementary Table 4 and Extended Data Fig. 1). From this model we derive estimates for the lensing magnification and time delay of each of the SN images, including two predicted future images (Table 1). The lens model predicts that the SN should appear in the fourth MRG-M0138 image in the year 2037 ± 2, demagnified with μ=0.4 ± 0.2. A fifth image will also appear at a still later date, located near the centre of the cluster and much more strongly demagnified, so it will not be easily observable. We anticipate that future lens modelling of the cluster will improve on these predictions primarily by exploring a wider range of mass models and incorporating more observational constraints (Supplementary Note: Future work). For example, our analysis adopted only a single form for the density profiles, and did not incorporate constraints from stellar kinematics or pixel-level surface brightness data from the multiply imaged systems. Although our LENSTOOL model does a good job of reproducing the morphology of the host galaxy images H1–H3, it does not reproduce image H4 as well (Supplementary Note: Host image morphology comparison).

If we can estimate the age of each SN image then we can derive direct observational constraints on the relative lensing time delays. For this goal, it is helpful to have a firm determination of the transient's class. Expected time delays and magnifications from the lens model exclude any of the various rapidly evolving and low-luminosity stellar transient classes, strongly suggesting that it is a SN. The first-order SN distinction remaining is between a type Ia SN—the explosion of a white dwarf star in a binary system— and a core collapse SN (CCSN)—the end-point of a star with mass > 10 $M_\odot$. The properties of the host galaxy can inform this classification because CCSN are limited to galaxies with young stellar populations. Limits on the specific star formation rate and age for this host, MRG-M0138, show it to be a massive but very quiescent and evolved galaxy, unlikely to retain any substantial population of high-mass stars⁷. On the basis of observed properties of the host

<footnote n="1">¹Department of Physics and Astronomy, University of South Carolina, Columbia, SC, USA. ²Cosmic Dawn Center (DAWN), Copenhagen, Denmark. ³Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark. ⁴Centre de Recherche Astrophysique de Lyon, ENS de Lyon, Université de Lyon, CNRS, Saint-Genis-Laval, France. ⁵Department of Physics, University of Connecticut, Storrs, CT, USA. ⁶Department of Astronomy, University of Massachusetts, Amherst, MA, USA. ✉e-mail: srodney@sc.edu; gabriel.brammer@nbi.ku.dk</footnote>
<figure><img src="image_4.png" /></figure>

**Fig. 1** | Overview of the MACS J0138.0-2155 cluster field and the AT 2016jka discovery. **a**, The wide-field view is 40" on a side, with tick marks indicating 10" intervals. Labels indicate the locations of the observed AT 2016jka images (SN1-SN3), the expected future images (SN4 and SN5) and the multiply imaged host galaxy (H1-H4). Labels 3.1-3.4 indicate the locations of a separate multiply imaged [O II] emitter at z = 0.77 used to help constrain the lensing potential. **b-i**, 4" cut-outs around the lensed SN images with 1" tick marks. **b-e**, The imaging from July 2016, where the SN was visible. **f-i**, The later imaging from July 2019, where the SN has faded away. Panels **a**, **e** and **i** include the location of the fourth image predicted to appear in -2037, with an ellipse demarcating the 68% probability region predicted from the LENSTOOL model, and **a** the expected location of the final and highly demagnified fifth image ('SN5'). The three-colour images are generated from the Wide Field Camera 3 infrared filters with the blue channel from F105W or F110W, green from F125W and red from F160W or F140W (as indicated in **d** and **h**). All panels use the late-epoch F125W imaging for the green channel; nevertheless, it is immediately clear that the SN2 image is substantially bluer than the other two, which helps to constrain the relative age of each SN image and the transient classification as a probable type Ia SN explosion.

<table><caption>**Table 1** | SN lensing observables</caption><thead><tr><th>Image</th><th>RA<sup>a</sup> (hh:mm:ss)</th><th>Dec.<sup>a</sup> (dd:mm:ss)</th><th>Age<sup>b</sup> (d)</th><th>Δt<sup>c</sup><sub>obs</sub> (d)</th><th>Δt<sup>c</sup><sub>pred</sub> (d)</th><th>μ</th></tr></thead><tbody><tr><td>SN1</td><td>01:38:03.63</td><td>-21:55:50.38</td><td>92<sup>+21</sup><sub>-19</sub></td><td>—</td><td>—</td><td>5 ± 1</td></tr><tr><td>SN2</td><td>01:38:02.96</td><td>-21:55:47.26</td><td>-24<sup>+16</sup><sub>-7</sub></td><td>114<sup>+28</sup><sub>-31</sub></td><td>82 ± 62</td><td>7 ± 3</td></tr><tr><td>SN3</td><td>01:38:02.42</td><td>-21:55:38.47</td><td>107<sup>+26</sup><sub>-21</sub></td><td>-17<sup>+19</sup><sub>-16</sub></td><td>-19 ± 34</td><td>3.9 ± 0.5</td></tr><tr><td>SN4</td><td>01:38:04.15 ± 0.36</td><td>-21:55:24.73 ± 0.43</td><td>—</td><td>—</td><td>7,742 ± 540<sup>d</sup></td><td>0.4 ± 0.2</td></tr><tr><td>SN5</td><td>01:38:03.77 ± 0.1</td><td>-21:55:31.74 ± 0.1</td><td>—</td><td>—</td><td>9,463 ± 770<sup>d</sup></td><td>&lt;0.01</td></tr></tbody></table>

galaxy alone, we find a 62–75% probability that AT 2016jka is a type Ia SN (Classification).

Adopting the lens model magnifications for the three observed SN images (Table 1), we can locate each SN image in colour–magnitude space (Fig. 2). After magnification correction, all three images are still brighter than expected for a type Ia SN, which may indicate that a lens-modelling degeneracy is at play. Nevertheless, the magnification-corrected AT 2016jka data are more consistent with the type Ia population than any CCSN subclass (Fig. 2a,b and Extended Data Fig. 2). AT 2016jka also demonstrates the expected evolution of a type Ia SN colour and brightness over ~100 d (Fig. 2c). By also including the model-predicted time delays, we can treat the three SN images effectively as three points on a common SN light curve, and we find $P(Ia)=94\%$ (Classification).

This composite light curve is shown in Fig. 3, with the best-fitting type Ia SN model. Extended Data Fig. 3 also shows a random draw of light curves from the Monte Carlo sampling for all three major SN subclasses. For the remainder of this analysis, we proceed under the assumption that AT 2016jka is indeed a type Ia SN. An improved classification could be achieved with spectroscopy and multiband photometry upon arrival of the fourth image. In that case the analysis that follows here could be revised to achieve similar results with a different underlying SN model.

The colour of a type Ia SN evolves substantially over its lifetime as the photosphere expands and cools, revealing different layers of the expanding shell and driving episodes of recombination¹¹. Since the phenomenon of gravitational lensing in general is achromatic, this colour evolution makes it possible to derive an age constraint

<footnote><sup>a</sup>Coordinates are given in the J2000.0 reference frame, as measured for images 1-3 and as predicted for images 4 and 5. Uncertainties in the predicted position of images 4 and 5 are in arc seconds.</footnote>

<footnote><sup>b</sup>The measured age of each SN image is the number of days relative to peak brightness, in observer-frame days. A negative value means the image was observed before peak brightness. <sup>c</sup>Time delays are reported relative to image SN1. <sup>d</sup>The combination of the models for the light curve evolution and the lens time delays implies a date of 2037 ± 2 for the peak brightness of image 4 and 2042 ± 4 for image 5 (which will probably be undetectable).</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 2** | Classification information for AT 2016jka on the basis of its position in colour-magnitude space. **a** The observed photometry for the three SN images is shown as open markers. Vertical dotted lines show the lens magnification corrections using our preferred lens model, with filled grey symbols indicating the corrected (demagnified) magnitude. Error bars indicate the observational and systematic uncertainty, including the range of alternative magnification corrections encompassed by lens model variations. Contours show the population distributions for normal SNe of type Ia (red), type Ib/Ic (gold) and type II (green), drawn with two contour levels enclosing 68% and 95% of each SN population. Each SN subclass was simulated at $z=1.95$, and samples from their expected light curves were drawn uniformly in time. **b** Marginalized distributions along the colour dimension for the three SN subclasses (using the same colour scheme). The simulated populations have been scaled according to the expected explosion rates in the SN host galaxy, based on its stellar-population properties, so the y-axis scaling is effectively a relative probability. Open symbols show the observed colours again, at arbitrary y position. **c** Zoomed-in view of the colour-magnitude space marked by the grey box in **a**. The evolution of a typical type Ia SN at $z=1.95$ is shown by a coloured line, with the line colour indicating SN age in observer-frame days relative to peak brightness. White diamonds correspond to the times labelled on the colour bar below. Grey shading shows the typical range of luminosities and colours observed for the type Ia SN population in the nearby universe. Although the magnification-corrected data are brighter than expected for most type Ia SNe, they are consistent both with the overall type Ia SN population and with the type Ia SN colour-magnitude versus time curve.</figcaption></figure>

that is largely independent of the lens model (Colour curve age constraints and Extended Data Figs. 4 and 5). Combining this information with magnification constraints from the lens modelling helps break parameter degeneracies, yielding measured delays in Table 1 (Light curve + lens model age constraints and Extended Data Figs. 6 and 7).

Using these measured time delays, we created a reconstructed form of the intrinsic light curve and colour curve of AT 2016jka, shown in Fig. 3. Remarkably, the ages of images 1 and 3 are constrained to better than ±20 d, despite having only a single epoch of photometric data. These uncertainties may be further reduced when the future fourth image is observed with high-precision, multi-epoch photometry. Such a light curve will pin down the intrinsic SN light curve parameters that are shared by all images, and break remaining parameter degeneracies. Improvements to the lens modelling will also be essential, to better estimate and minimize systematic biases that may arise from the necessary magnification and time-delay corrections.

AT 2016jka and other lensed supernovae (SNe) like it could eventually contribute to mapping the cosmic expansion history and

measuring the effects of dark energy, which appears to be driving an accelerating cosmic expansion rate$^{12,13}$. Recent investigations of the expansion rate of the universe (the Hubble-Lemaître constant, $H_0$) have found that measurements from the local universe are significantly different from the value inferred from measurements of the cosmic microwave background radiation$^{14,15}$. The community is actively attempting to resolve this ‘$H_0$ crisis’ by mitigating systematic uncertainties or discovering new physics from the early universe$^{16}$. Either resolution will require multiple independent cosmological probes. In recent years, lensed quasar time delays have provided a valuable independent tool for this, with seven high-precision measurements to date$^{17}$. As the sample of time-delay lenses grows to several dozen (including lensed SNe such as AT 2016jka), it is expected to deliver a measurement of $H_0$ with 1% precision$^{18}$.

Looking beyond the $H_0$ crisis, determining the nature of dark energy and how it may evolve over time is a primary goal for the large-scale cosmology experiments of the 2020s$^{19,20}$. A future sample with ~100 well measured lensing time delays would be a competitive tool for dark energy studies$^{4,18}$. Events such as AT 2016jka could be an important part of this time-delay cosmology sample, but to date there have been only two lensed SNe observed with multiple images. The first, SN Refsdal, was a peculiar type II SN whose image with the longest delay was missed$^{21,22}$. The second, SN 2016geu, was a type Ia SN with short delays that make high-precision time-delay measurements impossible$^{23,24}$. An earlier discovery of an unusually luminous SN was also shown to be a strongly lensed type Ia SN$^{25}$, though the multiple images were not resolved. This makes AT 2016jka just the third discovery of a lensed SN resolved into multiple images.

Future large-scale surveys such as the Vera C. Rubin Observatory and the Nancy Grace Roman Space Telescope will observe dozens to hundreds of lensed SNe over their mission lifetimes. The vast majority of these will be lensed by galaxy-scale deflectors$^{26,27}$ and thus will have substantially shorter delays (of order 10–100 d). Since it is the fractional time-delay uncertainty that propagates through to any time-delay distance measurement, the extraordinarily long time delays of cluster-lensed SNe such as AT 2016jka can deliver substantially better time-delay precision, with comparable observational cost (Supplementary Note: Future discoveries). In fact, the long time baseline of lensed SNe such as AT 2016jka effectively insures that their cosmological precision is not limited by time-delay measurement uncertainty. Cluster-scale lenses are more complex than galaxy-scale lenses, but they generally have ‘independent’ measurements of the magnification from several multiply imaged systems in the same field. Modelling cluster lenses is very different from galaxy lenses, so objects such as AT 2016jka can provide a valuable check on systematics for the larger sample of transients used in time-delay cosmology.

The first multiply imaged SN discovery, SN Refsdal, has shown that time-delay cosmography with a cluster-lensed SN is viable$^{22,28}$. However, the observing campaign for SN Refsdal was extraordinary, deploying more than 75 HST orbits over 3 yr. Considerable observational investment has also been required for high-precision time-delay measurement of lensed quasars, such as decade-long programmes$^{29}$ or daily monitoring for high-cadence light curves$^{30}$. In the case of AT 2016jka it will be possible to achieve similar time-delay precision over the 20 yr baseline with just a single imaging epoch as the anchor point. A sample of AT 2016jka-like events could be developed with regular monitoring of cluster-scale lenses, partnered with modest follow-up to characterize any lensed SN discovered (Supplementary Note: Future work).

HST observations enabled us to find this SN. We anticipate that HST may be deorbited and make its final plummet to Earth around the time of the reappearance of AT 2016jka, so we coin the name SN Requiem as an ode to the vast new discovery space that HST continues to unveil.
<figure><img src="image_6.png" /></figure>

**Fig. 3** | The reconstructed light curve and colour curve for AT 2016jka. **a,b,** The light curve (a) and the colour curve (b). The grey shaded region covers the 68% confidence interval (CI) of the best-fit type Ia model (using the SALT2 model$^{49}$), with the median model shown as a solid line. Observed photometric data are shown as coloured symbols. The x-dimension error bars on each data point represent the observed time-delay uncertainties (the 'Age' column of Table 1). The y-dimension uncertainties in both panels incorporate the photometric uncertainty. For the light curve (a), the y error bars also include uncertainty in the lens model magnification (Table 1). The colour measurements (y values, **b**) do not require lens model correction. Parameters for the best-fit SALT2 models in grey are shown in the legend, and were obtained using the joint posterior of the colour curve and light curve methods (Time-delay estimation and Extended Data Fig. 7). These fit parameters are $x_i$, the SALT2 'stretch' parameter; c, SALT2 colour parameter; $M_{B}$, absolute magnitude of a SN Ia.

## Methods

**Observations.** The observations used in this work are summarized in Supplementary Table 1. We processed all HST observations using the DrizzlePac software utilities$^{31}$, aligned to a common astrometric reference frame and resampled to 0.1 arcsec per pixel. We then identified isolated and unsaturated stars in each image and used them to create an effective point spread function with ×4 oversampling, using the photutils package from the Astropy software suite$^{32}$. To measure the SN photometry we followed two tracks. As our primary method we performed effective point spread function fitting directly on the F105W (Y band) and F160W (H band) images where the SN was apparent. This fitting allowed for a constant background flux to account for both the sky brightness and the background light of the cluster and host galaxy. As a second approach, we created 'pseudodifference images' by rescaling the F110W and F140W images collected in 2019 (in which the SN is not present). The transmission functions of the F110W and F140W filters are broader than those of F105W and F160W, and do not strictly overlap in wavelength. The optimal scaling factor to produce a clean subtraction therefore depends on the spectral energy distribution of the source. We set the scaling to 0.62 and 1.17 for F110W to F105W and F140W to F160W, respectively. These values produced visually clean subtractions of the SN host galaxy MRG0138—meaning that they minimize the residual flux in the pseudodifference images. We then performed effective point spread function fitting on the SN in each pseudodifference image, as before. Both sets of photometry agree to within 1 s.d. The reported photometry in Supplementary Table 2 is the measurements from the first method (collected directly from the unsubtracted images). Also reported in Supplementary Table 2 are flux densities and uncertainties measured within circular apertures of diameter 0.7" at the position of the SN3 image in the 2019 HST visits, where the SN is not detected.

**Very Large Telescope spectroscopy.** We make use of integral field spectroscopic data obtained on the cluster core of MRG0138 with Very Large Telescope (VLT)/ Multi-Unit Spectroscopic Explorer (MUSE), publicly available as part of the programme 0103.A-0777(A) (principal investigator Edge). Three exposures of 970 s each were taken with a small dithering offset and 90° rotations in between. This dataset was reduced and analysed using the MUSE data reduction pipeline v.2.7 (ref. 33) for basic calibration (bias, flat field, wavelength, line spread function, geometry) as well as flux calibration, sky subtraction and astrometry. We also make use of the self-calibration technique$^{34}$ to remove illumination systematics, specifically tuned for the case of crowded fields in the central region of galaxy clusters$^{15}$. The combined datacube is then processed through ZAP$^{35}$, which applies a principal component analysis technique to remove sky subtraction residuals. The final datacube covers the central $1 \times 1$ arcmin$^2$ around the cluster centre with $0.2^\circ \times 0.2^\circ \times 1.25$ Å$^3$ pixels.

We have extracted spectra for each HST-detected source and inspected them for redshift measurements. In addition, we have run the muselet software

(part of the MPDAF package$^{37}$, mpdaf.readthedocs.io) to search for line emitters not directly associated with HST sources$^{38,39}$. Apart from the lensed quiescent galaxy, we measured spectroscopic redshifts for cluster members and one ring-like background galaxy north of the brightest cluster galaxy (BCG) at $z=0.766$. Finally, we have measured the velocity dispersion of the BCG to be $390 \pm 10$ km s$^{-1}$. This measurement was not available for use in our blind lens modelling, but could provide a useful constraint for future lens model improvements (Supplementary Note: Future work).

**Lens modelling.** To model the mass distribution in the MACS J0138.0-2155 cluster core we use the latest version of LENSTOOL$^*$ (git-cral.univ-lyon1.fr/lenstool), which performs a Bayesian analysis with a Markov chain Monte Carlo sampler to estimate the best fit and uncertainty on each parameter of the mass distribution.

The strong-lensing constraints used are the locations of multiple images found in HST and VLT/MUSE. We group them into three systems: (1) the four images of the quiescent galaxy hosting MRG0138-SN, (2) the three observed images of AT 2016jka assumed to be at the same redshift and (3) the diffuse arc-like structure identified in HST and confirmed as an [O I]λ3727 emitter in the VLT/MUSE datacube (Observations). The image coordinates and redshifts used in the lens model are summarized in Supplementary Table 3.

The cluster mass modelling is performed similarly to that for other massive strong-lensing clusters observed with HST$^{40}$. The total mass distribution is parametrized as a combination of multiple dPIE (double pseudoisothermal elliptical) profiles describing both cluster-scale and galaxy-scale dark matter haloes. dPIEs are elliptical isothermal profiles with both a core radius and a cut radius, where the density flattens and drops, respectively. In the case of MRG0138 the mass distribution is dominated by a single mass concentration centred on its BCG. We therefore use a single cluster-scale halo at a fixed cut radius of 1 Mpc. We add a single galaxy-scale halo on each cluster member, where the shape parameters (halo centre, ellipticity) are fixed to their measured HST morphology and their core radius is negligible (fixed at 0.15 kpc). Extended Data Fig. 1 shows the locations of all components of the cluster model, including 32 cluster members. Cluster members were identified by the combination of red sequence selection (based on the F814W–F160W colour) and MUSE spectroscopy.

The majority of cluster members are elliptical galaxies selected from the red sequence, and to reduce the number of parameters we assume that they follow the scaling relations: $\sigma = \sigma^* (L^*/L_{\odot})^{(1/4)}$ for the velocity dispersion, and $r_{\text{cut}} = r_{\text{cut}}^* (L^*/L_{\odot})^{(1/2)}$ for the cut-off radius, assuming the Faber–Jackson relation and a constant mass-to-light ratio respectively. $\sigma^*$ and $r_{\text{cut}}^*$ are model parameters for a cluster member at the characteristic luminosity $L^*$. Following the discussion in ref. 41 we fix $\sigma^* = 158$ km s$^{-1}$ and $r_{\text{cut}}^* = 45$ kpc.

We individually optimize the $\sigma$ and $r_{\text{cut}}$ parameters for four specific galaxies that are not expected to follow the aforementioned scaling relations: the BCG and three perturbers P1 to P3. These perturbers are either blue gas-stripped galaxies infalling
into the cluster core, and/or located very close to the images of the AT 2016jka host galaxy, perturbing its apparent morphology with additional lensing. This choice of perturbers is similar to the ones used in the model in ref. ⁷.

LENSTOOL optimizes the parameters of the model by minimizing the overall root mean square dispersion between the predicted and observed locations of the multiple images. The best fit parameters of each mass component are provided in Supplementary Table 4. Uncertainties reported there are derived from the Markov chain Monte Carlo models, by sampling their posterior probability distribution.

We developed five lens model variants blindly (that is, without knowing the impact of each lens model variation on the transient classification or time-delay inferences). Model A was the first viable model developed; it did not include additional perturbers, and did not include the additional lensed background source at $z=0.7663$. In Model B we allowed the location of the main cluster dark matter halo to be free, with an offset from the reference position taken at the BCG centre. Model C allowed the same central position offset and also relaxed the constraints on the BCG ($\sigma$ and $r_{\text{crit}}$). Model D fixed the primary dark matter halo at the BCG centre, but still relaxed the constraints on the BCG $\sigma$ and $r_{\text{crit}}$. The final model, and the one selected as the preferred model before unblinding, is model E, which includes all four perturbers described above, and includes the additional background object at $z=0.7663$. Note that all of the lens model variants A–D would give a slightly larger magnification for all three SN images. This means that our estimated systematic uncertainties are one sided, as can be seen in Fig. 2 and Extended Data Fig. 2.

This model is then used to predict the magnification and time delays for the three observed images of AT 2016jka, as well as the locations of the fourth and fifth images, which are still to appear. These predictions are summarized in Table 1, but note that the magnification predictions for image SN4, as well as all predictions for SN5, should be treated with caution, as the lens models presented here have known shortcomings. For example, the best-fit velocity dispersion for the BCG is $700 \text{ km s}^{-1}$ (Supplementary Table 4), though we have measured this property from MUSE spectra to be $390 \pm 10 \text{ km s}^{-1}$. Further lens modelling is needed to incorporate such additional constraints and to fully quantify potential systematic biases (Supplementary Note: Future lens modelling).

The final LENSTOOL model E reproduces all multiple system positions with a root mean square of 0.15". The $1\sigma$ uncertainty from model E for the SN4 location is the ellipse in Fig. 1e,i. As the lens model reproduces the location of the SN images within a small uncertainty, these predictions are computed with LENSTOOL using the barycentre of all source positions corresponding to images SN1, SN2 and SN3 as the same reference source position. Comparing our lens modelling with the previous model of this cluster from ref. ⁷, we find that the magnification estimates are broadly consistent, though systematically lower (Supplementary Note: Comparison to previous lens modelling).

**Classification.** The lens modelled time delays between the images are ~100 observer-frame days, but we see that three images of the transient are visible simultaneously. From this we can infer that the visibility time of the transient in the $z=1.95$ rest frame must be at least ~30 d. Similarly, with expected magnifications in the vicinity of $\mu \sim 10$, the measured apparent magnitudes near 23 AB mag translate to a rest-frame absolute magnitude near $M_B \sim -19.5$ mag (too bright to be a nova, luminous blue variable or other low-luminosity stellar transient). Taken together, these indicators strongly suggest that the transient is a SN.

Although we have invoked the lens model in this analysis, we note that the inferences are not strongly dependent on the specific lens model predictions. To make the observed transient images consistent with a fast or low-luminosity transient, the time delays and/or magnifications would have to be changed by more than a factor of 2. In the analysis to follow, we will work under the assumption that AT 2016jka is a SN.

**SN subclassification on the basis of host galaxy.** With this transient identified as a SN, we now seek to identify the most likely SN type, under the assumption that it belongs to one of the three most common subclasses (Ia, II, Ib/c). We first use two methods that rely only on measured properties of the host galaxy to circumstantially infer the type. This inference is less strongly dependent on the lens model, helping to reduce any bias associated with a lens-model-dependent classification.

Although type Ia SNe are found in all types of galaxy, CCSN are limited to galaxies with relatively young stellar populations. We can therefore infer some information about the SN type using the observed host properties combined with knowledge of the relative rates of type Ia and CCSN in different stellar populations⁴⁻⁶. In the case of the host galaxy MRG0138, we have a very well constrained spectral energy distribution extending out to far-infrared wavelengths with Spitzer Infrared Array Camera data⁷,⁴³. From the spectral energy distribution fitting we derived the host galaxy's rest-frame $B-K$ colour and absolute magnitude, $M_K$, which serve as proxies for the stellar-population age and have been empirically calibrated with SN rates in the local universe⁴⁴. We adopt a lensing magnification correction using LENSTOOL model E to obtain $M_K$. The $B-K$ colour is not affected by the foreground lens. Using the galsnid method⁴⁵ we derive a 75% probability that AT 2016jka is of type Ia (Supplementary Table 5). We used host galaxy image 2 for this purpose, with the flux-weighted harmonic mean

$\mu = 8.3$ to derive $M_K$ (see Supplementary Table 6). We also evaluated the other host galaxy images and found no change in the resulting SN classification probability.

As an alternative host galaxy classification constraint, we use the derived properties of the host galaxy stellar population directly, rather than adopting colour and magnitude proxies. The MRG-M0138 galaxy has high mass ($\log_{10}(M/M_{\odot}) = 11.7$) but is a very quiescent galaxy, with a specific star formation rate of $\sim 10^{-11.3} \text{ yr}^{-1}$ and a stellar population that is well matched by an exponential star formation history with an age of 1.4 Gyr (ref. ⁷). The massive stars that end as CCSN explosions have main-sequence lifetimes of <40 Myr (ref. ⁶), making it unlikely that CCSN progenitors make up a large fraction of the MRG-M0138 stellar population—though the high total stellar mass makes it possible that pockets of young stars are present. We define classification probabilities based on the projected SN rate for each SN subclass, derived from the host galaxy's stellar mass and star formation rate⁴⁶. This yields a 62% probability that AT 2016jka is of type Ia (Supplementary Table 5).

**SN subclassification on the basis of SN photometry.** To improve the classification of the AT 2016jka subtype, we now bring in observed photometry of the SN itself, and again we adopt two methods. The first method uses only magnification information from the lens model, and the second uses both the modelled magnification and time-delay predictions. In both cases we adopt the stellar-population-based host galaxy classification probabilities as priors.

Extended Data Fig. 2 illustrates the first approach. After applying the magnification corrections, each of the three images of the SN is mapped to colour-magnitude space. We then treat each observed point (each SN image) separately, comparing its colour-magnitude location with a simulated population of un lensed SNe. Our simulation uses the sncosmo package⁷ to generate 10,000 SNe for each of the three principal SN subclasses (Ia, Ib/c, II), all at $z=1.95$. We then compute the number of simulated SNe within a rectangular region around each observed point. The width of this sampling region is set to three times the observed colour uncertainty, and the height is equal to the lens-modelling magnification uncertainty. The $\mu$ uncertainty used here includes an estimate of the systematic uncertainty for each SN image, derived from the spread of magnifications across lens model variants (similar to the methodology of ref. ⁷). Note, however, that the lens model variations evaluated here do not vary the assumption of the density profile, which can strongly affect magnifications (see also Supplementary Note: Comparison to previous lens modelling).

We take the number of simulated SNe for each type within the sampling region as an estimate of the likelihood that AT 2016jka belongs to that class. Note that in this case there is no need to apply a cut to the simulated sample to account for detectability, because the $5\sigma$ limiting magnitude of our HST observations is 26.5 AB mag, and after accounting for magnification of ~1.5 mag (Table 1) this becomes $m_{\text{lim}} \sim 28$ mag. This means that all the points shown in Extended Data Fig. 2 (and therefore all simulated SNe entering our classification counts) would be easily detectable in our HST imaging. Multiplying by the prior probabilities derived from host galaxy properties, we finally derive the probability that the SN is of type Ia as $P(\text{Ia}) = 0.92, 0.98$ and $0.95$ from the three SN images SN1, SN2 and SN3, respectively. Supplementary Table 5 reports the mean of our three classification probabilities for each subclass.

As a second photometric classification of this SN, we used the STARDUST2 Bayesian light curve classification tool¹⁸, which is also built on the underlying sncosmo framework. Here we adopt both the predicted magnifications and time delays from the best lens model, which allows us to put the photometry from the three images together as a composite 'light curve' and compare against simulated light curves. STARDUST2 uses the SALT2-extended model to represent type Ia SNe⁴⁹,⁵⁰ and a collection of 42 spectrophotometric time series templates to represent CCSN (27 type II and 15 type Ib/c). These CCSN templates comprise all of the templates developed for the Supernova Analysis software SNANA⁵¹, derived from the SN samples of the Sloan Digital Sky Survey⁵²⁻⁵⁴, Supernova Legacy Survey⁵⁵ and Carnegie Supernova Project⁵⁶⁻³⁸. With STARDUST2 we use a nested sampling algorithm to measure likelihoods over the SN simulation parameter space. Extended Data Fig. 3 shows the magnification- and time-delay-corrected photometry of AT 2016jka and a random sampling of light curve models from the sncosmo nested sampling algorithm employed by STARDUST2. Nested sampling is a Monte Carlo method that traverses the likelihood space in a manner that samples the Bayesian likelihood⁵⁹. These sample light curves and colour curves therefore give a visual representation of how well each SN subclass (Ia, II, Ib/c) can match the observed data. This figure shows that the limited photometric data can be reasonably well fitted by at least one model from any of these three subclasses. The density of curves in the top panels demonstrates that the type Ia SN model is consistently a good match to the data. However, the range of model parameters that allow such a fit to the data is much more limited for the heterogeneous CCSN types. To compute the posterior probability distribution we adopt priors for each of the three SN classes, again using the classification probabilities derived from the AT 2016jka host galaxy stellar-population properties (Supplementary Table 5). Marginalizing the posterior probability distributions over all free parameters, we find a 94% probability that AT 2016jka is of type Ia (Supplementary Table 5).

The combination of evidence from the derived host galaxy properties and SN photometry supports the conclusion that AT 2016jka is a type Ia SN with >90%
confidence. Although we have adopted some lens-modelling corrections for all of these methods, this conclusion is not sensitive to the choices we can reasonably make for modelling the lens. As shown in Extended Data Fig. 4, every LENSTOOL model we have evaluated locates the AT 2016jka demagnified position within the region of colour–magnitude space dominated by type Ia SNe. Dust is also not a confounding factor here. The simulations used in both SN-based classification methods include dust extinction at the source plane. If an appreciable screen of dust exists in the lens plane, this would have the effect of making the SN appear dimmer and more red, so correcting for this extinction would move the AT 2016jka points upward and leftward on Extended Data Fig. 2, which would not shift it into the regions occupied by type II and Ib/c SNe.

**Time-delay estimation.** We constrain the relative time delays of AT 2016jka by using two separate methods to estimate the age of the SN at each image during the single observed epoch. The preferred method of SN time-delay measurements involves measuring the time of peak brightness for the SN at each image by fitting the light curves, and taking the difference between each measurement as the relative time delay<sup>24,40,61</sup>. With only a single observed epoch, this method is impossible due to model parameter degeneracies, and we must rely on colour and brightness to constrain the age of each image of the SN. Such age estimates are sometimes referenced to the time of explosion, but in this case we use the observer's convention, setting age = 0 as the time of peak brightness in the rest-frame B band (λ ~ 4,500 Å). Each of these images stems from the same SN explosion, so the difference between the measured ages of the images is also a measure of the relative time delay.

In all of the light curve fitting exercises described below, we also fix $x_1 = 0$. This parameter defines the shape of the SN Ia light curve (the rate of decline in brightness). If we allow $x_1$ to be a free parameter, there is no useful constraint on it. It is highly degenerate with the time delay between the images, which is of course a free parameter in all the fits. As a check, we have also tried fixing $x_1$ to other values from -1 to +1, and the time-delay results change by less than 5 d, which is well within all error bars. Fixing $x_1$ in this way is comparable to the analysis that will be possible in the 2030s when the fourth SN image is observed. We expect that a SALT2 fit to a well sampled light curve from the fourth image will provide a tight constraint on $x_1$. This measured $x_1$ will then be propagated back as a fixed parameter (with small uncertainty) into revised fitting of SN images 1-3. Because of this, we do not incorporate the ±5 d systematic uncertainty in the time-delay errors reported in Table 1.

**Colour curve age constraints.** We first attempt to constrain the relative time delays using the colour of each observed image, which is independent of the lens model and possible because the phenomenon of gravitational lensing is intrinsically achromatic. One important caveat to this principle is that microlensing effects are not generally achromatic, because the microlensing caustics may cause differential magnification on the scale of the SN radius<sup>62–64</sup>. Hence, if the expanding SN shell has a colour gradient then microlensing may introduce spurious features in the observed colours of the SN<sup>65,66</sup>. Reference 62 found that such chromatic microlensing is most likely not present for lensed type Ia SNe in the period up to about 25 rest-frame days after explosion (~15 observer days after peak brightness for AT 2016jka). Only image 2 is probably in the achromatic microlensing phase, but ref. 62 found extremely small deviations in the rest-frame U–V colour curve due to microlensing at the 68% confidence interval, and up to a ~0.2–0.4 mag difference with 99% confidence. While such extreme microlensing could alter the results for images 1 and 3, it would not alter the measurement of image 2, as it is probably in the achromatic phase. Fortuitously, images 1 and 3 are minima, which are less susceptible to high deviations owing to microlensing when compared with image 2, which is a saddle.

We use version 2 of the SuperNova Time Delays (SNTD) package (publicly available at github.com/jpieler14/sntd) with documentation at sntd.readthedocs.io, which has several improvements over the original SNTD package<sup>60</sup>. The SNTD package employs a nested sampling algorithm within three separate methods to measure time delays, and is designed to fully utilize the information present in SN light curve templates<sup>49,50,67,68</sup> to reduce the impacts of microlensing and make more accurate measurements. We use the 'colour' method present in SNTD, which attempts to reconstruct the intrinsic colour curve using the SALT2 model as a template<sup>10</sup>. This method fits the age of each image simultaneously, while also varying the SN model parameters. This means we are finding a single set of SN model parameters to describe the intrinsic photometric evolution of the SN, and also finding the age (time from peak brightness) for each of the three images.

The colour curve constraints resulting from this process are shown in Extended Data Fig. 4. Joint and marginalised posterior distributions from SNTD for the SALT2 SN model parameters and the measured ages for each image are shown in Extended Data Fig. 5. In Extended Data Fig. 4, we can see that the measured colours intersect the model at two distinct locations for images 2 and 3 of AT 2016jka—meaning that there are two plausible ages. This results in a double-peaked posterior distribution in Extended Data Fig. 5. This is caused by a model parameter degeneracy that could be broken in a way independent of the lens model if a sufficiently precise colour curve of image 4 is obtained in the future.

**Light curve + lens model age constraints.** To break the age degeneracies in the colour-based constraints using only data available today, we need to use some information about the relative brightnesses of the AT 2016jka images. For this step we can no longer be independent of the lens models, as we must use the lens-model-predicted magnification values to demagnify the observed photometry for comparison with SN models (note, however, that we again do not use any time-delay information from the lens model for this method).

For the five lens models (A–E) described above, we correct the observed flux density of each image (in both F105W and F160W bands) using the predicted lensing magnification ($\mu$). Next we employ SNTD's 'series' method, as it is most effective for sparse sampling, to attempt a reconstruction of the intrinsic SN light curve<sup>60</sup>. Once again the ages of the images are constrained simultaneously, while also varying the SN model parameters. At this stage we adopt weak priors on the intrinsic type Ia SN luminosity<sup>69</sup> and type Ia SN colour<sup>70</sup> to help break degeneracies in the light curve model. Extended Data Fig. 6 shows the resulting light-curve-based constraints on the age of each SN image.

As our final method to incorporate the colour and brightness information together, we use the colour-based posterior probability distributions (Extended Data Fig. 4) as priors for the light-curve-based constraints. In this approach, we must use only a single photometric band for the light curve constraint, so that we are not 'double-counting' the colour information by simultaneously fitting to two bands together. We adopt F160W as the single band, since it is close to the rest-frame V band, where the SALT2 type Ia SN model is very well constrained. The joint posterior distributions from this method are shown in Extended Data Fig. 7.

## Data availability

All HST images used in this work are available from the Mikulski Archive for Space Telescopes (mast.stsci.edu). HST data from 2016 when the transient was active are from the programme HST-GO-14496 (archive.stsci.edu/proposal_search.php?id=14496&mission=hist). Data collected in 2019 are from the REQUIEM programme, HST-GO-15663 (archive.stsci.edu/proposal_search.php?id=15663&mission=hist). All VLT/MUSE spectroscopic data used in this work are available from the ESO Archive Science Portal (archive.eso.org/dataset/ADP.2019-10-07T18:14:24.762, archive.eso.org/dataset/ADP.2019-10-07T18:14:24.751 and archive.eso.org/dataset/ADP.2019-10-07T18:14:24.776). All derived data supporting the findings of this study (photometry, lens model inputs and so on) are available within this Letter and Supplementary Information. Source data are provided with this paper.

## Code availability

All software tools used in the analysis are publicly available, as indicated in the text. The software used for figure creation, including input data files, can be downloaded from github.com/gbrammer/mrg0138_supernova.

Received: 9 July 2020; Accepted: 14 July 2021;

Published online: 13 September 2021

## References

1. Einstein, A. Lens-like action of a star by the deviation of light in the gravitational field. *Science* **84**, 506–507 (1936).

2. Refsdal, S. On the possibility of determining Hubble's parameter and the masses of galaxies from the gravitational lens effect. *Mon. Not. R. Astron. Soc.* **128**, 307–310 (1964).

3. Holz, D. E. Seeing double: strong gravitational lensing of high-redshift supernovae. *Astrophys. J.* **556**, L71–L74 (2001).

4. Treu, T. & Marshall, P. J. Time delay cosmography. *Astron. Astrophys. Rev.* **24**, 11 (2016).

5. Brammer, G. et al. *Discovery of Probable Multiply Imaged Supernova in Archival HST Data* Transient Discovery Report 113893 (Transient Name Server, 2021).

6. Akhshik, M. et al. Recent star formation in a massive slowly quenched lensed quiescent galaxy at $z=1.88$. *Astrophys. J. Lett.* **907**, L8 (2021).

7. Newman, A. B., Belli, S., Ellis, R. S. & Patel, S. G. Resolving quiescent galaxies at $z>2$. II. Direct measures of rotational support. *Astrophys. J.* **862**, 126 (2018).

8. Ebeling, H., Edge, A. C. & Henry, J. P. MACS: a quest for the most massive galaxy clusters in the universe. *Astrophys. J.* **553**, 668–676 (2001).

9. Jullo, E. et al. A Bayesian approach to strong lensing modelling of galaxy clusters. *New J. Phys.* **9**, 447–447 (2007).

10. Kneib, J.-P. et al. LENSTOOL: a gravitational lensing software for modeling mass distribution of galaxies and clusters (strong and weak regime). *Astrophysics Source Code Library* ascl:1102.004 (2011).

11. Kasen, D. & Woosley, S. E. On the origin of the type Ia supernova width–luminosity relation. *Astrophys. J.* **656**, 661–665 (2007).

12. Riess, A. G. et al. Observational evidence from supernovae for an accelerating universe and a cosmological constant. *Astron. J.* **116**, 1009–1038 (1998).
13. Perlmutter, S. et al. Measurements of $\omega$ and $\lambda$ from 42 high-redshift supernovae. *Astrophys. J.* **517**, 565–586 (1999).

14. Riess, A. G., Casertano, S., Yuan, W., Macri, L. M. & Scolnic, D. Large Magellanic Cloud Cepheid standards provide a 1% foundation for the determination of the Hubble constant and stronger evidence for physics beyond ACDM. *Astrophys. J.* **876**, 85 (2019).

15. Planck Collaboration et al. Planck 2018 results. VI. Cosmological parameters. *Astron. Astrophys.* **641**, A6 (2020).

16. Verde, L., Treu, T. & Riess, A. G. Tensions between the early and late Universe. *Nat. Astron.* **3**, 891–895 (2019).

17. Birrer, S. et al. TDCOSMO. IV. Hierarchical time-delay cosmography—joint inference of the Hubble constant and galaxy density profiles. *Astron. Astrophys.* **643**, A165 (2020).

18. Birrer, S. & Treu, T. TDCOSMO. V. Strategies for precise and accurate measurements of the Hubble constant with strong lensing. *Astron. Astrophys.* **649**, A61 (2021).

19. Spergel, D. et al. Wide-Field Infrared Survey Telescope–Astrophysics Focused Telescope Assets WFIRST-AFTA 2015 report. Preprint at https://arxiv.org/abs/1503.03757 (2015).

20. Ivezic, Z. et al. LSST: from science drivers to reference design and anticipated data products. *Astrophys. J.* **873**, 111 (2019).

21. Kelly, P. L. et al. Multiple images of a highly magnified supernova formed by an early-type cluster galaxy lens. *Science* **347**, 1123–1126 (2015).

22. Kelly, P. L. et al. Déjà vu all over again: the reappearance of supernova Refsdal. *Astrophys. J.* **819**, L8 (2016).

23. Goobar, A. et al. iPTF16geu: a multiply imaged, gravitationally lensed type Ia supernova. *Science* **356**, 291–295 (2017).

24. Dhawan, S. et al. Magnification, dust and time-delay constraints from the first resolved strongly lensed type Ia supernova iPTF16geu. *Mon. Not. R. Astron. Soc.* **491**, 2639–2654 (2020).

25. Quimby, R. M. et al. Detection of the gravitational lens magnifying a type Ia supernova. *Science* **344**, 396–399 (2014).

26. Goldstein, D. A., Nugent, P. E. & Goobar, A. Rates and properties of supernovae strongly gravitationally lensed by elliptical galaxies in time-domain imaging surveys. *Astrophys. J. Suppl. Ser.* **243**, 6 (2019).

27. Pierel, J. D. R. et al. Projected cosmological constraints from strongly lensed supernovae with the Roman Space Telescope. *Astrophys. J.* **908**, 190–207 (2021).

28. Rodney, S. A. et al. SN Refsdal: photometry and time delay measurements of the first Einstein cross supernova. *Astrophys. J.* **820**, 50 (2016).

29. Millon, M. et al. COSMOGRAIL. XIX. Time delays in 18 strongly lensed quasars from 15 years of optical monitoring. *Astron. Astrophys.* **640**, A105 (2020).

30. Millon, M. et al. TDCOSMO. II. Six new time delays in lensed quasars from high-cadence monitoring at the MPIA 2.2m telescope. *Astron. Astrophys.* **642**, A193 (2020).

31. Gonzaga, S. et al. *The DrizzlePac Handbook* (Space Telescope Science Institute, 2012).

32. Astropy Collaboration et al. The Astropy Project: building an open-science project and status of the v2.0 core package. *Astron. J.* **156**, 123 (2018).

33. Weilbacher, P. M. et al. The data processing pipeline for the MUSE instrument. *Astron. Astrophys.* **641**, A28 (2020).

34. Bacon, R. et al. The MUSE Hubble Ultra Deep Field Survey: I. Survey description, data reduction, and source detection. *Astron. Astrophys.* **608**, A1 (2017).

35. Richard, J. et al. An atlas of MUSE observations towards twelve massive lensing clusters. *Astron. Astrophys.* **646**, A83 (2021).

36. Soto, K. T., Lilly, S. J., Bacon, R., Richard, J. & Conseil, S. ZAP—enhanced PCA sky subtraction for integral field spectroscopy. *Mon. Not. R. Astron. Soc.* **458**, 3210–3220 (2016).

37. Iqueras, L. et al. MPDAF—a Python package for the analysis of VLT/MUSE data. In *Astronomical Data Analysis Software and Systems XXVI: Astronomical Society of the Pacific Conference Series* Vol. 521 (eds Molinaro, M. et al.) 545 (Astronomical Society of the Pacific, 2019).

38. Mahler, G. et al. Strong-lensing analysis of A2744 with MUSE and Hubble Frontier Fields images. *Mon. Not. R. Astron. Soc.* **473**, 663–692 (2018).

39. Lagattuta, D. J. et al. Probing 3D structure with a large MUSE mosaic: extending the mass model of Frontier Field Abell 370. *Mon. Not. R. Astron. Soc.* **485**, 3738–3760 (2019).

40. Richard, J. et al. Mass and magnification maps for the Hubble Space Telescope Frontier Fields clusters: implications for high-redshift studies. *Mon. Not. R. Astron. Soc.* **444**, 268–289 (2014).

41. Richard, J. et al. LoCuSS: first results from strong-lensing analysis of 20 massive galaxy clusters at $z=0.2$. *Mon. Not. R. Astron. Soc.* **404**, 325–349 (2010).

42. Mannucci, F. et al. The supernova rate per unit mass. *Astron. Astrophys.* **433**, 807–814 (2005).

43. Newman, A. B., Belli, S., Ellis, R. S. & Patel, S. G. Resolving quiescent galaxies at $z>2$. I. Search for gravitationally lensed sources and characterization of their structure, stellar populations, and line emission. *Astrophys. J.* **862**, 125 (2018).

44. Foley, R. J. & Mandel, K. Classifying supernovae using only galaxy data. *Astrophys. J.* **778**, 167 (2013).

45. Smartt, S. J. Progenitors of core-collapse supernovae. *Annu. Rev. Astron. Astrophys.* **47**, 63–106 (2009).

46. Li, X., Hjorth, J. & Richard, J. The rates and time-delay distribution of multiply imaged supernovae behind lensing clusters. *J. Cosmol. Astropart. Phys.* **2012**, 015 (2012).

47. Barbary, K. et al. sncsmo/sncsmo: v1.4.0. Zenodo https://doi.org/10.5281/zenodo.168220 (2016).

48. Rodney, S. A. et al. Type Ia supernova rate measurements to redshift 2.5 from CANDELS: searching for prompt explosions in the early universe. *Astron. J.* **148**, 13 (2014).

49. Guy, J. et al. SALT2: using distant supernovae to improve the use of type Ia supernovae as distance indicators. *Astron. Astrophys.* **466**, 11–21 (2007).

50. Pierel, J. D. R. et al. Extending supernova spectral templates for next-generation space telescope observations. *Publ. Astron. Soc. Pac.* **130**, 114504 (2018).

51. Kessler, R. et al. SNANA: a public software package for supernova analysis. *Publ. Astron. Soc. Pac.* **121**, 1028–1035 (2009).

52. Frieman, J. A. et al. The Sloan Digital Sky Survey—II Supernova Survey: technical summary. *Astron. J.* **135**, 338–347 (2008).

53. Sako, M. et al. The Sloan Digital Sky Survey—II Supernova Survey: search algorithm and follow-up observations. *Astron. J.* **135**, 348–373 (2008).

54. D'Andrea, C. B. et al. Type II-P supernovae from the SDSS-II Supernova Survey and the Standardized Candle Method. *Astrophys. J.* **708**, 661–674 (2010).

55. Astier, P. et al. The Supernova Legacy Survey: measurement of $\Omega_{\mathrm{XP}}$, $\Omega_{\mathrm{A}}$ and $w$ from the first year data set. *Astron. Astrophys.* **447**, 31–48 (2006).

56. Hamuy, M. et al. The Carnegie Supernova Project: the low-redshift survey. *Publ. Astron. Soc. Pac.* **118**, 2–20 (2006).

57. Stritzinger, M. et al. The He-rich core-collapse supernova 2007Y: observations from x-ray to radio wavelengths. *Astrophys. J.* **696**, 713–728 (2009).

58. Morrell, N. I. Carnegie Supernova Project: spectroscopic observations of core collapse supernovae. In *Death of Massive Stars and Gamma-Ray Bursts: IAU Symposium* Vol. 279 (eds Roming, P. et al.) 361–362 (IAU, 2012).

59. Skilling, J. Nested sampling. *Am. Inst. Phys. Conf. Ser.* **735**, 395–405 (2004).

60. Pierel, J. D. R. & Roddum, S. Turning gravitationally lensed supernovae into cosmological probes. *Astrophys. J.* **876**, 107 (2019).

61. Huber, S. et al. Strongly lensed SNe Ia in the era of LSST: observing cadence for lens discoveries and time-delay measurements. *Astron. Astrophys.* **631**, A161 (2019).

62. Goldstein, D. A., Nugent, P. E., Kasen, D. N. & Collett, T. E. Precise time delays from strongly gravitationally lensed type Ia supernovae with chromatically microlensed images. *Astrophys. J.* **855**, 22 (2018).

63. Foxley-Marrable, M., Collett, T. E., Vernardos, G., Goldstein, D. A. & Bacon, D. The impact of microlensing on the standardization of strongly lensed type Ia supernovae. *Mon. Not. R. Astron. Soc.* **478**, 5081–5090 (2018).

64. Bonvin, V. et al. Impact of the 3D source geometry on time-delay measurements of lensed type-Ia supernovae. *Astron. Astrophys.* **621**, A55 (2019).

65. Kochanek, C. S. Quantitative interpretation of quasar microlensing light curves. *Astrophys. J.* **605**, 58–77 (2004).

66. Vernardos, G. A joint microlensing analysis of lensing mass and accretion disc models. *Mon. Not. R. Astron. Soc.* **480**, 4675–4683 (2018).

67. Hsiao, E. Y. et al. K-corrections and spectral templates of type Ia supernovae. *Astrophys. J.* **663**, 1187–1200 (2007).

68. Kessler, R. et al. Results from the Supernova Photometric Classification Challenge. *Publ. Astron. Soc. Pac.* **122**, 1415–1431 (2010).

69. Wang, X., Wang, L., Pain, R., Zhou, X. & Li, Z. Determination of the Hubble constant, the intrinsic scatter of luminosities of type Ia supernovae, and evidence for nonstandard dust in other galaxies. *Astrophys. J.* **645**, 488–505 (2006).

70. Mosher, J. et al. Cosmological parameter uncertainties from SALT-II type Ia supernova light curve models. *Astrophys. J.* **793**, 16 (2014).

Acknowledgements

We thank P. Kelly and L. Moustakas for helpful commentary on earlier drafts of this work. Data are based on observations made with the NASA/ESA HST, obtained from the data archive at the Space Telescope Science Institute, and on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 0103.A-0777(A). Support for this work was provided by NASA through grant numbers HST-GO-14622 (K.E.W.), HST-AR-15050 (J.D.R.P.), HST-GO-15663 (M.A.) and HST-GO-16264 (S.A.R.) from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. The Cosmic Dawn Center of Excellence is funded by the Danish National Research Foundation under grant no. 140. Support was provided by NASA Headquarters under the NASA Future Investigators in Earth and Space Science and Technology (FINESST) awards 80NSSC19K1418 (M.A.) and 80NSSC19K1418 (J.D.R.P.). K.E.W. acknowledges funding from the Alfred P. Sloan Foundation.
**Author contributions**

Conceptualization, S.A.R., G.B.B. and S.T.; methodology, S.A.R., J.D.R.P., J.R. and K.F.O.; investigation, G.B.B., J.R., S.T., M.A. and K.E.W.; writing—original draft, G.B.B. and S.T.; writing—review and editing, S.A.R., G.B.B., J.D.R.P., J.R., S.T., K.F.O., M.A. and K.E.W.; visualization, S.A.R., G.B.B., J.D.R.P., J.R. and K.F.O.; supervision, S.A.R., G.B.B., S.T. and K.E.W.; funding acquisition, M.A., K.E.W., J.D.R.P. and S.A.R.

**Competing interests**

The authors declare no competing interests.

**Additional information**

*Extended data* is available for this paper at https://doi.org/10.1038/s41550-021-01450-9.

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41550-021-01450-9.

**Correspondence and requests for materials** should be addressed to Steven A. Rodney or Gabriel B. Brammer.

**Peer review information** *Nature Astronomy* thanks Ariel Goobar, Anupreeta More and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2021
<figure><img src="image_7.png" /></figure>

**Extended Data Fig. 1** | Elements of the MRG0138 cluster lens model. The model comprises 37 potentials in total: the BCG (red), 32 cluster members (yellow), three perturbers (cyan), and the main cluster potential (pink). Labeled × symbols indicate the positions of the SN, host, and one additional multiply imaged galaxy with a secure redshift used as model constraints (Supplemental Table 3). The filters used to generate the color image are as in Fig. 1, and tick marks are separated by 10 arcsec.
<figure><img src="image_8.png" /></figure>

**Extended Data Fig. 2** | The position of AT 2016jka in color-magnitude space. Colored points show simulated photometry for normal SNe of Type Ia (red), Type Ib/Ic (gold), and Type II (green), with 10,000 simulated SNe in each subclass (not all apparent on this plot). Histograms above and below show the marginalized distributions that have been rescaled to represent posterior probability density functions. They are normalized to integrate to unity, then multiplied by the SN subclass priors based on the host galaxy stellar population (row b in Supplementary Table 2). Open markers show the observed photometry of the SN. Dotted vertical lines mark the magnification correction based on the preferred LENSTOOL model (model E, described in Methods: Lens Modeling). Closed markers show the resulting magnification-corrected photometry, with asymmetric error bars reflecting the systematic uncertainty derived from the five lens model variants. Horizontal error bars in the upper panel indicate the observed uncertainty in the SN color (not affected by lensing). The relevant SN photometry markers are repeated in the histogram side-panels with arbitrary vertical positions. All three SN images are located in regions of color-magnitude space that are expected to be dominated by Type Ia SNe.
<figure><img src="image_9.png" /></figure>

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

**Extended Data Fig. 3** | A representative set of light curve and color curve models from the STARDUST2 classification algorithm. Panels a-f show F105W and F160W, as indicated, plotting the model light curves in black and photometry as red markers. Panels g-i show the F105W-F160W color curves and color data. All data points as shown have been corrected for magnification and shifted in time using the preferred LENSTOOL model (model E, described in Methods: Lens Modeling). Plotted error bars include the measurement uncertainty and the lens modeling magnification uncertainty. Data points in the right column also include this magnification uncertainty, even though cluster-scale lensing is achromatic, because the STARDUST2 analysis was done on the light curve data, not the color data directly. In all panels the first data point is SN image 2, followed by image 1 and image 3. In each panel the black curves show 200 SN light curve models drawn at random from the nested sampling sequence of the STARDUST2 (sncosmo) classification.
<figure><img src="image_11.png" /></figure>

**Extended Data Fig. 4 | Color-based age constraints for AT 2016jka.** Constraints are shown separately for Image SN1 (panels a and b), image 2 (c and d), and image 3 (e and f), using the methodology described in Methods: Color Curve Age Constraints. Large lower panels (b, d and f) show the observations and model fits. Each magenta shaded region shows the 1σ range of the measured F105W-F160W color, which corresponds to a U-V color in the rest-frame. The model fits are shown as grey shaded regions, indicating the 68% confidence interval of the best-fit SALT2 color curve, with the median model shown as a solid line. The small upper panels (a, c and e) show the posterior for the age of each image from SNTD, using a prior on the SALT2 color parameter (c) based on known population characteristics of SNIa. The effect of adding this prior is slight, with no significant deviation from the best-fit value of $c = 0.02^{+0.04}_{-0.05}$.
<figure><img src="image_12.png" /></figure>

**Extended Data Fig. 5** | Marginalized and joint posterior distributions for the color curve age constraints measured in this analysis. Two-dimensional plots show MCMC sampling points as discrete dots, with contours for the high density regions, drawn at 0.5, 1, 1.5 and 2σ levels. Marginalized (1D) distributions are shown at the top of each column, with dashed vertical lines at the mean and ±1σ (marking 16%, 50% and 84% levels). We use a weak prior on the SALT2 color parameter (c), and set the SALT2 stretch parameter (x₁) to 0. This method is fully independent of lens modeling. The table in the upper right lists all priors, observations, and lens model information used for SN age estimates in this work. Only the highlighted components were used for the constraints shown here.
<figure><img src="image_13.png" /></figure>

**Extended Data Fig. 6** | **Light-curve-based age constraints for AT 2016jka.** Constraints are shown separately for Image SN1 (panels a and b), image 2 (c and d), and image 3 (e and f), using the methodology described in Methods: Light Curve Age Constraints. Large lower panels (b, d and f) show the observations and model fits. Each orange shaded region shows the $1\sigma$ range of the measured F160W magnitude after lens model correction (see Table 1), which corresponds to roughly V band in the rest-frame. The model fits are shown as grey shaded regions, indicating the 68% confidence interval of the best-fit SALT2 light curve, with the median model shown as a solid line. Small upper panels (a, c and e) show the posterior distributions from SNTD for the age of each image that is independent of the lens model (magenta, same as Extended Data Figure 7), using the preferred lens model E (light blue), and the combination of both methods (orange).
<figure><img src="image_14.png" /></figure>

**Extended Data Fig. 7** | Marginalized and joint posterior distributions for the final age constraints measured in this analysis. Two-dimensional plots show MCMC sampling points as discrete dots, with contours for the high density regions, drawn at 0.5, 1, 1.5 and 2σ levels. Marginalized (1D) distributions are shown at the top of each column, with dashed vertical lines at the mean and ±1σ (marking 16%, 50% and 84% levels). We use the color curve posterior as the prior for light curve fitting with lens model E, and include weak priors on the absolute magnitude of a SNla (M<sub>B</sub>) and the SALT2 color parameter (c), and set the SALT2 stretch parameter (x<sub>1</sub>) to 0.