© The Authors 2024

1 Introduction

Gravitationally lensed transients, such as quasars (QSOs) and supernovae (SNe), have emerged as novel powerful cosmic probes capable of delivering cosmological information independently of other probes (see, e.g., Treu et al. 2022 for a review). For instance, one can estimate the Hubble constant (H₀) using lensed transients without relying on additional datasets (Refsdal 1964; Refsdal & Bondi 1964; Saha et al. 2006; Oguri 2007; Bonvin et al. 2017; Wong et al. 2020; Birrer et al. 2020; Birrer & Treu 2021; Kelly et al. 2023; Grillo et al. 2024; Pascale et al. 2024); importantly, the value of H₀ remains quite insensitive to the chosen cosmological model. This capability can address the significant tension observed in H₀ calculations between early Universe probes (inverse distance ladder, Aghanim et al. 2020) and late Universe probes (distance ladder, Riess et al. 2022). With physically motivated lens mass models, the H0LiCOW team has further demonstrated that H₀ can be tightly constrained by combining constraints from multiple strongly lensed QSOs (Wong et al. 2020), thus underscoring the importance of having a large sample size.

Among variable sources, QSOs and SNe are currently the most prominent types, each possessing distinct advantages. For instance, lensed QSOs are more abundant and, as opposed to SNe, do not fade away, making them more reliable for time-delay cosmography at present. On the other hand, lensed SNe (LSNe) are poised to take the forefront in the coming decade due to a better understanding of their light curves, feasibility of improved stellar kinematic and SN host-galaxy measurements after the SN fades away (Ding et al. 2021; Suyu et al. 2024). Additionally, type Ia SNe serve as standard candles, offering direct luminosity distances to the sources, provided that magnification effects of microlensing (lensing by stars in the foreground lens galaxy) are mitigated (Foxley-Marrable et al. 2018; Weisenbach et al. 2021). However, a significant challenge stems from the rarity of LSNe, with only a handful of confirmed detections to date (Kelly et al. 2015; Goobar et al. 2017; More et al. 2017; Rodney et al. 2021; Chen et al. 2022; Goobar et al. 2023; Frye et al. 2023; Pierel et al. 2024). Furthermore, many LSNe will remain unresolved in the wide-field surveys due to their angular resolution being limited by seeing constraints (Quimby et al. 2014; Goldstein et al. 2018). Detecting these unresolved LSNe could notably enhance the sample size, thereby helping us to achieve precision cosmology. Therefore, the importance of detecting unresolved LSNe in time-domain wide-field surveys, such as the Legacy Survey of Space and Time (LSST) by the Vera C. Rubin Observatory, cannot be overstated.

Various techniques, such as the magnification method (Goldstein & Nugent 2017; Goobar et al. 2017; Wojtak et al. 2019), colour-magnitude diagram (Quimby et al. 2014; Arendse et al. 2023), and shifting centroid (Kodi Ramanah et al. 2022), can detect unresolved LSNe from photometry. As the observed light curve of an unresolved LSN is the sum of the individual image light curves, the former contains distinct features in its shape and can also be used to discern the lensed systems from the unlensed ones. This constitutes the central focus of our work¹. We developed several methods based on forward modelling for this purpose (Bag et al. 2021; Denissenya et al. 2021). In the present article, inspired by Denissenya & Linder (2022), we employ a simple 1D convolutional neural network (CNN) that leverages the shape of the light curves to detect lensed cases. We apply the CNNs to light-curve observations consistent with those achievable with LSST by Vera C. Rubin Observatory. LSST, which is deeper than the ongoing Zwicky Transient Facility (ZTF) survey (Bellm et al. 2019), is expected to observe approximately 20 000 square degrees over 10 years (LSST Science Collaboration 2009, 2017; The LSST Dark Energy Science Collaboration 2018), thereby notably increasing the potential for discovering new unresolved LSNe.

Our primary emphasis is on Type Ia supernovae (SNe Ia), given their pivotal role in cosmology as standardisable candles. We simulate mock LSST light curves using the baseline v3.2 observing strategy. The sample comprises joint light curves of unresolved LSNe Ia as positive examples and unlensed (i.e. intrinsic) SNe Ia light curves as negative examples. We employ two types of CNN: one for detection and the other intended to estimate time delays directly from the joint light curves in cases when a follow-up observation is not possible or conclusive. It is worth noting that this approach relies on the precise photometric classification of SNe Ia from other transients, for example using the algorithms proposed by Lochner et al. (2016); Alves et al. (2022), prior to applying the approach.

The article is organised as follows. In Section 2, we define the unresolved cases in LSST and compare their forecasted numbers and properties against the resolved ones. Next, in Section 3, we explain the process of simulating mock LSST light curves. Section 4 presents a brief overview of the method and data processing. Our main results are explained in Section 5. Finally, in Section 6, we conclude and discuss further scopes in this direction.

2 Unresolved SNe Ia in LSST

2.1 Unresolved light curves

The blended joint light curve of an unresolved LSN with N images is the sum of the individual image fluxes (neglecting microlensing for the moment), $$F(t)={\displaystyle {\displaystyle \sum }_{i\in \text{images}}^N}\left|{\tilde{\mu }}_i\right|ℱ\left(t-{\tilde{t}}_i\right),$$(1)

where ${\tilde{\mu }}_i$ and ${\tilde{t}}_i$ are the magnification (which can be positive or negative depending on the parity of an image) and time delay of the ith image. Note that the light curve of all the images can be described by the same function ℱ(t) (the unlensed light curve; here we ignore microlensing and differential dust extinction for simplicity) but with different magnifications and time delays. It is useful to recast the above expression in terms of the image that arrives the earliest, which we refer to as the ‘first’ image, $$F(t)={\displaystyle {\displaystyle \sum }_i^N}{\mu }_{i}f\left(t-\text{Δ}{t}_i\right),$$(2)

where $f(t)=\left|{\tilde{\mu }}_1\right|ℱ\left(t-{\tilde{t}}_1\right)$ is the light curve of the first image and $${\mu }_i\equiv \left|{\tilde{\mu }}_i/{\tilde{\mu }}_1\right|\text{\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}Δ}{t}_i\equiv {\tilde{t}}_i-{\tilde{t}}_1$$(3)

represent the magnification ratio and time delay of the i-th image with respect to the first image, respectively. By definition Δt_i, µ_i≥ 0 for all the images and naturally, µ₁ = 1 and Δt₁ = 0.

In this work, we only consider the light-curve information and aim to detect whether or not the system is lensed, that is, if the observed light curve stems from the blending of multiple images or corresponds to an unlensed SN. It is important to note that the task becomes challenging when dealing with extreme values of µ or small values of Δt. To understand this issue, let us first consider the double systems for simplicity. When the magnification ratio deviates significantly from unity (with Δt not being too high), one image dominates over the other; here, µ₂ ≪ 1 (µ₂ ≫ 1) implies that the second (first) image has negligible brightness compared to the first (second) image. On the other hand, if Δt₂ → 0, the two images coincide with each other. Therefore, in both of these conditions, it is mathematically impossible to detect the lensing nature in the unresolved light curve. This can be easily generalised for quads; if two images have extreme magnification ratios or very small time delay, then they cannot be discerned as two images from the unresolved light curve. In simpler terms, to detect lensed systems, it is essential to have at least one pair of images (which is the only pair in the case of a double) with µ ~ 𝒪(1) and Δt > a few days. Furthermore, to differentiate quads from doubles, multiple such pairs of images are required.

2.2 OM10 simulation for LSST

