GausSN: Bayesian Time-Delay Estimation for Strongly Lensed Supernovae

We present GausSN, a Bayesian semi-parametric Gaussian Process (GP) model for time-delay estimation with resolved systems of gravitationally lensed supernovae (glSNe). GausSN models the underlying light curve non-parametrically using a GP. Without assuming a template light curve for each SN type, GausSN fits for the time delays of all images using data in any number of wavelength filters simultaneously. We also introduce a novel time-varying magnification model to capture the effects of microlensing alongside time-delay estimation. In this analysis, we model the time-varying relative magnification as a sigmoid function, as well as a constant for comparison to existing time-delay estimation approaches. We demonstrate that GausSN provides robust time-delay estimates for simulations of glSNe from the Nancy Grace Roman Space Telescope and the Vera C. Rubin Observatory's Legacy Survey of Space and Time (Rubin-LSST). We find that up to 43.6% of time-delay estimates from Roman and 52.9% from Rubin-LSST have fractional errors of less than 5%. We then apply GausSN to SN Refsdal and find the time delay for the fifth image is consistent with the original analysis, regardless of microlensing treatment. Therefore, GausSN maintains the level of precision and accuracy achieved by existing time-delay extraction methods with fewer assumptions about the underlying shape of the light curve than template-based approaches, while incorporating microlensing into the statistical error budget. GausSN is scalable for time-delay cosmography analyses given current projections of glSNe discovery rates from Rubin-LSST and Roman.


INTRODUCTION
Strong lensing of a background variable source, such as a supernova (SN), by a foreground galaxy or galaxy cluster, can result in the appearance of multiple images of the source.Refsdal (1964) was the first to demonstrate that the time delay between the multiple images of a gravitationally lensed SN (glSN) images could be related to  0 when combined with a model of the lens mass distribution.A independent  0 estimate from strong lensing is motivated by the persistent 5 tension between the present-day expansion rate of the universe,  0 , measured using light from early-times (i.e. the Cosmic Microwave Background; Planck Collaboration et al. 2020) and from late-times (i.e. the distance ladder; Riess et al. 2022), known as the Hubble tension.If determined to be physical, this tension may be suggestive of new physics beyond our present model of the universe, ΛCDM (Mörtsell & Dhawan 2018;Di Valentino et al. 2021).As time-delay distances inferred from strong lensing are completely independent of the luminosity distances used in the distance ladder, time-delay cosmography represents a promising technique for crosschecking local measurements of  0 (see e.g.Treu et al. 2022;Suyu et al. 2023 for recent reviews).In this paper, we present GausSN: a ★ Contact e-mail: eeh55@cam.ac.ukBayesian semi-parametric approach for the extraction of time delays from glSN systems using Gaussian Processes (GPs).GausSN is a publicly available package that can be found at: https://github.com/erinhay/gaussn.
Because of the rarity of glSNe, the first cosmological constraints from strong lensing came from quasars (e.g.Keeton & Kochanek 1997;Wong et al. 2020;Birrer et al. 2020).However, SNe have several advantages over quasars for strong lensing analyses which have the potential to push the field further in precision and accuracy.Firstly, that SNe disappear after a few months allows for the study of the lens and host galaxy in greater detail.Improved knowledge of the lens helps to break the mass-sheet degeneracy (Falco et al. 1985), which was cited as the main source of uncertainty in the H0LiCOW estimate of  0 (Wong et al. 2020;Birrer et al. 2020).If the lensed object is a Type Ia SN (SN Ia), further leverage can be gained from the fact that SNe Ia are standardizable candles (Foxley-Marrable et al. 2018;Birrer et al. 2022).Secondly, the relatively simple light curves of SNe, in contrast to the stochastic variability of quasars, makes it easier to extract time delays from the systems.In particular, the achromatic expansion phase of SNe Ia, during which microlensing has been shown to cancel in color curves (Goldstein et al. 2018;Huber et al. 2021), may contribute to improved constraints on the microlensing affecting the system (Bonvin et al. 2019a).Microlens-ing, which can distort the shape of an image's light curve and lead to artificial shifts in the apparent peak of the light curve, is an important source of uncertainty that has proven difficult to disentangle from the underlying variability of the source (Foxley-Marrable et al. 2018;Pierel et al. 2023b;Weisenbach et al. 2024).Lastly, time-delay estimation with glSNe can also be done with a shorter time frame of observations because their variability is constrained to a period of a few weeks to months, compared to the years required to measure reliable time delays from lensed quasars (Bonvin et al. 2017(Bonvin et al. , 2018(Bonvin et al. , 2019b)).Of course, the rarity of SNe and short window in which they are active make glSNe more difficult to discover.
Therefore, it was not until November 2014 that the first resolved multiply-imaged SN, SN Refsdal, was discovered (Kelly et al. 2015).SN Refsdal, a peculiar Type II SN at  = 1.49(Kelly et al. 2016b), was discovered with four images in an Einstein cross due to lensing by an elliptical galaxy in the MACS J1149.6+2223galaxy cluster at  = 0.54.Unfortunately, the time delays from galaxy-scale lensing, originally estimated in (Rodney et al. 2016), were determined to be too short and imprecise to obtain a measurement of  0 .However, the host galaxy of SN Refsdal was also imaged multiple times.It was predicted that a 5 th image of SN Refsdal would later appear in another image of the host galaxy with a time delay of just under a year (Kelly et al. 2015;Oguri 2015;Sharon & Johnson 2015;Grillo et al. 2016;Diego et al. 2016;Jauzac et al. 2016;Kawamata et al. 2016;Treu et al. 2016).As predicted, SN Refsdal reappeared less than a year later in October 2015 (Kelly et al. 2016a).Using the time delay between SN Refsdal's fifth image and the other four images (Kelly et al. 2023b), the first measurement of  0 was published in May 2023 (Kelly et al. 2023a).Finding  0 = 64.8+4.4  −4.3 km s −1 Mpc −1a 6.0% precision estimate -SN Refsdal alone has demonstrated the importance of glSNe as precise local cosmological probes.
In the next decade, the Vera C. Rubin Observatory's Legacy Survey of Space and Time (Rubin-LSST) and Roman Space Telescope (Roman) are expected to discover tens to hundreds of glSNe, if we adopt the right search strategies (Huber et al. 2019;Wojtak et al. 2019;Pierel et al. 2021;Craig et al. 2021).With this data, it will be possible to reach percent level constraints on the cosmological parameters competitive with SH0ES and Planck measurements (Huber et al. 2019;Arendse et al. 2023).As samples of glSNe grow, so does their potential to resolve the Hubble tension with precise local measurements of  0 .
The relatively simple and well-understood variability of SNe lends itself to time-delay estimation techniques based around SN light curve templates -empirical models for the light curves or spectral energy distributions (SEDs) of SNe.Supernova Time Delays (sntd; Pierel & Rodney 2019) -which uses light curve templates or parametric models -has emerged as the standard for glSNe time-delay estimation due to its flexibility and accessibility.The package is built upon sncosmo (Barbary et al. 2023), which contains a diverse library of SN light curve templates.A variety of methods have been used to construct SN light curve templates, so now there are a selection of templates for each of the many types of SNe, e.g.Type Ia, Type Ib, or Type IIP.Of course, the use of templates therefore requires knowledge of the SN type, which is most reliably determined via classification of an observed SN spectrum (Filippenko 1997;Gal-Yam 2017).Also, the redshift of the SN is always required to use templates to shift the rest-frame template to the observer-frame in wavelength and time.Alternatively, flexible parametric models, such as the Bazin function (Bazin et al. 2009), have been created to fit SN light curves if spectroscopic classification is unavailable or ambiguous.Given a specified SN light curve template and redshift, sntd uses the nested sampling functionality available through sncosmo to yield posterior distributions for the light curve template parameters, the time delay, and the magnification.
However, uncertainties relating to template choice and the unknown redshift evolution of SN light curves, especially for types of SNe other than SNe Ia, may introduce systematic errors into the timedelay analysis that are difficult to quantify.While template-based approaches are motivated by the well-described nature of some SN light curves, a complementary method which is independent of these potential systematics is needed.GausSN is a template-independent approach for time-delay estimation that models the underlying light curve non-parametrically using a Gaussian Process (GP).This model is motivated by previous work on quasar time-delay estimation with GPs, in particular that of Tak et al. (2017), as well as that of Hojjati et al. (2013); Hojjati & Linder (2014); Meyer et al. (2023).
By modelling the underlying light curve non-parametrically with a GP, GausSN is able to fit the light curves of any type of SN without prior knowledge of the type or the redshift of the object.In addition, the model naturally fits for the time delays of all images, observed in any number of bands, simultaneously.By fitting all images in all filters together, we take advantage of the full information available from the data by leveraging the knowledge that for each band, the multiple images' time-series are time-shifted and magnified realization of the same underlying light curve.GausSN also implements a novel microlensing treatment which occurs alongside time-delay estimation for a one-step time-delay measurement.A time-varying magnification term is used to account for both macrolensing and microlensing simultaneously with the time-delay estimate.Because the time delay and microlensing are jointly inferred, GausSN does not require post-processing to account for a systematic error budget due to microlensing, but instead incorporates uncertainty from microlensing in the statistical uncertainty.Therefore, GausSN provides a coherent Bayesian time-delay estimate in one-step, without the need for post-processing to account for systematics from template choice or microlensing.
The peculiar nature of SN Refsdal's light curve exemplifies the need for a template-independent approach for time-delay estimation for glSNe, such as GausSN.Although most similar to SN 1987A (Woosley et al. 1988;Kelly et al. 2016b), SN Refsdal's light curves are not well described by any existing templates (although some theoretical models have been developed, e.g.Baklanov et al. 2021).While the light curves of SNe Ia are highly uniform, the shapes of the light curves of other types of SNe are less well studied and more highly varied from object to object.GPs are well suited for this application, as they are flexible and data driven, with minimal assumptions about the properties of the true underlying light curve.
Indeed, in the analysis of the time delays of SN Refsdal's five images, the Bayesian Gaussian Process method, upon which GausSN is based, was highly successful (Kelly et al. 2023b).Four methods, including sntd, used with a Bazin function as the parametric model because an SED template does not exist for SN Refsdal, and a custom piece-wise polynomial template approach, were tested on simulations of SN Refsdal-like light curves to determine the quality of each method's fits.In simulations, the Bayesian Gaussian Process method outperformed the parametric model-based approaches for extracting the time delays of the 2 best-sampled images S2 and S3 relative to the first image.On the other hand, sntd and the custom piecewise polynomial approach were more successful in estimating the time delay of the more poorly sampled images S4 and SX relative to the first image.Ultimately, the custom piece-wise polynomial approach was the most generally successful in fitting for time delays and magnifications for all the images in simulations.
GausSN provides coherent Bayesian time-delay estimates with complementary strengths to existing template-based time-delay estimation methods.As a template-independent approach, GausSN does not require prior knowledge of SN type or redshift, making the model quick and easy to implement.Furthermore, GausSN is not subject to potential sources of systematic uncertainties and biases from template choice; fine tuning of generic SN models, such as the Bazin function; or the need to construct custom models, which may not be feasible for larger samples of glSNe and makes validation of an analysis difficult.Finally, the model provides Bayesian timedelay estimates in one step, without the need for post-processing to account for systematic uncertainty from microlensing, by incorporating a time-varying magnification treatment.
In §2, we describe how we model glSN in flux-space with GPs.We describe the microlensing model in §3 and the sampling algorithms implemented in GausSN in §4.In §5, we show the fits to simulated doubly-imaged glSN light curves as expected to be seen from Rubin-LSST and Roman.We then apply GausSN to the five images of SN Refsdal in §6.In §7, we discuss the performance of GausSN compared to other leading time-delay extraction techniques, as well as comment on future work possible with GausSN.In §8, we conclude.

