Detecting strongly-lensed type Ia supernovae with LSST

Strongly-lensed supernovae are rare and valuable probes of cosmology and astrophysics. Upcoming wide-field time-domain surveys, such as the Vera C. Rubin Observatory's Legacy Survey of Space and Time (LSST), are expected to discover an order-of-magnitude more lensed supernovae than have previously been observed. In this work, we investigate the cosmological prospects of lensed type Ia supernovae (SNIa) in LSST by quantifying the expected annual number of detections, the impact of stellar microlensing, follow-up feasibility, and how to best separate lensed and unlensed SNIa. We simulate SNIa lensed by galaxies, using the current LSST baseline v3.0 cadence, and find an expected number of 44 lensed SNIa detections per year. Microlensing effects by stars in the lens galaxy are predicted to lower the lensed SNIa detections by $\sim 8 \%$. The lensed events can be separated from the unlensed ones by jointly considering their colours and peak magnitudes. We define a `gold sample' of $\sim 10$ lensed SNIa per year with time delay $>10$ days, $>5$ detections before light-curve peak, and sufficiently bright ($m_i<22.5$ mag) for follow-up observations. In three years of LSST operations, such a sample is expected to yield a $1.5\%$ measurement of the Hubble constant.


INTRODUCTION
When a supernova (SN) is positioned behind a massive galaxy or cluster, it can be gravitationally lensed to form multiple images.Such an event is a rare phenomenon that can give valuable insights into high-redshift SN physics, substructures in massive galaxies, and the cosmic expansion rate.An absolute distance measurement between the observer, lens and SN can be made using the arrival time delays between the appearance of the multiple images, combined with a model for the gravitational potential of the lens galaxy and line-of-sight structures (Refsdal 1964).This distance measure can be converted into the Hubble constant ( 0 ) -the present-day expansion rate of the Universe.The exact value of the Hubble constant is an unresolved question, with different techniques yielding different results.Measurements from the cosmic microwave background (CMB) radiation (Planck Collaboration et al. 2020) are in 5 tension with local observations from Cepheids and type Ia supernovae (SNIa) (Riess et al. 2021) and from gravitationally-lensed quasars (Wong et al. 2020).It is worth noting that several other local measurements agree with the CMB results, such as the Tip of the Red Giant Branch (TRGB), as measured by the Carnegie-Chicago Hubble Project (Freedman 2021), and the analysis of seven lensed quasars with less restrictive mass model priors by the TDCOSMO collaboration (Birrer et al. 2020).
Strongly-lensed SNe are promising probes for obtaining a measurement of the Hubble constant which is independent of the distance ladder and early Universe physics (Treu & Marshall 2016;Suyu et al. 2023).Lensed SNe are complementary to lensed quasars and have several advantages for time-domain cosmography.(1) Since they are found at lower redshifts, lensed SNe have generally shorter time delays than lensed quasars, and hence, require shorter observing campaigns.(2) While quasar light curves are stochastic, the light curves of SNe are more predictable and there are accurate light curve templates for several subtypes, like SNIa.This makes it easier to constrain their time delays with sparse observations.(3) When the SN fades away, follow-up observations of the lensed host galaxy light and stellar kinematics of the lensing galaxy can help to constrain the lens mass model and minimise the mass-sheet degeneracy (Falco et al. 1985;Gorenstein et al. 1988;Kolatt & Bartelmann 1998;Saha 2000;Oguri & Kawano 2003;Schneider & Sluse 2013;Birrer et al. 2021a) -a transformation of the lens potential and source plane coordinates that leaves the lensing observables unchanged.Obtaining detailed observations of the lens and host galaxy for lensed quasars is more challenging, as they are always on and often outshine the galaxies.( 4) If the lensed SN is a type Ia, their standardizable-candle nature makes them easier to identify when gravitationally magnified.Knowledge of their intrinsic brightness also helps to minimise the mass-sheet degeneracy, although recent studies show that microlensing effects from stars in the lensing galaxy complicate this by adding stochastic uncertainty to the SNIa luminosities (Dobler & Keeton 2006;Yahalomi et al. 2017;Foxley-Marrable et al. 2018;Weisenbach et al. 2021Weisenbach et al. , 2024)).
While the potential of lensed SNe for time-delay cosmography has been known for decades, these events have proved extremely difficult to detect.To date, only eight multiply-imaged lensed SNe have been discovered.Six of them were found in galaxy clusters: 'SN Refsdal' (Kelly et al. 2015), 'SN Requiem' (Rodney et al. 2021), 'AT2022riv' (Kelly et al. 2022), 'C22' (Chen et al. 2022), 'SN H0pe' (Frye et al. 2024), and SN Encore (Pierel et al. 2023).Furthermore, two SNIa have been discovered lensed by a single elliptical galaxy: 'iPTF16geu' (Goobar et al. 2017) and 'SN Zwicky' (Goobar et al. 2023).SN Refsdal and SN H0pe yielded Hubble constant measurements with uncertainties of ∼ 6% (Kelly et al. 2015) and ∼ 9% (Pascale et al. 2024), respectively, where the precision was primarily limited by the cluster mass models.While SN Encore is expected to enable an  0 measurement with ∼ 10% precision, the other lensed SNe had either insufficient data or too-short time delays to usefully constrain the Hubble constant.
The study of lensed SNe is currently at a turning point, as we will go from a handful of present discoveries to an order-of-magnitude increase with the next generation of telescopes such as the Vera C. Rubin Observatory (Ivezić et al. 2019) and the Nancy Grace Roman Space Telescope (e.g., Pierel et al. 2020).In particular, the Legacy Survey of Space and Time (LSST) to be conducted at the Rubin Observatory is predicted to discover several hundreds of lensed SNe per year according to studies based on limiting magnitude cuts (Wojtak et al. 2019;Goldstein & Nugent 2017;Oguri & Marshall 2010b).
In this work, we take into account the full LSST baseline v3.0 observing strategy to quantify the expected number of lensed SNIa detections per year and investigate the properties of the predicted sample.We focus on lensed SNIa because they are especially promising for cosmology by virtue of their standardizable-candle nature and predictable light curves.We examine the colours and apparent magnitudes of the simulated lensed SNIa sample and study how to best separate them from unlensed SNIa.Finally, we measure time delays from simulated light curves with only LSST data and show how to construct a 'gold sample' which is promising for follow-up observations.Our results include the effects of microlensing due to stars in the lens galaxy.Throughout this work, we assume a standard flat ΛCDM model with  0 = 67.8km s −1 Mpc −1 and Ω m = 0.308 (Planck Collaboration et al. 2016).
We outline our lensed SNIa simulation procedure in Sec. 2, the LSST observing strategy in Sec. 3, and the different aspects of detecting lensed SNIa in Sec. 4. Our results are presented in Sec. 5 and our discussion and conclusions in Sec. 6.