In this work, we concentrate on samples of type Ia SNe due to their significant cosmological utility, well-understood light curves, and ease of simulation. Using the Oguri and Marshall 2010 simulations (OM10 hereafter, Oguri & Marshall 2010) with minor updates (Oguri 2018), we forecast the number of lensed type Ia supernovae (LSNe Ia) that will remain unresolved in LSST. A key observation from Equation (1) is that (lensed but) unresolved sources will appear brighter than its individual images (if they were resolved) because light from different images can be accumulated for the unresolved cases. This in turn suggests that unresolved lensed sources could be visible at greater luminosity distances (or redshifts). Our primary focus is on the LSST i-band because LSNe Ia are at such redshifts that they are usually observed with the highest signal-to-noise ratio (S/N) in this filter (Oguri & Marshall 2010; Huber et al. 2022). Additionally, the LSST i-band presumably offers a more regular cadence, less affected by the moon phase, and typically provides better seeing conditions compared to shorter wavelengths. In the simulation, we allow all image separations and source brightness up to 30 mag in i-band to consider all possible unresolved LSNe Ia that could be observed in LSST.

We define the unresolved LSNe, detectable in LSST, with the following criteria: (i) image separation has to be smaller than a critical value above which one can resolve the images: θ < θ_crit and (ii) peak brightness (of the joint light curve) in i-band has to be brighter than a critical magnitude, m_crit. The 5σ depth (m₅) of LSST i-band is 23.3. However, we need the SN to be brighter than this limiting magnitude in order to have light-curve observations of sufficiently good quality to allow time-delay measurements with adequate accuracy. Therefore, we set the critical magnitude to be 0.7 magnitudes brighter (i.e. lower) than the limiting magnitude (Oguri & Marshall 2010; Dimitriadis et al. 2019), i.e. m_crit = 22.6.

On the other hand, we define the resolved systems as those that have at least one pair of images with angular separation larger than the critical value, θ > θ_crit, with both images being brighter than m_crit = 22.6². Figure 1 illustrates how the number of resolved (orange) and unresolved (blue) LSNe Ia varies with the critical angular separation (θ_crit) for one cumulative year of LSST observations. Therefore, this figure presents the theoretical prediction of the number of resolved and unresolved LSNe Ia expected to occur within a year (irrespective of the observing strategy) across the LSST footprint of 20 000 deg². Not all of these events will be observed due to seasonal gaps. Over its 10 year operation, LSST will, on average, cumulatively monitor each sky location for approximately 4 years. Therefore, for a realistic prediction, it is necessary to scale all the curves by a factor of approximately ~4.

The solid curves in Figure 1 assume only the macromagnifications, whereas the dashed curves are obtained by considering (rough) microlensing effects in the light curves³. For this range of θ_crit, we always obtain a few times more unresolved systems than resolved ones. Note that the sum of the number of unresolved and resolved cases is not constant, as the latter must meet stricter criteria. Specifically, for resolved cases, at least one pair of lensed images must be brighter than m_crit, while for unresolved cases, it is the combined flux that has to be brighter. In other words, when we lower θ_crit by a certain amount, some of the previously unresolved systems will now have θ > θ_crit and the number of unresolved cases will decrease; however, not all of them will be considered as resolved, resulting in smaller growth in the resolved case number.

The angular resolution of wide field surveys, such as LSST, is limited by seeing. For LSST, Oguri & Marshall (2010) adopt an angular resolution of θ_th = 0.5 arcsec − which is roughly two-thirds of the median θ_psf ≈ 0.7 arcsec − as the minimum separation that could be resolved in LSST. Following this, we also set θ_crit = θ_th = 0.5 arcsec for our estimations of resolved and unresolved cases. The final definitions of unresolved and resolved LSNe Ia are as follows:

• Resolved cases are those for which we have at least one pair of images with (i) angular separation of θ > 0.5 arcsec, and (ii) i-band magnitude of m < 22.6 for each of the two images.

• Unresolved cases are those where the angular separations between all the pairs are below 0.5 arcsec and (ii) the combined blended light curve has peak brightness m < 22.6 in i-band.

Let us first focus on the macro-lensing-based calculation; we return to the estimations considering microlensing below. Following the above definitions, we get 27.8 resolved and 71.1 unresolved LSNe Ia per cumulative year of LSST observation – these numbers are marked by filled circles on the respective solid curves in Figure 1. Therefore, we expect to see ~2.6 times as many unresolved LSNe Ia in LSST as resolved lensed cases.

Therefore, the detection of these unresolved cases is an important task.

In addition to introducing fluctuations in light curves, microlensing also either suppresses or boosts the brightness of images, while maintaining the same ensemble average as the macrolensing magnification. Intriguingly, microlensing tends to suppress brightness more frequently than boosting it (Goldstein et al. 2018; Arendse et al. 2023), resulting in fewer systems being brighter than the critical magnitude (m_crit) as compared to the estimations based on solely macrolensing magnification. This observation is evident as the dashed curves consistently lie below their macrolensed-only counterparts. However, since unresolved sources are permitted to be fainter than resolved sources (while both still need to be brighter than the limiting magnitude), the impact of microlensing effects on the number of unresolved lensed systems is less pronounced. In other words, summing up light curves of multiple images – each affected differently by microlensing – reduces the impact of microlensing on the peak brightness of the unresolved cases. Thus, the ratio of unresolved to resolved cases is further enhanced when one considers microlensing. For example, with θ_crit = 0.5, we get approximately 21.5 resolved and 68.2 unresolved systems per cumulative year, boosting the ratio to ~3.2. However, while exploring the comparison of resolved versus unresolved systems in terms of multiple statistics we restrict our analysis to systems obtained using macro-magnification only. Note that in the upcoming sections, where we aim to detect unresolved systems, we incorporate comprehensive microlensing effects into the light curves.

Now, let us return to the forecast with macro-magnificationbased calculation, which predicts 27.8 resolved and 71.1 unresolved LSNe Ia systems per effective year of LSST observation. We need to use the convergence and shear (κ, γ) values at the LSN image positions of the unresolved systems when simulating their joint light curves with the microlensing effect, as the microlensing maps depend on these (κ, γ) values. Figure 2 compares the 95% confidence levels (CLs) of the (κ,γ) distributions for the unresolved (blue contour) and resolved samples (orange contour). Note that the requirement of lower θ for the unresolved cases does not affect the (κ, γ) distribution as compared to that of the resolved cases. However, since unresolved sources are allowed to be fainter, they can accommodate lower magnifications and thus their (κ, γ) distribution extends towards higher values. We use the distribution for the unresolved samples when simulating their microlensed light curves in the following section. We present the source and lens redshift distributions of the unresolved systems in Appendix A, comparing them against those of the resolved cases.

Next, we compare the µ − Δt distribution for the resolved and unresolved cases in Figure 3. We remind the reader that we always sort the images chronologically and define the magnification ratios and time delays of the images with respect to the reference image, which arrives first in time (i.e. the first image). Consequently, a double system possesses only one set of µ and Δt, while a quad system entails three sets of µ and Δt. However, in the case of quads, the image with the largest time delay holds the greatest importance in the detection process. Therefore, we exclusively plot the µ and Δt of the largest time-delay image for each quad system in Figure 3.

The left panel of Figure 3 illustrates the 68 and 95% CLs of the µ − Δt distribution. The middle panel presents a 1D histogram of µ (specifically for the image with the largest time delay relative to the first image for quads). This histogram reveals a median µ of 0.40, with a range of (0.17,0.61) at the 68% percentile and (0.11,2.46) at the 95% percentile. Given that the first image, which acts as our reference, is usually the brightest for doubles and is brighter than the last image to appear for quads, the majority of the µ values fall below unity. The right panel displays the cumulative time-delay distribution, indicating how many systems have the largest time delay above a given t and comparing it for resolved (orange) and unresolved (blue) samples. (We emphasise that the number of systems presented here corresponds to one cumulative year of observations, ensuring a fair and straightforward comparison between resolved and unresolved systems.)

As anticipated, unresolved systems generally exhibit smaller time delays compared to their resolved counterparts, rendering them less suitable for time-delay cosmography. For the unresolved systems, the median gives Δt_max = 2.03 days. Nevertheless, if a mechanism can be devised to detect the majority of unresolved systems, regardless of their small time delays, their sheer abundance could contribute significantly to achieving precision cosmology. Although detecting these low-time- delay systems solely through light-curve data is challenging, the aim of this article is to provide a method to accomplish this.