MODELLING GLSNE WITH GAUSSIAN PROCESSES IN FLUX-SPACE
Strong lensing of a SN by a foreground galaxy or galaxy cluster results in the appearance of multiple images, which are copies of the same underlying light curve.These multiple images are observed with some shift in time and magnification or de-magnification relative to each another.We model the true underlying light curve in flux space,  (), as a draw from a Gaussian Process: where () is the mean function, which can depend on time, t, in the case of time series data, and  (,  ′ ) is the covariance (or kernel) function, which gives the covariance between  () and  ( ′ ).
Throughout this work, we use a zero mean function, such that () = 0.The zero mean function encourages the inferred function to go to zero flux outside of where there is data, which reflects the physical expectation that before the SN explosion and after the SN has faded, we expect an average of zero background-subtracted flux.For the covariance function, we choose the squared exponential kernel: where , the amplitude, and , the length scale, are the two hyperparameters controlling the kernel.This kernel enforces strong correlation between points nearby to one another, with weaker correlation as points further away are considered.It is also stationary, or invariant to overall time shifts, and produces smooth functions (Rasmussen & Williams 2006).The squared exponential kernel has proven to be well-suited for fitting SN light curve data, as it has been used to do so in previous work with success (Kim et al. 2013;Boone 2019;Vincenzi et al. 2019;Qu et al. 2021).
The GP is defined such that any vector  consisting of the evaluations of this continuous function,  (), at the finite set of points, , has a joint multivariate Gaussian prior distribution: where the elements of the covariance matrix,  (, ), are given by    =  (  ,   ) for   ,   ∈ .Note that we use  ∼ N ( , ) to denote a multivariate normal vector with mean  and covariance matrix .

Two Images in One Band
We first consider the simplest case of two images observed in one band, i.e. a filter covering a defined range of wavelengths.Taking image 1 to be the arbitrarily chosen reference image, the light curves of the two images are described by: where Δ is the relative time delay of image 2 compared to image 1 and  is the relative magnification of image 2 compared to image 1. Say we observe image 1 at a set of  1 times, t1 , and image 2 at a set of  2 times, t2 .We define f1 ( t1 ) and f2 ( t2 ), the observations of image 1 and 2 respectively, such that for the  th observation of image 1: and for the  th observation of image 2: where  1, and  2,  are the standard deviations of the measurement errors of image 1 and image 2, respectively.We assume that the measurement errors are independent.We re-scale the flux data and measurement standard deviations for each band and image by a single, constant factor determined by the range of fluxes over all images and bands.
We concatenate these two vectors to give the flux data vector, F, as: and the time data vector, T, as: In addition, we define the de-shifted time vector, TΔ , which depends on the time delay, Δ: To infer  () from the data F and T, we fit for a time delay, Δ, and magnification, , in addition to the two kernel hyperparameters.Therefore,  = ( , , Δ, ).
The marginal likelihood over  is given as: where we can think of   in four quadrants: Recall that for a GP, the covariance is given as Cov (  (),  ( ′ )) =  (,  ′ ).With Equations 5-8, the covariances in the first quadrant are: assuming there is no covariance between  1, and  1,  , as the physical process generating the light curve is independent of the measurement process.In the second quadrant, the covariances are: It follows that: and for the remaining two quadrants.Because our choice of kernel is stationary,  (t 2, − Δ, t2,  − Δ) =  (t 2, , t2,  ).
Based on this formulation, we can decompose   into three factors: where ⊙ represents a Hadamard product (i.e. the elementwise product).The first factor, , depends only on .Defining 1 as a matrix of ones,  is given as: where, again  1 is the number of observations of image 1 and  2 is the number of observations of image 2. The second factor, , depends on Δ and the kernel parameters and is defined as: where  ( t, t′ ) is the matrix whose (, ) th element is  (t  , t′  ).Finally, the diagonal matrix  contains the variances of the measurement errors, σ2 , along the diagonal.
With a specified prior on  and Equation 12, we can construct the posterior: which we can sample from and marginalize over to determine the probability distributions over the parameters.The dependence of ( | T) on T arises because we will restrict the time delay, Δ, to be within the range of observations.

Two Images in Two Bands
We will now consider the case with two images in two bands -band A and band B. We now have two functions we wish to model as a draw from a GP: and where the draws from the GP for each band are independent.While we choose to neglect correlations between the underlying light curves in different bands, these correlations can in principle be included, for example as in Hu & Tak (2020).Alternatively, the model could be used to fit color curves, which would encode correlations between bands, instead of the light curves in individual bands.Indeed, Pierel et al. (2021) found that fitting color curves with sntd gives better performance for time-delay estimation in simulations.
The light curves of the images are therefore defined by: We note that the magnification and time delay of the second image relative to the first is the same in band A as it is in band B. Therefore, there is still only one  and one Δ for which to fit, so it remains that  = ( , , Δ, ).It is possible to adopt unique GP kernel hyperparameters for each band, but we find that these extra parameters are unnecessary for the examples that follow in this paper.We construct the flux data vector, F, as: F = ( f1,A , f2,A , f1,B , f2,B ) ⊤ (25) and the time data vector, T, as: where f, are the observations of image  in band  at times t, .
We also define the de-shifted time vector, T , as: The marginal likelihood is: where, again,   = (  ⊙ ) + .
We make the simplifying assumption that there is no shared information between bands beyond the time delay and magnification.In other words, there is no covariance of either image in band A with either image in band B. With this in mind, we write  as: where 0 denotes a matrix of all zeros and  A and  B are defined as in Equation 20 for each band separately.We point out that there does not need to be the same number of observations of the images in band A as there are in band B. Therefore, the shape of 0 need not necessarily be square.Finally, the posterior is the same in this case as in Equation 21 with the new marginal likelihood terms defined for this data case.The case of two images in two bands can be generalized to fit multiple (> 2) images in multiple (> 2) bands.We fit simulated multi-band data as expected from Rubin-LSST and Roman in §5 and multi-image data from SN Refsdal in §6.

MICROLENSING
The effect of microlensing from stars and substructure in the lensing galaxy or galaxy cluster can additionally introduce additional timevarying magnification to individual images.Microlensing in systems of glSNe is extremely difficult to model because of complex layouts of stars and dark substructure in the foreground of lensing systems.The discrepancy between the observed and the model-predicted brightnesses of the four images from SN Zwicky without considerations of microlensing and the observed brightnesses of the images -up to 1.5 mag for one image -demonstrates the difficulty of modeling this effect and the significant impact which microlensing/substructure can have on a system (Pierel et al. 2023b;Goobar et al. 2023b).
Following Tak et al. (2017), we model microlensing as a timedependent extension of , so the previously constant  becomes ().Because we can only learn relative magnification effects from the light curve data available, we take all magnification to be relative to an arbitrarily chosen image 1.Therefore, Equations 5-8 become: where again  1, and  2,  are the standard deviations of the measurement errors of image 1 and image 2, respectively.2021) illustrate from microlensing maps examples of microlensing as a function of time.In these example microlensing curves, a period of constant magnification is typically interrupted by one change in brightness which occurs over a range of timescales, from less than a day to weeks.Physically, this shape corresponds to the SN crossing a microcaustic as the photosphere expands.After this change, the observed magnification appears to remain constant, as the small size of the SN and microlensing sub-structures means the photosphere most often crosses only one significant microcaustic.Therefore, we have chosen to model microlensing as a sigmoid function, given by: where  0 is the macrolensing effect,  1 is the scale of the microlensing effect,  is the rate of change in the microlensing effect, and  0 is the location of a change in microlensing.
We make the simplifying assumption of achromatic microlensing, or microlensing which is identical across bands for a given image.It is demonstrated by Goldstein et al. (2018); Huber et al. (2021) that this assumption is valid for the first three rest-frame weeks of lensed SN Ia, in which time microlensing is predicted to be largely achromatic.This model for the time delay and magnification therefore takes 5 parameters per image.The parameters we fit for in a doublyimaged system are now  = ( , , Δ,  0 ,  1 , ,  0 ).
Consider again the case of two images observed in one band, where image 1 is observed at a set of  1 times and image 2 is observed at a set of  2 times.For this system, the matrix  from Equation 19 is now: where 1 is a vector of ones of length  1 and  2 is a vector of length  2 whose  th element is given by  2,  = (t 2,  − Δ).This formulation generalizes to an arbitrary number of images observed in an arbitrary number of bands.We will refer to the model with a constant  as the "constant magnification" model and the model with a time-varying magnification term as the "sigmoid magnification" model.Figure 1 demonstrates how a time-varying magnification term can have a significant impact on the shape of a SN light curve, to the point that the inferred time delay could be significantly skewed if the glSN system is fit with a constant magnification model.A combined treatment of macrolensing and microlensing mitigates the potential source of bias that arises if these magnification effects are considered separately.We explore this effect in greater detail in Appendix A.
We emphasize that the choice of a single parameterization to describe the microlensing may introduce a source of systematic uncertainty, which should be carefully considered in a cosmological analysis.In §7, we discuss future work which will mitigate systematics that may arise from a choice of parameterization of microlensing.In addition, we test a second parameterization of microlensing -a flexible sinusoidal function -on a subset of the data used in the main analysis in Appendix B.

