Uncovering a population of gravitational lens galaxies with magnified standard candle SN Zwicky

Detecting gravitationally lensed supernovae is among the biggest challenges in astronomy. It involves a combination of two very rare phenomena: catching the transient signal of a stellar explosion in a distant galaxy and observing it through a nearly perfectly aligned foreground galaxy that deflects light towards the observer. Here we describe how high-cadence optical observations with the Zwicky Transient Facility, with its unparalleled large field of view, led to the detection of a multiply imaged type Ia supernova, SN Zwicky, also known as SN 2022qmx. Magnified nearly 25-fold, the system was found thanks to the standard candle nature of type Ia supernovae. High-spatial-resolution imaging with the Keck telescope resolved four images of the supernova with very small angular separation, corresponding to an Einstein radius of only θE = 0.167″ and almost identical arrival times. The small θE and faintness of the lensing galaxy are very unusual, highlighting the importance of supernovae to fully characterize the properties of galaxy-scale gravitational lenses, including the impact of galaxy substructures.

and faintness of the lensing galaxy is very unusual, highlighting the importance of supernovae to fully characterise the properties of galaxy-scale gravitational lenses, including the impact of galaxy substructures.

Main
Our understanding of gravitational lensing due to the curvature of space-time, and the analogy with the deflection of light in optics, dates back to the work by Einstein in 1936 [1].In this pioneering work he considered the case where both the lens and the magnified background source were stars in our Galaxy.Einstein concluded that the deflection angles were too small to be resolved with astronomical instruments.It was Zwicky [2] who one year later pointed out that if the source was extragalactic, entire galaxies or clusters of galaxies could be considered as gravitational deflectors.Hence, the image separation between multiple images of the source could be large enough to be resolved by astronomical facilities, as the size of the image separation scales with the lens mass and distance as the Einstein radius, , where M ⊙ is the mass of the Sun, M l and D l are the lensing mass and lens angular size distance, and D s and D ls are the distances from the observer to the source and between the lens and the source, respectively.