Fig. 1 Number of unresolved LSNe Ia (blue) is compared with the number of resolved cases (orange) for one cumulative year of LSST observation at different critical angular resolution (θcrit) above which the images can be resolved. The solid curves are calculated assuming the macro-magnifications only, whereas the dashed curves consider the microlensing effects. LSST can resolve images up to θth ≈ 0.5 arcsec ~ 2/3θpsf. For this threshold, we obtain approximately 71.1 unresolved (blue dot) and 27.8 resolved (orange dot) LSNe Ia based solely on macro-magnification calculation, indicating roughly 2.6 times more unresolved systems. The ratio increases to approximately 3.2 when we account for microlensing, as evident from the dashed curves. We emphasize that these numbers correspond to one effective year of LSST observation over 20 000 deg2 in i-band. LSST will observe a sky location cumulatively for ~4 years on average during its ten-year run. Therefore, it is essential to scale these figures by a factor of ~4 to get the realistic predictions.

Fig. 2 The 95% confidence level of the (κ, γ) distribution compared for the images of resolved systems (orange contour) and unresolved systems (blue contour). The greyscale represents the macro-magnification in log-scale, clipped between [0.1,10]. Notably, the (κ, γ) distribution for the unresolved systems extends towards higher values (compared to that of the resolved systems), accommodating lower magnifications. The crosses mark the (κ, γ) values of the images used to simulate the microlensed light curves, with green and red crosses distinguishing images in doubles and quads, respectively. The diagonal dashed line denotes the κ = γ line.

3 Simulating mock blended joint light curves

In this section, we detail the process of constructing the training, validation, and test sets comprised of the light curves of unresolved LSNe Ia – consistent with the OM10 simulations – and of unlensed SNe Ia, across multiple bands.

3.1 Microlensed light curves

In the previous section, we neglect microlensing when exploring the statistics of the forecasted unresolved LSNe Ia samples. Microlensing can affect the lensed light curves of individual images differently, and these effects may also differ across filters⁴. Therefore, we must take microlensing into account when training the models.

To calculate microlensed light curves of LSNe Ia, we need (i) specific intensity profiles, which we get from theoretical SN Ia models, and (ii) microlensing magnification maps. Details on the calculation of microlensed light curves are presented in Huber et al. (2019). As in Huber et al. (2019); Suyu et al. (2020); Huber et al. (2021, 2022), we use theoretical SN Ia models from ARTIS (Kromer & Sim 2009) and magnification maps from Gerlumph (Vernardos & Fluke 2014; Vernardos et al. 2014, 2015; Chan et al. 2021). For the microlensing maps, we assume a Salpeter initial mass function (IMF), as in Huber et al. (2022) and Huber & Suyu (2024), with a mean mass of the microlenses of 0.35 M_⊙ and a resolution of 20 000 × 20 000 pixels with a total size of 20 R_ein × 20 R_ein on the source plane.

Another important factor for generating microlensing maps is the smooth matter fraction (s ∈ [0,1]) that determines what fraction of the total matter is smooth while the rest is in the form of stars at an image location; lower s would lead to more severe microlensing modifications to the macrolensed light curve of an image. The saddle images, which often lie closer to the lens centres especially for doubles, tend to have lower smooth-matter fractions. Conversely, the minima images often lie further away from the lens centres and thus typically have higher values of s. Therefore, to maintain a balance between realism and simplicity, we assign s = 0.4 for all the saddle images and randomly choose either 0.6 or 0.8 for the minima images.

Finally, we simulate the microlensed light curves for a total of 250 images, stemming from randomly selected 49 doubles and 38 quads. The κ, γ values of these 250 images are shown by cross symbols in Figure 2; with green and red crosses indicating images of doubles and quads, respectively. As evident from the figure, the (κ, γ) values of the downsampled images adequately represent the overall distribution of these parameters for all unresolved systems. For a given image, we generate a microlensing map using its (κ, γ, s) values. Subsequently, we simulate 100 microlensing realisations corresponding to 100 random locations on the microlensing map. The selection of the SN Ia model can also play an important role. To maintain a model-agnostic approach, we should incorporate a mixture of different SN Ia models during the training process. Among the four theoretical models investigated in Huber et al. (2022), the merger model (Pakmor et al. 2011, 2012) exhibits the most diverse behaviour, the other three being W7 (carbon deflagration, Nomoto et al. 1984), N100 (delayed detonation, Röpke et al. 2012; Seitenzahl et al. 2013), and sub-Ch (sub-Chandrasekhar mass detonation, Sim et al. 2010). For simplicity, while still retaining diversity in the SN Ia models, we opt to randomly choose between the W7 and merger models when simulating the light curves for each system.

When we construct the blended light curve of a lensed unresolved SN, we randomly select a system from the corresponding sample (double or quad) with microlensing simulations available. We then randomly choose one microlensing realisation for each of the images and sum up the microlensed light curves separately in each filter.

Fig. 3 Distribution of the largest time delay (Δtmax) and the corresponding magnification ratio in all the unresolved systems is shown in the left panel; the lighter and darker shades of blue show the 68 and 95% credible regions, respectively. The middle panel presents the 1D marginalised distribution of the magnification ratio. The right panel shows the cumulative distribution of Δtmax, i.e. indicating the number of systems (in one cumulative year) that have Δtmax higher than a given value – resolved (orange) vs. unresolved (blue). It is evident that unresolved systems typically exhibit low time delays, with a median value of Δtmax of ~2.03 days and approximately 10% having values above 10 days.

3.2 Unlensed light curves

While simulating the light curves for an unlensed SN Ia, we randomly select one of the lensed systems for which the simulated light curves are available. Subsequently, we obtain the intrinsic light curve by dividing the macrolensed light curve of any of the images by the corresponding macro-magnification.

3.3 LSST cadence distribution and noise implementation

We implement a realistic LSST cadence using publicly available simulations from the Rubin Operations Simulator (OpSim), which emulates the field selection and image acquisition process of LSST over its 10 year duration (Delgado et al. 2014; Delgado & Reuter 2016; Naghib et al. 2019). The main survey mode of LSST is the Wide-Fast-Deep (WFD) programme, which will cover around 18 000 square degrees, corresponding to ~90% of the observing time. Following the latest recommendations from the Survey Cadence Optimization Committee⁵, the WFD survey will employ a rolling cadence, in which alternating areas of the sky will receive more frequent visits than others. As a result, objects discovered in the high-cadence region will have more observations and, hence, a better light-curve coverage. In our simulations, we consider LSNe Ia in the WFD region as well as in the remaining sky regions: the Galactic plane, polar regions, and the Deep Drilling Fields, encompassing a total footprint of 20000 square degrees. We implement the baseline v3.2 observing strategy to simulate our mock light curves. Among the six LSST filters, we consider only ‘r’, ‘i’, ‘z’, and ‘y’ bands, as the remaining two, namely ‘u’ and ‘𝑔’ bands, have much ‘lower’ cadences.

To sample the (joint) light curve according to the observing strategy, we randomly choose 7 110 000 sky locations (which corresponds to an oversampling of 100 000 times the number of unresolved SNe per cumulative year of LSST observations) within the LSST footprint. We use OpSimSummary (Biswas et al. 2020) to retrieve all observation times, filters, and 5σ depths corresponding to a specific sky location within the field of view of the telescope. For each sky location, we then randomly choose a time within the 10 year run of LSST that represents the epoch of SN explosion. Many of these SNe will remain unobserved due to season gaps and some will be observed only for a few epochs, leading to partially monitored light curves. To filter out the suitable systems that can have light curves observed for sufficient epochs, we apply the following criteria:

1. Number of observations of >10 within 120 days from the explosion.

2. First observation within 20 days of the explosion.

3. At least one observation after 40 days post-explosion.

