LensWatch. I. Resolved HST Observations and Constraints on the Strongly Lensed

Supernovae ( SNe ) that have been multiply imaged by gravitational lensing are rare and powerful probes for cosmology. Each detection is an opportunity to develop the critical tools and methodologies needed as the sample of lensed SNe increases by orders of magnitude with the upcoming Vera C. Rubin Observatory and Nancy Grace Roman Space Telescope. The latest such discovery is of the quadruply imaged Type Ia SN 2022qmx ( aka, “ SN Zwicky ” ) at z = 0.3544. SN Zwicky was discovered by the Zwicky Transient Facility in spatially unresolved data. Here we present follow-up Hubble Space Telescope observations of SN Zwicky, the ﬁ rst from the multicycle “ LensWatch ( www.lenswatch.org ) ” program. We measure photometry for each of the four images of SN Zwicky, which are resolved in three WFC3 / UVIS ﬁ lters ( F475W, F625W, and F814W ) but unresolved with WFC3 / IR F160W, and present an analysis of the lensing system using a variety of independent lens modeling

Facility (ZTF) in spatially unresolved data.Here we present follow-up Hubble Space Telescope observations of SN Zwicky, the first from the multi-cycle "LensWatch" program.We measure photometry for each of the four images of SN Zwicky, which are resolved in three WFC3/UVIS filters (F475W, F625W, F814W) but unresolved with WFC3/IR F160W, and present an analysis of the lensing system using a variety of independent lens modeling methods.We find consistency between lens model predicted time delays (≲ 1 day), and delays estimated with the single epoch of HST colors (≲ 3.5 days), including the uncertainty from chromatic microlensing (∼ 1-1.5 days).Our lens models converge to an Einstein radius of θ E = (0.168 +0.009  −0.005 ) ′′ , the smallest yet seen in a lensed SN system.The "standard candle" nature of SN Zwicky provides magnification estimates independent of the lens modeling that are brighter than predicted by ∼ 1.7 +0.8  −0.6 mag and ∼ 0.9 +0.8 −0.6 mag for two of the four images, suggesting significant microlensing and/or additional substructure beyond the flexibility of our image-position mass models.
1. INTRODUCTION Strong-gravitationally lensed supernovae (SNe) are rare events.In the strong lensing phenomenon, multiple images of a background source appear as light propagating along different paths are focused by a foreground galaxy or galaxy cluster.This requires a chance alignment along the line of sight between us the observers, the background source, and the foreground galaxy.Depending on the relative geometrical and gravitational potential differences of each path, the SN images typically appear delayed by hours to months (for galaxy-scale lenses) or years (for cluster-scale lenses).
Robust measurements of this "time delay" can constrain the Hubble constant (H 0 ) and the dark energy equation of state (e.g., w) in a single step (e.g., Refsdal 1964;Linder 2011;Paraficz & Hjorth 2009;Treu & Marshall 2016;Birrer et al. 2022b;Treu et al. 2022).Lensed SNe have several advantages relative to quasars, which have historically been used for this purpose (e.g., Vuissoz et al. 2008;Suyu et al. 2010;Tewes et al. 2013;Bonvin et al. 2017;Birrer et al. 2019;Bonvin et al. 2018Bonvin et al. , 2019b;;Wong et al. 2020): 1) SNe fade quickly, enabling predictive experiments on the delayed appearance of trailing images more accurate models of the lens and source, as the SN (or quasar) and host fluxes are otherwise highly blended (Ding et al. 2021), 2) SNe have predictable light curves, simplifying time-delay measurements and enabling SN progenitor system constraints, 3) the impact of microlensing is somewhat mitigated including a small (∼ 0.1 day) "microlensing time delay" (Tie & Kochanek 2018;Bonvin et al. 2019a) and less pronounced chromatic effects (Goldstein et al. 2018;Foxley-Marrable et al. 2018;Huber et al. 2019), though this can still be a significant source of uncertainty when time delays are ≲ 1 day (e.g., Goobar et al. 2017), and 4) time delay measurements for lensed SNe require much shorter observing campaigns than lensed quasars.
SNe of Type Ia (SNe Ia), those employed for decades as "standardizable candles" to measure cosmological parameters by way of luminosity distances and the cosmic distance ladder (e.g., Garnavich et al. 1998;Riess et al. 1998;Perlmutter et al. 1999;Scolnic et al. 2018;Brout et al. 2022), are particularly valuable when strongly lensed.In addition to having a well-understood model of light curve evolution (Hsiao et al. 2007;Guy et al. 2010;Saunders et al. 2018;Leget et al. 2020;Kenworthy et al. 2021;Pierel et al. 2022), their standardizable absolute brightness can provide additional leverage for lens modeling by limiting the uncertainty caused by the mass-sheet degeneracy (Falco et al. 1985;Kolatt & Bartelmann 1998;Holz 2001;Oguri & Kawano 2003;Patel et al. 2014;Nordin et al. 2014;Rodney et al. 2015;Xu et al. 2016;Birrer et al. 2022a), though only in cases where millilensing and microlensing are not extreme (see Goobar et al. 2017;Foxley-Marrable et al. 2018;Dhawan et al. 2019).The first such discovery was iPTF16geu (Goobar et al. 2017), which had image separations resolved using adaptive-optics (AO) and HST, with very short time delays of ∼ 0.25-1.5 days (Dhawan et al. 2019).Nevertheless, the detection and analysis of objects like iPTF16geu are critical to the future of lensed SN research as unresolved, galaxy-scale lenses are expected to be relatively common amongst lensed SN discoveries made with the Vera C. Rubin Observatory (Collett 2015;Goldstein et al. 2019;Wojtak et al. 2019).
LensWatch1 is a collaboration with the goal of finding gravitationally lensed SNe, both by monitoring active transient surveys (e.g., Fremling et al. 2020;Jones et al. 2021) and by way of targeted surveys (Craig et al. 2021).The collaboration maintains a Cycle 28 HST program2 , given longterm (3-cycle) target of opportunity (ToO) status due to the relatively low expected lensed SN rates.The program includes three ToO triggers (two non-disruptive, one disruptive), and was designed to provide the high-resolution followup imaging for a ground-based lensed SN discovery, which is critical for galaxy-scale multiply-imaged SNe due to their small image separations (e.g., Goobar et al. 2017).
A new multiply-imaged SN Ia was discovered in 2022 August by the Zwicky Transient Facility (ZTF; Fremling et al. 2020) 3 , subsequently classified and analyzed by Goobar et al. (2022, hereafter G22).The separate four images of this SN 2022qmx (aka "SN Zwicky") were spatially unresolved in ground-based imaging with separations of ≲ 0.3 ′′ .In order to provide reliable photometry and the data necessary for accurate lens modeling of the system, optical space-based observations are ideal.We therefore report the first observations and results of the LensWatch collaboration, which triggered HST GO program 16264 to schedule follow-up imaging of SN Zwicky.
This work is the first in a series of papers that utilize data from the LensWatch program.Section 2 gives an overview of SN Zwicky and presents the final HST observation characteristics including triggering, orbit design, and implementation.Section 3 details our lens modeling methodology and constraints on the lensing system, and our analysis of SN Zwicky (including photometry and measurements of time delays and magnifications) are reported in Section 4. We conclude with a discussion of implications of this new dataset, as well as future observation plans, in Section 5.