Sampling Algorithms
Within the publicly available GausSN framework, we implement three methods for sampling from the posterior distributions: • emcee: affine invariant MCMC ensemble sampler (Foreman-Mackey et al. 2013;Goodman & Weare 2010) • dynesty: nested sampling (Speagle 2020;Skilling 2004;Skilling 2006) • zeus: MCMC ensemble slice sampler (Karamanis & Beutler 2021;Karamanis et al. 2021) We choose to sample parameters with dynesty for the analyses in this paper.The nested sampling approach is best able to handle the multi-modal nature of the posteriors on the time delay and other magnification parameters.Therefore, for the most accurate timedelay estimates and associated uncertainties, we opt for dynesty.
Within dynesty, we use both uniform sampling with multiellipsoidal bounding (Feroz et al. 2009), for systems with fewer than 100 data points and fewer than 5 parameters, and random slice sampling, for systems with more than 100 data points or more than 5 parameters, with 500 live points.The slice sampling technique was developed by Neal (2003) and first implemented within the nested sampling framework by Handley et al. (2015a,b).We retain the default stopping criteria in dynesty, which is met when the remaining, unaccounted for evidence is less than a threshold which depends on the number of live points used (Speagle 2020).We note that the specifications of dynesty sampling, and all other sampling methods, can be easily adjusted within the GausSN parameter optimization function.
Given these specifications, GausSN performs as follows.For an object with 320 total data points, i.e. 40 observations of 2 images in 4 bands, a single evaluation of the marginal likelihood takes ∼ 2 ms using standard CPU resources on a desktop computer.The stopping criteria is typically met after 200,000-800,000 evaluations of the likelihood function, which roughly corresponds to 5,000-15,000 nested sampling iterations.Therefore, the nested sampling algorithm takes between 5 and 30 minutes to run for such an object.Depending on the quantity and quality of the data, sampling may take anywhere from 30 s to an hour for the simulated Rubin-LSST and Roman data.Given the rarity of glSNe, even in the best case scenarios for future observatories, such run times will not be a barrier for future applications of our model to real data.

Plotting the Fits
When considering example glSNe data in §5 and §6, we will plot draws of the fitted light curves to demonstrate the quality of the GausSN fits.Here, we describe the procedure for creating the fitted curves for a system with  images.First, we take a random draw of the time delay, Δ, and relative magnification, , from the joint posterior and de-shift and un-magnify the data from the ( − 1) trailing images according to the drawn parameters.Setting the mean and covariance of a multivariate normal distribution to the posterior predictive mean and covariance conditioned on the now overlapping data from all images, we take a realization of the fit.We plot a copy of this realization for each of the  images, where each copy of the realization is time-shifted and magnified back to the frame of the observed data, so that it overlaps with the raw light curve data of each image.We repeat this process for 200 samples of the joint posterior to get a representative visualization of the fitted light curve and uncertainty for each image.Therefore, the final plot shows the raw data for all  images with 200 realizations of the fitted light curve per image.
The procedure follows the same steps when considering a timevarying magnification term.In the first step, the data from each image will be un-magnified according to the drawn time-varying magnification function.When the realization of the light curve is being plotted, it is then put back into the observed data frame by magnifying according to the time-varying magnification term to match the original data.In this case, we will additionally plot realization of the magnification function for each posterior sample in a separate panel of the figure.The realizations of the magnification function alongside the data and fitted light curves allows for physical interpretations of the magnification term.

Data
Recently, Pierel et al. (2021) (hereafter P21) simulated 2.4 million glSN light curves -600,000 each of Type Ia, Ib/c, IIn, and IIPas expected from the Roman Space Telescope. 1 The cadence, depth, and detection thresholds for the simulations are based on the Roman SN survey "All-z" strategy described in Hounsell et al. (2018), with modifications made for more recent survey updates.Based on the current plans for the instrument, P21 predicts Roman will discover glSNe up to  = 4.
The simulation pipeline works as follows: for each SN type, a sample of 50 simulated galaxy-scale lenses are used to simulate 10 distinct glSN light curves with 2-4 images based on the structure of the lens.Each of these 500 systems are then subject to 100 iterations of microlensing for each of 12 microlensing maps, yielding 1200 variations of each glSN light curve.P21 assumes achromatic microlensing, so any individual image experiences the same microlensing effects across wavelength space.The 12 microlensing maps are based on different choices of stellar mass model, which vary the effective radius, initial mass function, and Sérsic index of the galaxy profile.Objects are required to pass a series of data cuts, which are described in P21.For simplicity, we consider only the doublyimaged glSNe from the P21 Roman simulations.The distributions of parameters for the P21 Roman simulations are shown in Figure 2.
We fit each object with both the constant magnification model and the sigmoid magnification model.For the constant magnification model, we use the following priors on the two kernel parameters, the time delay, and the magnification 2 : ∼ U (10, 40), (37) where  Δ is the difference between the times of the brightest observations of image 2 and 1 and   is the ratio of the brightest observation of image 2 relative to the brightest observation of image 1.We define t1 and t2 to be all the times of observation of image 1 and 2, respectively.The prior on Δ therefore requires that there always be overlap between data from image 1 and data from image 2. Although the scaled flux data has a maximum value of 1, we allow the prior on  to range from 0 to 5 so as not to overly restrict the amplitude of the underlying light curve.We adopt wide priors on all parameters to ensure the diversity of behavior present in the simulations is represented in our parameter space.With real glSN events, which will be fit on an individual basis because of their rarity, these priors can be adjusted to account for additional contextual information for each specific event.
For the sigmoid model, we set the following priors on the three 1 Publicly available at: https://dx.doi.org/10.17909/t9-k8w7-zk32. 2We use U (, ) to denote a uniform distribution over (, ) and T N ( ,  2 , , ) to denote a truncated normal distribution where  and  are the mean and variance of an untruncated normal distribution, which is then truncated on the left at location  and on the right at location .additional magnification parameters: with T as defined in Equation 10and TΔ as defined in Equation 11.
We compare the results from the two GausSN models to those from sntd presented in P21.While the results from P21 provide a helpful benchmark, these fits assume the correct redshift and the correct SN template, which may not always be well-known for real glSNe.However, the Roman objects were fit by P21 using sntd's "parallel" method, which fits each band separately, because it is less computationally expensive.The "series" method, which fits all bands simultaneously, and the "color" method, which fits color curves simultaneously with the lensing parameters, have been shown to be more robust.For this reason, these two methods are the only two that have been used in the analyses of real glSNe events with sntd (Kelly et al. 2023b;Pierel et al. 2023b).Therefore, the sntd results may be improved by using a different fitting method.