SIMULATING LENSED SNIA
In this work, we develop the publicly available Python code called lensed Supernova Simulator Tool (lensedSST 1 ), which simulates a sample of lensed SNIa light curves with the LSST baseline v3.0 observing strategy to perform our analysis.Additionally, we simulate a sample of unlensed SNIa as a "background" population, to compare their colours to lensed SNIa.Ancillary catalogue information such as the Einstein radius, SN image positions, time delays between images and magnifications of the lens systems are also saved and used in our work.In this section, we describe our assumptions in terms of the lens galaxy mass model, SNIa light curves, and microlensing simulations.

Lens galaxy mass profile assumptions
We employ the multi-purpose lens modelling package Lenstronomy2 (Birrer & Amara 2018;Birrer et al. 2021b) to generate lens galaxies with a power-law elliptical mass distribution (PEMD) to describe the projected surface mass density, or convergence : where  lens is the projected axis ratio of the lens,  lens corresponds to the logarithmic slope, and   denotes the Einstein radius.The coordinates (, ) are centred on the position of the lens centre, and rotated by the lens orientation angle  lens , such that the -axis is aligned with the major axis of the lens.We model the external shear from the line-of-sight structures with a shear modulus  ext and a shear angle  ext , which we assume to be uncorrelated from the lensing galaxy orientation.The adopted parameter distributions are given in Table 1.

Redshift distribution
Most of our input parameters are uncorrelated, except for the Einstein radius  E , lens redshift  lens and source redshift  src , which we sample from a kernel density estimation probability model fitted to ( E ,  lens ,  lens ) of gravitationally lensed SNIa generated in Wojtak et al. (2019).The aforementioned simulation assumes a population of lens galaxies with the velocity dispersion function derived from the Sloan Digital Sky Survey observations (Choi et al. 2007) and a model of the volumetric rate of SNIa fitted to measurements of the SNIa rate as a function of redshift (Rodney et al. 2014).The lensed SNIa rate is determined by drawing random realizations of lens galaxies and supernovae in a light cone and counting the events where strong lensing occurs.More details can be found in Wojtak et al. (2019).We impose an additional upper limit on the source redshift of  src < 1.5 to ensure that the SNe are not redshifted out of the filters.The resulting combinations of  lens ,  src and  E values are depicted in Fig. 1 and the projected 1D distributions can be found in Table 1.
For the background population of unlensed SNIa, we assume a volumetric redshift rate of (Dilday et al. 2008)   () = 2.5 • 10 −5 (1 + ) 1.5 Mpc −3 yr −1 . (2) Our sample of detected unlensed SNIa consists of objects with -band magnitude brighter than 24.0.For lensed SNIa we adopt stricter detection criteria, which are described in Sec.4.1.The resulting redshift distributions of detected lensed and unlensed SNIa and lens galaxies in LSST are displayed in Fig. 2. The Rubin Observatory will be able to discover unlensed SNe at redshifts up to  ∼ 1, resulting in a significant overlap in redshift space between lensed and unlensed SNIa.