OBSERVING WITH HST
As possible discovered lensed system configurations are highly variable, it is necessary to design a custom followup campaign for each new discovery.We therefore give an overview of the lensing system and SN characteristics for SN Zwicky, and then the subsequent observational choices made for the LensWatch HST ToO trigger.

The Multiply-Imaged SN Zwicky
The discovery, description of ground-based observations, and initial analysis of SN Zwicky are presented by G22.Briefly, the SN was discovered by ZTF at Palomar Observatory under the Bright Transient Survey (BTS; Fremling et al. 2020) on August 1, 2022 (MJD 59792).The SED Machine (SEDM) and Nordic Optical Telescope (NOT) provided spectroscopic classification of SN Zwicky as a Type Ia (SN Ia) at z = 0.35 and near maximum light on August 21-22, 2022 (MJD 59812-59813).Although the multiple images were not resolved by ZTF, the inferred absolute magnitude of SN Zwicky for this redshift was ∼ 3 magnitudes brighter than normal, suggesting the presence of strong gravitational lensing.G22 also obtained subsequent spectroscopic observations from the Keck observatory, Hobby-Eberly Telescope, and the Very Large Telescope (VLT), which led to a final SN redshift of z = 0.3544 and lensing galaxy redshift of z = 0.22615.The multiple images of SN Zwicky were first resolved with the Keck telescope Laser Guide Star aided Adaptive Optics (LGSAO) Near-IR Camera 2 (NIRC2) on September 15, 2022 (MJD 59837; see G22 for details).

ToO: Filter Choices & Orbit Design
3 https://www.wis-tns.org/object/2022qmx Roughly 12 days after the spectroscopic classification of SN Zwicky, we used a non-disruptive HST ToO trigger to obtain follow-up WFC3/UVIS and IR images of the lensing system.The average turnaround for a non-disruptive ToO trigger is ≳ 21 days, but close coordination with the HST scheduling team at STScI led to receiving our first images after 17 days on September 21, 2022 (MJD 59843), or ∼ 44 observer-frame (∼ 32 rest-frame) days post-discovery and ∼ 37 observer-frame (∼ 27 rest-frame) days after maximum brightness.
The anticipated image separations for a galaxy-scale lens of this mass and redshift are small enough that resolving the individual images with WFC3/IR (0.13 ′′ /pix), where the point-spread function (PSF) is severely undersampled, is unlikely.For the purposes of accurate photometry and lens modeling, the highest possible resolution imaging is required, and we therefore turned to WFC3/UVIS (0.04 ′′ /pix) to resolve the multiple images.We selected the F814W, F625W, and F475W filters to provide non-overlapping coverage across the full optical wavelength range (∼ 3,500-6,000 Å in the rest-frame; see Figure 1 and Table 1).Additionally, the ground-based follow-up campaign of SN Zwicky included (resolved) H-band Keck-AO imaging, and we therefore included (unresolved) WFC3/IR F160W observations to provide overall calibration and extra information about the lensing system.
The four filters were efficiently packed into a single orbit of observing using the 512 × 512 subarrays for both UVIS and IR imaging, even with three dithers per filter to reduce the impact of cosmic rays and provide optimal sampling of the (Figure 2).The four images of SN Zwicky were successfully resolved in the three UVIS filters, which provided a full-color image (Figure 3).In this section, we summarize the lens modeling analysis we carried out using the HST data presented in Section 2, leading to insights into the lensing galaxy mass (see Appendix).Given the very low number of identified strongly lensed SNe, this procedure has mainly been applied to strongly lensed quasars, e.g., by the H 0 Lenses in COSMOGRAIL's Wellspring (H0LiCOW) collaboration (e.g., Wong et al. 2020 Observer Wavelength (Å) Figure 1.The HST filters used to observe SN Zwicky, with restframe wavelength on the lower axis and observer-frame wavelength of the upper axis.The three bluer filters are from WFC3/UVIS, while F160W is from WFC3/IR.
For galaxy-scale lenses, the lens mass distribution is usually described by profiles such as the singular isothermal ellipsoid (SIE; Kormann et al. 1994) or the singular powerlaw elliptical mass distribution (SPEMD; Barkana 1998), in combination with an external shear component describing the influence from massive line-of-sight objects (McCully et al. 2017).These parameters can be constrained by the observed image positions alone (those measured in Section 4.1), and/or through the pixel intensities of the HST images.This requires a model of the lens light distribution, which is typically described by one or more stacked Sérsic profiles (De Vaucouleurs 1948;Sérsic 1963), as well as a model for the lensed SN represented by a PSF (described in Section 4.1).
Given these different potential methodologies, we used three independent software packages and five total methods to carry out independent analyses of SN Zwicky, which provides an examination of potential modeling systematics and allows us to marginalize over them (e.g., Ertl et al. 2022;Shajib et al. 2022).The three software packages are lfit_gui (Shu et al. 2016a), LENSTRONOMY (Birrer et al. 2015;Birrer & Amara 2018;Birrer et al. 2021) and the Gravitational Lens Efficient Explorer (GLEE; Suyu & Halkola 2010;Suyu et al. 2012;Ertl et al. 2022), and the five methods explored here are summarized by Table 2.
Using the first four models reported in Table 2, we found a significantly higher reduced-χ 2 when fitting both image positions and fluxes compared to fitting only image positions, indicating the presence of substructure and/or microlensing not captured by the lens modeling process.Additionally, the SALT+LS modeling method first performs PSF photometry to obtain individual image fluxes, which are fit directly with the SALT2.4model (Betoule et al. 2014).This initial fit is used to provide a prior on the magnifications given the standardizability of SN Zwicky as an SN Ia (see Section 4.2 for details), then the lens model parameters are constrained with the image positions.As is apparent in the following section, the SALT+LS model magnification estimates agree with the models constrained with positions only for images B and D but not A and C, further supporting our assumption that images A and C are impacted by additional factors.
As a result of the above, we rely only on the positions of the multiple images to constrain the initial four lens models.Our interpretation is that, pending updated difference image photometry, there is some combination of additional microlensing and/or millilensing impacting images A and C. We discuss this more in Section 5, and will wait for the upcoming template image to improve both our lens models and photometry.For the remainder of this work, we refer to the models used to constrain the lensing system parameters (lfit_gui, LS1, LS2, GLEE) as the "primary" models, and we refer to the SALT+LS model by name.