Analysis
On an individual object basis, the GP fits to the data and posterior distributions show that GausSN is effectively and accurately inferring the time delays from the glSNe systems.In Figure 3, we demonstrate the quality of the fit on the light curve level from the sigmoid magnification model.We plot the observed data with draws from the posterior predictive distribution.The bottom panel of the figure shows realizations of the magnification function for each posterior sample.Notably, because we are only able to constrain the relative magnification, the magnification function shows greater uncertainty around the beginning and end of the time series because these regions are where the light curve lacks overlapping data from the two images.
In Figure 4, we show the posteriors of the constant magnification fit.GausSN recovers Δ = 38.44 +1.23  −1.30 days, whcih captures the true time delay of 38.53 days well within the 68% credible interval (CI).As we don't necessarily expect the posteriors to be Gaussian, we calculate the 68% CI from the 16 th and 84 th percentiles of the posterior samples.The magnification is not well recovered, as GausSN finds  = 0.60 ± 0.02, when the true  = 0.24.We attribute this discrepancy to the presence of microlensing from lens substructure.Indeed, there is more uncertainty in the recovered magnification when fitting with the sigmoid magnification model.With this model, GausSN finds  0 = 0.59 +0.05 −0.30 , which is more consistent with the truth.The uncertainty on the time delay is however increased when fitting with the sigmoid magnification model, as expected, with Δ = 37.72 +1.76  −1.97 .This represents a 5.22% precision time-delay estimate with the sigmoid magnification model, up from 3.38% precision with the constant magnification model.Of course, the constant magnification model gives only a statistical uncertainty, so it does not take into consideration systematic uncertainty from microlensing.A separate systematic microlensing error would have to be accounted for in post-processing.Therefore, the time delay estimate with the sigmoid magnification model may, in the end, be more precise.
On a population level, GausSN performs well across all types of SNe and mass models.We note that the true parameters for the SN IIP SC2 and SC4 mass models were not available, so these results are excluded from the reported statistics.The left column of Figure 5 shows the distribution of Δ fit − Δ true for the sigmoid and constant magnification models.We also compare to the results from sntd, which are included in the P21 simulated Roman data release.With the constant magnification model we find 35.32% of time delays within ±1 day, 70.33% within ±3 days, and 83.73% within ±5 days.With the sigmoid magnification model, we find that 27.13% of time delays are recovered within ±1 day, 62.12% are recovered within ±3 days and 79.07%are recovered within ±5 days.For comparison, sntd recovers 31.95%within ±1 day, 66.25% within ±3 days, and 81.28% within ±5 days.
We also find that the distribution of Δ fit − Δ true is generally unbiased.The median of the distribution is -0.168 days for the sigmoid magnification model and -0.249 days for the constant magnification model.The diminished bias of the sigmoid magnification model relative to the constant model reiterates the need to account for microlensing alongside time-delay estimation.Future work on improved methods for marginalizing over a choice of microlensing parameterization to further reduce biases are discussed in §7.
We also consider the fractional error, |Δ fit − Δ true |/Δ true , on the time-delay estimates.This metric gives a better sense of closeness to the truth scaled by the size of the time delay.For the constant magnification model, we find that 32.55% of systems have fractional errors of less than 5% and for the sigmoid magnification model, 26.48% of systems have fractional errors of less than 5%.These results are comparable to the sntd results, which find 30.84% of systems to have a fractional error of less than 5%.
GausSN, therefore, is effective at estimating time delays close in absolute and fractional value to the true time delay.That the sigmoid magnification model is in general further from the truth according to this metric than the constant magnification model or sntd is not unexpected.As this model has more flexibility, often leading to more complex posteriors, the point estimate for the time delay may be further from the truth.The full posteriors give a more accurate sense of the time-delay estimate by a given model.In the right column of Figure 5, we show the CDF of the standardised error, (Δ fit − Δ true )/ Δ distribution for the sigmoid and constant GausSN magnification models, sntd, and a unit normal distribution for reference.If the uncertainties are well-calibrated, we expect that across the sample of SNe, 68% of Δ true will fall within the 68% CI and 95% of Δ true will fall within the 95% CI.Although the full posterior distributions from sntd are not available, we can compare these statistics to the fraction of SNe which have time delays recovered within 1 and 2 of the truth from the sntd point estimate and uncertainty.
For the constant magnification model, we find 54.46% in 68% CI and 81.31% in 95% CI.For the sigmoid magnification model, we find 54.65% of SNe with Δ true in 68% CI and 82.65% in 95% CI.For sntd, 50.92% of SNe have Δ true which fall in the 1 CI and 78.33% fall in the 2 CI.While these statistics do not meet the benchmark set by the normal distribution, they are comparable to, if not a slight improvement on, the results from sntd.Therefore, GausSN provides time-delay estimates that are as close, if not closer, in absolute value to the truth with relatively well-calibrated uncertainties compared A = 0.24 +0.04 0.04 to sntd, the leading time-delay estimation technique for glSNe.As there are some variations for different types of SNe, we report the above statistics broken down by type of SN for all three models in Table 1.
GausSN seems to have an improved calibration of uncertainties compared to sntd, which does not just arise from large uncertainties from GausSN.The mean uncertainty on the time delay is 2.90 days from the constant magnification model -comparable to the mean uncertainty of 3.02 days from sntd -and 3.81 days from the sigmoid magnification model.Furthermore, we can compare the fraction of glSNe with time delays measured to a certain level of precision, or which have | Δ /Δ fit | less than a threshold.We find that with the constant magnification model, 5.10% of fitted Roman objects are measured to 1% precision -or have | Δ /Δ fit | < 0.01 -and 25.50% are measured to 5% precision.For the sigmoid magnification model, 2.88% and 17.87% of objects are measured to 1% and 5% precision, respectively.
Finally, sntd measures 5.94% of glSNe to 1% precision and 25.85% to 5% precision.The median precision is 14.55% for the constant magnification model, 22.44% for the sigmoid magnification model, and 14.57% for sntd.The uncertainties from GausSN are still underestimated in general, though, possibly due to the single microlensing parameterization considered.Future work to reduce systematic uncertainties from microlensing is described in §7.
We note that 2.30% and 4.81% of glSNe have (|Δ fit − Δ true |)/ Δ > 5 for the sigmoid and constant magnification model, respectively.It is unsurprising that the more conservative sigmoid magnification Table 1.The results of GausSN fits to the Roman simulations by type of SN and model used.We compare the GausSN results to those of sntd, a template-based approach.While GausSN is agnostic to SN type and redshift, sntd has assumed the correct template of the correct SN type and redshift.Note that the results from sntd were not available for two mass models (SC4 and SC7) of SNe IIP.However, there is little variation in the results when further broken down by mass model, so we would not anticipate significant changes in the SNe IIP results if these mass models were included.The number of SNe included in total.The SNe fit by GausSN were selected at random from the full sample.While there are not equal numbers of each mass model fit by each model, we weight by the number from each mass model, so there is equal contribution from all mass models.b The fraction of SNe for which Δ true falls within the 68% and 95% credible intervals (CI) computed from the 16 th and 84 th percentiles of the posterior samples.Note that because the full posteriors are not available for the SNTD fits, we use the 1 and 2  intervals as a proxy for this value.
treatment has fewer catastrophic outliers compared to the constant magnification model.However, we still exceed the rate of 5 outliers consistent with a normal distribution.Concerningly, many of the outliers, particularly from the constant magnification fit, have fitted light curves that appear convincing based on visual inspection.It is likely that significant effects from microlensing, which may cause a time-varying magnification that is not well described by a sigmoid function, are at play in these simulations.Testing additional parameterizations of microlensing and using model comparison metrics will be necessary in analyses of real glSNe with significant microlensing effects.That said, these outliers will remain an issue for cosmological analyses of glSNe and uncertainties due to the presence of such outliers should be propagated through further analyses.
There are some additional outliers which arise from multi-modal posteriors on the time delay which are not well-described by the mean and standard deviation of the samples.These are of less concern for two reasons.Firstly, for many objects, a more informative prior on the time delay can resolve this issue.More careful tailoring of priors based on visual inspection of the light curves and other contextual information will be possible with real glSNe, due to their rarity, though it is difficult to do on the large scales needed for this analysis.Secondly, the full posterior, rather than a point estimate, should be used in the cosmological analysis of any glSN.There are also some outliers with very low SNR for one or both images, making any constraint on the time delay difficult.These two types of outliers are also expected at the low rates we see and are not of concern because of the ease with which they can be identified as needing more careful treatment.
Furthermore, such rates of outliers are not inconsistent with existing methods for time-delay estimation.For reference, 5.96% of sntd fits to the Roman simulations have GausSN has fewer outliers is consistent with the statistics reported in Table 1, which show the GausSN uncertainties are better calibrated.In addition, the methods presented in Kelly et al. (2023b) show similar rates of outliers.Therefore, GausSN again meets, if not exceeds, the benchmarks set by existing time-delay estimation methods.
Together, these results demonstrate that GausSN provides competitive time-delay estimates.When comparing to the results from sntd, the GausSN constant magnification model, which has the more similar magnification treatment to sntd, consistently retrieves time delays that are as close, if not closer depending on the type of SN, to the true time delay.The estimates also tend to have better calibrated uncertainties, with minimal differences in the percent of time delays estimated to 1% and 5% precision.Furthermore, GausSN is agnostic to the SN type and redshift, while sntd has assumed the correct template of the correct SN type and redshift to achieve the reported results.Because the Roman simulations were not generated from the GausSN forward model, our results demonstrate the ability of GausSN to generalize.
Including a time-varying magnification model with GausSN, as in the sigmoid magnification model, leads to an increase in the uncertainties on the time delay, as expected given our uncertainty in the effect of microlensing.Unsurprisingly, the fraction of SNe with time delays recovered to 1% and 5% precision is lower for this model with minimal decreases also seen in the time-delay recovery within 1, 3, and 5 days.However, the uncertainties are the best calibrated out of all three models, which makes sense given this model is the most conservative approach to microlensing considered in this work.The performance of this model suggests promise for the ability to constrain a time-varying magnification function with GausSN.We further comment on future work possible with time-varying magnification models in §7.

Constant Magnification Simulations
We apply this method to LSST-like simulations of doubly-imaged SNe.To achieve a realistic survey cadence and depth, we use the simulated survey strategy from the Rubin Operations Simulator (OpSim), which simulates the field selection and image acquisition process of Rubin-LSST over the 10-year duration of the planned survey (Delgado et al. 2014;Naghib et al. 2019).We use OpSimSummary (Biswas et al. 2019(Biswas et al. , 2020) ) to retrieve the times of observations of a specific point on the sky with the baseline v3.0 strategy from Rubin-LSST's Survey Cadence Optimization Committee (Rubin Observatory Survey Cadence Optimization Committee 2023).While baseline v3.3 was recently released, we do not expect the updates to have a significant impact on the results of this analysis, as the biggest change Table 2. Priors on the hyperparameters used when fitting the Rubin-LSST simulations.If a simulation model is specified, with "C" denoting the constant magnification simulations and "A" denoting the Arendse et al. (2023) simulations, then the prior applies to the objects simulated from that model.If a simulation model is not specified, then the prior on that hyperparameter is the same for both sets of simulations.Note that   , t1 , and t2 are defined as in Equations 39, and T and TΔ are defined as in Equations 10 and 11, respectively.