SNIa light curves
We model the SNIa as point sources, using synthetic light curves in the observer frame for their variability.The light curves are simulated using SNCosmo 3 (Barbary et al. 2016) and its in-built parametric light curve model SALT3 (Guy et al. 2007;Kenworthy et al. 2021), which takes as input an amplitude parameter  0 , stretch parameter  1 , and a colour parameter .We sample the  1 and  parameters from asymmetric Gaussian distributions that have been derived by Scolnic & Kessler (2016) for the Supernova Legacy Survey (Guy et al. 2010), the Sloan Digital Sky Survey (Sako et al. 2018), Pan-STARRS1 (Rest et al. 2014), and several low-redshift surveys (Hicken et al. 2012; 3 https://sncosmo.readthedocs.io/en/stable/z lens = 0.34 +0.18   Stritzinger et al. 2011).We compute the distance modulus  of each SNIa, based on its  1 and  parameters: Here,  0 is the expected absolute magnitude of a SNIa with  1 =  = 0 and  is the SN flux normalisation in magnitude units.We assume  0 = −19.43 in the -band, corresponding to a universe with Hubble constant  0 = 67.8km s −1 Mpc −1 . and  are the linear stretch and colour correction coefficients, as first found in Phillips (1993) and Tripp (1998) respectively, which specify the correlation of absolute magnitude with the stretch and colour parameters.We assume  = 0.14 and  = 3.1 (Scolnic & Kessler 2016).The resulting absolute magnitude values for each SNIa are used as input for SNCosmo to generate the corresponding unlensed light curves.
The final, lensed light curves are computed in the following way.
After drawing  lens ,  src and  E from the joint distribution from Wojtak et al. (2019) (Fig. 1), we sample random positions of the source and the remaining lens parameters from the distributions given in Table 1 until we find a system that is detectable by LSST.The criteria for what we consider as a 'detection' are described in Sec.4.1 for lensed SNIa.Using Lenstronomy, we obtain the image positions, time delays and magnifications for each lensed SNIa.The time delays and magnifications are applied to the unlensed SNIa light curve to obtain the final, lensed light curves.

Microlensing
Stars (and dark matter substructures) in the lens galaxy can give rise to additional gravitational lensing effects on top of the lens' macro magnification.Such microlensing effects from stars are typically able to magnify or demagnify the lensed SN images by approximately one magnitude.Since the resulting microlensing magnifications are not symmetric -some systems will be highly magnified while the majority will be slightly demagnified -their effects can change the number of lensed SNe that will pass our detection thresholds.We included microlensing in our simulations and investigated the resulting impact on the annual lensed SNIa detections.
To calculate microlensed light curves we follow the approach as described by Huber et al. (2019), where synthetic observables from theoretical SNIa models calculated via ARTIS (Kromer & Sim 2009) are combined with microlensing magnification maps.The maps are generated following Chan et al. (2021) and with software from GERLUMPH (Vernardos & Fluke 2014;Vernardos et al. 2014Vernardos et al. , 2015)).As in Huber et al. (2022) we create maps with a Salpeter initial mass function with a mean mass of the microlenses of ⟨⟩ = 0.35 M ⊙ , a resolution of 20000 × 20000 pixels and a total size of 20  E × 20  E .Here,  E corresponds to the physical Einstein radius of the microlenses at the source redshift and can be calculated via where  l ,  s and  ls are the angular diameter distances between the observer and the lens, the observer and the source, and the lens and the source, respectively.Further, we list in Table 2 the convergence , the shear  and the smooth matter fraction  ( = 1− * /, where  * is the convergence of the stellar component) for all magnification maps considered in this work.The specific realisations of ,  and  were chosen because they correspond to the most commonly occurring combinations amongst the simulated lensed SNIa.We normalise the microlensing magnification maps to have the same mean as the theoretical magnification predicted from the map's  and  values: (5) For each map we have 40000 microlensed spectra coming from 10000 random positions in the map and four theoretical SN models, the same as used by Suyu et al. (2020); Huber et al. (2021) and Huber et al. (2022).For all the maps listed in Table 2 we assumed a source redshift of 0.77 and a lens redshift of 0.32, which corresponds to the median values of the OM10 catalog (Oguri & Marshall 2010a) and defines the total size of the map  E .For our lensed SNe we are interested in  src between 0.0 and 1.4.To reduce the computational effort we grid the  src space in steps of 0.05.Given that the calculation of 10000 microlensed spectra for a single magnification map with a certain  E is on the order of a week, we approximate the microlensing contributions, as we now describe.For any source redshift of interest  src we use the microlensed spectra calculated for the source redshift of 0.77.We then rescale the spectra such that they correspond to  src in terms of absolute flux, wavelength and time after explosion.From the corrected spectra we can then calculate the exact light curves for  src following Huber et al. (2019), with the approximation that the total size of the microlensing map is the same as for the source redshift of 0.77.Using the same total size can slightly overestimate or underestimate the impact of microlensing, but Huber et al. (2021) tested different  E values, where no significant dependence between the strength of microlensing and  E was found.
For each simulated lensed SN, we compute the convergence, shear and smooth matter fraction at the position of the SN images and draw a random microlensing realisation from the magnification map with the closest ,  and  values.The local convergence and shear are calculated from the lens galaxy's mass model and the smooth matter fraction is obtained by approximating the stellar convergence  * at the image positions, for which we assume a spherical de Vaucouleurs profile (Dobler & Keeton 2006): where  = 7.67,  is the radius of the SN image position to the lens centre,  eff is the effective radius of the lens and  is a normalisation constant that is calibrated for each lens system such that the maximum  value is 1.The effective radius  eff of the lens is determined through a scaling relation between the radius in kpc h −1 and velocity dispersion  in km s −1 of elliptical galaxies (Hyde & Bernardi 2009): The above relations also assume spherical symmetry of the lens galaxy, but are only used to determine which ,  and  values are the best approximations for the image positions.
Finally, the microlensing contributions are obtained by drawing a random position in the chosen magnification map for each lensed SN image.Since the SN explosion models comprise different sizes at different wavelengths, the chromatic microlensing contributions are computed for each LSST filter and added to the simulated lensed SN light curves, as illustrated in Fig. 3.

THE VERA C. RUBIN OBSERVATORY
The Vera C. Rubin Observatory is a survey facility currently under construction on Cerro Pachón in Chile.It will host the Legacy Survey of Space and Time (LSST), a wide-field astronomical survey scheduled to start operations around 2025.The survey will take multi-colour  images and cover ∼20,000 square degrees of the sky in a ten-year period.Due to its depth and sky coverage, LSST is the most promising transient survey for observing gravitationally lensed SNe, with initial predicted numbers of several hundred discoveries a year (Goldstein & Nugent 2017;Goldstein et al. 2019;Wojtak et al. 2019;Oguri & Marshall 2010b).