Strongly lensed supernovae in the era of wide-field time-domain surveys
Besides the many observations of lensed galaxies and quasars, the feasibility to observe strong gravitational lensing of explosive transients in the distant universe has only been demonstrated in recent years, see [3][4][5] and references therein.{PS1-10afx was the first highly magnified SN Ia discovered.However, the lensing interpretation was made three years after the explosion [6,7], by then the SN it was too faint to resolve multiple images.Since then, the use of wide-field optical cameras in robotic telescopes at the Palomar observatory has led to notable breakthroughs.In [8] we reported the first discovery of a multiply-imaged SN Ia, iPTF16geu (SN 2016geu), by the intermediate Palomar Transient Factory (iPTF), a time-domain survey using a 7.3 sq.deg camera on the P48 (1.2-meter) telescope from 2013-2017.In 2018, a new camera was installed [9] with a field of view of 47 sq.deg.The project, known as the Zwicky Transient Facility (ZTF) [10,11], has been monitoring the northern sky with a two-to threeday cadence in at least two optical filters for the past four years [12].The very large sky coverage makes ZTF well-suited to search for rare phenomena, such as gravitational lensing of SNe.On the other hand, the distance (redshift) probed by ZTF is limited by the relatively small mirror of the telescope, light pollution, non-optimal atmospheric conditions, and only having three optical filters at the P48 telescope.Furthermore, ZTF typically obtains image quality (angular resolution) of 2 ′′ full-width half maximum (FWHM), and the camera has relatively large 1 ′′ pixels.Hence, in most instances, it is practically impossible to spatially resolve multiple-image systems with ZTF.Instead, the search for lensed sources makes use of the "standard candle" nature of Type Ia supernovae, i.e., they have nearly identical peak luminosity.These explosions are used as accurate distance estimators in cosmology, which led to the discovery of the accelerated expansion of the universe, see [13] and references therein.In addition to an imaging survey telescope, ZTF has access to a low spectral-resolution integral field spectrograph, SEDM [14], on the neighboring 1.5-meter telescope at Palomar (P60), used to spectroscopically classify about ten SNe every night as part of the Bright Transient Survey (BTS), where transients brighter than 19 mag are classified within timescales of a few days, aiming to obtain > 95% spectroscopic completeness to 18.5 mag or brighter [15].
Besides providing the classification of the transients, the SEDM spectrum is used to measure the SN redshift.

The discovery of SN Zwicky
Lensed system candidates are selected by ZTF for further spectroscopic screening when a SN Ia is found at a redshift above z = 0.2, where there is essentially negligible sensitivity for detection in BTS, unless the SN is greatly magnified by an intervening deflector.This was the case for "SN Zwicky" (a.k.a.ZTF22aaylnhq and SN 2022qmx), located at right ascension 17 h 35 m 44.32 s and declination 4 • 49 ′ 56.90 ′′ (J2000), where an SEDM spectrum from 2022 August 21 showed it to be a Type Ia SN at z = 0.35, as shown in Fig. 1.At that point we alerted the SN community of the discovery of a lensed SN Ia [16].
Spectroscopic observations following the time-evolution of the SN were carried out using multiple facilities, the 2.56-meter Nordic Optical Telescope at the Canary Islands, the Keck observatory in Hawaii, the 11-meter Hobby-Eberly Telescope at the McDonald Observatory in Texas and at ESO's 8-meter Very Large Telescopes (VLT) at the Paranal Observatory in Chile.In particular, multiple narrow emission and absorption lines of the SN host galaxy were found with LRIS/Keck and MUSE/VLT, refining the source redshift to z = 0.3544, as shown in the bottom panels of Fig. 1.As the SN faded and the foreground galaxy spectral energy distribution (SED) became more prominent, the Ca II doublet λλ3933, 3968 was found in absorption lines redshifted to z = 0.22615, the location of the deflecting galaxy.

Follow-up observations
The discovery from ZTF was also followed-up with high-spatial resolution instruments.Observations with laser guide star enhanced seeing at the Very Large Telescope with the HAWKI imaging camera in the near-IR, and optical spectrophotometry with MUSE, reduced the point spread function (PSF) width to about 0.4 ′′ .However, this was still not enough to resolve the system.On 2022 September 15, multiple images of the system were first resolved at near-IR wavelengths at the Keck telescope, using the Laser Guide Star aided Adaptive Optics (LGSAO) with the Near-IR Camera 2 (NIRC2) [17], yielding an image quality of 0.09 ′′ FWHM in the J-band centered at 1.2 µm, shown in Fig. 3, where the four SN images are labeled A-D.
On September 21, following our announcement of the discovery [16], a previously approved program aimed to target lensed SNe by the LensWatch collaboration resolved the multiple images of SN Zwicky using the optical filters F 475W , F 625W and F 814W (where the names correspond to the approximate location of the central wavelength in nanometers) on the UVIS/WFC3 Camera on the Hubble Space Telescope (HST) [18].A detailed description of the HST observations of SN Zwicky is presented in [19].

Results
Fig. 2 shows the unresolved photometric ground-based observations collected at P48 and the Liverpool telescope on the Canary Islands between 2022 August 1 and October 30.These were used to estimate the peak flux and lightcurve properties of the SN with the SALT2 lightcurve fitting tool [20], including corrections for lightcurve shape and colour excess given  The solid lines show the best fit SALT2 [20] model to the spatially unresolved data.The blue dashed lines indicate the expected lightcurves at z = 0.354 (without lensing) where the bands represent the standard deviation of the brightness distribution for SNe Ia.In order to fit the observed lightcurves, a brightness increase corresponding to 3.5 magnitudes is required.Also shown as dotted lines are the modeled individual lightcurves for the four SN images A-D.The flux ratios were measured from the NIRC2 data in Fig. 3. From these lightcurves we extract the time-delays between the images, all in units of days, ∆t AB = −0.4± 2.9, ∆t AC = −0.1 ± 2.3, and ∆t AD = −0.1 ± 2.7, as described in the Supplementary Material section.The shaded areas in the lower panels indicate the regions with data outside the phase range where the SALT2 model is defined, therefore, excluded from the lightcurve fit analysis.The error bars correspond to 1 standard deviation.
the SN redshift, as well as the extinction in the Milky-Way in the direction of the SN.Furthermore, the four resolved SN images were used to explore the possibility of additional reddening by dust in the lensing galaxy.Unaccounted dimming of light would lead to an underestimation of the lensing amplification.The HST and NIRC2 observations for each SN image were compared with the SN Ia spectral template from [21] allowing for possible differential dust extinction in the lens following the reddening law in [22].
No evidence for differential extinction between the different images was found.The lightcurve fit model included the four individual SN images, described by the SALT2 model with arbitrary time delays, but otherwise sharing the same lightcurve shape parameters, x 1 and c.The time delays were constrained by a prior on the image flux ratios from the NIRC2 observations shown in Fig. 3.We find ∆t AB = −0.4± 2.9, ∆t AC = −0.1 ± 2.3, and ∆t AD = −0.1 ± 2.7 (in units of days), where the indices A-D refer to the SN images in Fig. 3.The resulting lightcurve fit is shown in Fig. 2 and compared to the SALT2 model.The small time-delays are also consistent with the spectra of SN Zwicky shown in Fig. 1 being at a single SN phase.The fitted SN model parameters, i.e. lightcurve shape and colour from the SALT2 model are x 1 = 1.083 ± 0.094 and c = −0.007± 0.007.The lack of colour excess confirms that differential extinction is negligible.Since the lightcurve parameters errors do not include the model covariance, we conservatively add the SALT2 model error floor of σ(x 1 ) = 0.1 and σ(c) = 0.027 mag [23] in quadrature to the fit errors.Using the inferred apparent magnitude and the SALT2 parameters above, we find a total magnification of ∆m = −3.44 ± 0.14 mag, assuming standard cosmological parameters from [24] and restframe B-band SNe Ia peak absolute magnitude of −19.4 mag for the average SN Ia lightcurve width and colour, and intrinsic brightness scatter of 0.1 mag.Since the inferred stellar mass of the host galaxy is M ⋆ < ∼ 10 10 M ⊙ , mass-step corrections for the SN Ia absolute magnitude as suggested in [25] are not required.In summary, we find that including the four images, SN Zwicky is 23.7 ± 3.2 times brighter than the observed flux of normal SNe Ia at the same redshift, after applying colour and lightcurve shape corrections.
The Keck NIRC2 J-band image was used to obtain a lens model to account for the observed SN image positions, irrespectively of their fluxes.Assuming a singular isothermal ellipsoid [26,27] for the lens potential, we report an ellipticity ϵ e = 0.35 ± 0.01 and Einstein radius θ E = 0.1670 ± 0.0006 ′′ .The mass enclosed within the ellipse with semi-major axis 0.78 kpc and semi-minor axis 0.51 kpc is M = (7.82±0.06)•10 9M ⊙ .The lens model predictions for the time delays are in excellent agreement with the fitted values from the lighturves in Fig. 2. Further details regarding the lens modelling are presented in the Methods Section.
Interestingly, the individual image magnifications predicted for SN Zwicky by the smooth macro lens model are inconsistent with the observed flux ratios.According to the lensing model, the observed fluxes of the SN images A and C are factors of > 4 and > 2 too large, respectively, compared to images B and D. Given the small time-delays, this discrepancy cannot be accounted by different phases between the SN images.Other explanations need to be considered, e.g., excess magnification and demagnification from milli-and microlensing effects arising from stars and substructure in the lens galaxy [28,29].Since microlensing effects are capable of perturbing magnifications without altering image locations, these were also put forward to explain the observed flux ratios of iPTF16geu [30], displaying differences between the observed and modelled image flux ratios of similar magnitude.Probing microlensing in these central regions opens a new window to directly measure the central stellar initial mass function (IMF) [31] and test claims that the IMF may be heavier in the centres of galaxies [32].As detailed in Methods Section, the lack of time dependence in the image flux ratios, or anomalous variation in the unresolved lightcurves, gives a lower limit for the substructure masses of 0.02 M ⊙ , if the discrepancy from the smooth lens model is caused by microlensing.From the lack of further image splitting of the four individual SN images, we infer an upper limit for the substructure mass of 3 • 10 7 M ⊙ .
In order to check the impact of added macro lens model complexity, we have studied cases where the lens mass distribution is modelled with two matter components; one where the surface mass density follows the lens light distribution (with an arbitrary mass to light ratio, possibly interpreted as a baryonic mass component), and a second one introducing a dark matter halo with additional flexibility on density profiles.In spite of the added extra complexity, the quality of the fit to the SN image positions does not improve, and only induces shifts in the predicted flux ratios below 5 %, i.e., very small in comparison with the observed flux ratio anomalies.

SN Zwicky PS1-10afx
Lensed quasars BELLS lens galaxies SLACS lens galaxies SL2S lens galaxies Fig. 4 Stellar mass versus Einstein radius for lens galaxies discovered in galaxy surveys, demonstrating that SN Zwicky, iPTF16geu, and the unresolved lensed supernova PS1-10afx point to a poorly known population of small image separation lensing systems.Strongly lensed galaxy systems are represented by yellow diamonds for the BELLS sample [70], blue triangles for SLACS lenses [71] and purple squares for SL2S lenses [72].The green dots correspond to lensed quasars [73], of which the filled dots have been detected optically and the open dots through radio emission.The shaded grey contours show the 90% and 68% confidence levels for the full sample of 155 lensed galaxies and 45 lensed quasars.The lensed supernova data are presented as median values ± one standard deviation.For the unresolved lensed supernova PS1-10afx, only an upper limit of the Einstein radius is available [7].The stellar mass derivation for SN Zwicky is detailed in the Methods Section.
The demonstrated ability to discover multiply-imaged SNe makes it feasible to accomplish Refsdal's pioneering proposal [33] to use time-delays for strongly lensed SNe to measure the Hubble constant.This will require systems with time-delays of several days, i.e., longer than for SN Zwicky.The small physical scale of the lens probed by SN Zwicky, as well as iPTF16geu [8], make these SNe unique tracers for uncovering a population of systems that otherwise would remain undetected, as shown in Fig. 4, probing the mass distribution of the central, densest, regions of lensing galaxies.Both of these multiply-imaged SNe Ia were identified without preselections on e.g., association with bright galaxies or clusters, emphasizing the importance of untargeted surveys for unexpected discoveries.

Supernova survey and follow-up
The Zwicky Transient Facility at Palomar Observatory is monitoring the transient sky at optical wavelengths since 2018 [9][10][11]34].'SN Zwicky" (a.k.a.SN 2022qmx) [16] was discovered under the Bright Transient Survey (BTS) program [15].The first detection of the SN was in a ZTF g-band image from 2022 August 1.It was saved to the BTS as a SN candidate by on-duty scanners on August 3 [35] and subsequently assigned to the queue for spectroscopic follow-up with the SED Machine (SEDM) mounted on the Palomar 60-inch telescope (P60) under standard BTS protocols.The SEDM spectrum, obtained on August 21 [36], shows an excellent match to a normal Type Ia supernova at a redshift of z = 0.35 close to maximum light.The redshift and spectral classification were confirmed with a higher resolution spectrum obtained at the Nordic Optical Telescope on La Palma on the following night.We followed SN Zwicky up with P48 in the g and r bands.For our analysis we use the forced photometry as provided by the Image Processing and Analysis Center (IPAC) as detailed in [34].Observations in the SDSS g, r, i and z filters were taken with the IO:O optical imager on the Liverpool Telescope (LT) [37].The LT photometric data is processed with custom datareduction and image subtraction software (Taggart et al 2022 in prep.).Image subtraction is performed using PS1 references image.Images were stacked using SWarp to combine multiple exposures where required.The photometry is measured using a point spread function (PSF) fitting methodology relative to PS1 standards and is based on techniques in [38].