Hyperparameter Sim. Model
Prior has resulted in a decrease in the -band sensitivity and increase in the sensitivity of the other bands.While 10 years of data can be simulated, we choose to clip the data in time to just the time frame covered by the template to avoid artifacts in the light curves that might make fitting for the time delay easier.The objects are observed in as many of Rubin's  as possible at each object's simulated redshift given the template limitations in wavelength and the survey cadence.
The distribution for the absolute luminosities of CC SNe are taken from Grayling et al. (2023).The redshift distributions and absolute magnifications of glSNe as expected to be discovered by Rubin-LSST are approximated from Goldstein et al. (2019).For each simulated object, a redshift, magnification for the first image, and magnification for the second image are randomly assigned from these distributions.The time delays are randomly drawn from a normal distribution with a mean of 0 and standard deviation of 50 days.The distributions of simulated parameters are shown in Figure 2. Simulating realistic microlensing and dust extinction effects is beyond the scope of this work, so we choose not to include these effects in the simulations.We also do not consider correlations between the redshift, magnification, and time delay distributions.These simulations are intended to test a wide range of parameter space to probe the limits of the method.
The simulated objects must pass a series of data cuts to ensure the objects are high enough quality to firstly, produce an alert in the LSST data stream and secondly, fit for a time delay which may be reasonably informative.We require at least two observations of each image and the total number of observations to be greater than two times the number of bands in which there is data.The light curves must also have two data points with signal to noise ratio greater than 5 and at least two data points in any band within ±5 days of the time of peak magnitude for each image.We furthermore impose a requirement that  be less than 15.These quality cuts are conservative, as we expect to need higher quality data than the minimum required in our simulations for a competitive time-delay constraint.
For the Rubin-LSST simulations, we fit the objects using the priors given in Table 2. Again, the prior on the time delay requires that there always be overlap between data from image 1 and data from image 2. Although the simplified Rubin-LSST simulations created for this analysis do not include microlensing, we fit with both the constant magnification model, which is the model the objects were simulated from, and the sigmoid magnification model to test if any biases arise from using an incorrect magnification model to fit the data.Arendse et al. (2023) Recently, Arendse et al. (2023) released the publicly available "lensed Supernova Simulation Tool" (lensedSST3 ), which simulates glSN Ia with a realistic microlensing treatment.The lensedSST pipeline simulates glSNe Ia from the SALT3 model (Guy et al. 2007;Kenworthy et al. 2021) using sncosmo, where the cadence and observing conditions are drawn from OpSimSummary using the baseline v3.0 observing strategy.The redshift of each SN is drawn jointly with the lens redshift and the Einstein radius of the system to produce systems with realistic time delays and magnifications using the Lenstronomy lens modeling package (Birrer & Amara 2018;Birrer et al. 2021).In Figure 2, the distributions of simulated parameters are shown.

Realistic glSN Ia Simulations from
Microlensing, which is caused by the expanding photosphere of the SN, is accounted for using realistic SN Ia explosion models from ARTIS (Kromer & Sim 2009) and microlensing maps from GERLUMPH (Vernardos & Fluke 2014;Vernardos et al. 2014Vernardos et al. , 2015)).Given the microlensing parameters at the position of the images on the microlensing maps from Lenstronomy, a realization of microlensing is drawn for each wavelength filter and applied to the light curves.Finally, the simulated glSNe Ia are subject to selection effects, described in detail in Arendse et al. (2023), to determine whether an object would be detected via several detection pipelines.
The explosion models for other types of SNe are less well-studied, making the simulation of realistic chromatic microlensing for non-SNe Ia difficult.Pierel et al. (2021) assume a simple achromatic model for the expansion of the SN photosphere to produce realistic achromatic microlensing.The simulation of chromatic microlensing, however, is only possible for glSNe Ia.Therefore, in addition to the constant magnification simulations, which probe the sparse data regime for all types of glSNe, we use the lensedSST simulations to test GausSN's performance in the sparse data regime with microlensing for glSNe Ia.In this work, we continue to only focus on simulations of doubly imaged glSNe.
We use the same priors on the hyperparameters as the constant magnification Rubin-LSST simulations, as given in Table 2, with the exception of the time delay.In the lensedSST simulations, the leading image is assigned as image 1, so we restrict the time delay to be positive.In addition, we remove objects for which the two images do not overlap in phase space.This choice is motivated firstly by the need to set some general prior to fit the simulated objects in batches.The chosen prior would exclude the true value of the time delay for these objects, so we remove them.Secondly, GPs are data-driven, so their ability to extrapolate outside of where there is data is limited.If the multiple images of a glSN do not overlap in phase space, the GausSN fit may be unreliable because of this limitation of the model.We comment further on future work with GausSN on very sparsely sampled light curves in §7.As with the constant magnification simulations, we fit with both the constant magnification model and the sigmoid magnification model, both of which assume achromatic microlensing, to test the bias introduced by incorrect assumptions about the underlying microlensing.

Analysis
As with the Roman simulations, we first demonstrate the quality of GausSN's performance on an individual object -a glSN Ia from the constant magnification simulations at  = 1.39 and with Δ = 25.62 and  = 1.30.We demonstrate the quality of the GausSN fit at the data level in Figure 6.As in Figure 3, we plot draws of the posterior predictive distribution over the light curve data.Even with the sparse Rubin-LSST data, GausSN recovers a realistic fitted light curve.
In Figure 7, we show the joint posteriors for the four fitted parameters from the constant magnification fit in a corner plot.GausSN estimates Δ = 26.41+3.33  −3.29 days, a 12.61% precision estimate, and  = 1.37 +0.21  −0.17 .Although modeling this object is simpler than modeling the Roman objects, for example, because it is not subject to microlensing effects, the sparsity of the data clearly provides a challenge to the model.Even so, the time delay and magnification are still accurately recovered by GausSN.
On a population level for the constant magnification simulations, GausSN shows effectivenesss in recovering time delays that are close in absolute value to the true time delays.In the top row of Figure 8, we show the distribution of Δ fit − Δ true by type of SN for both the constant magnification fits and sigmoid magnification fits.With the constant magnification model, GausSN retrieves 40.78% of time delays within ±1 day of the truth, 69.20% within ±3 days, and 79.92% within ±5 days.We also find that 40.67% of systems have fractional errors of less than 5%.With the sigmoid magnification model, the performance is slightly reduced, but follow the same trends seen in the results on the Roman simulations.Notably, there does appear to be a slight bias in the median of the distribution of Δ fit − Δ true of -0.417 days that is introduced when fitting with the sigmoid magnification model.As discussed further in §7, it will be important to marginalize over the choice of microlensing parameterization to mitigate any potential bias introduced by a single microlensing function, as may be present in these results.On the whole, though, the GausSN time-A = 0.14 +0.03   delay estimate is consistently close to the truth, with only a slight and anticipated performance decrease when assuming an incorrect magnification model.
As with the Roman simulations, we check that the uncertainties on the time delays are well-calibrated by evaluating the distribution of (Δ fit − Δ true )/ Δ .We plot the CDF of this distribution in the bottom row of Figure 8.The fraction of SNe for which the true time delay is captured within the 68% and 95% CI, calculated from the 16 th and 84 th percentiles of the posterior samples, is 51.88% and 80.57%, respectively, for the constant magnification fits.With the sigmoid magnification model, 50.54% and 78.62% of true time delays are recovered within the 68% and 95% CI.The uncertainties are generally underestimated, as was seen in the results from the Roman simulations.We report these results broken down by type of SN in Table 3.
When fitting the lensedSST simulations, GausSN maintains a high quality performance, despite the greater complexity of these simulations due to their inclusion of realistic microlensing.Slight reductions are seen in the fraction of objects with time delays measured within 1 day -26.70% for the constant magnification fit and 23.23% for the sigmoid magnification fit.However, the fraction of objects with time delays measured within three and five days remain consistent at 58.67% and 72.70% for the constant magnification fit and 54.00% and 70.59% for the sigmoid magnification fit.There is a reduction in the number of objects with fractional errors of less than 5% on the time delay, but such behavior is expected given that microlensing increases uncertainty in time delay estimation.The distribution of (Δ fit − Δ true ) shows a slight bias, with a median of 0.106 days for the constant magnification fits and of -0.223 days for the sigmoid magnification fits.Interestingly, this offset in the median is reduced relative to the results on the constant magnification simulations.
Although the model tends to be further from the truth for the lensedSST, the uncertainties on the time delay accurately reflect the model's ability to recover the true time delay; with the constant magnification model, 61.54% and 88.54% are recovered within the 68% and 95% CI and with the sigmoid magnification model, 65.01% and 91.10% are recovered within the 68% and 95% CI.Further improvements on the calibration of uncertainties can be made through future work discussed in §7.Despite the challenge of sparsely sampled light curves and realistic microlensing effects, GausSN continues to provide precise and accurate time-delay estimates with reliable uncertainties when fitting the lensedSST simulations.
Additionally, with the current observing strategy, Arendse et al. (2023) predict that an appreciable fraction of glSNe discovered with Rubin-LSST will have few to no observations before light-curve peak for one or more images.In Table 4, we summarize the performance on GausSN on the subset of objects for which the peak in the light curve of one or both images is missed.Unsurprisingly, the recovered time delays for this subset of glSNe are on average further from the truth, as indicated by the reduced number of objects with fractional errors in the measured time delay of less than 5%.Similar drops in performance are also seen in the rate of time delays within 1, 3, and 5 days of the true time delay.Importantly, though, the fraction of objects with a true time delay that falls within the 68% and 95% CI remains equally well-calibrated, indicating an appropriate increase in uncertainty when estimating the time delays of these objects.
Only 5.73% and 7.87% of the constant magnification simulations have fitted time delays which are more than 5 away from the true value for the constant magnification and sigmoid magnification models, respectively.Based on visual inspection of the fits, it is very easy to tell when there has been a catastrophic failure with GausSN for the Rubin-LSST simulations.For the majority of these objects, the posteriors show that the nested sampling chains have run up against the prior bounds, either for the GP length scale or the time delay, suggesting that the prior was not well-chosen for the data.Real glSNe will not be fit in bulk, as the simulated objects are for this analysis, so again more informative priors can be selected for objects on an individual basis to mitigate many of these extreme outliers.Indeed, fitting a few of these objects with more informative priors enabled more accurate time-delay recovery.Therefore, these outliers are not of great concern, as they could be easily flagged and refit with more informative priors to get reliable time-delay estimates.
As with the Roman simulations, there are also several outliers with multi-modal posteriors for which the mean and standard deviation of the samples are simply a poor summary of the posterior.Again, fitting with more informative priors easily solves this problem.Furthermore, we re-emphasize that the full posterior should be used in the cosmological analysis of any glSN.These types of outliers are, therefore, expected and not as concerning.
A similarly small fraction of lensedSST objects -2.11% for the constant magnification fit and 0.90% for the sigmoid magnification fit -have fitted time delays which are more than 5 away from the true value.Outliers appear to be dominated by objects that experience extreme microlensing effects, which appear to alter the location of the peak in the light curve.Because the sigmoid microlensing model is not flexible enough to capture all possible effects of microlensing, these outliers are unsurprising.Future work to account for uncertainty in the microlensing model choice is discussed in §7 and would reduce the frequency of these fits with very poorly calibrated uncertainties.Unfortunately, these outliers are difficult to identify based on visual inspection of the fitted light curves.
The results in Table 3 are roughly consistent with those reported for the Roman simulations.However, the mean uncertainty on the time-delay estimate from GausSN for the constant magnification Table 3.The results of GausSN fits to the Rubin-LSST simulations by type of SN.
Sim.The model from which the data was simulated, where "C" denotes the constant magnification model simulations from this work and "A" denotes fully realistic glSN Ia simulations from Arendse et al. (2023).b The GausSN model that was used to fit the data, where again "C" denotes the constant magnification model and "S" denotes the sigmoid magnification model.c The fraction of SNe for which Δ true falls within the 68% and 95% credible intervals (CI) computed from the 16 th and 84 th percentiles of the posterior samples.simulations is 7.53 days for the constant magnification model and 10.22 days for the sigmoid magnification model.For the lensedSST simulations, the mean uncertainty on the time-delay estimate is 5.02 days for the constant magnification model and 5.84 days for the sigmoid magnification model.Relative to the Roman simulations, the time-delay estimates for the Rubin-LSST simulations are more uncertain in general.Given that GausSN is a data driven approach, and the Rubin-LSST data is more sparse than the Roman data, this in-crease in uncertainty is expected.Most importantly, the uncertainties remain equally well-calibrated across two very different sets of data.
Finally, we note that there are differences in the constant magnification and lensedSST simulations beyond the microlensing treatment, which may also have an effect on the results.The lensedSST pipeline retrieves time delays and magnifications from physical modeling of the lens given the source redshift and lens geometry.Therefore, these parameters are highly correlated in these simulations.On the other hand, the constant magnification simulations from this work independently sample from the distributions of redshift, time delay, and magnification.Therefore, the constant magnification simulations may contain parameter combinations that would be unlikely to occur in reality.The regions of parameter space which the lensedSST simulations fall into appear to lend themselves better to fitting in the GausSN model, as seen in the improved calibration of the time delay uncertainties.The distributions of source redshifts, time delays, and relative magnifications for these two sets of simulations are shown in Figure 2.However, the number of objects with fractional errors on the time delay estimate of less than 5% is reduced for the lensedSST simulations, which we attribute to microlensing uncertainty.On the whole, GausSN shows strength in fitting Rubin-LSST simulations, despite the sparse nature of the light curves.