LSST observing strategy
LSST will operate several survey modes.The main programme, comprising ∼ 90% of observing time, will be the Wide-Fast-Deep (WFD) survey, consisting of an area of 18,000 square degrees.The other major survey programmes include the Galactic plane, polar regions and  the Deep Drilling Fields (DDFs); the latter will be observed with deeper coverage and higher cadence.Following the latest recommendations from the Survey Cadence Optimization Committee 4 , the WFD survey is expected to proceed using a rolling cadence, in which certain areas of the WFD footprint will be assigned more frequent visits, with the focus of increased visits "rolling" over time.This improves the light curve sampling of the objects discovered in those high-cadence areas, the active regions.The drawback is that since the sky coverage is not homogeneous in any given period, there is 4 https://pstn-055.lsst.io/ a greater chance of missing rare events if they occur in an undersampled area, the so-called background region.For unlensed SNe, it has been shown by Alves et al. (2023) that the active region yields a 25% improvement in type classification performance relative to the background region.One of the goals of this work is to determine the impact of the rolling cadence on lensed SNIa.Huber et al. (2019) investigated the effects of several LSST observing strategies on lensed SNIa discoveries, but this previous work considered earlier versions of the observing strategy that differed significantly from the current implementation of the rolling cadence.
We implement the baseline v3.0 observing strategy, which adopts a half-sky rolling cadence with a ∼ 0.9 rolling weight (corresponding to the background regions receiving only 10% of the standard number of visits, and the active regions the rest).While the baseline v3.3 cadence was recently released, we do not expect it to significantly alter our results.The biggest change with respect to v3.0 is an updated mirror coating which results in decreased -band sensitivity and increased sensitivity in the  bands.Since lensed SNIa are red and our analysis does not consider -band detections, the change from v3.0 to v3.3 is minimal, and only expected to result in slightly more lensed SNIa detections.A sky map with the observations corresponding to the baseline v3.0 strategy is depicted in Fig. 4. The rolling cadence begins roughly 1.5 years after the start of the survey, to first allow for a complete season of uniform observations.This gives us an ideal opportunity to compare the effects of the rolling cadence on the discovery of lensed SNIa.For the annual lensed SN detections computed in this work, we only consider observations in the WFD and the DDFs, since those do not suffer from severe dust extinction found e.g. in the Galactic plane regions.An example of a lensed SNIa with baseline v3.0 WFD cadence observations is depicted in Fig. 3.

Simulating LSST observations
In order to simulate observations of lensed SNe at the catalogue level with sufficient information to define useful metrics, we need to find the set of times at which SNe at particular locations are observed, along with the observational metadata required to estimate the uncertainty with which the SN flux will be measured.This information can be accessed through the Rubin Operations Simulator (OpSim), which simulates the field selection and image acquisition process of LSST over the 10-year duration of the planned survey (Delgado & Reuter 2016;Delgado et al. 2014;Naghib et al. 2019).Detailed information about each simulated pointing of the telescope is stored in an output data product in the form of a 'sqlite' database, where each pointing forms a row in the database.Using OpSimSummary (Biswas et al. 2020), we find all the observations (and associated metadata) that include the position of a given SN within the field of view of the telescope.
In order to separate the Galactic plane region and the WFD and DDF surveys, we use a threshold based on the number of visits a sky location has received after 10 years of LSST observations.This serves as a proxy for OpSim's distinction between the Galactic plane and WFD regions.We assume that regions with  visits,10yr < 400 belong to the Galactic plane and polar regions,  visits,10yr > 1000 are the DDFs, and all remaining sky locations are assigned to the WFD.Within the WFD, we distinguish between observations taken during the non-rolling phase (year 0-1.5;MJD < 60768) and the first rolling period (year 1.5-3; 60768 < MJD < 61325).
From the OpSim database, we obtain the observing times, filters, mean 5-sigma depth ( 5 ), and point-spread functions (psf).We compute the 1 noise on the SN flux (which contains contributions from the sky brightness, the airmass, the atmosphere, and the psf) from the 5-sigma depth in the following way: where ZP corresponds to the instrument zero-point for a given band: 28.38, 28.16, 27.85, 27.46, 26.68 . (10)

DETECTING LENSED SNIA
In this section, we describe our methods for calculating the number of lensed SNIa detections and the properties of the detected sample.We examine the simulated lensed SNIa based on their colours, magnitudes, time-delay measurements, and prospects for follow-up.

Annual lensed SNIa detections
Studies that predict the number of lensed SN discoveries generally take into account two distinct detection methods.The first is the image multiplicity method, which looks for multiple resolved images of the lensed SN (Oguri & Marshall 2010b), and the second is the magnification method, which looks for objects that appear significantly brighter than a typical SN at the redshift of the lens galaxy (which acts as the apparent host galaxy) (Goldstein & Nugent 2017).For the latter method, the lensed SN images do not need to be resolved.We build our estimates of the lensed SN discoveries upon the results from Wojtak et al. (2019).They combine the image multiplicity and the magnification method and predict that LSST will discover around 89 lensed SNIa per year, assuming a 0.2 mag buffer above 5 average limiting magnitudes.For comparison, the predicted rate of unlensed SNe Ia, after quality cuts for cosmological utility, is 104000 for the ten-year sample.(The LSST Dark Energy Science Collaboration et al. 2018) The number of lensed SN detections from Wojtak et al. ( 2019) considers whether a lensed SN passes the image multiplicity and magnification cuts based on its full light curve information.There is no observation cadence information included in the predictions, which would alter the results if a lensed SNe will occur between observing seasons or when sufficient observations are missing around the peak.Here, we update the annual lensed SNIa detections taking into account the cadence from the baseline v3.0 observing strategy.For each simulated lensed SN, we draw a random observation sequence and assess whether the object still passes the detection cuts for the magnification and image multiplicity method.
The criteria for being "detected" with both methods are the following:

Image multiplicity method
• The maximum image separation  max is larger than 0.5 ′′ and smaller than 4 ′′ ; • For doubles: the flux ratio between the images (measured at peak) is between 0.1 and 10; • At least three or two images are detected (signal-to-noise ratio > 5) for quads (systems with four images) and doubles (systems with two images), respectively.

Magnification method
• The apparent magnitude of the unresolved lensed SN images should be brighter than a typical SNIa at the lens redshift at peak: with  X () the apparent magnitude of the transient in band X at time , ⟨ X ⟩( peak ) the absolute magnitude of a standard SNIa in band X at peak,  the distance modulus,  XX the K-correction, and Δ the magnitude gap (adopted here to be −0.7 for consistency with Wojtak et al. ( 2019) and Goldstein & Nugent ( 2017)); • The combined flux of the unresolved data points should be above the detection threshold (signal-to-noise ratio > 5).