Lightcurve fit, magnification and Time Delay Inference
We used the publicly available, python-based software sntd [39] for inferring the restframe B-peak magnitude, the lightcurve shape and colour SN Ia SALT2 parameters [23] and the time-delays between the SN images.Data points with ≥ 3σ detections from the g and r filters from P48 and the g, r, i, z filters from the Liverpool telescope were included in the fit.We adopted an iterative procedure in two steps.First, all the data was used.Next, only data points in the range where the SALT2 model is defined, i.e. −20 to +50 days were kept for the second iteration.The final lightcurve fit parameters were derived from data in this phase range.
We estimated the time delays, i.e., the relative phase between the SN images, by fitting the unresolved ground-based flux data with the model which includes the flux contributions from the four sets SN lightcurves, F j (t, λ), each one with their own time of B-band max, t j 0 .The fit is constrained by imposing a prior on the image ratios at the date of the Keck/NIRC2 observations, reported in Supplementary Table 3.We let the lightcurve fit parameters vary withing the ranges x 1 = [−3.0,3.0], c = [−0.3,0.3], but assume them to be the same for all four images.This is an excellent approximation in the absence of appreciable differential reddening in the lensing galaxy, confirmed by the result of the fit, c = −0.01±0.01,consistent with no colour excess.
Galactic extinction was included in the model, adopting E(B − V ) MW = 0.1558 mag, based on the extinction maps in [40].We used the wavelength dependence from the dust from [22] and the measured mean value of total to the Galactic selection extinction ratio, R V = 3.1.
From the location of the fitted restframe B-band peak luminosity in the SALT2 model for the summed fluxes we inferred the time of maximum for SN Zwicky, t 0 = 59808.67± 0.19, corresponding to August 17, 2022.The total lensing magnification, µ = 24.3 ± 2.7 was calculated after standardising the SN Ia fitted B-band peak magnitude with the standard SALT2 lightcurve shape and colour correction parameters, (α, β) = (0.14, 3.1).Throughout the analysis of SN Zwicky, we adopted a flat ΛCDM cosmological model with H 0 = 67.4km s −1 Mpc −1 and Ω M = 0.315 [24].
The uncertainties we report for the time delays account for parameter degeneracies in the the fit.Corner plots with posteriors for the relative time delays of B, C, and D with respect to A are illustrated in Supplementary Fig. 1.The unresolved data used for this analysis will be available on WiseRep after paper acceptance.It should be emphasised that the lightcurve and time-delay fits were done completely independently from the lens modelling.