APPLICATION TO SN REFSDAL
We apply GausSN to the five images of SN Refsdal to demonstrate the model's performance on real data, as well as on an object with more than two images.Following the publication of Kelly et al. (2023a) and Kelly et al. (2023b), the data was made public 4 .The five images of SN Refsdal, images S2-S4 and SX, were observed by the Hubble Space Telescope (HST) in three filters -F105W, F125W, and F160W.As the data in the F105W filter is limited and of worse quality relative to the other two bands, we follow Kelly et al. (2023b) and exclude this data from the time-delay fit.
We first fit SN Refsdal with the constant magnification model.The priors on the hyperparameters of this fit are given in Table 5.We note that this analysis is not blinded in the way that the analysis in Kelly et al. (2023b) was.The corner plot for the fit is shown in Figure 9.We measure the time delay of the fifth image, SX, which is the most cosmologically interesting due to its long time delay, to be Δ 1, = 377.98+5.10  −4.98 days -a 1.35% precision measurement.In addition, our results are in very good agreement with those reported by Kelly et al. (2023b) for all five images.As shown in Figure 9, the time-delay estimates from GausSN fall within the 16 th and 84 th percentiles of the posterior distributions on these parameters given in Kelly et al. (2023b).The fitted light curves are shown as an inset in Figure 9.
As noted in Kelly et al. (2023b), there is evidence for microlensing affecting the light curves, particularly those of images S2, S4, and SX.Therefore, we additionally fit SN Refsdal with the sigmoid magnification function, with the priors on the additional hyperparameters given in Table 5.With the sigmoid magnification function, the time delay for SX, Δ 1, = 376.46+5.37  −5.29 days, is still highly consistent with the previous results.There is only a slight reduction in the precision, 4 Available at https://dx.doi.org/10.3847/1538-4357/ac4ccb.Table 5. Priors on the hyperparameters describing the SN Refsdal system.If an image is specified for a prior, the hyperparameter is taken as relative to image 1.If an image is not specified, then the prior on that hyperparameter is the same for every image.

Hyperparameter
Image Prior with the time delay now measured to 1.43% precision.The time delay of image S3 is also still in very good agreement with Kelly et al. (2023b).

𝐴
For images S2 and S4, the results are more discrepant from those in the constant magnification fit to the data.With the sigmoid magnification model, GausSN estimates Δ 1,2 = −1.95+2.53  −2.54 days and Δ 1,4 = −0.40+5.45  −6.02 days.On the other hand, Kelly et al. (2023b) reports Δ 1,2 = 9.7 +4.4  −3.7 days and Δ 1,4 = 19.44 +8.0 −4.9 days.We emphasize, though, that the results from the individual methods used in the original analysis show additional dispersion in their time-delay estimates that is not well-represented in the values reported above.Therefore, the discrepancy seen in the Kelly et al. (2023b) time-delay estimates and the time-delay estimates reported from the sigmoid magnification model fit in this analysis is not unexpected or of great concern.Most importantly, the results for image SX remain consistent.
We show realizations of the sigmoid magnification function using draws from the posteriors, as shown in Figure 10.As in Kelly et al. (2023b), we find considerable uncertainty about the effect of microlensing on image SX.There appears to be evidence for microlensing affecting each of images S2-S4, as well.At early times, there is uncertainty in the relative magnification because of the lack of overlapping data from multiple images.However, even around peak, where there is data from other images, the microlensing effect is uncertain.For image S3, the microlensing effect switches on in the tail, where the light curve appears to begin rising again.As this behavior is not seen in the tails of image S2 or S4, it is attributed to microlensing.The magnification from microlensing does not change the shape of the light curve around peak, so it has a minimal effect on the time delay.
Images S2 and S4 show signs of microlensing affecting the beginnings of the light curves.This behavior suggests that images S2 and S4 may exhibit a faster rise than the other images.Interestingly, the time delays of these two images are more discrepant with the Kelly et al. (2023b) results.The rise of the light curves is well-observed for all five images, so we expect that a considerable amount of time-delay constraining power comes from this region of the light curves.Therefore, microlensing affecting images S2 and S4 in this time frame may result in the more significant shift in the time delay from the constant magnification model to the sigmoid magnification model.There is clearly uncertainty in the time-delay fit due to microlensing effects on the light curves, highlighting the importance of accounting for microlensing alongside the time-delay estimate.For comparison, we show the 68% CI for the time delays from the 16 th and 84 th percentiles reported in Kelly et al. (2023b) as the shaded red region, with the reported point estimate shown by the red line.We do not show the shaded region for the magnifications because they cover all of the parameter space shown in each 1D posterior plot.We use the 68% CI as a proxy for the full posterior distributions from Kelly et al. (2023b), which is not publicly available.The time delays from GausSN are all in good agreement with Kelly et al. (2023b).Inset: The flux data and fitted light curves of SN Refsdal's five images from the constant magnification model.As described in detail in §4.2, the red fitted curves are time-shifted and magnified copies of the blue fitted curves, which have been conditioned on the data from both images.

DISCUSSION
There now exist several methods for extracting time delays from lensed quasars and glSNe.A range of time-delay estimation methods have been used in analyses of lensed quasars.A number of optimization schemes have been paired with spline and GP models of the underlying light curve shape, due to the stochastic nature of the quasar light curves.Among these methods is the widely popular spline-based method PyCS (Tewes et al. 2013;Millon et al. 2020), which was used in the H0LiCOW 2.4% precision  0 estimate (Wong et al. 2020) and adapted for use on SN Refsdal (Kelly et al. 2023b).On the other hand, template-based approaches, such as sntd, have The realizations of the fitted light curves are created according to the procedure described in §4.2.We find that there is considerable uncertainty regarding the effect of microlensing on image SX, as also found in the analysis from Kelly et al. (2023b).GausSN more confidently constraints the relative magnification affecting images S2-S4.There is evidence for microlensing affecting the beginnings of the light curves for images S2 and S4 and the tail of the light curve for image S3.
been favored in recent analyses of glSNe (Pierel et al. 2023b;Kelly et al. 2023b).We consider below a few advantages to the GausSN approach, which motivate its use alongside other methods in analyses of glSNe discovered in the future, as well as point out where improvements on GausSN can be made.