Colours and magnitudes of lensed SNIa
SNe affected by strong gravitational lensing are expected to look different from unlensed SNe in several ways.Fig. 2 shows that in general, lensed SNe will be found at higher redshifts than unlensed ones and hence, they will be observed as redder and more slowly evolving.Additionally, for lensed and unlensed SNIa at the same redshifts, the lensed SNIa will appear brighter because of the gravitational lensing magnification.
In this analysis, we investigate which observables are best suited to separate the populations of lensed and unlensed SNIa in LSST data.We aim to investigate optimal selection criteria based on the brightness and colours of lensed SN candidates.We measure light curve properties in the observer frame from the simulated sample of lensed and unlensed SNIa as observed with LSST.The resulting observables are apparent magnitudes for each LSST filter (,,,,) and all colour combinations ( − ,  − ,  − ,  − ,  − ,  − ,  − ,  − ,  −  and  − ), at different epochs at the light curve peak, which is determined in the following way.A polynomial fit is performed on every light curve of the sample to find the peak time from the filters with the best detection cadence, which mostly corresponds to the  or  bands.We use the same polynomial fits to obtain the expected apparent magnitudes at the given epochs to compute the colours.Error bars from the detections are considered in the fits and propagated into uncertainties on the measured magnitudes and colours.
The redshift distributions for unlensed and lensed SNIa observed with LSST are largely overlapping, as shown in Fig. 2. Due to this, using the apparent magnitude and colours alone will not serve as good metrics to separate lensed from unlensed SNIa, as is also illustrated in Fig. A1 and A2.Therefore, we chose to investigate cuts based on all combinations of colours versus apparent magnitudes.In Quimby et al. (2014), such a colour-magnitude cut in the  and -band is shown to be very promising for distinguishing lensed SNIa.We aim to devise linear cuts in colour-magnitude space for all LSST bands that exclude most of the unlensed events while preserving the lensed ones.Our method consists of obtaining the 90% contour for the unlensed SNe in each colour-magnitude space and fitting a linear function to this contour to extract a simple linear cut.

Time-delay measurements
For a fraction of the lensed SN discoveries, LSST will be able to resolve the individual images.Here, we compute the fraction of those systems for which we will be able to infer the time delay precisely using LSST data only.Even for events that are on the cusp of being resolvable and with variable seeing, a single epoch with sufficient seeing to resolve the images will enable the extraction of the full light curves using forced photometry.To find the objects with the best time-delay measurements, we limit ourselves to systems with an angular separation of   > 0.5".We use the most commonly implemented model for the SED of an SNIa, SALT3 (Guy et al. 2010;Betoule et al. 2014;Kenworthy et al. 2021) as implemented in SNCosmo (Barbary et al. 2016) and fit each light curve with a common stretch and colour parameter.Since our aim is to infer what fraction of objects have an accurately estimated time delay, we do not use a simultaneous inference for extinction and magnification like for iPTF16geu (Dhawan et al. 2020) or SN Zwicky (Goobar et al. 2023).We assume the SN redshift will be known from spectroscopic followup observations, which can be carried out after the SN has faded away.The difference in returned  0 values provides the predicted time delay between the images, Δ  ,pred .We classify a system as having a "good" time delay when Δ  ,pred − Δ  ,true /Δ  ,true < 5%.An example case is shown in Fig. 5.We only use LSST data for Δ measurements and hence our results are a conservative estimate which will improve further with follow-up observations, especially if the follow-up also resolves the individual images.

Gold sample for follow-up observations
In order to conduct timely follow-up observations of the lensed SNe, we should find them early on in their evolution.After the initial detection of a lensed SN candidate in LSST, spectroscopic followup observations are needed to verify its lensed nature.A spectrum will reveal the SN type and its redshift, which for SNIa will identify the objects that are magnified by strong gravitational lensing.From simulated detected lensed SNe, we construct a 'gold' sample that satisfies the following criteria: •  premax > 5 in at least two filters; •   < 22.5 mag; • Δ > 10 days, with  premax the number of detections with signal-to-noise ratio > 3 before the SN peak and   the apparent -band magnitude at peak.We use SALT3, as implemented in SNCosmo (Barbary et al. 2016) to fit the light curves of the simulated lensed SNIa and infer the time of peak.To trigger spectroscopic follow-up observations we require that the object is detected in at least two filters and has a minimum of five observations before the inferred time of maximum.This is because it allows for spectroscopic follow-up when the lensed SN is close to its brightest, while still having ample time for scheduling highresolution follow-up.In addition, we apply the constraint that the lensed SNe should be bright enough to get a classification spectrum with a 4-m class telescope, e.g.4MOST (Swann et al. 2019) or e.g. the New Technology Telescope (Snodgrass et al. 2008) with exposure times ≲ 1hr, corresponding to a brightness of   < 22.5 mag.Alternatively, with shorter exposure times, the spectra can be  obtained with instruments on 8m class telescopes, e.g. the Gemini Multi-Object Spectrographs (GMOS; Crampton et al. 2000).