When a SN satisfies the above criteria in a filter, we consider that the system has ‘good’ observation in that particular filter. Therefore, we have four sets of well-observed systems in the four LSST filters that we consider: r, i, z, and y. Various statistics for these sets, derived from the observation strategy, are presented in Table 1 in two parts (each with four columns corresponding to four bands). The left part considers the four sets of good systems individually, whereas the right part presents the same statistics but for the subset of the systems that are well-observed in all four bands, that is, for the intersection of the four sets. For example, the second row presents the percentage of the unresolved systems that will be well observed in different bands in the left part, whereas the right part informs us that 6.2% will be well observed in all four bands. The union of good systems covers 18.8% of the total sample. The third row shows the mean number of observed epochs (up to 120 days from the explosion) in different bands. The next two rows present the distributions of the cadence and the largest gap in observation epochs, respectively, in terms of the mean and median of the respective quantities. The last row indicates the mean 5σ depth in different bands; these numbers should remain the same for the intersection set as well.

Overall, we find that r and i bands provide better cadence distributions and smaller maximum gaps in observational epochs, as compared to the other two remaining bands. Therefore, the number of ‘good systems’ is higher in these two (r and i) bands. Fortunately, SNe Ia exhibit higher S/N in these two bands due to their higher values of 5σ depth (m₅), making them ideal for this study.

3.4 Uncertainties

The Baseline v3.2 observing strategy provides the 5σ depth (m₅) in each filter at every epoch. The bottom row of Table 1 presents the mean m₅ in each filter. Following LSST Science Collaboration (2009), we simulate the uncertainties on the magnitudes as $${\sigma }_{\text{tot}}^2={\sigma }_{\text{sys}}^2+{\sigma }_{\text{rand}}^2,$$(4)

where we set the systematic noise to σ_sys = 0.005 (the upper limit) for all filters, the random noise is given by $${\sigma }_{\text{rand}}^2=(0.04-\gamma )x+\gamma x^2,\text{\hspace{0.17em}\hspace{0.17em}where\hspace{0.17em}\hspace{0.17em}}x={10}^{m-m_5}\text{,}$$(5)

and the filter-dependent quantity γ = 0.039 for r and i bands and 0.040 for z and y bands; see LSST Science Collaboration (2009) for details.

Note that when the true (or theoretical) magnitude is much larger than m₅, that is, the SN light curve is too faint compared to the 5σ depth, the σ_rand can be quite high. In these cases, a simple normal distribution can often incorrectly produce very bright observations (significantly brighter than the m₅ itself). To avoid this, we generate the observed magnitude by randomly choosing a value m ~ Uniform(m₅, m_true) (see Huber et al. 2022); however this is of little importance as these are typically the faintest parts of the SN light curves. The left panel of Fig. 4 presents the unresolved light curves of a typical quad system in four bands as an example.

Fig. 4 (a) Blended light curves of a typical unresolved system are shown in the left panel across four LSST bands. We adopt LSST-like cadence and noise to generate mock light curves (as detailed in Sect. 3.3). Data processing begins by smoothing these light curves using different smoothing scales. The grey curves, delineated with dashed, dash-dotted, and dotted line styles, represent the smoothed light curves corresponding to smoothing scales of 8, 13, and 16 days, respectively. Subsequently, we convert the smoothed magnitudes into fluxes and sample the light curves on a predefined set of epochs, kept fixed for each band. For simplicity, we opt for 100 epochs with 1 day spacing, starting from the initial observation epoch (set as t = 0), for all the bands. We then normalise every light curve separately by scaling the fluxes in the range [0,1]. (b) Finally, for each smoothing scale, we stack the normalised flux-versus-time light curves corresponding to the four bands side by side, thereby recasting the data into an array of 400 in length. Subsequently, we arrange three such arrays – corresponding to the three choices of smoothing scale – into three channels. Consequently, the final data size for each system becomes (400,3), as depicted in the right panel. This sample of data, of fixed size, is then fed into the CNNs.

4 Methodology: One-dimensional CNN

Inspired by Denissenya & Linder (2022), we also employ 1D CNNs to detect the lensed cases from the light curve data.

4.1 Data processing

In reality, we do not know the explosion epoch of any LSN. To avoid providing any information regarding the explosion time to the network, we set the time to zero at the epoch of the first observation of each light curve. The SNe Ia light curve sample from LSST will contain a different number of observations per system, as well as different observation epochs. Nevertheless, to be able to employ a simple 1D CNN on the mock observed light curves we need the input to the networks to have the same dimensions.

To achieve this, in addition to considering observation noise, we adopt the following data processing.

• We first smooth the observed magnitude versus time light curves using the iterative smoothing process of Shafieloo et al. (2006), which takes the uncertainties into account. Then we convert the smoothed magnitudes into the fluxes (in arbitrary units) to get flux versus time light curves. Next we interpolate (cubic spline) the smoothed light curves and extrapolate them at the tails if necessary to have the same data length; for example, N_j for the j-th filter. Next we normalise the light curves by scaling all the fluxes from 0 to 1 (for all systems, separately in all filters). At first, we want to test the method on long light curves by empirically choosing N_j = 100 for all the four filters, which means we use lightcurve data until 100 days from the first observation for all the systems across all bands. Later in the article, we reduce N_j while trying to detect the lensed system early enough so that the SNe Ia could still be visible in the follow up.

• We then concatenate (i.e. stack side by side) the normalised smoothed light curves for different bands (if available) to have an 1D array of $N_{\text{tot}}={\displaystyle {\sum }_{j\in \text{filters}}^4{N}_j}$ data points comprising light curves for all the four bands. Note that we test the classification results using both stacked multi-band light curves and single-band observations.

• Furthermore, we use three different smoothing scales − 8,13, and 16 days – and pass the 1D processed data (each having N_tot length) to the CNNs in three different channels. This choice of the smoothing scale is consistent with the cadence distribution (considering the mean value and the gaps) and is kept the same across all filter⁶. The smoothed magnitude versus time light curves, corresponding to the three smoothing scales, are shown in each of the four panels of Figure 4a. We expect each channel to provide features at three different timescales, if any, in the smoothed light curves.

• Therefore, at the end, our input data eventually take the shape (N_tot, 3). Figure 4b illustrates how we arranged the input data for the example quad system. The three rows show the three 1D input arrays, each having a length of N_tot = 400, corresponding to three different smoothing scales, δ = 8,13, and 16 days.

In addition to accounting for different numbers of observational epochs for different systems, this choice of data preprocessing handles observation noise, which would influence the final smoothed light curves at different smoothing scales.

4.2 Training and test sets

Figure 3 shows the distribution of maximum time delay and the corresponding magnification ratio for the unresolved systems in LSST. However, instead of using these distributions, we use simple uniform distributions for µ and Δt to construct balanced training sets. As explained in Section 2, low values of time delays (Δt → 0) or extreme values of magnification ratios (µ ≫ 1 or µ ≪ 1) make our goal of detecting the lensed systems using blended joint light curves challenging.

At first, we set a reasonable limit on the magnification ratio, µ ~ U[1/3,3.0]. However, to test the working limits of this approach in terms of the minimum feasible time delay, we constructed three different controlled sets, each with a different minimum Δt: (i) Δt ∈ U[5.0,24.0] for set A, (ii) Δt ∈ U[4.0,24.0] for set B, and (iii) Δt ∈ U[3.0,24.0] for set C. These sets are summarised in Table 2. We also created a separate test set (along with the validation set) for each training set with the same distributions of µ and Δt for initial tests. We keep a ratio of 0.5 : 0.25 : 0.25 between the population of unlensed, doubles, and quads, respectively, in all training and test sets used for classification.

Note that even when the macro-magnification ratio, µ, is non-extreme (e.g. around unity), microlensing can stochastically make it extreme in some cases. Inclusion of these joint light curves (comprising extreme values of µ) in the training set can confuse the training process, as the light curves can be very similar to an unlensed case. This can eventually result in lower accuracy in classification (as well as time-delay estimations). Therefore, we always scale the microlensed light curves of the images of a system such that their µ match the input values (drawn from a uniform distribution), and the i-band peak brightness matches the corresponding value in the OM10 simulation. By scaling with a constant factor, we preserve the microlensing effects (i.e. those that lead to wiggles) on each light curve. However, although microlensing should roughly preserve the ensemble average of µ, we can expect some cases to have extreme µ in reality due to microlensing. We anticipate that these systems will simply be classified as unlensed cases by the model, which is trained on a set that excludes extreme µ.