Lens Model Constraints
Each of the primary lens modeling methods described in Section 3.1 was used to independently constrain the magnifications and time delays for each SN image, as well as the Einstein radius (θ E ) of the lensing galaxy.These results are summarized in Table 3, where ∆t iA refers to the relative delay between the i th image and image A (i.e., t i − t A ), and we also report the "Final" combined measurement for each parameter.This final value is the equally-weighted average of each lens modeling result, and the uncertainty is a combination of the standard deviation of the primary model results and the statistical scatter in the average posterior distribution.As the models are equally weighted to obtain our final values for each parameter, we show normalized posterior distributions simply for visual comparison in Figures 4-6.Note that the models from Table 2 with fewer parameters also have narrower posterior distributions.We also include the results of the SALT+LS model, which reveals the impact of including information about the SN Ia standardizability.The additional specific lensing model parameters measured by each method, as well as more details about the modeling processes, are given in the Appendix.
While we expect improved constraints following the LensWatch template image scheduled for 2023 as we will be able to disentangle the SN and lensing galaxy flux more reliably, the level of agreement between the primary modeling methods gives us confidence in the final constraints and uncertainties.We also find good agreement between these key parameters with the modeling of G22, which used only the resolved near-IR Keck data.We note that our measured Einstein radius of θ E = (0.168 +0.009 −0.005 ) ′′ is the smallest detected value for a multiply-imaged SN thus far, and corresponds to a lens mass of ∼ 8 × 10 9 M ⊙ .The similar lensed SN iPTF16geu had a measured θ E = 0.29 ′′ (More et al. 2017;Mörtsell et al. 2020), with time delays of ∼ 1 day (Dhawan et al. 2019).Here we predict time delays of ∼ 0.2-0.5 days for each of the images of SN Zwicky, which is well below the predicted time-delay measurement uncertainty for even a  The SALT+LS model also includes constraints from the Type Ia absolute magnitude (see Section 3.1).
resolved and well-sampled lensed SN Ia due to the impacts of microlensing (e.g., Pierel et al. 2021;Huber et al. 2022).
We also note that the SALT+LS method, which uniquely uses the measured photometry to infer a standardized absolute magnitude measurement of SN Zwicky and sets a prior on the image magnifications (see Section 3.1), is generally in good agreement with other methods apart from the predicted magnifications for images A and D and a slightly lower θ E .The method also significantly reduced the plausible model parameter space (see Appendix), which lends weight to the claims that SN Ia standardization can significantly improve lens modeling efforts when microlensing is minimal.By implementing models that did not include this extra step alongside SALT+LS, the relative agreement (or disagreement) between methods was a useful indicator of additional substructure/microlensing beyond the primary lens modeling flexibility.Due to the compact nature of the lensing system and difficulty in disentangling the SN and lens galaxy flux, an identical "template" epoch has been scheduled for ∼ 6-12 months after the first, once the SN has long faded.This will provide more precise measurements and constraints for SN Zwicky and the lensing system in general.In the meantime, we have used PSF photometry to optimally measure the brightness of each SN image.HST photometry for SN Zwicky was measured using PSF photometry on the WFC3/UVIS "FLC" images, which are individual exposures that have been bias-subtracted, darksubtracted, and flat-fielded but not yet corrected for geometric distortion.The UVIS data processing also includes a charge transfer efficiency (CTE) correction, which results in FLC images instead of the FLT images used for WFC3/IR.The WFC3/UVIS2 pixel area map (PAM) for the corresponding subarray was also applied to each exposure to correct for pixel area variations across the images4 .In most cases, the individual exposures for each filter are "drizzled" together to create a single image (e.g., Figure 3).Here, we primarily are concerned with precisely measuring the position and brightness of each SN image, which (without a template image) requires accurate fitting of a PSF model.Drizzled images can introduce inconsistencies into the modeling of a PSF, and so we restrict ourselves to the FLCs to preserve the PSF structure.We use the standard HST PSF models5 to represent the PSF, which also take into account spatial variation across the detector.
For each UVIS filter, we implement a Bayesian nested sampling routine6 to simultaneously constrain the (common) SN flux and relative position in all three FLCs for all four SN images.Each PSF was fit to the multiple SN images within a 5 × 5 pixel square in an attempt to limit the contamination of both the lensing galaxy (as we assume a constant background in the fitted region) and other SN images.The PSF full-width at half-maximum (FWHM) for WFC3/UVIS is < 2 pixels, so this PSF size should include ∼ 99% of the total SN flux and not be contaminated by significant flux from the other images (each ≳ 5 pixels away).
The final measured flux is the integral of each full fitted PSF model, which is 101×101 pixels and large enough to approximately contain all of the SN flux.These corrected fluxes were converted to AB magnitudes using the time-dependent inverse sensitivity and filter pivot wavelengths provided with each data file.The final measured magnitudes and colors are reported in Table 4.