Template Choice Uncertainty
SNe are known to have some level of intrinsic variation within each type and a number of methods have been used to extract templates from SN light curves.As a result, there are several templates available for each SN types which one may chose to fit with using a templatebased approach for time-delay estimation.However, the effect of fitting for a time delay with an incorrectly chosen template has yet to be explored.This unknown systematic contributes unaccounted for uncertainties to the time-delay estimate, which may be difficult to propagate through to the uncertainty in a motivated way.Even for types of SNe which exhibit a great diversity of light curves, which may be difficult to fit with a template-based approach, GausSN maintains optimal performance.In addition, template-based approaches depend on reliable redshift estimates to shift the template into the correct frame of reference.In the era of Rubin-LSST and Roman, spectroscopic follow-up will be limited, so spectroscopic redshifts will likely not be available for every object.Therefore, GausSN, a method which is agnostic to source redshift, is an important complement to methods which depend on redshifts.
Furthermore, lensing enables us to probe SNe at higher redshifts than would otherwise be observable by existing telescopes.The redshift evolution of SNe is an open question in the field, with an unaccounted for evolution potentially contributing to systematics in cosmological analyses (Nicolas et al. 2021).Therefore, the use of templates based on SNe light curves in the nearby universe similarly introduces a potential source of unaccounted for systematic uncertainty in the time-delay estimate.Especially because templates tend to be based on data from very-nearby SNe to ensure it is high quality, the extremity of this difference will be at its peak for glSNe expected to be discovered by Rubin-LSST, up to  = 2, and by Roman, up to  = 4.
Being completely independent of templates, GausSN is not subject to this systematic and therefore provides complete Bayesian uncertainties on time-delay estimates without further post-processing to account for such an uncertainty.Given glSNe are attractive for having minimal systematics relative to the distance ladder, keeping sources of systematics to a minimum is important in these analyses.
Finally, if the SN light curve shape is peculiar or without a good template, as was the case with SN Refsdal, template-based approaches must turn to generic templates, such as the Bazin function (Bazin et al. 2009;Karpenka et al. 2013;Villar et al. 2019).The flexibility of generic SN templates can make fitting the templates to the data tedious and sensitive to fine tuning.Alternatively, a custom template may be created for the system.While this approach was possible for SN Refsdal (Kelly et al. 2023b), it will not be feasible in the future as larger numbers of glSNe are discovered.Furthermore, it is difficult to validate an analysis which is based off a custom-made template.GausSN is capable of fitting any object, regardless of its spectroscopic classification, without constructing or fine-tuning a template for the data beforehand.
However, in the very sparse data regime, a template-based approach may be preferable if it is possible to make stronger assumptions about the shape of the underlying light curve to better constrain the time delay.While GausSN performs well on Roman simulations and Rubin-LSST simulations, we cannot predict performance on only a few epochs of data per image, especially if these observations have little overlap in SN light curve phase, because of the limited ability of a GP to extrapolate between sparse observations.As demonstrated in Table 4, the performance of GausSN is reduced for objects in which pre-peak data is missing for one or both images.In future work, we plan to incorporate SN light curve templates as the mean function of the GP.This functionality would provide a more informative prior over the functions describing the light curve, while still allowing for residual variations which may arise from incorrect template choice or a redshift evolution of SN light curves.

Microlensing
Another advantage to GausSN is the treatment of microlensing.The wide ranging nature of the effect which microlensing from stars and substructure can have on a light curve is extremely difficult to constrain (Pierel et al. 2023b).For this reason, existing time-delay techniques consider only an additional uncertainty due to microlensing rather than treating it alongside the time-delay estimate.By treating macro-and micro-lensing together as one time-varying magnification term, GausSN marginalizes over only the relative magnification realizations which are consistent with the data.Therefore, the uncertainties from microlensing are propagated through the measurement in a fully Bayesian way.
We adopt the sigmoid magnification model because it resembles the microlensing curves shown in Foxley-Marrable et al. ( 2018) and Pierel et al. (2023b) and is a convenient fitting function which can be interpreted in a qualitative sense.However, this parametric form is not directly derived from physical models of the microlensing of an expanding photosphere.Furthermore, this treatment does not take advantage of the knowledge of the microlensing magnification distribution.As the photosphere expands, the microlensing effect tends to become less strong because the microcaustics are being averaged over a larger area.Additional parametric models which take advantage of this information can and should be tested within the GausSN framework.
It would also be possible to implement a more flexible and expressive non-parametric microlensing representation within GausSN.Other time-delay estimation methods have used non-parameteric approaches such as splines (Tewes et al. 2013) or simulations from a Gaussian process (Pierel & Rodney 2019) to quantify the effect of microlensing.Within GausSN such a representation could be incorporated seamlessly into the Bayesian model.A natural kernel choice in a GP representation of microlensing could be the non-stationary Gibbs (1997) kernel (see Revsbech et al. 2018, for an example).This would allow for microlensing functions with a short burst of fast magnification change (corresponding to the moment of caustic crossing) and a tendency towards something flatter away from this.
However, any one parameterization of microlensing will not be able to capture the full diversity of microlensing curves possible.In future work, we intend to explore different parametric and nonparametric magnification models, including models which do not assume achromatic microlensing.In addition, we plan to develop methods for Bayesian model averaging (e.g.Nixon et al. 2023, for a recent application in astronomy) to marginalize over multiple possible microlensing models when estimating time delays.

Fully Bayesian Treatment
The GausSN framework provides a fully Bayesian treatment of timedelay estimation, but this is just one component of the  0 estimate.A model of the gravitational potential of the lens is necessary to estimate a time-delay distance, which can be used to constrain  0 (see e.g.Treu et al. 2022;Suyu et al. 2023).Of course, the time delay and lensing potential are intrinsically linked, meaning there is shared information between these two methods that is not taken advantage of by considering them separately.In the future, a fully Bayesian estimate of the time-delay distance from a combined lens modeling and time-delay estimation method would provide the tightest possible constraints on an  0 estimate.
The most outstanding progress in that direction thus far in the field came from Birrer et al. (2020), which used a hierarchical Bayesian model to jointly infer  0 and galaxy density profiles.GausSN represents a step towards this goal, as well.The GausSN posteriors over the time delay and time-varying magnification parameters could be incorporated into a larger hierarchical Bayesian model for  0 to best utilize all the information available.

CONCLUSION
We have demonstrated that precise and reliable time-delay estimates can be obtained for resolved glSN systems using GausSN, a Bayesian semi-parametric Gaussian Process model.This approach has a number of complementary attributes to existing approaches.In particular, GausSN: • easily fits glSNe of any type, regardless of whether it has a well-understood light curve shape, without fine-tuning or redshift information, • is not sensitive to systematic effects introduced by template selection, such as a potential redshift evolution of SNe light curves or fitting with a template that is not representative of the true underlying light curve shape, • accounts for microlensing alongside time-delay fitting with a modular implementation, which enables any number of microlensing models to be used, • provides fully Bayesian uncertainties for time-delay estimates.
GausSN has proven successful on simulations of glSNe as expected from Rubin-LSST and Roman.The Roman simulations from P21 are highly realistic, with careful treatment of survey cadence and depth, dust extinction, and microlensing.Without assuming that the SN type is known and a valid template is available, and without knowledge of the redshift, GausSN maintains a similar performance to the template-based approach sntd.While the results from sntd are based on the assumption of the correct SN-type, template and redshift, GausSN is entirely agnostic to this information.Furthermore, these simulations were not generated from the GausSN forward model, so our results demonstrate the ability of GausSN to generalize.
The simplified Rubin-LSST simulations, created for this analysis using templates from sncosmo, emulate the cadence and depth of the survey.While these simulations do not take into account effects such as dust extinction and microlensing, they test GausSN's performance on sparsely sampled data.GausSN is successful in predicting time delays which are very close to the true time delay, but the uncertainties on the fits are underestimated in the case of some catastrophic failures.Given the current and projected rates of discovery of glSNe, we expect to have the capacity to inspect the fits individually and choose object-specific priors, which we expect will mitigate any highly inaccurate outliers.
We finally demonstrate the performance of GausSN on the five images of SN Refsdal.Using GausSN, we find that the time delay of image SX is highly consistent with the results from Kelly et al. (2023b), regardless of whether a constant or sigmoid magnification model is used.The time delays of images S2-S4 show more dispersion around the values estimated in Kelly et al. (2023b), but not so much to be concerning, especially given the lessened cosmological interest in these images.As in the original analysis, there is considerable uncertainty in the effect of microlensing on image SX.However, even with the uncertainty due to microlensing accounted for, we are still able to estimate the time delay to 1.43% precisiona level of precision competitive with the estimates reported in Kelly et al. (2023b).GausSN, by fitting for the time delay alongside a time-varying magnification term using the data from all five images simultaneously, provides competitive time-delay constraints using real glSN data, without compromising on fit quality.
There remain many possibilities for extensions of this work to address i) restrictions on the kernel which reflect known physical characteristics of the system, such as constraining the kernel to positive flux only, ii) a test of additional treatments of microlensing, which ideally would be based on physically-derived parameters, iii) the incorporation of additional spectroscopic or unresolved light curve information, iv) a treatment of the covariance between light curve shapes in different bands, as we know photometric data to be integrals of slices of a continuous SED in time-and wavelength-space, v) incorporation of dust effects from the host and lens into the model, and vi) development of quantitative model checking approaches to assess the reliability of time-delay estimates.We leave addressing these issues to future work.
A percent-level  0 measurement will require more than precise and accurate time-delay estimation.Optimistically, even if Rubin-LSST and Roman on their own provide the necessary photometric data for time day estimation measurements, additional data will be needed from other instruments to enable the complex lens modeling required for an  0 measurement.Fast and synergistic follow-up of glSNe from a variety of instruments will be needed to complete the full cosmological analysis of these objects.
Despite these challenges, the recent increase in the rate of discovery of glSNe and first measurement of  0 from SN Refsdal suggests promise for this young technique.A measurement of  0 to percent level precision from glSNe would be an important new piece of evidence for efforts to determine whether the Hubble tension is driven by systematics or new physics.With upcoming instruments, such as Rubin-LSST and Roman, and precise statistical analysis tools for time-delay extraction, such as GausSN, glSNe show promise to reach this benchmark in the coming decade.The joint posteriors for this fit are shown in Figure A3.The time delay is estimated to be Δ = −0.21+0.40 −0.51 , which is both closer to the truth than the constant magnification fit and captures the truth within the 1 confidence interval.The estimated uncertainties are larger, but this more appropriately reflects our level of uncertainty due to the difficulty of constraining time-varying magnification from the light curves of the images alone.
Furthermore, we can use this example to check that the magnification curve is well-recovered by GausSN.The inset plot in Figure A3 shows 1,000 realizations of the magnification function based on draws from the posterior distribution of the sigmoid magnification fit to the data.The true magnification function, over-plotted in the figure, is consistent with the fitted magnification curves, demonstrating qualitatively the ability of GausSN to recover the true underlying relative magnification.