Spectroscopic follow-up
The first classification spectrum of SN Zwicky was obtained with the IFU on the SEDM on 2022 August 21.The data were reduced using a custom IFU pipeline developed for the instrument [41,42].Flux calibration and correction of telluric bands were done using a standard star taken at a similar airmass.The details of all spectroscopic observations are listed in Table ??.
We obtained three epochs of spectroscopy between 21 August and 11 September 2022 with the Alhambra Faint Object Spectrograph and Camera (ALFOSC) on the 2.56 m Nordic Optical Telescope (NOT) at the Observatorio del Roque de los Muchachos on La Palma (Spain), (see http://www.not.iac.es/instruments/alfosc for further information).Observations were taken using grism 4, providing wavelength coverage over most of the optical spectral range (typically 3700-9600 Å).Reduction and calibration were performed using PypeIt version 1.8.1 [43,44].
We obtained three epochs of spectroscopy between 22 August and 19 October 2022 with the Low Resolution Imaging Spectrometer (LRIS) on the Keck I 10 m telescope [45].All spectra were reduced and extracted with LPipe [46].

Adaptive Optics observations
Observations with the VLT

HAWKI observations and data reduction:
We obtained 7 epochs of near-IR in Y JHK between 23 August and 30 September 2022 with the High Acuity Wide field K-band Imager (HAWK-I) [47][48][49] at the ESO Very Large Telescope at the Paranal Observatory (Chile).All observations, except of the first epoch, were performed with the ground layer adaptive optics offered by the GALACSI module [50][51][52] to improve the image quality.The first three were observed in Y JH filters and the following 4 in Y JHK s filters.For the first three epochs we exposed for 3 × 60 s in the Y and J bands and 6 × 60 s in H band.For the fourth epoch we also observe for 10 × 60 s in the K s band.To account for the brightness evolution, we exposed for 10 × 100 s and 6 × 60 s for epochs five and six.For the final epoch we also increased H-band exposure times to 10 × 60 s.For the HAWK-I observations we used offsets of (−115 ′′ , 115 ′′ ) to place the target on the optimal detector chip.
The data used in our work has been reduced using the HAWK-I pipeline version 2.4.11 and the Reflex environment [53].The data reduction included subtracting bias and flat fielding.The world coordinate system was calibrated against stars from GAIA.

MUSE observations and data reduction:
We obtained 4 epochs of integral-field spectroscopy between 24 August and 30 September 2022 with the Multi-Unit Spectroscopic Explorer (MUSE) [54] at the ESO Very Large Telescope at the Paranal Observatory (Chile).Each pointing has an approximately 1 ′ × 1 ′ field of view with spatial sampling of 0.2 ′′ /pixel and covers the wavelength range from 4700 to 9300 Å with a spectral resolution of 1800 to 3600.All observations were performed with the ground layer adaptive optics offered by the GALACSI module [50,52] to improve the image quality.The integration time of each epoch was 1800 s.
The data used in our work has been reduced using the MUSE pipeline version 2.8.7 [55] and the Reflex environment [53].The data reduction included subtracting bias, flat fielding, wavelength calibration and flux calibration against spectrophotometric standard stars.Afterwards, we improved the sky-subtraction with the Zurich Atmosphere Purge (ZAP) [56] module in the ESO MUSE pipeline.The world coordinate system was calibrated against stars from GAIA.

Laser Guide Star Adaptive Optics imaging from Keck
The NIRC2 J-band observations consist of 9 images in a 5-point dither pattern based on a 2 ′′ × 2 ′′ grid size, to facilitate sky background subtraction.At the first dither location we obtained 2 images; one co-added 600s exposure and one 200s exposure.At the second location we obtained one 200s exposure.At each of the third, fourth and fifth dither locations, we obtained two 200s exposures.To correct for flat-fielding and bias we acquired a set of ten bias frames (flat lamp off) and ten dome flat (flat lamp on) frames.Sky background and dark current was removed as part of the sky subtraction, which utilized a different sky map for each dither position, created by median combining the frames from all other dither positions, excluding the current dither position.The final combined image was created by aligning each dither position to each other using the centroid of the brightest SN image (Image A), and median combining.The reduction was carried out using custom python scripts.The NIRC2 J-band data (Fig. 2) provides the highest resolution (FWHM 0.086 ′′ ) image in our dataset of the system where the four SN images are visible.

Modelling of the lens galaxy
The Keck NIRC2 J-band image was used to model the lens galaxy in terms of its Einstein radius θ E , semi-minor to semi-major axis ratio q (or ellipticity ϵ e = 1 − q) and orientation angle ϕ.The mass profile used in our analysis is a singular isothermal ellipsoid (SIE) [26,27] : where κ corresponds to the convergence (i.e. the dimensionless projected surface mass density) and the coordinates (x, y) are centred on the position of the lens centre and rotated counterclockwise by ϕ.The projected mass M inside an isodensity contour of the SIE is given by: [27] where D l , D s , and D ls are the angular diameter distances between the observer and the lens, the observer and the source, and the lens and the source, respectively.In order to calculate these distances, we assumed a flat ΛCDM cosmology with H 0 = 67.4km s −1 Mpc −1 and Ω M = 0.315 [24].
In addition to the lens mass model, we included light models for the lens galaxy and SN host galaxy in the form of elliptical Srsic profiles: where I e is the intensity at the half-light radius R e , b n = 1.9992n − 0.3271 [57] and with q S the axis ratio of the Srsic profile.The SN images were modelled as point sources.
We used a Moffat PSF with power index 2.94 and FWHM 0.091 ′′ to model the full image.
We simultaneously reconstructed the lens mass model, SN images, and the surface brightness distributions of the lens and host galaxy.The lens mass model is constrained only by the positions of the lensed SN images and not by their fluxes, since the latter may be considerably affected by substructures, such as stars, in the lensing galaxy.Our model contains 13 nonlinear free parameters: θ E , ϕ, q, x, y for the lens mass model, R e , n, ϕ S , q S , x S , y S for the lens light model and x SN , y SN for the SN position in the source plane.Our results are obtained using LENSTRONOMY (https://lenstronomy.readthedocs.io/en/latest/),an open-source python package that uses forward modeling to reconstruct strong gravitational lenses [57].The result of the fit and comparison with the observations is shown in Supplemental Fig. 2. As a crosscheck, we also modeled the HST photometry data and redid the analysis with LENSMODEL [58,59], which produced consistent results.The resulting best-fit values for the lens mass and light profiles are summarised in Supplemental Table 2. Additionally, we derived the gravitational mass within the isodensity contour given by the critical line of the lens (with radius θ E = 0.167 arcsec, corresponding to 0.628 kpc) to be M = (7.82± 0.06) • 10 9 M ⊙ .
Supplemental Table 3 displays the observed time delays and individual fractional flux contributions from each SN image as detailed in the subsection Lightcurve fit and Time Delay Inference (t obs and f obs ) compared to the predictions from the lens model (t mod and f mod ).Here, time delays are given with respect to image A, e.g.t i ≡ t i −t A .The observed fractional flux ratios are measured from the Keck J-band image after subtracting the lens galaxy light, and the model predictions, f mod , are computed from the magnifications predicted by the lens model, f i ≡ µ i / j µ j .In addition to the uncertainties obtained from the J-band image analysis, we make a conservative error estimate by also including the scatter in f obs and f mod obtained when modelling the system using data from the HST optical filters F 475W , F 625W and F 814W , as well as the Keck near-IR J-band data.Using this approach, we also take into account possible error contributions from uncertainties in dust extinction, lens galaxy subtraction, and lens mass modelling.
The observed individual image magnifications can be obtained from f obs by multiplying the individual fractional fluxes with the total observed SN magnification of µ tot obs = 24.3±2.7.The total lens model image magnification, µ tot mod , is sensitive to the lens mass slope (which is unconstrained by the imaging data), such that flatter halos predicts a larger magnification.For an isothermal lens, µ tot mod = 14.9 ± 0.9, indicating a flatter halo profile (see also [8]).However, the predicted flux ratios remain approximately constant, which means that to first approximation, we can multiply the derived f mod with an arbitrary µ tot mod .In contrast to the predicted individual fractional flux contributions, the observed flux is dominated by images A and C. In order to match the observed values, substructure lensing is needed to additionally magnify images A and C with factors of > 4 and > 2, respectively, compared to any additional substructure (de)magnifications of images B and D.
In order to check the impact of added macro lens model complexity, we have studied cases where the lens mass distribution is modelled with two matter components; one where the surface mass density follows the lens light distribution (with an arbitrary mass to light ratio, possibly interpreted as a baryonic mass component), and a second one introducing a dark matter halo with additional flexibility on density profiles.In spite of the added extra complexity, the quality of the fit to the SN image positions does not improve, and only induces shifts in the predicted flux ratios below 5 %, i.e., very small in comparison with the observed flux ratio anomalies.
We do not detect any time dependence in the flux ratios between the Keck and HST images, observed just one month after lightcurve peak, six days apart, nor any other anomalous variation in the unresolved ∼ 80 days long lightcurve shown in Fig. 2. Comparing the expected size of SN photosphere with the Einstein radii of compact objects in the lensing galaxy, we infer that if the discrepancy from the smooth lens model is caused by microlensing, the deflectors must exceed 0.02 M ⊙ .
With the aim of distinguishing between lensing effects by stars or larger substructures, we investigated the upper limit of image splitting for the brightest image, A. We approximated the maximum Einstein radius of a large structure at the position of image A by putting an upper limit on the difference in the full width at half maximum between PSF A and PSF BCD .Using equation 2, we inferred a 95% confidence upper limit on the substructure's mass within its Einstein radius.Combined with the lower mass limit derived from the lack of flux ratio variability, this constrains the mass of the substructure deflector in the line-of-sight to image A to 0.02 < M/M ⊙ < 3 • 10 7 .
We conclude that, similarly to the case of iPTF16geu, a smooth lens density fails to explain the individual image magnitudes and additional sub-structure lensing is needed.Since the observed properties of lens systems to first order only depend on the integrated mass within the images and/or the surface mass density of the lens at the image positions or in the annulus between the images [60], small image separation systems provide a unique probe of the central regions of gravitational lensing galaxies.The probability for lens substructures, such as stars, to accommodate the needed (de)magnifications for a range of lens density slopes will be investigated in a future lens modeling paper.

Lens galaxy photometry
We retrieved science-ready coadded images from the Panoramic Survey Telescope and Rapid Response System (PS) DR1 [61].Using a circular aperture with a radius of 1.1 ′′ , we obtain the following apparent magnitudes of the lens galaxy: g = 22.09 ± 0.09, r = 20.71± 0.02, i = 20.14 ± 0.02, z = 19.84 ± 0.02, and y = 19.63 ± 0.05 (all errors are of statistical nature).We model the spectral energy distribution with the software package prospector version 1.1 [62].Prospector uses the Flexible Stellar Population Synthesis (FSPS) code [63] to generate the underlying physical model and python-fsps [64] to interface with FSPS in python.The FSPS code also accounts for the contribution from the diffuse gas (e.g., H II regions) based on the Cloudy models from [65].Furthermore, we assumed a Chabrier initial mass function [66] and approximated the star formation history (SFH) by a linearly increasing SFH at early times followed by an exponential decline at late times (functional form t × exp (−t/τ )).The model was attenuated with the [67] model.
The best fit galaxy model points to a moderately massive galaxy with stellar mass log M ⋆ /M ⊙ = 10.1 ± 0.3.The other parameters, such as star-formation rate, attenuation and age are poorly constrained and we report their values only for the sake of completeness: SFR = 15 +21 −15 M ⊙ yr −1 , E(B − V ) star = 0.6 +0.1 −0.4 mag, Age = 1.6 +5.4 −1.2 Gyr.We acknowledge that the aperture might not encircle the entire galaxy and, therefore, we might underestimate the stellar mass of the lens.Changing the radius of the measurement aperture to 4 ′′ increases the mass by ≈ 0.3 dex.This large aperture also includes the contribution of the SN host galaxy and therefore overestimate the stellar mass of the lens galaxy.Nonetheless, this upper bound does not alter our conclusion about the compact nature of the lensing galaxy, relative to other lensing systems.

Host galaxy photometry
Numerous studies have shown that the peak absolute magnitudes of SNe Ia depend on their host galaxy masses, see e.g., [25].Although the host and the lens galaxy are well separated in the HST images, both galaxies have diffuse emission extending well beyond the separation of two galaxies.To first order, we can remove the contribution of the lens by subtracting the lens images predicted by LENSTRONOMY from the HST images.Using a circular aperture with a 1 ′′ radius and appropriate zeropoints from HST, we measure for the host: 22.31±0.11,21.81± 0.06 and 21.13 ± 0.06 mag in F 475W , F 625W and F 814W , respectively.Fitting the SED with prospector with the same assumptions as in the previous Section gives a galaxy stellar mass of log M ⋆ /M ⊙ = 9.6 +0.4 −0.3 .A closer inspection of the subtracted HST images reveals that the 1 ′′ aperture does not encircle the total emission of the host.Increasing the aperture radius to 1.5 ′′ encircles almost the entire host galaxy (F 475W = 21.80 ± 0.09, F 652W = 21.38 ± 0.04, F 814W = 20.74 ± 0.03), but it also includes non-negligible contribution from residuals of the lens galaxy.The galaxy mass increases marginally to log M/M ⊙ = 9.7 +0.5 −0.4 but is consistent with the previous measurement.This suggest that a mass-step correction is not required.To get a more robust estimate of the galaxy mass, NIR observations after the SN has faded are required.IA: data analysis; TC, CL: manuscript contribution; EB, JB,MG, MK, SK, AM, JDN, JN, KN, RR, MR, BR, GS, AT, AW: ZTF observations; YS, RS: Keck spectroscopy; JV, JCW: follow-up observations; JP; HST observations.J.R., Very Large Telescope observations.

Fig. 1
Fig. 1 Spectroscopic identification of SN Zwicky as a Type Ia supernova and redshift measurements of the SN host galaxy and the intervening lensing galaxy.The SN spectra obtained with the P60, NOT and Keck telescopes (black lines) are best fit by a normal SN Ia spectral template.The green line shows a comparison with the nearby Type Ia SN 2012cg [69] at a similar rest-frame phase, redshifted to z = 0.354.The SN flux peaked on August 17, 2022.Bottom panels show a zoomed-in view of a VLT/MUSE spectrum from 2022 September 30, displaying narrow absorption and emission lines, from which precise redshifts of the lens (z = 0.2262) and host (z = 0.3544) galaxies were determined.[O II], Ca II H & K, Na I D, Hα, [N II] and [S II] lines can be seen at the rest-frame of the lens (blue lines) and host (red lines) galaxies.

Fig. 3
Fig. 3 Image of the field of SN Zwicky using pre-explosion g and r−band images from ZTF.The top right panel shows a composite image using VLT MUSE (g/r) and HAWK-I J-band data.The bottom panel shows the Keck/NIRC2 J-band image where the four SN images are resolved.The top left panel shows a 2 ′ ×2 ′ section of the ZTF g and r-band pre-SN images (FWHM 2.3 ′′ ), centered on the location of SN Zwicky.Top right panel shows a zoomed-in composite image of SN Zwicky using AO enhanced VLT MUSE and HAWK-I observations on September 2 and 4 (FWHM ∼ 0.4 ′′ ).Bottom panel shows a 2 ′′ × 2 ′′ portion of the Keck/NIRC2 LGSAO J-band image (FWHM 0.09 ′′ ), resolving the four multiple images of SN Zwicky (labelled A, B, C, D).The blue dashed ellipse shows the critical line of the lens, corresponding to the inferred Einstein radius θ E = 0.167 ′′ (0.6 kpc at z = 0.2262), enclosing the lens M = (7.82± 0.06) • 10 9 M ⊙ .The host galaxy nucleus is located 1.4 ′′ to the northeast of the lens, implying that SN Zwicky exploded at a projected distance of 7 kpc from the centre of its host galaxy.

Fig. 5 Fig. 6 A
Fig. 5 Corner plot showing the posterior distribution for the relative time-delay of the individual images B, C, and D to A,the fitted supernova parameters x 1 and c, as well as the reconstructed flux ratios at peak from the fit with the relevant correlations between each parameter.The 1D plots show the mean value with the 1 sigma uncertainty.The 2D plots show the 1 and 2 sigma contours.
Multi-color lightcurve of SN Zwicky showing that the supernova is 3.5 magnitudes brighter than an unlensed SN at the same redshift.The magnitudes are measured with respect to time of maximum light (Modified Julian Date 59808.67,corresponding to August 17, 2022) in ZTF g and r and Liverpool Telescope griz filters.