Single Epoch Time Delays and Magnifications from HST Photometry
We use the single epoch of HST photometry from Table 4 to constrain the time delays and magnifications for the multiple images of SN Zwicky in the manner of Rodney et al. (2021).As measuring the difference in time of peak brightness for each image directly (e.g., Rodney et al. 2016) is not possible with a single epoch, we instead constrain the age of each SN image given a single light curve model.The relative age difference for each image is also a measure of the time delay, though we note this method is only possible because we have a reliable model for the light (and color) curve evolution as SN Zwicky is of Type Ia.As we use some information from the unresolved light curve in G22, we also fit the data with the commonly used Spectral Adaptive Lightcurve Template (SALT2; Guy et al. 2010), which provides a simple parameterization of the Type Ia SN normalization (x 0 ),  We fit the photometry of the multiple images simultaneously using the SNTD software package (Pierel & Rodney 2019), where we also include the known effects of Milky Way dust (E(B − V ) = 0.16, R V = 3.1) based on the maps of Schlafly & Finkbeiner (2011) and extinction curve of Fitzpatrick (1999).Additionally, we use the simulations of Goldstein et al. (2018) to estimate the additional uncertainty introduced by chromatic microlensing.Using the time of peak estimate from G22 our observations are ∼ 46 days post explosion, which corresponds to ∼ 0.05, 0.05, and 0.11mag of additional color uncertainty in rest-frame U −B (∼F475W−F625W), B −V (∼F625W−F814W), and U − V (∼F475W−F814W) respectively (95% confidence; see Goldstein et al. 2018, figure 5).We add these uncertainties in quadrature to the color uncertainties shown in Table 4 for the fitting process.
We follow the methods outlined by Rodney et al. ( 2021) to measure time delays for SN Zwicky, which performed a similar analysis with the single epoch of SN Requiem.This process uses the SN color curves to constrain the time delay (with the SNTD "Color" method), and then fits for relative magnifications (with the SNTD "Series" method).Unlike the analysis of SN Requiem, an unresolved light curve exists for SN Zwicky and in G22 was analyzed to give t pk = 59808.6,c = 0.005, x 1 = 1.16.While our single epoch of photometry should constrain the color parameter, it will be unable to constrain the x 1 parameter and there will be significant degeneracies between time delays and t pk (as seen in Rodney et al. 2021).We therefore allow the t pk parameter, which here describes the time of peak for image A (see Figure 3 for naming convention), to vary only within fifteen days of 59808.6.We also fix x 1 to the parameter derived by G22 (x 1 = 1.16), mainly to ensure an accurate light curve standardization.We repeated the fitting first following the choice in Rodney et al. (2021) to set x 1 = 0 and second allowing x 1 to vary within 3σ of the value measured in G22, and found these varied the time delays by ≲ 0.5 days, well within our measurement error bars.We also checked the difference in measured time delays when fixing the value of t pk to 59808.6 and found a difference of ≲ 0.5 days.Finally, we note that the additional uncertainty added due to chromatic microlensing changes the measured time delay by ≲ 0.5 days as well, but increases the measurement uncertainties by ∼ 1-1.5 days.
SNTD finds a common value for c amongst all SN Zwicky images while varying the time delays of images B-D relative to the value of t pk , which describes image A, within relative bounds of [−15, 15] days.Figure 8 shows the measured colors and time delays for SN Zwicky with the bestfit SALT2 model overlaid.While all colors are used in the fit simultaneously, the photometric and model uncertainties mean F625W-F814W provides the most constraining power, followed by the colors that include the rest-frame ultraviolet.After these time delays have been measured with the SNTD Color method, we fix all best-fit parameters and use the SNTD Series method to estimate magnification ratios for images B-D (within bounds of [0.01, 100]) relative to image A. The fitting procedure for the light curve parameters is summarized in Table 5.As mentioned above, SNTD measures an overall normalization (x 0 ) and relative magnifications, and so we convert the combination of x 0 and magnification ratios to absolute magnitudes by assuming SN Zwicky is a perfect standardizable candle.Specifically, we apply light curve corrections based on Table 5 for stretch (x 1 = 1.16, with luminosity coefficient α = 0.14) and color (c = 0.03 with a luminosity coefficient of β = 3.1) in the manner of Scolnic et al. (2018) to obtain absolute magnitude estimates.We then compare the distance modulus of each image to the value predicted by a flat ΛCDM model (with H 0 = 70 kms −1 Mpc −1 , Ω m = 0.3) for an average SN Ia (M B = −19.36,Richardson et al. 2014) at z = 0.3544, which results in a measure of the absolute magnifications.We combine the statistical uncertainties on each measured magnification with a systematic uncertainty based on the intrinsic scatter of SN Ia absolute magnitudes (0.1 mag; Scolnic et al. 2018).The measured time delays and magnifications (with subscript "meas") are shown in Table 6 compared with lens model-predicted values from Section 3.2 (with subscript "pred").The posterior distributions for all parameters fit with SNTD (using the conversions listed above) are shown in Figures 9 and 10.While the relative time delay uncertainties are too large to provide a useful direct cosmological constraint, these results are a valuable check on our lens modeling predictions.The agreement also supports the plausibility of measuring time delays in a single epoch, at least when the lensed SN is of Type Ia and there is some constraint on the overall explosion date.