Annual lensed SNIa detections
We simulate a set of 5, 000 doubly imaged lensed SNe and 5, 000 quadruply imaged SNe that pass either the image multiplicity or the magnification method as described in Sec.4.1, where the size of the sample is chosen such that it is large enough for our statistical analysis.These numbers are subsequently scaled with the predicted lensed SNIa rates from Wojtak et al. (2019): 64 doubles and 25 quads per year.We then compute what fraction of the simulated lensed SNIa remains detectable when the baseline v3.0 observing strategy is applied.We find that 46% of doubles and 70% of quads remain from the full simulated sample.We also assess the impact of the rolling cadence on the annual lensed SNIa detections, by separating the sample into objects that are detected in the first 1.5 years of the survey (MJD < 60768) and objects discovered in years 1.5 to 3 (60768 < MJD < 61325).
Fig. 6 shows the number of detections with signal-to-noise ratio > 5 ( detect ) per lensed SN system for the rolling and non-rolling cadence.The distribution for a non-rolling baseline v3.0 cadence peaks around 20 detections per lensed SN, while the rolling cadence has a large tail towards systems with higher numbers of detections.The background regions (corresponding to the dark areas in Fig. 4) acquire a lower number of detections, while the active regions (light areas in Fig. 4) receive a higher  detect .As a result, the non-rolling cadence scans a larger area of the sky with a medium cadence and therefore discovers more lensed SNe, while the rolling-cadence provides more detections for the systems it discovers.This effect is illustrated in Fig. 7, which shows the expected number of lensed SNIa per year with  detect above a certain threshold.The non-rolling cadence will discover more lensed SNe up to  detect < 20, while the rolling cadence will find a larger sample with well-sampled light curves ( detect > 20).Nevertheless, we note that the differences are relatively small.We also compute the number of lensed SNIa that fall in the Deep Drilling Fields (DDFs) in our simulation, since those objects will be observed with a much higher cadence and better depth.However, we find that only ∼ 0.2 lensed SNIa per year are expected to be in the DDFs, which is not surprising given the small area covered relative to WFD.
Our findings are summarised in Table 3, which contains the predicted annual number of lensed SNIa detections for a non-rolling versus a rolling cadence.Our results are consistent with a recent study by Sainz de Murieta et al. (2023), which predicts the number of lensed SNIa detected with the magnification method for an approximate LSST survey strategy.

Microlensing impact
While all results presented so far included the effects of microlensing, we also generate each lensed SN light curve both with and without microlensing in order to clearly quantify the microlensing impact for each system.We distinguish three scenarios, in which microlensing effects (i) do not change the detectability of the lensed SN; (ii) make a detected lensed SN undetectable; (iii) make an undetected lensed SN detectable.Fig. 8 investigates these three scenarios.It shows the difference in apparent peak -band magnitude for the 5,000 simulated doublyimaged SNe.The red dots are the systems that become undetectable because of microlensing, while the green dots are the ones that have become detectable due to the microlensing magnifications.The sum of these effects is that we detect a handful fewer lensed SNe; Table 3 presents that we go from a total annual number of 48 (54) without microlensing to 44 (50) with microlensing for a rolling (non-rolling) cadence.For the events with longer time delays than 10 days, we predict to find 29 (31) without microlensing and 26 (30) with microlensing for a rolling (non-rolling) cadence.This weak effect of detecting fewer objects when microlensing is included in the simulations can be Table 3.The annual expected number of discovered lensed SNIa (doubles, quads, and the total number) in the baseline v3.0 observing strategy, with separate predictions for the non-rolling and rolling cadence.The rows list the predicted numbers of lensed SNIa that are detected without and with microlensing (see 5.2), are detected with the magnification (mag.)method and image multiplicity (im.mult.)method, have time delays larger than 10 days, pass the colour-magnitude cut (described in 5.3), and are in the gold sample for follow-up (5.5).
understood when looking at the projected 1D distributions of Fig. 8, which shows that the majority of events will be slightly demagnified due to microlensing, while a rare few will be highly magnified.Note also the asymmetry in the distributions.The first image is a minimum of the time delay and can never be demagnified -microlensing can therefore only affect it subject to this constraint.The second image is a saddlepoint with no such restriction on its magnification and can therefore be more severely influenced by microlensing (Schechter & Wambsganss 2002).

Colour and magnitudes of lensed and unlensed SNIa
For each of the simulated lensed and unlensed SNIa, we calculate the apparent magnitudes at peak in every band, following the procedure outlined in Sec.4.2.The best separation between lensed and unlensed SNIa is achieved with the  −  peak colour versus observed apparent -band peak magnitude, which is shown in Fig. 9.Other colour and magnitude combinations are included in Appendix A. Due to their higher redshift distributions, lensed SNe are expected to appear redder than unlensed ones.However, since LSST will also  detect unlensed SNe at high redshifts (see Fig. 2), this difference is less pronounced in LSST than in precursor surveys such as ZTF.The overlap in redshift constitutes a potential difficulty when it comes to distinguishing lensed SNe from unlensed ones in LSST.Nevertheless, Fig. 9 demonstrates that we can achieve a better separation by combining colours with apparent magnitudes, since lensed SNe (especially the quads) are magnified and hence brighter than unlensed ones at the same redshifts.We also see a few very red lensed SNe at high redshifts where unlensed SNe are not visible anymore with Rubin (redshifts ⪆ 1).
We investigate each colour and magnitude combination at multiple epochs and obtain the following linear cuts using the method described in Sec.We calculate the percentage of the lensed SNIa and the percentage of background contaminants from the simulated unlensed population that pass the colour-magnitude cuts.We find that 41% (43%) of the lensed doubles (quads) pass at least one of the mentioned colour magnitude cuts.2-3% of the unlensed SNIa would also pass this colour-magnitude cut.We find that the parameters combination that best separates lensed from unlensed SN for LSST is  −  peak colour versus observed apparent -band peak magnitude, which keeps 23% (26%) of the doubles (quads) with only 0.7% of the unlensed sample.28% (32%) of the lensed doubles (quads) would pass  − peak colour versus observed apparent -band peak magnitude, but with a almost 2% of the unlensed events.Since unlensed SNe outnumber lensed ones by a factor of ∼ 10 3 we want to keep the false-positive rate low.The combination with the lowest contaminants is  −  peak colour versus observed apparent -band peak magnitude, with only 0.1% unlensed events, but due to the poorer cadence and depth for the filter, requiring detections in  means we only preserve 15% (20%) of doubles (quads).As detailed in Table 3, this corresponds to around 20 lensed SNIa a year that pass one of the colour-magnitude cuts.While colour and magnitude cuts alone are not enough to distinguish completely between lensed and unlensed SNIa, this analysis shows that they can be useful additional tools to inform us about lensed SNIa candidates in LSST.