While predicting the time delays for a test set, if we use the model trained on the corresponding training set with the same Δt interval, say Δt ∈ (Δt_min, Δt_max), we will inevitably face bias in the results at the Δt boundaries. In this case, the model tends to over(under)-predict for smaller(larger) time delays. This is simply because the model will not learn to predict values lower than Δt_min and higher than Δt_max. To overcome this apparent bias issue, we train the model on systems with Δt ~ U[0.5,30.0] days (set D, which covers a larger range of Δt than sets A, B, and C) to predict the time delays for the test sets − A, B, and C. In this training set for time-delay estimation, we keep the populations of doubles and quads the same. Note that the test sets also share the same ratio.

Finally, to assess how our models – which are trained on the above sets – perform in realistic scenarios, we test them on a set with LSST-like unresolved systems (set LSST) without any restrictions on µ and Δt_true.

4.3 Network architectures

As mentioned above, this work focuses on two distinct goals: (i) classification, that is, discerning the lensed systems from the unlensed ones (if possible classifying them into unlensed, doubles, and quads) and (ii) estimation of the time delays. Therefore, we employ two types of 1D CNN designed for the above two purposes: classification and time-delay regression. However, for simplicity, we keep the architecture of the two types of network quite similar until the last layer – starting with two 1D convolution layers (each is followed by a max-pooling layer) followed by two dense layers⁷. The last layer is adapted according to the purpose of classification or regression. We emphasise that we only use simple network architectures in this work; however, a completely different (e.g. long-short term memory or LSTM) network can result in better results. Exploring all possible network architectures is beyond the scope of this work, the aim of which is simply to demonstrate the overall feasibility of classifying and analysing unresolved SN light curves with CNNs, but this could be an interesting avenue to explore in a follow-up project.

5 Results

5.1 Classification using long light curves

Initially, we use the long light curves, until 100 days from the first observation in each band.

5.1.1 3 × 3 classification: unlensed versus doubles versus quads

First we tried to classify the three types of systems – unlensed, doubles, and quads. The results are summarised in Table 3 for the three sets (A, B, and C in Table 2) with minimum ∆t being 5.0, 4.0, and 3.0 days from the top. The first row of a confusion matrix provides the percentage of the true unlensed systems that are classified as unlensed, double, and quads in the three respective columns; similarly, the next two rows provide the classification percentages of true doubles and quads. In each cell, we provide the result for all band data, considering only r-band, and only i-band data. As expected, the classification results are the best when the multi-band data are used; however, r or i band only data also result in reasonably good classification matrices.

As the lensing probability is low (as we expect only one lensed system out of ~ 10³–10⁴ randomly chosen SNe la; Oguri & Marshall 2010), we need to control the false-positive rate (FPR) very well. For example, if the lensing probability is 0.01%, a FPR of 10% gives roughly 1000 false alarms per one true detection, assuming the perfect true-positive rate (TPR) of 100%. We find that for set A (with ∆t_min = 5.0 days), the networks can classify with a FPR of as low as ~l.6%. However, the FPR increases to ~4.8% for set C, where ∆t_min is 3 days. The results mentioned above appears slight inferior as compared to those presented in Denissenya & Linder (2022), due to differences in cadence, uncertainties, and the distribution of (µ, ∆t) between LSST observation conditions considered here and the setup employed in Denissenya & Linder (2022).

Table 3 (the four bottom-right cells) clearly shows that the network confuses between doubles and quads very often. This confusion arises because we do not keep any minimum separation between the time delays of the later three images (with respect to the first image) for the quads. As a consequence, the joint light curve of a quad with three images with similar time delays, say ∆t′ from the first image, will be similar to a double with the time delay ∆t′ but with an aggregated magnification. Therefore, the network is expected to confuse between certain quads and doubles. Also, owing to this confusion, the FPRs in Table 3 exhibit a variation of ~1-2% (additional variation on the FPR, not percentage of FPR) due to stochasticity in the network training, that is, with different network initialisations. Therefore, it is more logical to classify the systems as unlensed versus lensed.

5.1.2 2 × 2 classification: lensed versus unlensed

Here we consider that all doubles and quads have the same truth label – ‘lensed’. The classification results are presented in Tables 4, 5 and 6 for the three sets – A, B, and C. The top number (in boldface) in each cell corresponds to the scenario with complete four-band light-curve data. The subsequent four numbers are obtained using single-band light-curve data, as labelled. The TPR remains at ≳90%, with an FPR of around 3%, for all three sets when we consider the four-band data. However, excluding the lower time-delay systems (going from set C to A) leads to enhanced overall accuracy, as expected.

Note that, as opposed to the y-band, the three bands r, i, and z individually produce well-controlled FPR of approximately 4-7%, while maintaining acceptable TPR of approximately 80-92%. However, it is evident that i-band light curves provide the most reliable classification because of the improved cadence and S/N. Therefore, even with single-band light-curve data from any of these three bands, we can achieve reasonably accurate classification.

Fig. 5 Classification results for the light curves truncated at different tcut values. True-positive rates (TPR; solid curves) and true-negative rates (TNR=l-FPR; dotted curves) are plotted for four different choices of tcut (note that the far-right points corresponds to the results already described in Section 5.1 and Tables 4–6). The results for sets A, B, and C are shown with blue, orange, and green colours, respectively. As we truncate the light curves at lower epochs, both TPR and TNR decrease, and therefore FPR increases.

5.2 Early detection

In the previous section, we show that we could classify the systems well, but we had to use light curves until 100 days from the first observation in each band. However, our aim is to detect SNe early enough – preferably in the rising phase of the light curve – so that the follow up can resolve them and subsequently measure the time delays. However, the nature of our method limits how early we can detect. In this section, we explore the minimum length of light-curve observations required after the first detection to confidently determine that the SN is lensed. We refer to this duration as t_cut. (The previous section provides results for t_cut = 100 days.) Therefore, we cut the light curves at an earlier epoch, t_cut days from the epoch of first observation. We keep t_cut the same for all bands⁸ and try to perform the classification separately for different choices of t_cut, such as at = 50, 40, 30, 20, and 10 days, to check how early we can detect the lensed systems with sufficient accuracy. For our goal of triggering a follow up, we only need the unlensed versus lensed classification, and do not require the 3 × 3 classification (unlensed vs double vs quad), which is more difficult.

Instead of the confusion matrix for individual sets, here we present the results in Figure 5, which shows true positive (solid curves) and true negative (dotted curves) detection rates as a function of t_cut. Indeed, the classification accuracy remains reasonably high until t_cut ~ 20 days, with a TPR of above 82%, while maintaining an FPR of below 9% for the most difficult set C. Therefore, this method shows the potential for detecting the lensed cases early enough to trigger a follow up.

5.3 Classifying LSST set

Finally, we test the models – which are trained on the controlled sets A, B, and C – on an LSST-like sample of unresolved systems (i.e. without any restrictions on µ or ∆t). We focus on unlensed vs lensed classifications to avoid any confusion in detection of doubles and quads. The FPRs of these models are already given in Tables 4–6 (top right element); for example, the FPR for these models is ~3% when data from all four bands are used. This exercise unambiguously determines the fraction of unresolved LSNe Ia in LSST that these models can detect.

For simplicity, here we limit our analysis to light curves observed over a period of 100 days, specifically focusing on systems that have been well-monitored across all four LSST bands – r, i, z, and y. The top part of Table 7 summarises what fractions of the unresolved systems are recovered by the model trained on sets A, B, and C (i.e. TPR as all the test samples are lensed). The first three rows present the TPR considering all systems, only doubles, and only quads, respectively. For all three sets, roughly ~30% of all the systems are recovered. However, the recovery rate slightly improves left to right (from set A to C) when the training set includes lower time-delay systems, which is expected because most of the LSST samples have ∆t < 2 days⁹. Training on set C(Δt_min = 3 days) produces the best recovery rate: 31.0%.

Interestingly, for all three sets, doubles are recovered at a better rate than the quads, because the former typically have larger time delays.

