Detecting strongly lensed type Ia supernovae with LSST

ABSTRACT

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 lensing galaxy are predicted to lower the lensed SNIa detections by ∼8 per cent. The lensed events can be separated from the unlensed ones by jointly considering their colours and peak magnitudes. We define a ‘gold sample’ of ∼10 lensed SNIa per year with time delay >10 d, >5 detections before light curve peak, and sufficiently bright (m_i < 22.5 mag) for follow-up observations. In 3 yr of LSST operations, such a sample is expected to yield a 1.5 per cent measurement of the Hubble constant.

gravitational lensing: strong, methods: statistical, transients: supernovae

1 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, sub-structures 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 (H₀) – 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 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 Time-Delay COSMOgraphy (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 minimize the mass-sheet degeneracy (Falco, Gorenstein & Shapiro 1985; Gorenstein, Falco & Shapiro 1988; Kolatt & Bartelmann 1998; Saha 2000; Oguri & Kawano 2003; Schneider & Sluse 2013; Birrer, Dhawan & Shajib 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 minimize 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, Schechter & Wambsganss 2017; Foxley-Marrable et al. 2018; Weisenbach, Schechter & Pontula 2021; Weisenbach 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 per cent (Kelly et al. 2015) and ∼9 per cent (Pascale et al. 2024), respectively, where the precision was primarily limited by the cluster mass models. While SN Encore is expected to enable an H₀ measurement with ∼10 per cent 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 (Oguri & Marshall 2010; Goldstein & Nugent 2017; Wojtak, Hjorth & Gall 2019).

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 Lambda cold dark matter (ΛCDM) model with |$H_0 = 67.8 \ \textrm {km} \, \textrm {s}^{-1} \textrm {Mpc}^{-1}$| and Ω_m = 0.308 (Planck Collaboration 2016).

We outline our lensed SNIa simulation procedure in Section 2, the LSST observing strategy in Section 3, and the different aspects of detecting lensed SNIa in Section 4. Our results are presented in Section 5 and our discussion and conclusions in Section 6.

2 SIMULATING LENSED SNIA

In this work, we develop the publicly available python code called Lensed Supernova Simulator Tool (lensedsst¹), 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.

2.1 Lens galaxy mass profile assumptions

We employ the multipurpose lens modelling package lenstronomy² (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 κ:

$$\begin{eqnarray} \kappa (x, y) = \frac{3-\gamma _{\rm lens}}{2} \left(\frac{\theta _E}{\sqrt{q_\textrm {lens} x^2 + y^2/q_\textrm {lens}}} \right)^{\gamma _{\rm lens}-1}, \end{eqnarray}$$

(1)

where q_lens is the projected axis ratio of the lens, γ_lens corresponds to the logarithmic slope, and θ_E denotes the Einstein radius. The coordinates (x, y) are centred on the position of the lens centre, and rotated by the lens orientation angle ϕ_lens, such that the x-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.

Table 1.

Parameter distributions for lensed SNIa. The distribution of input parameters employed in the simulation pipeline to generate lensed SNIa light curves. |$\mathcal {N}(\mu , \sigma)$| indicates a normal distribution with mean μ and standard deviation σ, |$\mathcal {N}^{\mathcal {S}}(a, \mu , \sigma)$| denotes a skewed normal distribution with skewness parameter a, while |$\mathcal {U}(x, y)$| represents a uniform distribution with bounds x and y. The skewed normal distributions for z_lens, z_src, and θ_E are one-dimensional representations of the full joint distribution from Wojtak et al. (2019) depicted in Fig. 1.

Parameter	Distribution
Hubble constant	|$H_0 = 67.8\ \textrm {km}\, \textrm {s}^{-1} \textrm {Mpc}^{-1}$|
Lens redshift	|$z_{\textrm {lens}} \sim \mathcal {N}^{ \mathcal {S}}(3.88, 0.13, 0.36)$|
Lensed source redshift	|$z_{\textrm {src}} \sim \mathcal {N}^{ \mathcal {S}}(3.22, 0.53, 0.55)$|
Unlensed source rate	|$r_{v} (z) = 2.5 \cdot 10^{-5} (1+z)^{1.5} \, \textrm {Mpc}^{-3}\textrm {yr}^{-1}$|
Source position (doubles)	|$x_\textrm {src}, y_\textrm {src} \sim \mathcal {U}(-\theta _{\textrm {E}}, \theta _{\textrm {E}})$|
Source position (quads)	|$x_\textrm {src}, y_\textrm {src} \sim \mathcal {U}(-0.4\theta _{\textrm {E}}, 0.4\theta _{\textrm {E}})$|
Lens galaxy
Elliptical power-law mass profile
Lens centre (arcsec)	xlens, ylens ≡ (0, 0)
Einstein radius (arcsec)	|$\theta _{\textrm {E}} \sim \mathcal {N}^{ \mathcal {S}}(5.45, 0.14, 0.63)$|
Power-law slope	|$\gamma _{\textrm {lens}} \sim \mathcal {N}(2.0, 0.2)$|
Axis ratio	|$q_{\textrm {lens}} \sim \mathcal {N}(0.7, 0.15)$|
Orientation angle (rad)	|$\phi _{\textrm {lens}} \sim \mathcal {U}(-\pi /2, \pi /2)$|
Environment
External shear modulus	|$\gamma _{\textrm {ext}} \ \sim \mathcal {U}(0, 0.05)$|
Light curve
Stretch	|$x_1 \ \sim \mathcal {N}^{ \mathcal {S}}(-8.24, 1.23, 1.67)$|
Colour	|$c \ \sim \mathcal {N}^{ \mathcal {S}}(2.48, -0.089, 0.12)$|
Absolute magnitude	|$M_{\mathrm{B}} \ \sim \mathcal {N}(-19.43, 0.12)$|
Milky Way extinction	|$E(B-V) \ \sim \mathcal {U}(0, 0.2)$|

2.2 Redshift distribution

Most of our input parameters are uncorrelated, except for the Einstein radius θ_E, lens redshift z_lens, and source redshift z_src, which we sample from a kernel density estimation probability model fitted to (θ_E, z_lens, z_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, Park & Vogeley 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 z_src < 1.5 to ensure that the SNe are not redshifted out of the filters. The resulting combinations of z_lens, z_src, and θ_E values are depicted in Fig. 1 and the projected one-dimensional distributions can be found in Table 1.

Figure 1.

The joint distribution of the lens redshift (z_lens), source redshift (z_src), and Einstein radius (θ_E) used to simulate the lensed SNIa systems. The z_lens, z_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 z_src < 1.5.

For the background population of unlensed SNIa, we assume a volumetric redshift rate of (Dilday et al. 2008)

$$\begin{eqnarray} r_{v} (z) = 2.5 \times 10^{-5} (1+z)^{1.5} \ \textrm {Mpc}^{-3} \,\, \textrm {yr}^{-1}. \end{eqnarray}$$

(2)

Our sample of detected unlensed SNIa consists of objects with i-band magnitude brighter than 24.0. For lensed SNIa we adopt stricter detection criteria, which are described in Section 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 z ∼ 1, resulting in a significant overlap in redshift space between lensed and unlensed SNIa.

Figure 2.

Normalized 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 Section 2.1 for unlensed SNIa and in Section 4.1 for lensed SNIa have been applied.

2.3 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³ (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 x₀, stretch parameter x₁, and a colour parameter c. We sample the x₁ and c 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), Panoramic Survey Telescope and Rapid Response System (Pan-STARRS1) (Rest et al. 2014), and several low-redshift surveys (Stritzinger et al. 2011; Hicken et al. 2012). We compute the distance modulus μ of each SNIa, based on its x₁ and c parameters:

$$\begin{eqnarray} &\mu = {m} + \alpha x_1-\beta c-\mathcal {N}(M_0, 0.12). \end{eqnarray}$$

(3)

Here, M₀ is the expected absolute magnitude of a SNIa with x₁ = c = 0 and m is the SN flux normalization in magnitude units. We assume M₀ = −19.43 in the B band, corresponding to a universe with Hubble constant |$H_0 = 67.8 \ \textrm {km} \, \textrm {s}^{-1} \textrm {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 z_lens, z_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 Section 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.

2.4 Microlensing

Stars (and dark matter sub-structures) 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 1 mag. 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. 2014, 2015). As in Huber et al. (2022) we create maps with a Salpeter initial mass function with a mean mass of the microlenses of 〈M〉 = 0.35 M_⊙, a resolution of 20000 × 20 000 pixels and a total size of 20 R_E × 20 R_E. Here, R_E corresponds to the physical Einstein radius of the microlenses at the source redshift and can be calculated via

$$\begin{eqnarray} R_{\rm E} = \sqrt{\frac{4G \langle M \rangle }{c^2} \frac{D_{\rm s} D_{\rm ls}}{D_{\rm l}}}, \end{eqnarray}$$

(4)

where D_l, D_s, and D_ls are the angular diameter distances between the observer and the lens, the observer and the source, and the lens and the source, respectively. Further, we list in Table 2 the convergence κ, the shear γ, and the smooth matter fraction s (s = 1 − κ_*/κ, where κ_* is the convergence of the stellar component) for all magnification maps considered in this work. The specific realizations of κ, γ, and s were chosen because they correspond to the most commonly occurring combinations amongst the simulated lensed SNIa. We normalize the microlensing magnification maps to have the same mean as the theoretical magnification predicted from the map’s κ and γ values:

$$\begin{eqnarray} \mu = \frac{1}{(1-\kappa ^2)-\gamma ^2}. \end{eqnarray}$$

(5)

For each map, we have 40 000 microlensed spectra coming from 10 000 random positions in the map and four theoretical SN models, the same as used by Suyu et al. (2020) and Huber et al. (2021, 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 catalogue (Oguri & Marshall 2010) and defines the total size of the map R_E. For our lensed SNe we are interested in z_src between 0.0 and 1.4. To reduce the computational effort we grid the z_src space in steps of 0.05. Given that the calculation of 10 000 microlensed spectra for a single magnification map with a certain R_E is on the order of a week, we approximate the microlensing contributions, as we now describe. For any source redshift of interest z_src, we use the microlensed spectra calculated for the source redshift of 0.77. We then rescale the spectra such that they correspond to z_src in terms of absolute flux, wavelength, and time after explosion. From the corrected spectra, we can then calculate the exact light curves for z_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 R_E values, where no significant dependence between the strength of microlensing and R_E was found.

Table 2.

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

Convergence (κ)	Shear (γ)	Smooth matter fraction (s)
0.362	0.342	0.443
0.655	0.669	0.443
0.655	0.952	0.443
0.956	0.669	0.443
0.956	0.952	0.443
0.362	0.342	0.616
0.655	0.669	0.616
0.655	0.952	0.616
0.956	0.669	0.616
0.956	0.952	0.616
0.362	0.342	0.790
0.655	0.669	0.790
0.655	0.952	0.790
0.956	0.669	0.790
0.956	0.952	0.790
0.362	0.280	0.910

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 realization from the magnification map with the closest κ, γ, and s 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):

$$\begin{eqnarray} \kappa _{*} (r) = A e^{-k(r/R_{\textrm {eff}})^{1/4}}, \end{eqnarray}$$

(6)

where k = 7.67, r is the radius of the SN image position to the lens centre, R_eff is the effective radius of the lens, and A is a normalization constant that is calibrated for each lens system such that the maximum s value is 1. The effective radius R_eff of the lens is determined through a scaling relation between the radius in kpc h⁻¹ and velocity dispersion σ in km s⁻¹ of elliptical galaxies (Hyde & Bernardi 2009):

$$\begin{eqnarray} &{\log _{10}(R_{\rm eff})} = 2.46 - 2.79 \log _{10}(\sigma) + 0.84 \log _{10}(\sigma)^2 , \end{eqnarray}$$

(7)

$$\begin{eqnarray} \textrm {with~} \ \sigma ^2 = \frac{c^2 \, \theta _{\rm E} \, D_{\rm s}}{4 \pi \, D_{\rm ls}}. \end{eqnarray}$$

(8)

The above relations also assume spherical symmetry of the lens galaxy, but are only used to determine which κ, γ, and s 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.

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 i-band light curves in the observer frame (dashed curves without microlensing and solid curves with microlensing). The 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 z_lens = 0.1, z_src = 0.52, θ_E = 1.26 arcsec, Δt = 22 d.

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 multicolour ugrizy images and cover ∼20 000 deg² of the sky in a 10 yr 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 (Oguri & Marshall 2010; Goldstein & Nugent 2017; Goldstein, Nugent & Goobar 2019; Wojtak et al. 2019).

3.1 LSST observing strategy

LSST will operate several survey modes. The main programme, comprising ∼90  per cent of observing time, will be the Wide-Fast-Deep (WFD) survey, consisting of an area of 18 000 deg². 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⁴, 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 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 per cent 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 per cent 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 u-band sensitivity and increased sensitivity in the grizy bands. Since lensed SNIa are red and our analysis does not consider u-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 yr 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.

Figure 4.

LSST survey footprint in equatorial coordinates showing the number of visits (N_visits) conducted per sky location. The large background area corresponds to the full LSST footprint, including the Galactic plane and polar regions. The darker region represents the WFD area (∼18 000 deg²) and the light small patches are the DDFs. In the first 1.5 yr 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 N_visits and active regions with higher N_visits.

3.2 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-yr duration of the planned survey (Delgado et al. 2014; Delgado & Reuter 2016; 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’ data base, where each pointing forms a row in the data base. 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 yr of LSST observations. This serves as a proxy for opsim’s distinction between the Galactic plane and WFD regions. We assume that regions with N_visits,10yr < 400 belong to the Galactic plane and polar regions, N_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 data base, we obtain the observing times, filters, mean 5σ depth (m₅), 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σ depth in the following way:

$$\begin{eqnarray} {\sigma _{f_{\rm SN}}} = \frac{10^{(\textrm {ZP}-{m_5})/2.5}}{5}, \end{eqnarray}$$

(9)

where ZP corresponds to the instrument zero-point for a given band: 28.38, 28.16, 27.85, 27.46, 26.68, respectively for the g, r, i, z, y bands.⁵ 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 f_SN are combined with the SN flux noise to calculate the error on the observed magnitude:

$$\begin{eqnarray} \left| \frac{-2.5 {\sigma _{f_{\rm SN}}}}{f_{\rm SN} \textrm {ln}(10)} \right|. \end{eqnarray}$$

(10)

4 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.

4.1 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 2010), 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 104 000 for the 10 yr sample (The LSST Dark Energy Science Collaboration 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 arcsec and smaller than 4 arcsec.

• 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:

  $$\begin{eqnarray} m_{\rm X} (t) \lt \langle M_{\rm X} \rangle (t_{\rm peak}) + \mu (z_{\rm lens}) + K_{\rm XX}(z_{\rm lens}, t_{\rm peak}) + \Delta m, \end{eqnarray}$$

  (11)

  with m_X(t) the apparent magnitude of the transient in band X at time t, 〈M_X〉(t_peak) the absolute magnitude of a standard SNIa in band X at peak, μ the distance modulus, K_XX the K-correction, and Δm 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).

4.2 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 (g, r, i, z, y) and all colour combinations (g − r, g − i, g − z, g − y, r − i, r − z, r − y, i − z, i − y, and z − y), 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 r or i 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 Figs 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 r and i 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 per cent contour for the unlensed SNe in each colour–magnitude space and fitting a linear function to this contour to extract a simple linear cut.

4.3 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 θ_E > 0.5 arcsec. 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 follow-up observations, which can be carried out after the SN has faded away. The difference in returned t₀ values provides the predicted time delay between the images, Δ_t,pred. We classify a system as having a ‘good’ time delay when |$\left| \Delta _{t, \textrm {pred}}-\Delta _{t, \textrm {true}} \right|/\Delta _{t, \textrm {true}} \, \lt \, 5~{{\ \rm per\ cent}}$|⁠. An example case is shown in Fig. 5. We only use LSST data for Δt 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.

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 r, i, z, and y bands and the solid curves show the salt3 fits used to infer the time delay. This object is at z_src = 0.464, lensed by a deflector galaxy at z_lens = 0.142, an input time delay of 11.14 d and a recovered time-delay estimate of 12.23 d. The Einstein radius of the system is θ_E = 0.63 so the images are treated as resolved.

4.4 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 follow-up 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:

• N_premax > 5 in at least two filters,

• m_i < 22.5 mag,

• Δt > 10 d,

with N_premax the number of detections with signal-to-noise ratio > 3 before the SN peak and m_i the apparent i-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 high-resolution 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. 4-metre Multi-Object Spectroscopic Telescope (4MOST) (Swann et al. 2019) or e.g. the New Technology Telescope (Snodgrass et al. 2008) with exposure times ≲1 h, corresponding to a brightness of m_i < 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).

5 RESULTS

5.1 Annual lensed SNIa detections

We simulate a set of 5000 doubly imaged lensed SNe and 5000 quadruply imaged SNe that pass either the image multiplicity or the magnification method as described in Section 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 per cent of doubles and 70 per cent 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 yr 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 (N_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 N_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 N_detect above a certain threshold. The non-rolling cadence will discover more lensed SNe up to N_detect < 20, while the rolling cadence will find a larger sample with well-sampled light curves (N_detect > 20). Nevertheless, we note that the differences are relatively small. We also compute the number of lensed SNIa that fall in the 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.

Figure 6.

Number of detections (N_detect) per lensed SNIa for the non-rolling cadence (first 1.5 yr 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 N_detect and with large N_detect. The mean number of observations is 15 (17) for a non-rolling (rolling) cadence.

Figure 7.

Annual expected number of lensed SNIa with N_detect above a certain threshold, for observing strategies without a rolling cadence (blue/upper line on the left) and with rolling cadence (orange/lower line on the left). 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.

Our findings are summarized 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.

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 Section 5.2), are detected with the magnification (mag.) method and image multiplicity (im. mult.) method, have time delays larger than 10 d, pass the colour–magnitude cut (described in Section 5.3), and are in the gold sample for follow-up (Section 5.5).

Doubles	Non-rolling	Rolling
Detected without microlensing	36	31
Detected	32	27
Detected with mag. method	13	12
Detected with im. mult. method	24	19
With Δt > 10 d	22	18
Pass colour–mag cut	13	11
Gold sample	8	6
Quads
Detected without microlensing	18	17
Detected	17	16
Detected with mag. method	13	13
Detected with im. mult. method	13	11
With Δt > 10 d	8	8
Pass colour–mag cut	7	7
Gold sample	5	4
Total
Detected without microlensing	54	48
Detected	50	44
Detected with mag. method	26	25
Detected with im. mult. method	37	31
With Δt > 10	30	26
Pass colour–mag cut	22	19
Gold sample	13	10

5.2 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

• do not change the detectability of the lensed SN;

• make a detected lensed SN undetectable; and

• make an undetected lensed SN detectable.

Fig. 8 investigates these three scenarios. It shows the difference in apparent peak i-band magnitude for the 5000 simulated doubly imaged 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 d, 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 understood when looking at the projected one-dimensional 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).

Figure 8.

Impact of microlensing on the lensed SNIa detections. The scatterplot shows the difference in apparent peak i-band magnitude (Δm_i) due to microlensing for image one (first occurring) and two of the 5000 simulated doubly imaged SNe. The light grey background points correspond to lensed SNe whose detection is not impacted by microlensing, the red points (mostly in upper right quadrant) are the systems that have become undetectable by microlensing, and the green ones (in upper left and lower quadrants) have become detectable by a microlensing magnification boost. The projected one-dimesional distributions show the microlensing effect on the apparent magnitude per lensed SN image.

5.3 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 Section 4.2. The best separation between lensed and unlensed SNIa is achieved with the r − z peak colour versus observed apparent z-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 the Zwicky Transient Facility (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).

Figure 9.

Peak colours and magnitudes of the detected lensed SNIa (doubles and quads) and unlensed SNIa. The r − z colours and z-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 |$r-z \, \gt\, 0.88 \, m_z-19.1$|⁠.

We investigate each colour and magnitude combination at multiple epochs and obtain the following linear cuts using the method described in Section 4.2:

$$\begin{eqnarray} g-r &\gt & 0.84 m_ r -16.9 \\g-i &\gt & 0.69 m_ i -13.5 \\g-z &\gt & 1.41 m_ z -29.4 \\g-y &\gt & 1.14 m_ y -23.6 \\r-i &\gt & 0.44 m_ i -9.3 \\r-z &\gt & 0.88 m_ z -19.1 \\r-y &\gt & 1.0 m_ y -21.7 \\i-z &\gt & 0.58 m_ z -12.6 \\i-y &\gt & 0.82 m_ y -17.9 \\z-y &\gt & 0.15 m_ y -2.9 \, . \end{eqnarray}$$

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 per cent (43 per cent) of the lensed doubles (quads) pass at least one of the mentioned colour–magnitude cuts. 2–3 per cent 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 r − z peak colour versus observed apparent z-band peak magnitude, which keeps 23 per cent (26 per cent) of the doubles (quads) with only 0.7 per cent of the unlensed sample. 28 per cent (32 per cent) of the lensed doubles (quads) would pass r − i peak colour versus observed apparent i-band peak magnitude, but with a almost 2 per cent of the unlensed events. Since unlensed SNe outnumber lensed ones by a factor of ∼10³ we want to keep the false-positive rate low. The combination with the lowest contaminants is r − y peak colour versus observed apparent y-band peak magnitude, with only 0.1 per cent unlensed events, but due to the poorer cadence and depth for the y-filter, requiring detections in y means we only preserve 15 per cent (20 per cent) 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.

5.4 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 x₁ and c is unsuccessful. Less than 2 per cent of our detected lensed SNIa sample allows for a Δt measurement with accuracy better than 5 per cent. Additionally, in Hayes et al. (2023) our lensed SNIa simulations are fitted with Gaussn⁶, a Bayesian semiparametric Gaussian Process time-delay estimation model. Gaussn recovers for |$\sim 13~{{\ \rm per\ cent}}$| of the resolved double systems the time delays within 5 per cent accuracy, corresponding to 6 per cent of the total detected sample. These results emphasize 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 Δt measurements with follow-up observations is discussed in Section 5.5.

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 Section 4.4, we find that 25 per cent 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.

Figure 10.

Normalized 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 (Δt_max), apparent i-band magnitude at peak (m_i), and the number of detections before peak (N_premax). The non-shaded regions indicate the objects that satisfy the conditions for the ‘gold’ sample.

To estimate the cosmological prospects of such a gold sample, we assume the availability, for each system in the sample, of ground-based 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 follow-up, will be ∼0.1–0.5 d, 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 per cent. For the lens mass model, we assume an uncertainty of 7 per cent. 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 per cent in H₀ for each system in the ‘gold’ sample. Consequently, we would need 30 lensed SNIa to reduce the uncertainty to 1.5 per cent in H₀, corresponding to ∼3 yr of LSST observations. Furthermore, we note that lensed SNe with shorter time delays (e.g. 5 < Δt < 10 d) will also contribute to improving the precision, even though the individual uncertainties per system would be greater than 8.6 per cent. In that case, we could reach the expected precision in a shorter duration of the LSST survey.

6 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 (2010) 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 |$\sim 40~{{\ \rm per\ cent}}$| 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 per cent of the truth value is expected to be measured for only a small fraction, |$\sim 2~{{\ \rm per\ cent}}$| of the systems. Hence, it is important to assess the feasibility of time-delay and H₀ 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 d to enable time delay cosmography. Assuming uncertainties of 8.6 per cent in H₀ per object, this sample is expected to enable a Hubble constant measurement of 1.5 per cent precision in 3 yr of LSST observations.

Our results only focus on SNIa; the expected number of lensed core-collapse SNe is likely even higher (Oguri & Marshall 2010; Goldstein & Nugent 2017; Wojtak et al. 2019). 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. 2024) 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.

ACKNOWLEDGEMENTS

This paper has undergone internal review in the LSST Dark Energy Science Collaboration. The authors would like to thank Simon Birrer, Philippe Gris, and Peter Nugent for their helpful comments and reviews. We are also grateful to Henk Arendse for data base advice, Christian Setzer for help with simulations, and Edvard Mörtsell, D’Arcy Kenworthy, Richard Kessler, and Luke Weisenbach for useful discussions.

Author contributions are listed below.

NA: conceptualization, methodology, software (lensed SN simulations), formal analysis, writing (original draft; review and editing), visualization;

SD: methodology, software, validation, formal analysis (time-delay measurements and gold sample), writing (original draft), visualization;

ASC: methodology, software, validation, formal analysis (colour–magnitude diagrams), writing (original draft), visualization;

HVP: conceptualization, validation, writing (review and editing), supervision, funding acquisition;

AG: validation, writing (review), funding acquisition;

RW: software, validation, writing (review);

CA: validation;

RB: software (OpSim Summary), writing (original draft);

SH: software (microlensing simulations), writing (original draft);

SB: writing (collaboration internal review).

This work has been enabled by support from the research project grant ‘Understanding the Dynamic Universe’ funded by the Knut and Alice Wallenberg Foundation under Dnr KAW 2018.0067. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 101018897 CosmicExplorer). The work of HVP was partially supported by the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine. This work was partially enabled by funding from the UCL Cosmoparticle Initiative. SD acknowledges support from the Marie Curie Individual Fellowship under grant ID 890695 and a Junior Research Fellowship at Lucy Cavendish College. SB thanks Stony Brook University for support. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 771776).