DISCUSSION
We have presented the first analysis of the LensWatch collaboration, which includes the only space-based observations of the gravitationally lensed and quadruply-imaged SN Zwicky.The images are resolved with HST WFC3/UVIS (PSF FWHM ∼ 0.07 ′′ ) but not with WFC3/IR (PSF FWHM ∼ 0.15 ′′ ).We have measured photometry for each SN image in three optical HST filters, and we use the resulting colors (with an additional uncertainty due to chromatic microlensing) to infer time delays (0.30We have also carried out an analysis of the lensing system using five distinct methodologies, of which we combine four primary methods to obtain our best constraints on the magnifications, time delays, and Einstein radius (θ E ) for the lensing system.We infer the smallest Einstein radius yet seen in a lensed SN (θ E = (0.168 +0.009 −0.005 ) ′′ ) and short time delays of ≲ 1 day, consistent with our color curve fitting results and G22.For SN images B and D, we find consistent magnification predictions across all of our lens models that are in good agreement with the measured values.However, our lens models are unable to fully explain the observed fluxes for images A and C, and we see a significant discrepancy between measured and predicted magnification (∼ 1.7 +0.8 −0.6 mag and ∼ 0.9 +0.8 −0.6 mag, respectively).We resort to variations from microlensing and/or millilensing (Metcalf & Zhao 2002;Foxley-Marrable et al. 2018;Goldstein et al. 2018;Hsueh et al. 2018;Huber et al. 2019) to explain this inconsistency.We do see some evidence for differential dust extinction and/or chromatic microlensing across the four images of SN Zwicky with relative differences in measured F475W-F814W of up to ∼ 0.18 ± 0.08 mag, and based on the work of Goldstein et al. (2018) our HST epoch is ∼ 3 rest-frame weeks outside of the "achromatic phase" of SN Ia microlensing where the impact on optical colors is expected to be ≲ 0.1 mag (95% confidence).We include the additional predicted uncertainty due to microlensing in our fitting.However, the F475W filter was most discrepant with the SALT2 fitting suggesting a possible systematic error in the photometry, and regardless this differential extinction is insufficient to explain the discrepancies we observe in images A and C relative to B (the largest color difference).We therefore estimate the additional (roughly achromatic) magnification to be ∼ 1 magnitude in both images, which is significant but well within expectations for average galaxy-scale microlensing (Pierel et al. 2021) or millilensing (Metcalf & Zhao 2002).While it is suggested that saddle images are more likely to be demagnified by microlensing (e.g., Schechter & Wambsganss 2002), we expect the method of detection for SN Zwicky would be significantly biased toward microlensing events with high magnifications.A template epoch is already scheduled for after SN Zwicky will have faded, which will drastically simplify and improve the photometry and lens modeling processes and provide more stringent constraints.
Although the multiply-imaged SN population is still small, this decade it is expected to grow by orders of magnitude with the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST; Ivezic et al. 2019) and Nancy Grace Roman Space Telescope (Pierel et al. 2021).Though both telescopes are expected to find a large number of spatially resolved lensing systems, unresolved lensed SNe discovered with Rubin will still be common, and a dedicated follow-up campaign from space similar to LensWatch will be necessary to provide accurate photometry and lens modeling for such systems.SN Zwicky is the first lensed SN analysis presented by the LensWatch program, and is an excellent example of the coordination that will be required for upcoming lensed SN cosmology efforts.The unresolved, ground-based discovery with ZTF and subsequent follow-up with HST is a glimpse at likely future discoveries with Rubin and follow-up with Roman (or HST), a strategy that can be extremely fruitful for the field of gravitationally lensed SN cosmology.

Individual Lens Modeling Results:
A. MODELING WITH GLEE A.1.Lens model parameterization We modeled SN Zwicky with an automated modeling pipeline that is based on the modeling software GLEE (Ertl et al. 2022;Suyu & Halkola 2010;Suyu et al. 2012), where we adopt the SPEMD (Barkana 1998) profile whose dimensionless surface mass density (or convergence) is given by where (x m , y m ) is the lens mass centroid, q m is the axis ratio of the elliptical mass distribution, θ E the Einstein radius, and γ is the power-law slope.The mass distribution is then rotated by the position angle ϕ m , where an elliptical mass distribution with an angle of ϕ m = 0 corresponds to elongation along the y-axis (after converting to the conventional definition of position angle).
The external shear strength is described by γ ext = γ 2 ext,1 + γ 2 ext,2 , with γ ext,1 and γ ext,2 the components of the shear.For a shear position angle ϕ ext = 0, the system is sheared along the x-direction.

A.2. Results based on SN image positions
First, we model the light of the lens galaxy with two Sérsic profiles, and the light of multiple lensed SN images by fitting a PSF model constructed from multiple stars in the field of the drizzled data.Ertl et al. (2022) showed that for lensed quasars we can achieve astrometric accuracy of 2 milli-arcseconds (mas) from the surface brightness (SB) fit, by comparing the modeled image positions to those measured by the Gaia satellite.We use our SN image positions (from PSF fitting) to constrain the mass parameters, since we did not find (and do not expect) any substantial lensed arc light (from the SN host galaxy) in the modeling residuals of the three UVIS bands.We show the results of our SB fit in Fig. 11.The measured astrometric positions of the four SN images in all three modeled bands are summarized in Tab. 7, and the lens light properties (based on the first and second brightness moments of the modeled lens light distribution that is a combination of the 2 Sérsics) in the F814W filter is in Tab. 8.The positions in the F475W and F625W bands are aligned with the F814W coordinate frame.For each band, we use the image positions reported in Tab.7 and adopt an uncertainty on the image positions of 4 mas to constrain the lens mass parameters.The 4 mas is an estimate based on the astrometric accuracy of 2 mas (Ertl et al. 2022) and to account for substructure lensing, which can perturb the image positions at the few mas level, as shown by Chen et al. (2007).We impose uniform priors on all 8 lens mass parameters that are tabulated in the leftmost column of Tab. 9. We report the mass model results from the F814W band because we achieved the lowest image-position χ 2 in this band.The results in this band are consistent with those of the other 2 bands, and also with the results of our light fit of the lens galaxy.We do not combine the constraints from all 3 bands in one single model, because positional uncertainties due to astrometric perturbations from e.g.substructures in the mass distributions are dominant and the measured SN image positions of the individual bands are thus not completely independent.Our final lens mass and shear parameters are presented in Tab. 9. Comparing to the modeled lens light in Tab. 8, the lens galaxy mass profile agrees well with the light in terms of centroid, axis ratio and position angle.In Tab. 10, we present the convergence κ and total shear strength γ tot at the (modeled) image positions.