The bottom part presents the median time delay (∆t_med) of the systems that are recovered by the models trained on different sets. Again the three rows here consider all systems, doubles, and quads separately. For comparison, the input LSST sample has ∆t_med = 2.03, 2.37, 1.10 days for all systems, doubles, and quads, respectively. It is clear that the recovered systems have higher time delays as compared to the total input sample.

Figure 6 explicitly shows the recovery rate of the systems across different ∆t bins. The results for models trained on the three controlled sets, shown in three different colours, exhibit a generally consistent pattern with each other. However, the model trained on set C demonstrates slightly superior performance on low time-delay systems. As anticipated from Table 7, the recovery rate overall improves for systems with higher time delays. For example, for systems with ∆t < 3 days, the recovery rate ranges between 16 and 20%, whereas it exceeds ∼60% for systems with ∆t between 9 and 12 days. Moreover, approximately 80-90% of the systems with ∆t > 15 days could be recovered. We note that the recovery rate for the LSST set, even for higher time-delay systems, is typically slightly lower than in the control sets. This difference arises because LSST samples include systems with µ smaller than the lower limit on the control sets. These systems are generally more challenging to detect solely based on the shape of the light curves.

Fig. 6 Recovery rate for the unresolved systems in set LSST – which have good multi-band light curves available – divided into eight timedelay bins. Three colours – blue, orange, and green – represent the results obtained using the model trained on the control sets A, B, and C, respectively.

5.4 Time-delay estimation

Ideally, our goal is to detect the lensed systems using the light curves early enough (not too long after the the peak of the first image) so that follow-up observations can resolve the multiple images (and thus determine if they are doubles or quads) and measure the time delays. However, as we need long light-curve measurements to confidently detect the lensed cases, the timedelay measurements might not be possible for many systems. In such cases, the ability to measure time delays directly from the light curves becomes valuable, and that is the focus of this section.

We train separate networks for predicting time delays in double and quad samples. This assumes proper classification of the system (as a double or quad), either through follow-up observations or a near-perfect classification process, before directing it to the corresponding model. As mentioned above, both networks share the same architecture until the last layer, which has three nodes for quads and one node for doubles, employing a linear activation (or no activation) function.

Figure 7 compares the estimated time delays against the true values for the three sets – A, B, and C – placed in three panels. In each panel, individual ∆t estimations are represented by small dots, with orange and blue dots corresponding to time delays in doubles and quads, respectively. The green dashed line signifies ∆t_pred = ∆t_true, and it is evident that the predicted time delays roughly align with this. This indicates the feasibility of estimating time delays solely from joint light curves, even with LSST-like cadence and noise.

We further segment the time delays into ten equally spaced bins, separately for doubles and quads. For each bin, we calculate the median of the predicted time delays and the 68% confidence level around it, depicted by error bars in green for doubles and red for quads. As anticipated, predictions are more accurate for higher time delays.

It is worth noting that, even after training on a set with a broader ∆t range (set D), a slight over-prediction is observed at the lower end of ∆t, which is particularly visible in the right panel of Figure 7. This is somewhat inevitable as, when approaching very small time delays, predictions become biased due to the inability to accommodate negative ∆t.

Importantly, predictions are notably more accurate for double systems compared to quads. Double systems, having only one time delay, exhibit better predictions, while quads, with three time delays, pose a greater challenge to the network. Unlike Denissenya & Linder (2022), we do not impose a minimum time-delay separation between successive images, guided by the distribution of time delays in unresolved systems. This choice makes time-delay estimations for quads considerably more challenging than for doubles.

In each panel, we also provide the standard deviation (σ) of the residuals and the mean of their absolute values (⟨ϵ⟩), separately for the double and quad samples. Across all three sets, we observe that σ falls within the range of 3.1–3.2 days, and ⟨ϵ⟩ is approximately 2.4 days for quads. For doubles, σ is around 2.5 days, with ⟨ϵ⟩ at approximately 1.8 days. In summary, based on the controlled sets, we show in this section that indeed timedelay estimation is possible from only the joint light curves for the controlled sets.

In our final assessment, we test the model on the LSST set, where the majority of systems exhibit low values of ∆t and µ (the median values of ∆t_max and corresponding µ are 2.0 days and 0.4, respectively). The predicted time delays are plotted against the ground truths in Figure 8. The estimated time delays for the quads (blue dots) deviate substantially from the true values, reflecting the challenge in accurate prediction, since quads often entail very small time delays (with a median of ∆t_max ≈ 1.1 day).

Consequently, we abstain from presenting statistics across different ∆t bins for the quads. On the other hand, the predictions for doubles (orange dots) fare somewhat better. Nonetheless, the model exhibits a tendency to overpredict for low time-delay systems, as expected from the above analysis. Even for moderate time-delay systems, a notable bias is observed, with more values underpredicted than otherwise. This bias is visually evident from the error bars, which depict the median and 68% confidence level around it across different ∆t bins. Notably, the error bars for the LSST dataset are somewhat larger than those for the controlled sets shown in Figure 7. This discrepancy can be attributed to the distribution of µ in the LSST dataset, which significantly falls below the lower limit of training set D. Overall, this evaluation underscores the inherent difficulty in estimating time delays solely from the shape of joint light curves of unresolved systems in LSST, as these systems often exhibit low time delays and magnification ratios.

Fig. 7 Time-delay estimates from blended joint light curves compared with truths for sets A, B, and C. All predictions come from the model trained on set D – which covers a larger ∆t range ([0.5, 30] days) – in order to minimise apparent bias at ∆t boundaries. Orange and blue dots in the plots represent individual predictions for double and quad systems, respectively. The standard deviation of residuals (σ) and the mean of their absolute values (⟨ϵ⟩) are quoted in each plot separately for doubles and quads. Samples are split into ten bins, and the 68% CL around the median is shown by error bars, again separately for doubles (black) and quads (red); the two kinds of error bars are plotted with a slight offset along the x-axis for better visualisation. It is evident that time-delay predictions are more accurate for doubles compared to quads, and, as expected, for higher time delays.

Fig. 8 Time-delay estimations for the set LSST vs. ground truths. Statistics of ∆tpred for the quads (depicted by blue dots) at different bins are omitted, as the ∆t values for quads are typically very low, and we lack sufficient numbers at higher ∆t. Therefore, we present the median and 68% confidence level of ∆tpred at different ∆ttrue bins only for the doubles.

6 Summary and discussion

Gravitationally lensed supernovae (LSNe) hold significant promise as vital cosmic probes in the near future. The primary challenge of their rarity can be overcome by identifying the abundant unresolved LSNe anticipated in wide-field surveys. In this article, we show how we can leverage the lensing features in the observed light curves to discern between unresolved LSNe Ia and their unlensed counterparts in LSST using a straightforward 1D CNN approach. Joint light curves of unresolved SNe Ia comprise the positive examples, while the negative examples are represented by the intrinsic light curves of (unlensed) SNe Ia. To simulate LSST-like mock light curves, we employ the cadence distribution and 5σ depth from the baseline v3.2 observing strategy. Details of the simulation process are provided in Section 3. Out of six LSST filters, we consider four, namely r, i, z, and y, because of their better cadences.

Based on OM10 simulations, we show that LSST will observe approximately three times as many unresolved LSNe Ia as resolved ones. However, for nearly half of these systems, the maximum time delay (∆t) is expected to be less than 2 days. It is crucial to avoid the limit as ∆t approaches zero or extreme magnification ratios (µ) in the training process, as in both scenarios, the distinctive lensing signature within the joint light curves diminishes, rendering it nearly impossible for any method to distinguish lensed cases based solely on light-curve shapes. To ensure robustness in our training process, we exclude instances with low values of ∆t and µ. Initially, we focus on controlled datasets, denoted A, B, and C. These datasets adhere to specific criteria: both ∆t and µ are uniformly sampled, with µ following a uniform distribution in the range [1/3, 3] in all sets, and the minimum allowed ∆t is 5, 4, and 3 days, respectively, for sets A, B, and C. The outcomes for these controlled sets are provided below.

