{'fraction_non_alphanumeric': 0.03595049892839859, 'fraction_numerical': 0.027861636669508905, 'mean_word_length': 4.529306829765545, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
LensWatch: I. Resolved HST Observations and Constraints on the Strongly-Lensed Type Ia Supernova 2022qmx ("SN Zwicky")

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"; Goobar et al. 2022) at z = 0.3544. SN Zwicky was discovered by the Zwicky Transient 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 (www.lenswatch.org). 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 produce an analysis of the lensing system using a variety of independent lens modeling methods. We find consistency between time delays estimated with the single epoch of HST photometry and the lens model predictions constrained through the multiple image positions, with both inferring time delays of<1 day. Our lens models converge to an Einstein radius of (0.168+0.009-0.005)", the smallest yet seen in a lensed SN. The"standard candle"nature of SN Zwicky provides magnification estimates independent of the lens modeling that are brighter by ~1.5 mag and ~0.8 mag for two of the four images, suggesting significant microlensing and/or additional substructure beyond the flexibility of our image-position mass models.


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 appear de-layed 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 (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;Tewes et al. 2013;Bonvin et al. 2017;Birrer et al. 2019;Bonvin et al. 2018Bonvin et al. , 2019bWong 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 are 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 stan-dardizable 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 adaptiveoptics (AO) and HST, with very short time delays. 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).
LensWatch 1 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 program 2 , 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 nondisruptive, one disruptive), and was designed to provide the high-resolution follow-up 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, here-after 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 obtain 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 follow-up 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
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.
For a galaxy-scale lens of this mass and redshift, the anticipated image separations are small enough that resolving the individual images with WFC3/IR (0.13 /pix) 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 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 pointspread function (PSF; 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 COS-MOGRAIL's Wellspring (H0LiCOW) collaboration (e.g., Wong et al. 2020).
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 power-law 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 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 cap- tured 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.4 model (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 as-sumption 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 de-lays 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 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.

HST Photometry
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, dark-subtracted, and flatfielded 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 images 4 . 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 models 5 to rep-resent the PSF, which also take into account spatial variation across the detector.
For each UVIS filter, we implement a Bayesian nested sampling routine 6 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  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).
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  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.
We fit (including model covariances) the data with the commonly used Spectral Adaptive Lightcurve Template (SALT3; Kenworthy et al. 2021), which provides a simple parameterization of the Type Ia SN normalization (x 0 ), shape or stretch (x 1 ), and color (c) used for light curve standardization. The remaining SALT3 parameters are a time of peak B-band brightness (t pk ) and the SN redshift. We fit the photometry of the multiple images simultaneously using the SNTD software package  (2011) and extinction curve of Fitzpatrick (1999). We employ the SNTD "Series" method, which is optimal for a sparse light curve and was used to measure time delays for the single epoch of SN Requiem (Rodney et al. 2021).
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 three 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 also repeated the fitting following the choice in Rodney et al. (2021) to set x 1 = 0, and found this varied the time delays by 0.5 days, well within error bars. We also checked the difference in measured time delays when fixing the value of t pk to 59808.6 and found a similar difference of 0.5 days.
SNTD finds a common value for c amongst all SN Zwicky images while varying the time delays (magnification ratios) of images B-D relative to the value of t pk (x 0 ), which describes image A. We set the (relative) time delay bounds to [−5, 5] and the magnification ratio bounds to [0.01, 10]. The fitting procedure for the light curve parameters is summarized in Table 5. We note excellent agreement be- tween our measurement of c and that of G22 using the full unresolved light curve. As mentioned above, SNTD measures an overall normalization (x 0 ) and relative magnifications as it is not specifically designed for lensed SNe Ia, 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.02 with a luminosity coefficient of β = 3.1) in the manner of Scolnic et al. (2018)   . 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.
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 Figure 8. 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 (0.04 /pix) but not with WFC3/IR (0.13 /pix). We have measured photometry for each SN image in three optical HST filters, which we use to infer time delays (∼ 0.24-1.58 days) and magnification ratios using the SNTD software package. Leveraging the fact that SN Zwicky is of Type Ia and therefore has a standardizable light curve, we apply a fiducial light curve standardization and obtain absolute magnification estimates of 7.75, 3.29, 6.26, and 3.38 for images A-D, respectively. We 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 light 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.5 and 0.8 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 dif-ferential dust extinction across the four images of SN Zwicky with relative differences in measured F475W-F814W of up to ∼ 0.3 mag, and based on the work of Goldstein et al. (2018) our HST epoch is just outside of the "achromatic phase" of SN Ia microlensing where the impact on optical colors is expected to be 0.1 mag (99% confidence). However, the F475W filter was most discrepant with the SALT3 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. The t pk parameter is given relative to the G22 value of 59808.6.
tially 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.

Acknowledgements
We would like to thank David Jones, Nao Suzuki and Taylor Hoyt for discussions helpful to the analysis in this work. This paper is based in part on observations with the NASA/ESA Hubble Space Telescope obtained from the Mikulski Archive for Space Telescopes at STScI; support was provided to JDRP and ME through program HST-  (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. 9. 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. 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.
Parameter Description position angle 155 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. 10, magnifications in Fig. 11, and predicted time delays in Fig. 12 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. 10.
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, 2017Marques-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 Portal 7 . The PSF models constructed in the GLEE approach are used. These photometric fitting results are shown in Figure 13. 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-to-noise 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 14. As shown in Figure 15 1σ. We find a clear anti-correlation between the Einstein radius and axis ratio (Figure 14), which is also observed in other lens modelling methods (e.g. Figure 19 and Figure 21).
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 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 ( 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.   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 builtin 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 16. 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

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 18, and the posterior computed model parameters for the same filter are shown in Figure 19.

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 20). 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 21. We find a somewhat steep mass profile slope for the lensing galaxy: γ = 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 22). The total predicted magnification is 17.73 (Table 3). Compared with the expectation of 24.3 ± 2.7 from G22, it is smaller by 2.5 σ. The model predictions for the brightnesses of images A and C are within ∼ 6% but are higher by just under 30% for the two fainter images (B and D). It appears that without taking into account microlensing and/or differential dust extinction, this is the best compromise the model can achieve.
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 (γ 1.75). In contrast, using the SN Ia brightness prior, we find γ 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 based on SN Ia brightness. 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 The gray scale background shows the convergence (κ) map. This image is shown in a slightly different orientation relative to most of the figures in this paper, but is the same as Figure 17. The brightest image, A, is still at the very top, and the other images are in the same alphabetical order clockwise from this reference point. model, the optimization can be distorted in a way to compensate for these effects. On the other hand, if the total magnification estimation of 20.68 in Table 6 is correct, our prediction is only lower by ∼ 16%, not far from the intrinsic scatter of SN Ia luminosity. Once follow-up HST observations are completed after the SN has faded and improved photometry has been obtained, we will revisit this model.