A.3. Impact on results due to flux constraints
To investigate the bias to higher magnification for images A and C, we include fluxes, which were obtained from the light fit to the SN images (listed in Tab.7), in our model.The flux amplitude had typical uncertainties of ∼1%.We find that models based on image positions and fluxes try to fit to the fluxes at the expense of the poorer image position recovery, so the model cannot fit well to both positions and fluxes.
The image positions are close to a critical curve, so small shifts lead to large change in magnification.This is especially evident for the model where we use flux uncertainties from the SB fit.Imposing higher flux uncertainties (either 10% or 20% of the modeled flux value) leads to a lower image position χ 2 and a higher magnification χ 2 and brings the modeling results closer to the models where we used only image positions.We show the impact of including fluxes in our model by plotting the distribution of θ E in Fig. 12, magnifications in Fig. 13, and predicted time delays in Fig. 14 for the 4 different model classes.Since the models with flux constraints do not fit well to both the image positions and fluxes (with χ 2 ≳ 10), these models result in underestimated mass parameter uncertainties, as indicated by the narrower distributions for the blue, red and green models in Fig. 12. B. MODELLING WITH LFIT_GUI lfit_gui is a lens modeling software introduced by Shu et al. (2016a), which has been applied to about 150 strong-lens systems (Shu et al. 2016a(Shu et al. ,b, 2017;;Marques-Chaves et al. 2017;Wang et al. 2017;Shu et al. 2018;Marques-Chaves et al. 2020).In order to maintain an independent analysis, in the lfit_gui approach, positions and fluxes of the four SN images in the three optical filter bands are independently measured by fitting a photometric model consisting of two concentric Sersic components, four PSF components, and a constant component to the drizzled data downloaded from the Barbara A. Mikulski Archive for Space Telescopes Portal7 .The PSF models constructed in the GLEE approach are used.These photometric fitting results are shown in Figure 15.Overall speaking, the photometric model considered is able to reproduce the main structures in the data.Some residuals are seen at the lensed SN positions, which are primarily caused by PSF mismatches.The measured positions and fluxes of the four SN images are summarised in Table 11.The positional uncertainties are clearly correlated with the signal-tonoise ratios (S/Ns) of the lensed SN images.As a reulst.they are the smallest in F814W (≈ 1mas) and the highest in F475W (≈ 4mas).In general, the measured positions of the four lensed SN images agree well across the three bands.The largest differences are seen in the relative x positions of images C and D between F475W and F625W, which are about 1.8σ.The measured photometry for the four lensed SN images are found to be systematically brighter than the measurements in Table 4.The differences are typically within 0.15 mag in F625W and F814W and become 0.3-0.6 mag in F475W.We think this is likely related to the different treatments of the lens galaxy light.It affects photometry in the F475W the most because the brightness contrast between lensed SN images and the lens galaxy is the smallest.
In terms of lens modeling, the lfit_gui approach considered an SIE lens model, the convergence of which follows the profile defined in Equation A1 but with γ fixed to 2, and used the measured positions of the four SN images to constrain the five SIE parameters (as well as the source position) in the three bands separately.The sampling was done using the EnsembleSampler from the emcee package assuming uniform priors with sufficiently wide ranges for all the seven free parameters.The maximum a posteriori (MAP) estimation and marginalised posterior distribution for the three key SIE parameters, i.e.Einstein radius, axis ratio, and position angle, are reported in Table 12, and the posterior probability density distributions (PDFs) are provided in Figure 16.As shown in Figure 17, this lens model well reproduces the four lensed SN positions.The root mean square of the differences between the predicted and observed image positions is 0.002 ′′ , 0.0004 ′′ , and 0.0001 ′′ in F475W, F625W, and F814W respectively.The tightest constraints on the lens model parameters are obtained in F814W that has the smallest positional uncertainties and the posterior PDFs are the broadest in F475W.Nevertheless, the lens model parameters are generally consistent within 1σ.We find a clear anti-correlation between the Einstein radius and axis ratio (Figure 16), which is also observed in other lens modelling methods (e.g. Figure 21 and Figure 23).
We use the Bayesian information criterion (BIC) to combine results from the three bands, which is defined as where p is the number of free parameters (i.e.7), n is the number of constraints (i.e.8), and L is the maximum likelihood of a model.The BIC values are 15.9757, 14.9495, and 14.6745 in the F475W, F625W, and F814W.Weighting the results in F475W  the combined results suggest that the Einstein radius is 0.1664 +0.0010 −0.0019 arcsec, the axis ratio is 0.632 +0.027 −0.019 , and the position angle is 69.0 +0.3 −0.3 degrees (i.e. the angle between the major axis of the lens surface mass density distribution and the x axis, measured counterclockwise).The total lensing mass (within the ellipse that corresponds to κ = 1) is thus estimated to be 7.47 +0.09  −0.17 × 10 9 M ⊙ .The predicted magnifications, time delays, and convergence/shear values for the four lensed SN images are reported in Table 13 (and also in Table 3).We note that, strictly speaking, the BIC values can only be compared and combined when different models are constrained by the same data set, which does not apply to our results from three diferent bands.Nevertheless, the adopted weighting scheme is equivalent to weighting by the likelihoods, which is still a sensible treatment.Table 13.Predicted magnifications, relative time delays, and κ/γ values for the four lensed SN images by the lfit_gui approach (using the MAP estimation for individual bands and BIC-weighted average for the combined result).For the SIE model, κ and γ values are always the same.
C. MODELING WITH LENSTRONOMY LENSTRONOMY is a multi-purpose, open-source, community-lead, ASTROPY-affiliated gravitational lensing and image modeling package (Birrer & Amara 2018;Birrer et al. 2021) 8 .LENSTRONOMY supports a large variety of lens models and surface brightness profiles, as well as multiple numerical options to treat point sources.The modularity of LENSTRONOMY supports imaging modeling as well as catalogue-based model fitting.LENSTRONOMY has been applied for time-delay cosmography of lensed quasars (Birrer et al. 2019;Shajib et al. 2020) and lensed SNe as well as a variety of other lens modeling and image analysis applications (Gilman et al. 2020;Shajib et al. 2021;Schmidt et al. 2022).

C.1. The LS1 Method
We described the mass profile of the lens galaxy by a singular isothermal ellipsoid (Kormann et al. 1994), where the convergence is given by κ Here, θ E , q m , and (x m , y m ) are defined similarly as for equation A1.In order to fit the lens galaxy light profile, we stacked two Sérsic profiles (De Vaucouleurs 1948;Sérsic 1963): where I e is the intensity at the half-light radius R e .The constant b n is equal to 1.9992n − 0.3271 (Birrer & Amara 2018), and R ≡ q S x 2 + y 2 /q S with q S being the axis ratio of the Sérsic profile.The supernova images were fitted as point sources on the image planes with a PSF model.We initiated the model fitting with the PSF model constructed by the GLEE team and then further improved the PSF model using a built-in feature in LENSTRONOMY's that minimizes the residuals between the observed and reconstructed image around the supernova positions (Shajib et al. 2019).The comparison between the initial PSF model in the F814W band and the final reconstructed one is illustrated in Figure 18.Additionally, we adopted a circular region around the lensing system for likelihood computation to avoid the boundary effect of the PSF convolution in the evaluated likelihood function.
We fitted the pixel-level data from the three optical HST bands in a joint likelihood.The uncertainties on the model parameters were obtained from a Markov Chain Monte Carlo (MCMC) sampling.The flux ratios of the supernova images were not included in our lens model, because they failed to provide a good fit to the data and increased the reduced χ 2 from 1.17 to 5.59 (for the F814W filter).The reconstructed image model, source, convergence, and magnification model using the best-fit parameters from the converged MCMC chain are shown in Figure 19.Our measured κ and γ parameters at each image location are given in Table 14.

C.2. The LS2 Method
The "LS2" team used the catalog-data modeling functionality of the LENSTRONOMY software package, using the positions and positional uncertainties and redshifts reported in this paper.We adopted an elliptical power law mass profile (Tessore & Metcalf 2015) plus external shear, and allowed all parameters to vary.Models were computed for each of the three WFC3/UVIS filters, F475W, F625W, and F814W.The MCMC parameter sampling for F475W is shown in Figure 20, the posterior computed model parameters for the same filter are shown in Figure 21, and κ, γ results are given in Table 15.

C.3. The SALT+LS Method
As a proof of concept, in lens modeling we use the expected SN Ia brightness as prior with a broad standard deviation of 0.3 mag.First, using SNCosmo (Barbary et al. 2016), we fit SALT2.4 for the publicly available ZTF photometric data in g and r bands.We obtain x 1 = 1.11 ± 0.43, c = −0.071± 0.029, with maximum light at MJD = 59808.54± 0.43 (Figure 22).These values are in good agreement with the best-fit SALT parameters from G22.Note that G22 also used data from the Liverpool Telescope, which provided additional observations in griz bands.To standardize the SN Ia brightness, we adopt peak B magnitude of M B = −19.05,α = −0.141,and β = 3.101 (Betoule et al. 2014).We combine the three HST optical bands by averaging the best-fit PSF centroids and adding the fluxes.We assume a flat ΛCDM universe with H 0 = 70 km/s/Mpc and Ω M = 0.3.Our model consists of an elliptical power law (EPL) main lens and external shear.Our loss function combines the summed squares of the difference in the delensed image positions and summed squares of the difference between the observed and the model predicted fluxes.The relative weight between the flux and position terms is adjusted to achieve the best overall fit.The results are summarized in Figure 23.We find a somewhat steep mass profile slope for the lensing galaxy: γ EPL = 2.50.As G22 pointed out, with such a small θ E , this system is in a regime of lensing galaxies that have seldom been studied before.Such systems can be used to probe the density profile at sub-kpc scales within the lensing galaxy core.The image positions and model predictions agree to better than 0.023 ′′ (Figure 24).The magnifications from this model for images B & D (Table 3) are in fairly good agreement with the corresponding expectations in Table 6, whereas for images A & C, the model predictions are lower by < 2 σ (see also Figure 24).It appears that without taking into account microlensing and/or differential dust extinction, this is the best compromise the model can achieve.The total predicted magnification is 17.73.Compared with the expectation of 24.3 ± 2.7 from G22, it is smaller by 2.5 σ.
We find that without using SN Ia brightness as prior, it is possible to find models with acceptable predictions for both image positions and flux values, but they tend to have a much shallower slope (γ EPL ≲ 1.75).In contrast, using the SN Ia brightness prior, we find γ EPL to be consistently ≳ 2.5, whether we use single band data or combine the different bands.We also note that if the host galaxy identification in G22 is correct, this system is possibly in a unique situation in that the core of the host galaxy is not multiply-imaged.This makes the modeling of this system especially challenging: we cannot separately perform lens modeling using the lensed host galaxy in contrast to the other small-θ E lensed SN Ia (Dhawan et al. 2019).
We now briefly compare the SALT+LS model with the other four models presented in this paper that only use image positions.With regard to θ E , when fluxes are taken into consideration by the GLEE team, they have found acceptable models with θ E in agreement with the SALT+LS best-fit value.We further note that the time delay predictions from the SALT+LS model are in good agreement with those from the GLEE model that has a similar θ E .With regard to magnification: 1) The total magnification from the these four models ranges from 9.2 to 15.7 (Table 3).The SALT+LS model predicts a magnification of 17.43, higher than these four models.2) Whereas the SALT+LS model predicts the magnifications for the brighter two images, A and C, to be higher than those of B and D, the other four models predict the opposite.And yet, as mentioned before, even the total magnification from the SALT+LS model is ∼ 30% lower than the expected magnification from G22 or from Table 6 based on SN Ia brightness.The (κ, γ) values at the locations of the images are, in order of A-D: (0.25, 0.86), (0.22, 0.56), (0.26, 0.88), and (0.21, 0.60).
Given that this appears to be a fairly normal SN Ia, it is possible that microlensing and/or differential dust extinction (likely to be small, given the small color differences for the four images shown in Table 4) have played a significant role.If so, for the SALT+LS model, the optimization can be distorted in a way to compensate for these effects.Thus, for the uncertainties for κ and γ, we report the largest uncertainties for the four images, which are 0.024 and 0.023 respectively (Table 16).Once follow-up HST observations are completed after the SN has faded and improved photometry has been obtained, we will revisit this model.