• We observe that the classification between unlensed, double, and quad can maintain a well-controlled false-positive rate (FPR). However, there is notable confusion while distinguishing some quads from doubles. This confusion arises due to the fact that joint light curves of certain quads may closely resemble those of doubles, particularly when the separations between time delays are not significant enough.

• Hence, the classification between unlensed and lensed cases appears more meaningful. For systems with ‘well-observed’ light curves across all four bands, we achieve a recovery rate of ∼90% with a FPR of ∼3% for all three sets. Even when considering only a single-band light curve observed among the filters r, i, and z, the true-positive rate (TPR) remains above 80%, with the FPR staying at ∼4-7%. Consequently, the confident detection of unresolved lensed systems is feasible even with just a single-band light curve.

• Early detection plays a pivotal role in initiating follow-up observations. Our demonstration illustrates that even with light curves truncated at 20 days from the initial observation epoch, the reliable detection of lensed cases is feasible with FPR ≲ 8% for all three sets. This underscores ample opportunity for triggering follow-up observations aimed at first confirming and then precisely measuring the time delays.

• In scenarios where follow-up observations are not possible, we can still estimate the time delays directly from the joint light curves provided the systems are accurately classified as doubles or quads. Time-delay estimations are notably more precise for doubles, as they possess a unique time delay¹⁰.

Finally, after training these models on controlled sets (A, B, and C), we test their performance on the full LSST sample (which has no restriction on ∆t or µ). We manage to recover nearly 30% of these systems provided their light curves are well observed in the four bands. Based on the testing conducted on the controlled sets, we are already aware that the FPR would remain around 3% for these cases. As expected, the recovery rate is higher for larger time-delay systems. For instance, we recover approximately 60% and 80-90% of the systems with time delays ranging between 9 and 12 days, and greater than 15 days, respectively. Subsequently, we apply the model trained on set D to the set LSST for time-delay estimation. The result reveals challenges in achieving accurate time-delay measurements, characterised by biases and significant uncertainties, even for doubles. This difficulty stems from the predominance of low ∆t and µ values in the LSST unresolved samples. Consequently, relying solely on joint light curve data for time-delay estimations may not yield reliable results.

Note that our method is specifically applicable to SN Ia light curves for distinguishing between lensed and unlensed SNe Ia. Therefore, this approach relies on the accurate identification of SN Ia events beforehand. Current algorithms, such as those of Lochner et al. (2016); Alves et al. (2022), demonstrate a high level of precision (as high as ∼90%, depending on the observing scenario) in detecting SNe Ia among other transients based on photometry. However, some contamination remains, and our method should be adapted to the other types of transients. Additionally, we need to ensure that unresolved LSNe Ia are detected by these algorithms, which have been exclusively developed for unlensed SNe Ia. Given that unresolved LSNe Ia are expected to exhibit short time delays (most with ≲2 days), it is likely that the majority of them will be identified. Addressing these challenges will require rigorous studies beyond the scope of the current work, and is left to future works. These issues are not only relevant to our method but also to other approaches, such as the magnification method, for detecting unresolved SNe Ia in large time-domain surveys like LSST. Therefore, this future work has broader implications for methods relying on photometry alone.

Overall, our approach of detecting the unresolved LSNe Ia solely based on the shape of the blended light curves shows great potential for the upcoming LSST survey. Consequently, the recovery rate (of all unresolved systems) is fundamentally limited by two factors. Firstly, we need well sampled light curves. Therefore, we follow the set of ‘good observation’ criteria, as described in Section 3.3, to select the well-monitored systems and leave out those that could not be sufficiently monitored, for example because of season gaps. As these criteria are mostly ad hoc, we have the flexibility to relax them, allowing the inclusion of a larger fraction of unresolved systems. This adjustment will enable us to explore the limit on these criteria, without significantly sacrificing the accuracy of the results.

Secondly, as expected, this method struggles to detect systems with low time delays, which is significant considering that most unresolved systems in LSST will likely have low time delays. Fortunately, lower time-delay systems often exhibit higher total magnification, where the traditional magnification method (and possibly the colour–magnitude diagram too) is expected to perform better. On the contrary, this method performs well for systems with time delays exceeding 5–6 days. Therefore, this approach complements the magnification method.

It is important to note that although our FPR can be as low as ≲3% under optimal conditions, it can still significantly affect the precision (the ratio of true positives among all positive outcomes) due to the inherently low lensing probability in reality. For example, for a lensing probability of 10⁻⁴ , an FPR of 3%, together with our TPR of ∼30% for the LSST dataset, we obtain a precision of ~1/1000, giving rise to about 1000 false alarms for every true detection. Thus, achieving a lower FPR is crucial. A reduction to an FPR of ~0.1% under the same conditions would significantly decrease the false alarm rate to a more manageable value of 30.

Another promising way to significantly improve the classification efficiency is by incorporating imaging and/or spectroscopic data alongside the light curves. For instance, in lensed cases, we anticipate the presence of the lensing galaxy near the SN. Moreover, in many instances, the host galaxy might be evident in the reference image. Additionally, the shape of an unresolved SN will not necessarily appear point-like; rather, it could be elongated in the difference imaging, possibly with its centroid shifting over time as different images dominate the observed light. Furthermore, for type Ia SNe, estimating the magnification is feasible as they can serve as standardisable candles. This allows us to verify whether the observed brightness aligns well with the photo-z estimation of the lens. Incorporating this supplementary information should lead to a better classification efficiency with lower FPR. Another potential path for improvement lies in the network architecture. Instead of employing a CNN, which necessitates fixed-length input data, we can use a long short-term memory (LSTM) network. LSTM networks have the advantage of being able to handle variable-length input sequences and can also accommodate flux uncertainties directly as inputs. These aspects will be further explored in subsequent research endeavors.

It is worth noting that unresolved SNe can be observed up to higher luminosity distances, presenting an opportunity to study the properties of these phenomena at high redshifts. For instance, they can be used to investigate whether or not there is any evolution in the light curves of SNe Ia.

In summary, this study showcases the viability of identifying unresolved LSNe Ia through light-curve observations with LSST using deep-learning techniques. Additionally, we illustrate the potential for extracting time delays solely from these unresolved light curves, enabling the prompt utilisation of such systems delivered by LSST, which may be overlooked by most other detection methods. Although the achieved FPR is below 3% in the optimal scenario, further exploration with more complex network architectures and inclusion of additional datasets may aid in further reducing this rate.

Data availability

The data products (including the catalog of unresolved LSNe Ia, microlensed light curves) are available at Zenedo DOI 10.5281/zenodo.13644602 The relevant codes (including ipython notebooks) for analyzing the unresolved systems and constructing the training and test sets comprising their blended joint light curves are available at https://github.com/deltasata/Unresolved_LSNeIa_in_LSST

Acknowledgements

We thank Masamune Oguri for providing the updated OM10 simulation code and for offering valuable comments. We also thank Taeho Ryu and Kaitlyn Szekerczes for their assistances in setting up the code. S.B. acknowledges the insightful discussions with Georgios Vernardos, Mikhail Denissenya and Wuhyun Sohn on various aspects of this project. S.B. acknowledges the funding provided by the Alexander von Humboldt Foundation. SHS thanks the Max Planck Society for support through the Max Planck Fellowship. 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). This research is supported in part by the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. S.S. has received funding from the European Union’s Horizon 2022 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 101105167 - FASTIDIoUS. A.S. would like to acknowledge the support by National Research Foundation of Korea NRF- 2021M3F7A1082056, and the support of the Korea Institute for Advanced Study (KIAS) grant funded by the government of Korea.

Appendix A Source and lens redshift of the unresolved systems in LSST

In Figure A.1 we compare between the resolved and unresolved LSNe Ia in LSST in terms of their source and lens redshifts in the left and right panel respectively. Unresolved systems, being observable at greater distances, can accommodate larger source redshifts compared to resolved ones. This presents a unique opportunity to investigate the properties of SNe at high redshifts, such as potential evolution in the light curve of SNe Ia. Understanding such evolution is crucial for cosmology, as Type Ia SNe serve as fundamental tools for studying the universe’s expansion. On the other hand, the need for low angular separations in unresolved systems, compared to resolved ones, requires smaller Einstein radii. This can be achieved either by having higher lens redshifts or lower lens masses for the same set of sources. The right panel of Figure A.1 explicitly shows that unresolved systems tend to have higher lens redshifts compared to resolved systems.