The DESC acknowledges ongoing support from the Institut National de Physique Nucléaire et de Physique des Particules in France; the Science & Technology Facilities Council in the United Kingdom; and the Department of Energy, the National Science Foundation, and the LSST Corporation in the United States. DESC uses resources of the IN2P3 Computing Center (CC-IN2P3 – Lyon/Villeurbanne – France) funded by the Centre National de la Recherche Scientifique; the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231; STFC DiRAC HPC Facilities, funded by UK BEIS National E-infrastructure capital grants; and the UK particle physics grid, supported by the GridPP Collaboration. This work was performed in part under DOE Contract DE-AC02-76SF00515.

Software:astropy (Astropy Collaboration 2013, 2018, 2022), jupyter (Kluyver et al. 2016), matplotlib (Hunter 2007; Caswell et al. 2020), numpy (Harris et al. 2020), pandas (Wes McKinney 2010; Team 2023), pickle (Van Rossum 2020), scipy (Virtanen et al. 2020), seaborn (Waskom et al. 2020), sqlite (Hipp 2020), chainconsumer (Hinton 2016), lenstronomy (Birrer & Amara 2018; Birrer et al. 2021b), sncosmo (Barbary et al. 2016), salt3 (Guy et al. 2007; Kenworthy et al. 2021), opsimsummary (Biswas et al. 2020).