APPENDIX B: FLEXIBLE SINUSOIDAL MICROLENSING PARAMETERIZATION
A single parameterization of microlensing will not be able to capture the full diversity of effects microlensing can have on glSNe.Therefore, it will be important to test multiple parameterizations of microlensing and marginalize over this choice with Bayesian model averaging to minimize systematics.In this appendix, we provide a brief exploration of a flexible sinusoidal microlensing parameterization.We define the time-varying magnification function as: where  0 controls the macro-scale magnification,  1 controls the amplitude of the time-varying microlensing magnification effect, T controls the timescale of the microlensing effect, and  0 controls the centering of the microlensing effect.We adopt the following priors on the 7 fitted parameters of the system: where T is given by Equation 10, TΔ is given by Equation 11, and   , t1 , and t2 are defined as in Equations 39.While the microlensing signal is not expected to be sinusoidal, the smooth and well-behaved curves produced by this function in the regime of timescales of greater than 20 days make them suitable for this application.When allowing for a time-varying magnification effect, the time-delay recovery is less biased than when using a constant magnification model (see Figure A2).The additional flexibility of a time-varying magnification treatment leads to a more accurate posterior distribution for Δ, which better reflects the uncertainty in the estimate due to microlensing.Inset: Realizations of the magnification function based on 1,000 draws from the joint posterior for the time delay and magnification parameters, with the true underlying magnification function over-plotted in red.The fitted magnification function is consistent with the true function, demonstrating GausSN's ability to constrain the shape of the time-varying magnification acting on the second image.
We present the results on a subset of the Rubin-LSST simulations in Table B1.The results are consistent with those reported for the other magnification models tested, indicating that no strong biases are introduced by this choice of magnification model.With the flexible sinusoidal model, the true time delays fall within the 68% and 95% CI of the posteriors for NN% and NN% of objects simulated from the constant magnification model, respectively.For the lensedSST simulations, NN% and NN% of objects have true time delays that fall within the 68% and 95% CI of the posteriors.These values are consistent with what is seen in results from the sigmoid magnification model.Furthermore, the median of the distribution of Δ fit − Δ true is -

Figure 1 .
Figure 1.The effect that a 10% increase in brightness due to a time-varying magnification term over 5 days can have on the shape of a SN light curve.The upper panel shows the -band flux of an unlensed SN Ia light curve from the Hsiao et al. (2007) template (blue, dashed line) and a microlensed version of the same underlying light curve (red, solid line).We emphasize that the true time delay is Δ = 0, although the peak in the microlensed light curve appears to be later because of the time-varying magnification term.The bottom panel shows the time-varying relative magnification affecting the second image,  ( ).

Figure 2 .
Figure 2. The distributions of  source , Δ, and  for the the three simulated glSNe samples.

FluxFigure 3 .
Figure 3.An example Roman-like simulated glSN system fitted with the sigmoid magnification model.The first image is shown in blue and the second in red.This object is at  = 1.93 and has time delay, Δ = 38.53days and relative magnification,  = 0.24.At this redshift, the Roman  -, -, and bands roughly cover 3000 − 6000 Å in the rest frame, or the -, -, and -band wavelength regimes.The top three panels show the fitted flux data and GP fit and the bottom panel of the figure shows the fitted sigmoid magnification, ( ), plotted according to the procedure described in §4.2.The red fitted curves are time-shifted and magnified copies of the blue fitted curves, which have been conditioned on the data from both images.

Figure 4 .
Figure 4.The corner plot for the constant magnification fit to the glSN data shown in Figure3.The constant magnification model has four parameters: the kernel amplitude, , the kernel length scale, , the time delay, Δ, and the relative magnification, .The three dashed lines in each of the 1D posterior distribution plots shows the 16th, 50th, and 84th percentiles from left to right, respectively, and the red line shows the true time delay and relative magnification.The time delay is well recovered and to 3.38% precision.The magnification is not as well recovered, potentially indicating additional magnification effects from microlensing that are not captured by the constant magnification model.

Figure 5 .
Figure 5. Left: The difference between the measured time delay and the true time delay for the Roman simulations by type of SN.Right: The CDF of the standardised error, (Δ fit − Δ true )/ Δ , distribution by type of SN.If the error bars are well calibrated, we expect this distribution to resemble that of a Gaussian CDF, shown by the gray line.The solid colored line shows the results from the GausSN sigmoid magnification model, the dashed colored line shows that of the GausSN constant magnification model, and the solid black line shows the results from sntd.

Figure 6 .
Figure6.The flux data and fitted light curve for a simulated Rubin-LSST doubly imaged glSN using the constant magnification GausSN model.This object is a glSN Ia at  = 1.39 with Δ = 25.62 and  = 1.30.As described in §4.2, the red fitted curves are time-shifted and magnified copies of the blue fitted curves, which have been conditioned on the data from both images.

Figure 7 .
Figure7.The corner plot for the fit to an example Rubin-LSSTlike light curve at  = 1.39.From left to right, the parameters fit are , the kernel amplitude, , the kernel length scale, Δ, the time delay, and , the magnification.Δ and  are retrieved within 1 of the truth.

Figure 8 .
Figure 8. Top: The distribution of Δ fit − Δ true by type of SN for the Rubin-LSST simulations.The solid lines show the results from the constant magnification model and the dashed lines show the results from the sigmoid magnification model.The results from fitting the constant magnification simulations are shown in color and the results from fitting the Arendse et al. (2023) glSN Ia simulations are shown in black.We take the mean of the posterior distribution for Δ to be Δ fit .Bottom: The CDF of the standardized error, (Δ fit − Δ true )/ Δ , distribution for the Rubin-LSST simulations by type of SN.The CDF of a normal distribution is shown in gray for reference.

Figure 9 .
Figure9.The posterior distribution from GausSN over parameters from the constant magnification model fit to SN Refsdal's five images.For comparison, we show the 68% CI for the time delays from the 16 th and 84 th percentiles reported inKelly et al. (2023b) as the shaded red region, with the reported point estimate shown by the red line.We do not show the shaded region for the magnifications because they cover all of the parameter space shown in each 1D posterior plot.We use the 68% CI as a proxy for the full posterior distributions fromKelly et al. (2023b), which is not publicly available.The time delays from GausSN are all in good agreement withKelly et al. (2023b).Inset: The flux data and fitted light curves of SN Refsdal's five images from the constant magnification model.As described in detail in §4.2, the red fitted curves are time-shifted and magnified copies of the blue fitted curves, which have been conditioned on the data from both images.

Figure 10 .
Figure10.The sigmoid magnification fit to SN Refsdal's images S2-SX in the F125W filter (top row) and the F160W filter (middle row) with realizations of the magnification function for each of images S2-SX for 100 draws from the posterior (bottom row).The realizations of the fitted light curves are created according to the procedure described in §4.2.We find that there is considerable uncertainty regarding the effect of microlensing on image SX, as also found in the analysis fromKelly et al. (2023b).GausSN more confidently constraints the relative magnification affecting images S2-S4.There is evidence for microlensing affecting the beginnings of the light curves for images S2 and S4 and the tail of the light curve for image S3.

Figure A1 .
Figure A1.The example data drawn from the underlying light curves shown in the upper panel of Figure 1, where the second image is subject to a time varying magnification function shown in the bottom panel of Figure 1.The light curves are "observed" with a 2 day cadence in Rubin-LSST's  filters.The true time delay is Δ = 0.

Figure A2 .
Figure A2.The joint posteriors for the four parameters in the constant magnification model, with the true time delay shown in red.The recovered time delay of Δ = 0.40 ± 0.10 is 4 from the truth.This bias results from the presence of a time-varying magnification function, which the constant magnification model cannot accommodate.

Figure A3 .
Figure A3.The joint posterior for the mock data shown in FigureA1for the seven parameters from the sigmoid magnification model.The true values of the time delay and magnification parameters are shown in red.When allowing for a time-varying magnification effect, the time-delay recovery is less biased than when using a constant magnification model (see FigureA2).The additional flexibility of a time-varying magnification treatment leads to a more accurate posterior distribution for Δ, which better reflects the uncertainty in the estimate due to microlensing.Inset: Realizations of the magnification function based on 1,000 draws from the joint posterior for the time delay and magnification parameters, with the true underlying magnification function over-plotted in red.The fitted magnification function is consistent with the true function, demonstrating GausSN's ability to constrain the shape of the time-varying magnification acting on the second image.
Vincenzi et al. (2019)ted using templates from sncosmo(Barbary et al. 2023).Included in the sample of simulated objects are SNe Ia, SNe Ib, SNe Ic, SNe IIb, SNe IIn, and SNe IIP.We use theHsiao et al. (2007)template to simulate SNe Ia, thePierel et al. (2018)templates to simulate SNe IIP, and theVincenzi et al. (2019)templates to simulate SNe Ib, SNe Ic, SNe IIb, and SNe IIn.Note that some templates are excluded due to the presence of unphysical artifacts in the light curves or poor NIR fits.

Table 4 .
The results from GausSN on the subset of Rubin-LSST data for which there are no observations before peak for one or both images.