Figure A.2 provides a detailed comparison of lens galaxies for unresolved and resolved systems. The left panel shows the distribution of velocity dispersion (σ), which serves as a proxy for the mass of elliptical galaxies. As anticipated, unresolved systems are associated with lower-mass galaxies. The right panel presents the 2D distribution of lens galaxies in the σ vs z plane, aiming at distinguishing unresolved (blue) and resolved (orange) systems. Darker and lighter shades represent the 68% and 95% credible regions, respectively. This panel indicates that unresolved systems are likely to be associated with either lower- mass galaxies or similarly massive galaxies at higher redshifts compared to resolved systems, though some overlap exists.

Fig. A.1 The unresolved and resolved LSNe Ia in LSST have been compared in terms of their distributions of source redshift (left panel) and lens redshift (right panel). Since unresolved systems can be observed up to larger luminosity distances, they can accommodate larger source redshifts as compared to the resolved ones. Conversely, their smaller angular separations necessitate larger lens redshifts compared to the resolved systems.

Fig. A.2 The left panel compares the distribution of the velocity dispersion (σ) in the lens galaxies for unresolved (blue) and resolved (orange) systems. The right panel shows 2D distribution of lens galaxies in the velocity dispersion vs redshift plane, separately for unresolved (blue) and resolved (orange) systems. Darker and lighter shades correspond to the 68% and 95% confidence intervals, respectively. Since velocity dispersion serves as a proxy for the mass of elliptical galaxies, the combined insights from both panels indicate that lenses in unresolved systems are either lower-mass galaxies, similarly massive galaxies at higher redshifts, or a mix of both factors.

References

All Tables

All Figures

Fig. 1 Number of unresolved LSNe Ia (blue) is compared with the number of resolved cases (orange) for one cumulative year of LSST observation at different critical angular resolution (θcrit) above which the images can be resolved. The solid curves are calculated assuming the macro-magnifications only, whereas the dashed curves consider the microlensing effects. LSST can resolve images up to θth ≈ 0.5 arcsec ~ 2/3θpsf. For this threshold, we obtain approximately 71.1 unresolved (blue dot) and 27.8 resolved (orange dot) LSNe Ia based solely on macro-magnification calculation, indicating roughly 2.6 times more unresolved systems. The ratio increases to approximately 3.2 when we account for microlensing, as evident from the dashed curves. We emphasize that these numbers correspond to one effective year of LSST observation over 20 000 deg2 in i-band. LSST will observe a sky location cumulatively for ~4 years on average during its ten-year run. Therefore, it is essential to scale these figures by a factor of ~4 to get the realistic predictions.

In the text

Fig. 2 The 95% confidence level of the (κ, γ) distribution compared for the images of resolved systems (orange contour) and unresolved systems (blue contour). The greyscale represents the macro-magnification in log-scale, clipped between [0.1,10]. Notably, the (κ, γ) distribution for the unresolved systems extends towards higher values (compared to that of the resolved systems), accommodating lower magnifications. The crosses mark the (κ, γ) values of the images used to simulate the microlensed light curves, with green and red crosses distinguishing images in doubles and quads, respectively. The diagonal dashed line denotes the κ = γ line.

In the text

Fig. 3 Distribution of the largest time delay (Δtmax) and the corresponding magnification ratio in all the unresolved systems is shown in the left panel; the lighter and darker shades of blue show the 68 and 95% credible regions, respectively. The middle panel presents the 1D marginalised distribution of the magnification ratio. The right panel shows the cumulative distribution of Δtmax, i.e. indicating the number of systems (in one cumulative year) that have Δtmax higher than a given value – resolved (orange) vs. unresolved (blue). It is evident that unresolved systems typically exhibit low time delays, with a median value of Δtmax of ~2.03 days and approximately 10% having values above 10 days.

In the text

Fig. 4 (a) Blended light curves of a typical unresolved system are shown in the left panel across four LSST bands. We adopt LSST-like cadence and noise to generate mock light curves (as detailed in Sect. 3.3). Data processing begins by smoothing these light curves using different smoothing scales. The grey curves, delineated with dashed, dash-dotted, and dotted line styles, represent the smoothed light curves corresponding to smoothing scales of 8, 13, and 16 days, respectively. Subsequently, we convert the smoothed magnitudes into fluxes and sample the light curves on a predefined set of epochs, kept fixed for each band. For simplicity, we opt for 100 epochs with 1 day spacing, starting from the initial observation epoch (set as t = 0), for all the bands. We then normalise every light curve separately by scaling the fluxes in the range [0,1]. (b) Finally, for each smoothing scale, we stack the normalised flux-versus-time light curves corresponding to the four bands side by side, thereby recasting the data into an array of 400 in length. Subsequently, we arrange three such arrays – corresponding to the three choices of smoothing scale – into three channels. Consequently, the final data size for each system becomes (400,3), as depicted in the right panel. This sample of data, of fixed size, is then fed into the CNNs.

In the text

Fig. 5 Classification results for the light curves truncated at different tcut values. True-positive rates (TPR; solid curves) and true-negative rates (TNR=l-FPR; dotted curves) are plotted for four different choices of tcut (note that the far-right points corresponds to the results already described in Section 5.1 and Tables 4–6). The results for sets A, B, and C are shown with blue, orange, and green colours, respectively. As we truncate the light curves at lower epochs, both TPR and TNR decrease, and therefore FPR increases.

In the text

	Fig. 6 Recovery rate for the unresolved systems in set LSST – which have good multi-band light curves available – divided into eight timedelay bins. Three colours – blue, orange, and green – represent the results obtained using the model trained on the control sets A, B, and C, respectively.
In the text

Fig. 7 Time-delay estimates from blended joint light curves compared with truths for sets A, B, and C. All predictions come from the model trained on set D – which covers a larger ∆t range ([0.5, 30] days) – in order to minimise apparent bias at ∆t boundaries. Orange and blue dots in the plots represent individual predictions for double and quad systems, respectively. The standard deviation of residuals (σ) and the mean of their absolute values (⟨ϵ⟩) are quoted in each plot separately for doubles and quads. Samples are split into ten bins, and the 68% CL around the median is shown by error bars, again separately for doubles (black) and quads (red); the two kinds of error bars are plotted with a slight offset along the x-axis for better visualisation. It is evident that time-delay predictions are more accurate for doubles compared to quads, and, as expected, for higher time delays.

In the text

	Fig. 8 Time-delay estimations for the set LSST vs. ground truths. Statistics of ∆tpred for the quads (depicted by blue dots) at different bins are omitted, as the ∆t values for quads are typically very low, and we lack sufficient numbers at higher ∆t. Therefore, we present the median and 68% confidence level of ∆tpred at different ∆ttrue bins only for the doubles.
In the text

Fig. A.1 The unresolved and resolved LSNe Ia in LSST have been compared in terms of their distributions of source redshift (left panel) and lens redshift (right panel). Since unresolved systems can be observed up to larger luminosity distances, they can accommodate larger source redshifts as compared to the resolved ones. Conversely, their smaller angular separations necessitate larger lens redshifts compared to the resolved systems.

In the text

Fig. A.2 The left panel compares the distribution of the velocity dispersion (σ) in the lens galaxies for unresolved (blue) and resolved (orange) systems. The right panel shows 2D distribution of lens galaxies in the velocity dispersion vs redshift plane, separately for unresolved (blue) and resolved (orange) systems. Darker and lighter shades correspond to the 68% and 95% confidence intervals, respectively. Since velocity dispersion serves as a proxy for the mass of elliptical galaxies, the combined insights from both panels indicate that lenses in unresolved systems are either lower-mass galaxies, similarly massive galaxies at higher redshifts, or a mix of both factors.

In the text