DATA AVAILABILITY

The simulation code lensed Supernova Simulator Tool (lensedsst) and lensed SNIa catalogues are publicly available at https://github.com/Nikki1510/lensed_supernova_simulator_tool.

Footnotes

1

https://github.com/Nikki1510/lensed_supernova_simulator_tool

2

https://lenstronomy.readthedocs.io/en/latest/

3

https://sncosmo.readthedocs.io/en/stable/

4

https://pstn-055.lsst.io/

5

https://smtn-002.lsst.io/#change-record

6

https://github.com/erinhay/gaussn

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 Section 4.2. Fig. A1 shows the one-dimensional distributions of peak apparent magnitudes in the g, r, i, and z 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 r − z zcolours versus z-band magnitudes, but we combine all colour–magnitude combinations for the best result.

Figure A1.

Peak apparent magnitudes in the g, r, i, z, and y bands for simulated lensed and unlensed SNIa.

Figure A2.

Colours (g − r, g − i, g − z, g − y, r − i, r − z, r − y, i − z, i − y, and z − y) for simulated lensed and unlensed SNIa.

Figure A3.

Colour–magnitude cuts described in Section 5.3 used (together with r − z versus z shown in Fig. 9) to calculate the total number of lensed events that are distinguishable from unlensed SNIa.