Time-delay measurements from LSST-only data
For a small fraction of the simulated lensed SNIa, we find that we can extract a useful time-delay measurement using only LSST data.An example of such an object is shown in Fig. 5.However, for most cases the light curve has a sparse sampling such that a SALT3 fit with constrained  1 and  is unsuccessful.Less than 2% of our detected lensed SNIa sample allows for a Δ measurement with accuracy better than 5%.Additionally, in Hayes et al. ( 2023) our lensed SNIa simulations are fitted with GausSN6 , a Bayesian semi-parametric Gaussian Process time-delay estimation model.GausSN recovers for ∼ 13% of the resolved double systems the time delays within 5% accuracy, corresponding to 6% of the total detected sample.These results emphasise the importance of follow-up observations to improve the quality of the time-delay measurements, in line with the conclusions from Huber et al. ( 2019).The number of lensed SNIa a year that qualifies for accurate Δ measurements with follow-up observations is discussed in Sec.5.5.

Gold sample and cosmological prospects
Fig. 10 shows the early detections, peak magnitude, and time delay distributions of the detected sample.When applying the cuts described in Sec.4.4, we find that 25% of the detected lensed SNIa belong to the 'gold' sample.This corresponds to roughly 10 systems per year, as outlined in Table 3, for which high-quality follow-up observations and precision cosmology measurements are expected to be feasible.
To estimate the cosmological prospects of such a gold sample, we assume the availability, for each system in the sample, of groundbased follow-up observations to sample the light curve well, and high-resolution imaging and spatially-resolved spectroscopy to constrain both the time delays and the lens mass model.We expect that the uncertainty in the inferred time delay, with high-resolution followup, will be ∼ 0.1-0.5 days, as inferred from local SN Ia samples (e.g., Johansson et al. 2021) corresponding to a conservative upper limit on the time-delay error of 5%.For the lens mass model, we assume an uncertainty of 7%.Since SNe are explosive transients, we can obtain post-explosion images to cross-check the lens model, and hence, reduce the uncertainty in the final lens mass model estimate.Ding et al. (2021) predicts a factor 4 improvement in the lens model for lensed SNe over lensed quasars, although the difference will be smaller for observations at higher angular resolution than the Hubble Space Telescope.Observations with Integral Field Units (IFUs) can measure the stellar kinematics and help to further break the mass sheet degeneracy (Birrer & Treu 2021).We note, furthermore, that systems which do not have very precisely measured time delays (e.g., iPTF16geu, SN Zwicky Dhawan et al. 2020;Goobar et al. 2023) can still be important for reducing uncertainties on the mass modelling, via a precisely measured model independent estimate of the lensing magnification (Birrer et al. 2021a).Combining the uncertainties from the time-delay measurement and the mass modelling, we obtain a precision of 8.6% in  0 for each system in the 'gold' sample.Consequently, we would need 30 lensed SNIa to reduce the uncertainty to 1.5% in  0 , corresponding to ∼ 3 years of LSST observations.Furthermore, we note that lensed SNe with shorter time delays (e.g. 5 < Δ < 10 days) will also contribute to improving the precision, even though the individual uncertainties per system would be greater than 8.6%.In that case, we could reach the expected precision in a shorter duration of the LSST survey.

DISCUSSION AND CONCLUSIONS
In this work, we studied the detectability of lensed SNIa in LSST.We have investigated the impact of the LSST baseline v3.0 observing strategy and of microlensing on the predicted annual lensed SNIa detections.Similar to Wojtak et al. (2019), we consider both the magnification method as used in Goldstein & Nugent (2017) and image multiplicity method as used in Oguri & Marshall (2010b) for detecting lensed SNIa.Compared to the aforementioned simulations, we go from limiting magnitude cuts to fully simulated LSST observations including weather effects and stellar microlensing.The LSST observing strategy is expected to proceed using a rolling cadence, in which certain areas of the WFD footprint will be assigned more frequent visits than others.The expected yearly number of lensed SNIa is higher for a non-rolling cadence (50 events) than for a rolling cadence (44 events), but the difference does not appear to be detrimental to the lensed SNe science case.Microlensing effects from stars in the lens galaxy result in a handful fewer detected lensed SNe per year.We found that ∼ 40% of lensed SNIa detected in LSST will stand out from unlensed SNIa with simple linear cuts in colour and peak magnitude.Using only LSST data, a time delay within 5% of the truth value is expected to be measured for only a small fraction, ∼ 2% of the systems.Hence it is important to assess the feasibility of time-delay and  0 measurements from follow-up observations.We have determined a set of detectability criteria that will allow for timely follow-up and cosmological inference.Our results predict ∼ 10 lensed SNIa per year that will have sufficient early detections and will be sufficiently bright for follow-up observations, while also having time delays larger than 10 days to enable time delay cosmography.Assuming uncertainties of 8.6% in  0 per object, this sample is expected to enable a Hubble constant measurement of 1.5% precision in three years of LSST observations.Our results only focus on SNIa; the expected number of lensed core-collapse SNe is likely even higher (Wojtak et al. 2019;Goldstein & Nugent 2017;Oguri & Marshall 2010b).Future work is needed to investigate the cosmological prospects of strongly-lensed core-collapse SNe, which display more intrinsic variation than type Ias but could be used efficiently for spectroscopic time-delay measurements (see e.g.Bayer et al. (2021) for type IIP SNe).Additionally, in this work we have not investigated other sources of background contamination than unlensed SNIa.Other future studies include testing the resolvable separation by injecting lensed SNe in the LSST difference imaging analysis (DIA) pipeline (e.g.Liu et al., in prep.) and testing the pixel-to-cosmology performance for lensed SNe in an analysis similar to Sánchez et al. (2022).
We have shown that lensed SNIa discovered with LSST will be excellent precision probes of cosmology.Crucial to the success of this programme is the availability and coordination of follow-up resources, both monitoring of light curves to measure time delays, and high-resolution imaging data and spatially resolved kinematics to constrain the lens mass model.With the selection cuts applied, our work shows that the sample will be sufficiently bright to measure the present-day expansion rate of the Universe with high precision.