Figure 2 .Figure 3 .
Figure 2. The layout of the orbit used for these HST observations.The dither sections (white) include other overheads such as filter changing.

Figure 4 .
Figure 4. Normalized posterior distributions for θE of the four primary lens models, and their combined constraint (dark grey mean, light grey uncertainty).The SALT+LS model constraint is also shown for comparison (red slashed).

Figure 5 .
Figure5.Normalized posterior distributions of the four primary lens models for the four image magnifications, and their combined constraints (dark grey mean, light grey uncertainty).The SALT+LS model constraints are also shown for comparison (red slashed).Note that images A,C have negative magnification and B,D have positive magnification, indicating all lens models agree on the parity of each image in addition to the absolute value of the magnification.

Figure 6 .
Figure 6.Normalized posterior distributions of the four primary lens models for the time delays relative to image A, and their combined constraints (dark grey mean, light grey uncertainty).The SALT+LS model constraints are also shown for comparison (red slashed).

F814WFigure 7 .
Figure 7. Results (right) of subtracting the best-fit PSF models from a single FLC for each WFC3/UVIS filter (left).Some residuals remain, particularly for the brightest images, but our planned template epoch will significantly improve the measured fluxes.

Figure 8 .
Figure 8. Measurements of time delay and color for each image of SN Zwicky, with image A the reference image (i.e., ∆t = 0).The vertical error bars are the photometric precision based on the work in Section 4.1 combined with an additional microlensing uncertainty , while the horizontal error bars are the 16 th and 84 th quantiles for the time delay posterior of images B-D (columns 2 and 3-5 in Figure 9), and t pk for image A. The grey shaded region is the best-fit SALT2 model from the SNTD Color method.

Figure 9 .Figure 10 .
Figure 9. Posterior distributions of the SNTD Color method fitting of HST photometry.The dashed vertical lines correspond to the distribution 16 th , th , and 84 th quantiles and the solid red lines show the final lens model predicted time delays and magnifications (those of A and C are off of the plot, see Section 3.1).The t pk parameter is given relative to the G22 value of 59808.6.

Figure 11 .
Figure 11.GLEE: Surface brightness fitting with GLEE in the three HST filters, as shown in the different rows.From left to right: observed image, model, and normalized residuals after modeling the light of the lens galaxy and the four SN images.

Figure 12 .
Figure 12.Exploration of the impact of flux constraints on the lens models with GLEE: Distribution of Einstein radii for the 4 different model classes.The flux constraints are based on the measured SN image amplitudes in Tab. 7, with typical uncertainties of ∼1%, unless boosted to 10% or 20% as indicated in the legend.

Figure 13 .
Figure 13.Similar to Fig. 12 from GLEE but for the SN image magnifications from the 4 different model classes.

Figure 14 .Figure 15 .
Figure 14.Similar to Fig. 12 from GLEE but for the distribution of predicted time delays from the 4 different model classes.

Figure 16 .
Figure16.Posterior PDFs for the key SIE model parameter in the lfit_gui approach.The blue, green, and red contours and histograms correspond to results inferred from the F475W, F625W, and F814W data, and the black contours and histograms correspond to the combined results.

Figure 18 .
Figure 18.The initial PSF model for the F814W band (left panel), the reconstruction after iterative improvement by LENSTRONOMY (middle panel), and the difference between the two (right panel).

Figure 19 .
Figure 19.The observed HST F814W-band data compared to the reconstructed model.Upper panel from left to right: the observed image, the reconstructed light intensity, and the normalized residuals.The circular mask illustrates the image region used in computing the likelihood.Lower panel from left to right: the reconstructed source, convergence (projected surface mass density), and magnification model with the four supernova images.In the reconstructed source panel, the blue star marks the unlensed position of the supernova, and no light from the host galaxy is detected above the noise level.

Figure 20 .
Figure 20.Parameter sampling for the LS2 Lenstronomy method, for the F475W filter position measurements.

Figure 21 .
Figure 21.The full corner diagram plot for the computed posterior estimates for the LS2 lens modeling fit, for the F475W filter position measurements.

Figure 22 .
Figure 22.ZTF g and r band light curves and best-fit SALT2.4 parameters.

Figure 23 .Figure 24 .
Figure 23.A corner plot of the posterior samples for the parameters of the SALT+LS lens model.

Table 2 .
Summary of lens modeling methodologies. a

Table 3 .
Lens modeling constraints on key parameters.
a See appendix for κ and γ results for each lens model.

Table 4 .
Photometry and colors measured for each image of SN Zwicky in AB magnitudes.

Table 5 .
Summary of the SALT2 parameters used in fitting SN Zwicky.

Table 6 .
Measured time delays and magnifications compared to the predictions from lens models from Section 3.

Table 7 .
GLEE: Astrometry and brightness of SN images -best-fit SN image positions and amplitudes from surface brightness fitting.

Table 8 .
GLEE: Centroid, axis-ratio, and position angle of the lens light computed from the second brightness moments of the two Sersic profiles in our best-fit model of the F814W filter.

Table 9 .
GLEE: Modeled mass and shear parameters in the F814W band.We present the median and 1σ uncertainties.Position angles are reported as East of North.

Table 10 .
GLEE: Convergence κ and total shear strength γtot at the (modeled) image positions.

Table 11 .
Relative positions (with respect to Image A that has the highest S/Ns in all three bands and thus smallest positional uncertainties) and AB magnitudes measured by the lfit_gui approach in F475W, F625W, and F814W respectively.A constant 0.1 mag uncertainty is assumed.

Table 12 .
Key SIE model parameters from maximum a posteriori (MAP) estimation and marginalised posterior distributions in the lfit_gui approach.
a Note that κ = γ for an SIE model.

Table 16 .
Best-fit κ/γ values for the SALT+LS lens model.