APPENDIX A: COLOUR AND MAGNITUDE CUTS
Here, we provide all results from the colour-magnitude investigation of simulated lensed and unlensed SNIa as described in Sec.4.2.Fig. A1 shows the 1D distributions of peak apparent magnitudes in the , , , -bands and Fig. A2 the colours at peak.The joint colour-magnitude diagrams are displayed in Fig. A3.Our results predict that the colour-magnitude combination with the strongest separation between lensed and unlensed SNIa is  −  colours versus -band magnitudes, but we combine all colour-magnitude combinations for the best result.9) to calculate the total number of lensed events that are distinguishable from unlensed SNIa.

Figure 1 .
Figure1.The joint distribution of the lens redshift ( lens ), source redshift ( src ) and Einstein radius ( E ) used to simulate the lensed SNIa systems.The  lens ,  src and  E combinations correspond to galaxy-source configurations where strong lensing occurs(Wojtak et al. 2019).The sample includes all lensed SNIa at redshift  src < 1.5.

Figure 2 .
Figure 2. Normalised redshift distributions of simulated lens galaxies, lensed SNIa and unlensed SNIa that will be detectable with the Rubin Observatory.The Figure shows the observed populations, after the detection criteria described in Sec.2.1 for unlensed SNIa and in Sec.4.1 for lensed SNIa have been applied.

Figure 3 .
Figure 3. Simulated observations of a doubly imaged SNIa in LSST.For each lensed SN image, the figure shows the microlensing magnifications maps and the corresponding -band light curves in the observer frame (dashed curves without microlensing and solid curves with microlensing).The coloured markers correspond to LSST observations in the active region of the WFD survey with the baseline v3.0 cadence.The properties of this simulated lensed SN system are  lens = 0.1,  src = 0.52,  E = 1.26 ′′ , Δ = 22 days.

Figure 4 .
Figure 4. LSST survey footprint in equatorial coordinates showing the number of visits ( visits ) conducted per sky location.The light gray area corresponds to the full LSST footprint, including the Galactic plane and polar regions.The coloured region represents the WFD area (∼ 18, 000 sq. deg.) and the light yellow patches are the DDFs.In the first 1.5 years of the survey (left panel) the sky coverage will be homogeneous, after which the rolling cadence will start (right panel).Rolling results in background regions with lower  visits and active regions with higher  visits .
, respectively for the , , , ,  bands 5 .The SN flux for each image is perturbed by drawing a new flux value from a normal distribution with mean of the model flux and width equal to the sky noise.Then, the new flux values  SN 5 https://smtn-002.lsst.io/#change-recordare combined with the SN flux noise to calculate the error on the observed magnitude: −2.5 •   SN  SN • ln(10)

Figure 5 .
Figure 5. Example of a lensed SNIa with a robust time-delay inference from LSST data only.The markers correspond to LSST observations taken in the  , , ,  bands and the solid curves show the SALT3 fits used to infer the time delay.This object is at  src = 0.464, lensed by a deflector galaxy at  lens = 0.142, an input time delay of 11.14 days and a recovered time-delay estimate of 12.23 days.The Einstein radius of the system is  E = 0.63 so the images are treated as resolved.

Figure 6 .
Figure 6.Number of detections ( detect ) per lensed SNIa for the non-rolling cadence (first 1.5 years of the survey) compared to the rolling cadence (years 1.5 -3 of the survey).The distributions show that compared to the non-rolling cadence, the rolling cadence produces both more systems with small  detect and with large  detect .The mean number of observations is 15 (17) for a non-rolling (rolling) cadence.

Figure 7 .
Figure 7. Annual expected number of lensed SNIa with  detect above a certain threshold, for observing strategies without a rolling cadence (blue line) and with rolling cadence (orange line).The figure shows that although non-rolling cadences are expected to discover a larger overall number of lensed SNe, rolling cadences will provide more lensed SNe with a high number of detections.

Figure 8 .
Figure8.Impact of microlensing on the lensed SNIa detections.The scatterplot shows the difference in apparent peak -band magnitude (Δ  ) due to microlensing for image one (first occurring) and two of the 5,000 simulated doubly-imaged SNe.The gray points correspond to lensed SNe whose detection is not impacted by microlensing, the red points are the systems that have become undetectable by microlensing, and the green ones have become detectable by a microlensing magnification boost.The projected 1D distributions show the microlensing effect on the apparent magnitude per lensed SN image.

Figure 9 .
Figure 9. Peak colours and magnitudes of the detected lensed SNIa (doubles and quads) and unlensed SNIa.The  −  colours and -band magnitudes are able to separate the populations because lensed SNe are brighter than unlensed ones at similar redshifts.The dashed black line shows a simple linear separation cut of  −  > 0.88   − 19.1.

Figure 10 .
Figure10.Normalised distributions showing the properties of doubles and quads from the detected lensed SNIa sample.From left to right, the panels show the total magnification ( tot ), maximum time delay between the images (Δ max ), apparent -band magnitude at peak ( i ), and the number of detections before peak ( premax ).The non-shaded regions indicate the objects that satisfy the conditions for the 'gold' sample.

Figure A3 .
Figure A3.Colour-magnitude cuts described in Sec.5.3 used (together with  −  versus  shown in Fig.9) to calculate the total number of lensed events that are distinguishable from unlensed SNIa.

Table 2 .
Combinations of the convergence (), shear (), and smooth matter fraction () used to simulate the microlensing contributions to the lensed SN light curves.