11institutetext: Technical University of Munich, TUM School of Natural
Sciences, Physics Department, James-Franck-Straße 1, 85748 Garching, Germany
22institutetext: Max-Planck-Institut für Astrophysik, Karl-Schwarzschild
Straße 1, 85748 Garching, Germany 33institutetext: Academia Sinica Institute
of Astronomy and Astrophysics (ASIAA), 11F of ASMAB, No. 1, Section 4,
Roosevelt Road, Taipei 10617, Taiwan 44institutetext: Oskar Klein Centre,
Department of Physics, Stockholm University, SE-106 91 Stockholm, Sweden
55institutetext: Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France
66institutetext: Lawrence Berkeley National Laboratory, 1 Cyclotron Road,
Berkeley, CA 94720, USA 77institutetext: Berkeley Center for Cosmological
Physics, University of California, Berkeley, CA 94720, USA 88institutetext:
Korea Astronomy and Space Science Institute (KASI), 776 Daedeok-daero,
Yuseong-gu, Daejeon 34055, Korea 99institutetext: KASI Campus, University of
Science and Technology, 217 Gajeong-ro, Yuseong-gu, Daejeon 34113, Korea
1010institutetext: Inter-University Centre for Astronomy and Astrophysics,
Post Bag 4, Ganeshkhind, Pune 411007, India 1111institutetext: Kavli Institute
for the Physics and Mathematics of the Universe (IPMU), 5-1-5 Kashiwanoha,
Kashiwa-shi, Chiba 277-8583, Japan 1212institutetext: Dipartimento di Fisica,
Università degli Studi di Milano, via Celoria 16, I-20133 Milano, Italy
1313institutetext: INAF - IASF Milano, via A. Corti 12, I-20133 Milano, Italy

# Detecting unresolved lensed SNe Ia in LSST using blended light curves

Satadru Bag<EMAIL_ADDRESS>Simon Huber 1122 Sherry H. Suyu 112233
Nikki Arendse 44 Irham Taufik Andika 1122 Raoul Cañameras 221155 Alex Kim 66
Eric Linder 6677 Kushal Lodha 8899 Alejandra Melo 1122 Anupreeta More 10101111
Stefan Schuldt 12121313 Arman Shafieloo8,9

###### Abstract

Context.  Strong-gravitationally lensed supernovae (LSNe) are promising probes
for providing absolute distance measurements using gravitational lens time
delays. Spatially unresolved LSNe offer an opportunity to enhance the sample
size for precision cosmology. We predict that there will be approximately $3$
times more unresolved than resolved LSNe Ia in the Legacy Survey of Space and
Time (LSST) by the Rubin Observatory. In this article, we explore the
feasibility of detecting unresolved LSNe Ia from the shape of the observed
blended light curves using deep learning techniques, and we find that $\sim
30\%$ can be detected with a simple 1D CNN using well-sampled $rizy$-band
light curves (with a false-positive rate of $\sim 3\%$). Even when the light
curve is well-observed in only a single band among $r$, $i$, and $z$,
detection is still possible with false-positive rates ranging from $\sim
4-7\%$, depending on the band. Furthermore, we demonstrate that these
unresolved cases can be detected at an early stage using light curves up to
$\sim 20$ days from the first observation, with well-controlled false-positive
rates, providing ample opportunities for triggering follow-up observations.
Additionally, we demonstrate the feasibility of time-delay estimations using
solely LSST-like data of unresolved light curves, particularly for doubles,
when excluding systems with low time delay and magnification ratio. However,
the abundance of such systems among those unresolved in LSST poses a
significant challenge. This approach holds potential utility for upcoming
wide-field surveys, and overall results could significantly improve with
enhanced cadence and depth in the future surveys.

Aims.

Methods.

Results.

###### Key Words.:

Gravitational lensing: strong, micro – methods: data analysis – supernovae:
Type Ia supernova

##  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_{0}$) 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_{0}$ remains fairly insensitive to the chosen cosmological
model. This capability can address the significant tension observed in $H_{0}$
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_{0}$ can be tightly constrained by combining constraints
from multiple strongly lensed QSOs (Wong et al., 2020), thus underscoring the
importance of 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 do not fade away, unlike SNe, 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 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
confirmed detection 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
aiding in achieving 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 lensed SNe from photometry. Since
the observed light curve of an unresolved LSN will be 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. 111Recently, similar
efforts have been undertaken to detect unresolved lensed QSOs from light curve
observations (Shu et al., 2021; Bag et al., 2022, 2023). However, it’s
important to note that the techniques employed are entirely different due to
the distinct nature of QSO light curves. We developed several methods based on
forward modeling for this purpose (Bag et al., 2021; Denissenya et al., 2021).
In this 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 on light curves observations
consistent 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 a 10-year span
(LSST Science Collaboration et al., 2009, 2017; The LSST Dark Energy Science
Collaboration et al., 2018), thereby making the potential for discovering new
unresolved LSNe notably significant.

Our primary emphasis is on Type Ia supernovae (SNe Ia), given their pivotal
role in cosmology as standardizable candles. We simulate mock LSST light
curves using the baseline v3.2 observing strategy. The sample comprises joint
light curves of unresolved lensed SNe Ia as positive examples and unlensed
(i.e. intrinsic) SNe Ia light curves as negative examples. We employ two types
of CNNs: 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’s worth noting that this approach relies on the
precise photometric classification of SNe Ia from other transients, e.g.,
using the algorithms proposed by Lochner et al. (2016); Alves et al. (2022),
prior to applying.

The article is organized 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 lensed SN with $N$ images is
the sum of the individual image fluxes (neglecting microlensing for the
moment),

$F(t)=\sum_{i\in\rm
images}^{N}|\tilde{\mu}_{i}|\mathcal{F}(t-\tilde{t}_{i})\;,$ (1)

where $\tilde{\mu}_{i}$ and $\tilde{t}_{i}$ are the magnification (which can
be positive or negative depending upon the parity of an image) and time delay
of the $i$th image. Note that the light curve of all the images can be
described by the same function $\mathcal{F}(t)$ (the unlensed light curve;
here we ignored 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)=\sum_{i}^{N}\mu_{i}f(t-\Delta t_{i})\;,$ (2)

where $f(t)=|\tilde{\mu}_{1}|\mathcal{F}(t-\tilde{t}_{1})$ is the light curve
of the first image and

$\mu_{i}\equiv|\tilde{\mu}_{i}/\tilde{\mu}_{1}|\leavevmode\nobreak\
\leavevmode\nobreak\ {\rm and}\leavevmode\nobreak\ \leavevmode\nobreak\ \Delta
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. By definition $\Delta t_{i},\leavevmode\nobreak\
\mu_{i}\geq 0$ for all the images and naturally, $\mu_{1}=1$ and $\Delta
t_{1}=0$.

In this work, we only consider the light curve information and aim for
detecting if the system is lensed or not, i.e. 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 $\mu$ or small values of $\Delta t$. To understand this
issue let us first consider the double systems for simplicity. When the
magnification ratio deviates significantly from unity (with $\Delta t$ not
being too high), one image dominates over the other – $\mu_{2}\ll 1$
($\mu_{2}\gg 1$) implies that the second (first) image has negligible
brightness compared to the first (second) image. On the other hand, if $\Delta
t_{2}\to 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
$\mu\sim\mathcal{O}(1)$ and $\Delta t>$ a few days. Furthermore, to
differentiate quads from doubles, multiple such pairs of images are required.

###  2.2 OM10 simulation for LSST

Figure 1: Number of unresolved lensed SNe Ia (blue) is compared with the
number of resolved cases (orange) for one cumulative year of LSST observation
at different critical angular resolution ($\theta_{\rm 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 $\theta_{\rm th}\approx
0.5\leavevmode\nobreak\ \text{arcsec}\sim 2/3\theta_{\rm psf}$. For this
threshold, we obtain approximately $71.1$ unresolved (blue dot) and $27.8$
resolved (orange dot) lensed SNe 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 $\sim 4$ years on the average during
its 10 year long run. Therefore, it is essential to scale these figures by a
factor of $\sim 4$ to get the realistic predictions.

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 since LSNe Ia are at such redshifts that they are
usually observed with 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 lensed SNe Ia that
could be observed in LSST.

We define the unresolved lensed SNe, 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: $\theta<\theta_{\rm crit}$ and (ii) peak
brightness (of the joint light curve) in $i$-band has to be brighter than the
critical magnitude, $m_{\rm crit}$. The 5$\sigma$ depth ($m_{5}$) of LSST
$i$-band is $23.3$. However, we need the SN to be brighter than this limiting
magnitude in order to have sufficiently good quality light curve observations
that is essential for the time-delay measurements with adequate accuracy.
Therefore, we set the critical magnitude to be $0.7$ magnitude brighter (i.e.
lower) than the limiting magnitude (Oguri & Marshall, 2010; Dimitriadis et
al., 2019), i.e. $m_{\rm crit}=22.6$ .

On the other hand, we define the resolved systems which have at least one pair
of images with angular separation larger than the critical value,
$\theta>\theta_{\rm crit}$ and both images are brighter than $m_{\rm
crit}=22.6$ 222To have an apples to apples comparison, we don’t care if it’s a
double or a quad at this moment; this can be established by a follow up..
Figure 1 illustrates how the number of resolved (orange) and unresolved (blue)
lensed SNe Ia varies with the critical angular separation ($\theta_{\rm
crit}$) for one cumulative year of LSST observation. Therefore, this figure
presents the theoretical prediction of 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$ deg2. 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 $\sim 4$.

The solid curves in Figure 1 assume only the macro magnifications whereas the
dashed curves are obtained by considering (rough) microlensing effects in the
light curves 333OM10 simulation for 10 years produces too many lensed SNe Ia
for it to be computationally feasible to calculate microlensing for all.
Therefore, we take a representative subset of systems to generate the
microlensing magnification maps, and for each system in the OM10, we pick the
most similar system in this representative subset to approximate the
microlensing effect. For each system, we overpopulate it by a factor of 100
and take 100 random positions on its corresponding microlensing magnification
maps to scale the fluxes of the multiple images in the system. After repeating
this procedure for all systems in OM10, we renormalize by the factor of 100..
For this range of $\theta_{\rm crit}$, we always get 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_{\rm crit}$, while for unresolved
cases, it’s the combined flux that has to be brighter. In other words, when we
lower $\theta_{\rm crit}$ by a certain amount, some of the previously
unresolved systems will now have $\theta>\theta_{\rm crit}$ and the number of
unresolved cases will decrease. But all of them won’t 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) adopts angular resolution
$\theta_{\rm th}=0.5$ arcsec, which is roughly $2/3$ of the median
$\theta_{\rm psf}\approx 0.7$ arcsec, as the minimum separation that could be
resolved in LSST. Following this, we also set $\theta_{\rm crit}=\theta_{\rm
th}=0.5$ arcsec for our estimations of resolved and unresolved cases. Final
definitions of unresolved and resolved lensed SNe Ia are given below,

  * •

resolved cases: if any pair of images have (i) angular separation $\theta>0.5$
arcsec, and (ii) $i$-band magnitude $m<22.6$ for each of the two images,

  * •

unresolved cases: if the angular separation between all the pairs is 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 will get back to
the estimations considering microlensing shortly. Following the above
definitions, we get $27.8$ resolved and $71.1$ unresolved lensed SNe 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 that
$\sim 2.6$ times more unresolved lensed SNe Ia will be seen in LSST than
resolved lensed cases. Therefore, 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 macro lensing magnification. Intriguingly,
microlensing tends to suppress brightness more frequently than boosting
(Goldstein et al., 2018; Arendse et al., 2023), resulting in fewer systems
being brighter than the critical magnitude ($m_{\rm crit}$) as compared to the
estimations based on solely macro magnification. This observation is evident
as the dashed curves consistently lie below their macro-lensed-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 unresolved cases. Thus the ratio of
unresolved to resolved cases is further enhanced when one considers
microlensing. E.g., with $\theta_{\rm crit}=0.5$, we get approximately $21.5$
resolved and $68.2$ unresolved systems per cumulative year, boosting the ratio
to $\sim 3.2$. However, while exploring the comparison of resolved vs
unresolved systems in terms of multiple statistics we only stick 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 go back to the forecast with macro-magnification-only based
calculation that predicts $27.8$ resolved and $71.1$ unresolved lensed SNe Ia
systems per effective year of LSST observation. Figure 2 compares the $95\%$
CLs of the $(\kappa,\gamma)$ distributions for the unresolved (blue contour)
and resolved samples (orange contour). Note that the requirement of lower
$\theta$ for the unresolved cases does not affect the $(\kappa,\gamma)$
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 $(\kappa,\gamma)$ distribution extends towards
higher values. We need to use this distribution when simulating the
microlensed light curves in the next 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 $\mu-\Delta 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 $\mu$ and $\Delta t$, while a quad system entails three sets of $\mu$ and
$\Delta t$. However, in the case of quads, the image with the largest time
delay holds most significant importance in the detection process. Therefore,
we exclusively plot the $\mu$ and $\Delta 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
$\mu-\Delta t$ distribution. In the middle panel, 1D histogram of $\mu$
(specifically for the image with the largest time delay relative to the first
image for quads) is presented. This histogram reveals a median $\mu$ 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 brighter than the
last appearing image for quads, the majority of the $\mu$ 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
emphasize that the number of systems presented here corresponds to one
cumulative year of observation, 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 $\Delta
t_{\rm 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, this article aims to provide a method to
accomplish that.

Figure 2: The $95\%$ confidence level of the $(\kappa,\gamma)$ distribution is
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 $(\kappa,\gamma)$
distribution for the unresolved systems extends towards higher values
(compared to that of the resolved systems), accommodating lower
magnifications. The crosses mark the $(\kappa,\gamma)$ 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 $\kappa=\gamma$ line.

Figure 3: Distribution of the largest time delay ($\Delta t_{\rm max}$) 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 $\Delta t_{\rm max}$, i.e. how many systems (in
one cumulative year) have $\Delta t_{\rm max}$ 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 $\Delta t_{\rm max}$
being $\sim 2.03$ days and approximately $10\%$ having values above $10$ days.

##  3 Simulating mock blended joint light curves

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

###  3.1 Microlensed light curves

In the previous section, we neglected microlensing when exploring the
statistics of the forecasted unresolved lensed SNe Ia samples. Microlensing
can affect the lensed light curves of individual images differently, and these
effects may also differ across filters 444However, multiple studies (e.g.
Goldstein et al., 2018; Huber et al., 2021) have shown that there is an
achromatic phase in the microlensing effects up to roughly three rest-frame
weeks from the explosion.. Therefore, we must take microlensing into the
account when training the models.

To calculate microlensed light curves of lensed SNe 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). We use as in Huber et al.
(2019); Suyu et al. (2020); Huber et al. (2021, 2022) 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 as in Huber et al. (2022) and Huber & Suyu
(2024) a Salpeter initial mass function (IMF) with a mean mass of the
microlenses of $0.35M_{\odot}$ and a resolution of 20000 $\times$ 20000 pixels
with a total size of 20 $R_{\mathrm{ein}}$ $\times$ 20 $R_{\mathrm{ein}}$ on
the source plane.

Another important factor for generating microlensing maps is the smooth matter
fraction ($s\in[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 macro lensed light
curve of an image. The saddle images, that 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
$\kappa,\gamma$ values of these $250$ images are shown by cross marks in
Figure 2; with green and red crosses indicating images of doubles and quads,
respectively. It is evident that they fairly cover the $(\kappa,\gamma)$
distribution of all unresolved systems. For a given image, we generate a
microlensing map using its $(\kappa,\gamma,s)$ values. Subsequently, we
simulate $100$ microlensing realizations corresponding to $100$ random
locations on the microlensing map. The selection of SN Ia model can also play
an important role. To maintain a model-agnostic approach, we should
incorporate a mix 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 behavior compared
to the other three, namely 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.

###  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

|  systems individually well-observed in 4 bands (4 sets) |  systems well-observed in all 4 bands (intersection of the 4 sets)
---|---|---
bands | $r$ | $i$ | $z$ | $y$ | $r$ | $i$ | $z$ | $y$
good systems |  $13.0\%$ |  $14.0\%$ |  $12.2\%$ |  $9.8\%$ | 6.2 %
avg no. of epochs |  $21.4$ |  $21.6$ |  $20.8$ |  $19.9$ | $24.9$ | $26.1$ | $24.2$ | $22.3$
cadence mean & median (days) |  $7.3\leavevmode\nobreak\ \&\leavevmode\nobreak\ 4.0$ |  $7.6\leavevmode\nobreak\ \&\leavevmode\nobreak\ 5.0$ |  $8.4\leavevmode\nobreak\ \&\leavevmode\nobreak\ 5.7$ |  $11.6\leavevmode\nobreak\ \&\leavevmode\nobreak\ 5.8$ | $6.7\leavevmode\nobreak\ \&\leavevmode\nobreak\ 4.0$ | $6.7\leavevmode\nobreak\ \&\leavevmode\nobreak\ 4.8$ | $7.6\leavevmode\nobreak\ \&\leavevmode\nobreak\ 5.0$ | $10.8\leavevmode\nobreak\ \&\leavevmode\nobreak\ 5.0$
largest gap mean & median (days) |  $28.3\leavevmode\nobreak\ \&\leavevmode\nobreak\ 24.9$ |  $27.9\leavevmode\nobreak\ \&\leavevmode\nobreak\ 23.9$ |  $29.1\leavevmode\nobreak\ \&\leavevmode\nobreak\ 25.8$ |  $33.3\leavevmode\nobreak\ \&\leavevmode\nobreak\ 28.9$ | $24.5\leavevmode\nobreak\ \&\leavevmode\nobreak\ 23.0$ | $23.9\leavevmode\nobreak\ \&\leavevmode\nobreak\ 21.9$ | $25.9\leavevmode\nobreak\ \&\leavevmode\nobreak\ 23.9$ | $31.6\leavevmode\nobreak\ \&\leavevmode\nobreak\ 28.8$
mean $m_{5}$ | $23.9$ | $23.3$ | $22.7$ | $21.9$ | same as in the individual 4 sets corresponding to 4 bands

Table 1:  Statistics from observing strategy baseline v3.2 in the four LSST
filters – $r$, $i$, $z$, and $y$. We have a set of ‘good systems’ that meet
the criteria for ‘good observation’, as described in Section 3.3, for each of
these four filters. The table is divided into two parts, each with four
columns. The left part presents various statistics, derived from the observing
strategy, considering these four sets of systems well-observed in the four
bands. The right half presents the same statistics but for the systems well-
observed in all four bands, i.e., for the intersection of the four sets. For
example, the left part of the second row shows the fraction of unresolved
systems that would be well-observed in different bands, whereas the right part
indicates that $6.2\%$ of the systems will be well-observed in all four bands.
Note that, in comparison, $\sim 40\%$ of the total resolved systems will be
detected due to the season gaps in baseline observing strategy (Huber et al.,
2019). The next three rows present the average number of observation epochs,
mean $\&$ median of cadence and the largest gap in it, respectively, in the
four bands considering the four sets of systems individually in the left part
and their intersection in the right part. The last row provides the mean of
5$\sigma$ depth ($m_{5}$) in the four bands; these numbers should remain the
same for the intersection part as well.

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 $\sim 90\%$ of observing time.
Following the latest recommendations from the Survey Cadence Optimization
Committee555https://pstn-055.lsst.io/, 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 lensed SNe 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 20,000 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 ‘$g$’ bands, have much “lower” cadence.

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 observation) within the LSST footprint. We use OpSimSummary
(Biswas et al., 2020) to retrieve all observation times, filters, and
5$\sigma$ depths corresponding to a specific sky location within the field of
view of the telescope. Now for each sky location we 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. 1.

no. of observation $>10$ within $120$ days from the explosion,

  2. 2.

the first observation has to be within $20$ days from explosion,

  3. 3.

there has to be at least one observation after $40$ days from 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 4 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, i.e.
for the intersection of the 4 sets. For example, the second row presents what
percentage of the unresolved systems will be well-observed in different bands
in the left part whereas the right part informs 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 no. of observed epochs (till $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 mean and median of the respective quantities. The
last row indicates the mean 5$\sigma$ depth in different bands, these numbers
should remain the same for the intersection set as well.

All in all, 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$\sigma$ depth
($m_{5}$), making them ideal for this study.

###  3.4 Uncertainties

Baseline v3.2 observing strategy provides the 5$\sigma$ depth ($m_{5}$) in
each filter at every epoch. The bottom row of Table 1 presents the mean
$m_{5}$ in each filter. Following LSST Science Collaboration et al. (2009), we
simulate the uncertainties on the magnitudes as

$\sigma_{\rm tot}^{2}=\sigma_{\rm sys}^{2}+\sigma_{\rm rand}^{2}$ (4)

where we set the systematic noise $\sigma_{\rm sys}=0.005$ (the upper limit)
for all filters and the random noise is given by

$\sigma_{\rm rand}^{2}=(0.04-\gamma)x+\gamma x^{2}\leavevmode\nobreak\
,\leavevmode\nobreak\ {\rm where}\leavevmode\nobreak\ x=10^{m-m_{5}}$ (5)

and the filter dependent quantity $\gamma=0.039$ ($r$ and $i$ bands), $=0.040$
($z$ and $y$ bands); see LSST Science Collaboration et al. (2009) for details.

Note that when the true (or theoretical) magnitude is much larger than the
$m_{5}$, i.e. the SN light curve is too faint compared to the 5$\sigma$ depth,
the $\sigma_{\rm rand}$ can be quite high. In these cases, a simple normal
distribution can often incorrectly produces very bright observation
(significantly brighter than the $m_{5}$ itself). To avoid this, we generate
the observed magnitude by randomly choosing a value $m\sim{\rm
Uniform}(m_{5},m_{\rm true})$ (see Huber et al. (2022)); however this does not
matter much as these are typically the faintest part 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.

(a) Light curves in 4 bands

(b) 1D input to the CNN in 3 channels (shown in vertical panels)

Figure 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 Sec. 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 normalize every light curve separately by scaling the fluxes in the
range $[0,1]$. (b) Finally, for each smoothing scale, we stack the normalized
flux versus time light curves corresponding to the four bands side by side,
thereby recasting the data into an array of length $400$. Subsequently, we
arrange three such arrays, corresponding to the three choices of smoothing
scales, into three channels. Consequently, the final data size for each system
becomes $(400,3)$, as depicted in the right panel. This data, of fixed size,
is then fed into the CNNs.

##  4 Methodology: 1d CNN

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

###  4.1 Data prepossessing

In reality, we don’t know the explosion epoch of any lensed SN. To avoid
providing any hint for 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 vs 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 unit) to get flux vs 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, e.g. $N_{j}$ for the
$j$-th filter. Next we normalize 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 4 filters; that means we use light curve data till
$100$ days from the first observation for all the systems across all bands.
Later, we reduce $N_{j}$ while trying to detect the lensed system early enough
so that the SNe Ia could be still visible in the follow up.

  * •

Then we concatenate (i.e. stack side by side) the normalized smoothed light
curves for different bands (if available) to have an 1D array of $N_{\rm
tot}=\sum_{j\in\rm filters}^{4}N_{j}$ data points comprising light curves for
all the 4 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_{\rm 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 filters. 666The selection of the smoothing scales is
important because larger scales can wipe out small-scale features, if present,
while smaller scales may overfit and struggle with long gaps in the
observation epochs. We have tested multiple combinations of smoothing scales
and found this choice to perform the best for the LSST cadence. The smoothed
magnitude vs time light curves, corresponding to the three smoothing scales,
are shown in each of the four panels of Figure 4(a). We expect that each
channel will provide features at three different time scales, if any, in the
smoothed light curves.

  * •

Therefore, at the end, our input data end up having the shape $(N_{\rm
tot},3)$. Figure 4(b) illustrates how we arranged the input data for the
example quad system. The three rows show the three 1D input arrays, each
having the length of $N_{\rm tot}=400$, corresponding to three different
smoothing scale, $\delta=8,13,16$ days.

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

###  4.2 Training and test sets

set name | $\Delta t_{\rm true}$ |  use in classification |  use in $\Delta t$ estimation
---|---|---|---
set A | $U[5.0,24.0]$ | train + test | test only
set B | $U[4.0,24.0]$
set C | $U[3.0,24.0]$
set D | $U[0.5,30.0]$ | – | train only
set LSST | LSST like | test only | test only
Table 2: Description of different sets used for training and testing. Note
that, except for set LSST, in rest of the controlled sets we use $\mu\sim
U[1/3,3.0]$ to have a balanced training. Set D is only used to train model for
predicting time delays. Set LSST represents the whole sample of unresolved
systems in LSST without any cuts being imposed. We use this set to test our
models trained on different sets depending on our purpose of classification or
time-delay estimation. Note that we randomly choose between W7 and merger
models while generating the light curves for each system (lensed or unlensed)
in all the 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 $\mu$ and
$\Delta t$ to construct balanced training sets. As explained in Section 2, low
values of time delays $(\Delta t\to 0)$ or extreme values of magnification
ratios ($\mu\gg 1$ or $\mu\ll 1$) makes our goal of detecting the lensed
systems using blended joint light curves challenging.

At first, we set a reasonable limit on the magnification ratio, $\mu\sim
U[1/3,3.0]$. However, to test the limit on minimum time delay this approach
can work on, we construct three different controlled sets, each with a
different minimum $\Delta t$: (i) $\Delta t\in U[5.0,24.0]$ for set A, (ii)
$\Delta t\in U[4.0,24.0]$ for set B and (iii) $\Delta t\in U[3.0,24.0]$ for
set C. These sets are summarized in Table 2. We also create a separate test
set (along with the validation set) for each training set with the same
distributions of $\mu$ and $\Delta t$ for initial testings. 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, $\mu$, is non-extreme (e.g.
around unity), microlensing can stochastically make it extreme in some cases.
Inclusion of these joint light curves (comprising extreme $\mu$’s) 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 $\mu$’s
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 put wiggles) on each light curve. However, although microlensing
should roughly preserve the ensemble average of $\mu$, we can expect some
cases to have extreme $\mu$’s 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 $\mu$’s.

While predicting the time delays for a test set, if we use the model trained
on the corresponding training set with the same $\Delta t$ interval, say
$\Delta t\in(\Delta t_{\rm min},\Delta t_{\rm max})$, we will inevitably face
bias in the results at the $\Delta t$ boundaries. In this case the model tends
to over (under) predict for smaller (larger) time delays. This is simply
because the model won’t learn predicting lower than $\Delta t_{\rm min}$ and
higher than $\Delta t_{\rm max}$ values. To overcome this apparent bias
issues, we train the model on systems with $\Delta t\sim U[0.5,30.0]$ days
(set D, which covers larger range of $\Delta t$ than that of set A, B, 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,
performs in realistic scenarios, we test them on a set with LSST like
unresolved systems (set LSST) without any restrictions on $\mu$ and $\Delta
t_{\rm true}$.

###  4.3 Network architectures

As mentioned before, this work focuses on two distinct goals – (i)
classification: discern the lensed systems from the unlensed ones (if possible
classify between unlensed, doubles and quads) and (ii) estimation of the time
delays. Therefore, we employ two types of 1D CNNs 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 2 dense layers.777When we try to optimize the
architecture parameters (such as the number of convolution layers, filter
sizes and numbers and the number of dense layers) using keras.tuner keeping
the overall structure the same, we could only find minute improvements.
Therefore, for simplicity, in this article we keep the architecture parameters
the same in two types of networks until the last layer. The last layer is
adapted according to the purpose of classification or regression. We emphasize
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 demonstrating the overall feasibility of classifying and
analysing unresolved SN light curves with CNNs, but it could be a nice follow
up project.

##  5 Results

###  5.1 Classification using long light curves

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

####  5.1.1 $3\times 3$ classification: unlensed vs doubles vs quads

| Predictions
---|---
| unlensed | double | quad
Truth | unlensed |  $\mathbf{all:98.4\%}$ $r$: $93.8\%$ $i$: $96.8\%$ |  $\mathbf{all:1.3\%}$ $r$: $4.7\%$ $i$: $2.6\%$ |  $\mathbf{all:0.3\%}$ $r$: $1.5\%$ $i$: $0.6\%$
double |  $\mathbf{all:14.0\%}$ $r$: $25.8\%$ $i$: $20.2\%$ |  $\mathbf{all:70.2\%}$ $r$: $54.5\%$ $i$: $55.1\%$ |  $\mathbf{all:15.8\%}$ $r$: $19.7\%$ $i$: $24.7\%$
quad |  $\mathbf{all:2.4\%}$ $r$: $9.2\%$ $i$: $8.2\%$ |  $\mathbf{all:9.9\%}$ $r$: $18.5\%$ $i$: $19.5\%$ |  $\mathbf{all:87.7\%}$ $r$: $72.3\%$ $i$: $72.3\%$

set A: $\Delta t_{\rm true}\sim U[5.0,24.0]\leavevmode\nobreak\ {\rm
and}\leavevmode\nobreak\ \mu_{\rm true}\sim U[1/3,3.0]$

| Predictions
---|---
| unlensed | double | quad
Truth | unlensed |  $\mathbf{all:97.0\%}$ $r$: $93.9\%$ $i$: $96.6\%$ |  $\mathbf{all:2.6\%}$ $r$: $4.2\%$ $i$: $2.7\%$ |  $\mathbf{all:0.4\%}$ $r$: $1.9\%$ $i$: $0.7\%$
double |  $\mathbf{all:14.0\%}$ $r$: $28.3\%$ $i$: $21.9\%$ |  $\mathbf{all:69.6\%}$ $r$: $52.8\%$ $i$: $52.0\%$ |  $\mathbf{all:16.4\%}$ $r$: $18.9\%$ $i$: $26.1\%$
quad |  $\mathbf{all:3.4\%}$ $r$: $10.3\%$ $i$: $9.6\%$ |  $\mathbf{all:11.1\%}$ $r$: $19.5\%$ $i$: $19.5\%$ |  $\mathbf{all:85.5\%}$ $r$: $70.2\%$ $i$: $70.9\%$

set B: $\Delta t_{\rm true}\sim U[4.0,24.0]\leavevmode\nobreak\ {\rm
and}\leavevmode\nobreak\ \mu_{\rm true}\sim U[1/3,3.0]$

| Predictions
---|---
| unlensed | double | quad
Truth | unlensed |  $\mathbf{all:95.2\%}$ $r$: $93.4\%$ $i$: $96.0\%$ |  $\mathbf{all:3.9\%}$ $r$ : $4.7\%$ $i$: $2.8\%$ |  $\mathbf{all:0.9\%}$ $r$: $1.9\%$ $i$: $1.2\%$
double |  $\mathbf{all:16.8\%}$ $r$: $30.0\%$ $i$: $24.2\%$ |  $\mathbf{all:65.0\%}$ $r$: $49.1\%$ $i$: $47.9\%$ |  $\mathbf{all:18.2\%}$ $r$: $20.9\%$ $i$: $27.9\%$
quad |  $\mathbf{all:3.6\%}$ $r$: $12.0\%$ $i$: $9.8\%$ |  $\mathbf{all:11.8\%}$ $r$: $17.4\%$ $i$: $17.6\%$ |  $\mathbf{all:84.6\%}$ $r$: $70.6\%$ $i$: $72.6\%$

set C: $\Delta t_{\rm true}\sim U[3.0,24.0]\leavevmode\nobreak\ {\rm
and}\leavevmode\nobreak\ \mu_{\rm true}\sim U[1/3,3.0]$

Table 3: Confusion matrices for the classification of unlensed, double and
quad systems. The results obtained from multi-band ($rizy$) light curves are
labeled as ’all,’ while those from individual $r$ and $i$ bands are labeled
accordingly.

| Predictions
---|---
| unlensed | lensed
Truth | unlensed |  $\mathbf{all:97.1\%}$ $r$: $92.9\%$ $i$: $95.8\%$ $z$: $93.5\%$ $y$: $86.1\%$ |  $\mathbf{all:2.9\%}$ $r$: $7.1\%$ $i$: $4.2\%$ $z$: $6.5\%$ $y$: $15.9\%$
lensed |  $\mathbf{all:7.7\%}$ $r$: $16.6\%$ $i$: $12.5\%$ $z$: $17.7\%$ $y$: $25.1\%$ |  $\mathbf{all:92.3\%}$ $r$: $83.4\%$ $i$: $87.5\%$ $z$: $82.3\%$ $y$: $74.9\%$
Table 4: set A: $\Delta t_{\rm true}\sim U[5.0,24.0]\leavevmode\nobreak\ {\rm
and}\leavevmode\nobreak\ \mu_{\rm true}\sim U[1/3,3.0]$

First we try 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 $\Delta t$ being $5.0$, $4.0$ and $3.0$ days from the
top. The first row of a confusion matrix provides what percentage of the true
unlensed systems 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 at the best when the multi-band data are used,
however, $r$ or $i$ band only data also results in reasonably good
classification matrices.

Since the lensing probability is low (as we expect only one lensed systems out
of $10^{4}-10^{5}$ randomly chosen SNe Ia), we need to control the false-
positive rate (FPR) very well. E.g. if the lensing probability is as good as
$0.01\%$, a FPR of $10\%$ gives roughly $1000$ false alarms per one true
detection. We find that for the set A (with $\Delta t_{\rm min}=5.0$ days),
the networks can classify with a FPR as low as $\sim 1.6\%$. However, FPR
increase to $\sim 4.8\%$ for the set C where $\Delta t_{\rm min}$ is $3$ days.
The results mentioned above show a slight inferiority compared to those
presented in Denissenya & Linder (2022), primarily attributed to differences
in cadence, uncertainties, and the distribution of $(\mu,\Delta t)$ between
LSST observation conditions considered here and the setup employed in
Denissenya & Linder (2022).

We readily notice (in the four bottom-right cells) of Table 3 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 having similar
time delays, say $\Delta t^{\prime}$ from the first image, will be similar to
a double with the time delay $\Delta t^{\prime}$ but with an aggregated
magnification. Therefore, it is expected that the network confuses between
certain quads and doubles. Also, owing to this confusion, the FPRs in Table 3
exhibit a variation of $\sim$$1$-$2\%$ (additional variation on the FPR, not
percentage of FPR) due to stochasticity in the network training, i.e., with
different network initializations. Therefore, it makes more sense to classify
the systems as unlensed vs lensed.

####  5.1.2 $2\times 2$ classification: lensed vs unlensed

Here we consider all doubles and quads having 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 a single-band light curve data as labeled. The
true-positive rate (TPR) stays $\gtrsim 90\%$, with the FPR 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 except for the $y$-band, the rest of 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
better cadence and signal to noise ratio. Therefore, even with single-band
light curve data from any of these three bands, we can achieve reasonably
accurate classification.

| Predictions
---|---
| unlensed | lensed
Truth | unlensed |  $\mathbf{all:96.8\%}$ $r$: $91.6\%$ $i$: $95.2\%$ $z$: $92.6\%$ $y$: $83.8\%$ |  $\mathbf{all:3.2\%}$ $r$: $8.4\%$ $i$: $4.8\%$ $z$: $7.4\%$ $y$: $16.2\%$
lensed |  $\mathbf{all:8.9\%}$ $r$: $17.2\%$ $i$: $13.9\%$ $z$: $18.9\%$ $y$: $24.8\%$ |  $\mathbf{all:91.1\%}$ $r$: $82.8\%$ $i$: $86.1\%$ $z$: $81.1\%$ $y$: $75.2\%$
Table 5: set B: $\Delta t_{\rm true}\sim U[4.0,24.0]\leavevmode\nobreak\ {\rm and}\leavevmode\nobreak\ \mu_{\rm true}\sim U[1/3,3.0]$ | Predictions
---|---
| unlensed | lensed
Truth | unlensed |  $\mathbf{all:96.7\%}$ $r$: $91.5\%$ $i$: $95.5\%$ $z$: $92.7\%$ $y$: $84.3\%$ |  $\mathbf{all:3.3\%}$ $r$: $8.5\%$ $i$: $4.5\%$ $z$: $7.3\%$ $y$: $15.7\%$
lensed |  $\mathbf{all:10.7\%}$ $r$: $19.1\%$ $i$: $15.8\%$ $z$: $21.0\%$ $y$: $26.6\%$ |  $\mathbf{all:89.3\%}$ $r$: $80.9\%$ $i$: $84.2\%$ $z$: $79.0\%$ $y$: $73.4\%$
Table 6: set C: $\Delta t_{\rm true}\sim U[3.0,24.0]\leavevmode\nobreak\ {\rm
and}\leavevmode\nobreak\ \mu_{\rm true}\sim U[1/3,3.0]$

###  5.2 Early detection

Figure 5: Classification results for the light curves truncated at different
$t_{\rm cut}$ values: true-positive rates (TPR; solid curves) and true-
negative rates (TNR=$1-$FPR; dotted curves) are plotted for 4 different
choices of $t_{\rm cut}$ (note that the far right points corresponds to the
results already described in Section 5.1 and Tables 4-6). The results for set
A, B and C are shown by the blue, orange and green colours respectively. As we
truncate the light curves at lower epochs, both TPR and TNR decrease;
therefore FPR increases.

In the previous section, we could classify the systems well but we had to use
light curves till $100$ days from the first observation in each band. However,
we want to detect them early enough, preferably in the light curve rising
phase, so that the follow up can resolve it 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_{\rm cut}$. (Previous section provides
results for $t_{\rm cut}=100$ days.) Therefore, we cut the light curves at an
earlier epoch, $t_{\rm cut}$ days from the epoch of first observation. We keep
$t_{\rm cut}$ the same for all bands 888Cutting the long processed light
curves (as in Figure 4(b)) at $t_{\rm cut}$ and feeding it to the CNN would
give a hint of the peak location which might actually be beyond $t_{\rm cut}$.
Therefore, for a given $t_{\rm cut}$, we always smooth the mag vs time light
curves considering data up to the first observation after $t_{\rm cut}$ and
then cut the smoothed light curve at $t_{\rm cut}$. Finally we normalise the
flux vs time light curves and feed it to the CNN.. We try to perform the
classification separately for different choices of $t_{\rm cut}$, e.g. 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 vs lensed classification, instead of the $3\times 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_{\rm cut}$. Indeed, the
classification accuracy remains reasonably well till $t_{\rm cut}\sim 20$ days
with TPR above $82\%$ while maintaining a FPR below $9\%$ for the most
difficult set C. Therefore, this method shows the potential for detecting the
lensed cases early enough for triggering a follow up.

###  5.3 Classifying LSST set

| model trained on
---|---
| set A | set B | set C
recovery | all | $28.0\%$ | $28.1\%$ | $30.2\%$
double | $29.5\%$ | $29.2\%$ | $31.3\%$
quad | $21.9\%$ | $23.3\%$ | $25.9\%$
$\Delta t_{\rm med}$ of recovered systems | all | $4.5$ | $4.4$ | $4.1$
double | $4.8$ | $4.7$ | $4.6$
quad | $2.6$ | $2.1$ | $1.6$
Table 7: The top part presents the percentage of mock LSST unresolved systems
recovered by models trained on different sets – A, B and C. The three sub-rows
show the recovery rate for all recovered sample and then divided it into
doubles and quads. The bottom part shows the median time delays of the
recovered systems, again divided into three sets – all, doubles and quads from
top to bottom. For comparison, the input LSST sample has $\Delta t_{\rm
med}=2.0,2.4,1.1$ days respectively for the full sample, doubles and quads. It
is evident that the models can recover the doubles better than the quads since
the former systems have typically larger time delays. Figure 6: Recovery rate
for the unresolved systems in set LSST, which have good multi-band light
curves available, divided into $8$ time-delay bins. Three colours – blue,
orange and green – represent the results obtained using the model trained on
the control set A, B and C, respectively.

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 $\mu$ and $\Delta t$). We focus on unlensed vs lensed classification to
avoid any confusion in detection of doubles and quads. The FPR of these models
are already given in Tables 4–6 (top right element); e.g. FPR for these models
is $\sim 3\%$ when all four band data are used. This exercise unambiguously
determines what fraction of unresolved LSNe Ia in LSST 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 set A, B and C (i.e. TPR as all the test samples are
lensed). The first three rows presents the TPR considering all systems, only
doubles and only quads respectively. For all three sets, roughly $\sim 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 as expected since most of the LSST samples have $\Delta
t\lesssim 2$ days.999When tested with models trained on a few other sets where
lower $\Delta t$ and/or $\mu$ systems are included, we found minimal
improvements in the recovery rate but at an expense of greater FPR. Therefore,
we do not report the results for these additional sets. Training on set C
($\Delta t_{\rm 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 since the former typically have larger time delays.

The bottom part presents the median time delay ($\Delta t_{\rm med}$) of the
systems that are recovered by the models trained on different sets. Again
three rows here consider all systems, doubles and quads separately. For
comparison, the input LSST sample has $\Delta t_{\rm
med}=2.03,\leavevmode\nobreak\ 2.37,\leavevmode\nobreak\ 1.10$ days for all
systems, doubles and quads. 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
$\Delta t$ bins. The results for models trained on three controlled sets,
shown by three 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. E.g., for
systems with $\Delta t<3$ days, the recovery rate ranges between $16-20\%$,
whereas it exceeds $\sim 60\%$ for systems with $\Delta t$ between $9$ and
$12$ days. Moreover, approximately $80-90\%$ of the systems with $\Delta t>15$
days could be recovered. 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 $\mu$
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.

###  5.4 Time-delay estimation

(a) set A

(b) set B

(c) set C

Figure 7: Time-delay estimates from blended joint light curves were compared
with truths for sets A, B, and C. All predictions comes from the model trained
on set D, which covers a larger $\Delta t$ range ($[0.5,30]$ days), to
minimize apparent bias at $\Delta t$ boundaries. Orange and blue dots in the
plots represent individual predictions for double and quad systems,
respectively. The standard deviation of residuals ($\sigma$) and mean of their
absolute values ($\left\langle\epsilon\right\rangle$) are quoted in each plot,
separately for doubles and quads. Samples are split into 10 bins, and the
$68\%$ CL around the median is shown by error bars, again separately for
doubles (black) and quads (red); the two kind of errorbars are plotted with a
slight offset along x-axis for a better visualization. It is evident that
time-delay predictions are more accurate for doubles compared to quads or at
higher time delays, as expected. Figure 8: Time-delay estimations for the set
LSST are compared against the ground truths. Statistics of $\Delta
t_{\text{pred}}$ for the quads (depicted by blue dots) at different bins are
omitted, as the $\Delta t$ values for quads are typically very low, and we
lack sufficient numbers at higher $\Delta t$. Therefore, we present the median
and $68\%$ confidence level of $\Delta t_{\text{pred}}$ at different $\Delta
t_{\text{true}}$ bins only for the doubles.

Ideally, our goal is to detect the lensed systems using the light curves early
enough (not too far from the peak of the first image) so that the follow ups
can resolve the multiple images (and thus determine if it’s a double or a
quad) and measure the time delays. However, since we need long light curve
measurements to confidently detect the lensed cases, the time-delay
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 previously
mentioned, both networks share the same architecture until the last layer,
which has three nodes for quads and one node for doubles, employing 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
$\Delta 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 $\Delta t_{\rm pred}=\Delta t_{\rm true}$, and it
is evident that the predicted time delays roughly align with this line. 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 10 equispaced 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’s worth noting that, even after training on a set with a broader $\Delta t$
range (set D), a slight over-prediction is observed at the lower end of
$\Delta t$, particularly visible in the right panel of Figure 7. This is
somewhat inevitable as, when approaching very small time delays, predictions
becomes biased due to the inability to accommodate negative $\Delta 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 for
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 ($\sigma$) of the
residuals and the mean of their absolute values
($\left\langle\epsilon\right\rangle$), separately for the double and quad
samples. Across all three sets, we observe that $\sigma$ falls within the
range of $3.1$ to $3.2$ days, and $\left\langle\epsilon\right\rangle$ is
approximately $2.4$ days for quads. For doubles, $\sigma$ is around $2.5$
days, with $\left\langle\epsilon\right\rangle$ at approximately $1.8$ days. In
summary, based on the controlled sets, we show in this section that indeed
time-delay 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 $\Delta t$ and $\mu$ (the median values of
$\Delta t_{\text{max}}$ and corresponding $\mu$ 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) deviates
substantially from the true values, reflecting the challenge in accurate
prediction, since quads often entail very small time delays (with a median of
$\Delta t_{\text{max}}\approx 1.1$ day) . Consequently, we abstain from
presenting statistics across different $\Delta 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 $\Delta 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 $\mu$ 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.

##  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. This article leverages the lensing features in the observed light
curves to discern between unresolved LSNe Ia and their unlensed counterparts
in LSST, utilizing 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$\sigma$ 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$, due to their better cadences.

Based on OM10 simulations we show that LSST will observe approximately three
times more unresolved LSNe Ia than the resolved ones. However, for nearly half
of these systems, the maximum time delay ($\Delta t$) is expected to be less
than $2$ days. It is crucial to avoid the limit as $\Delta t$ approaches zero
or extreme magnification ratios ($\mu$) 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 $\Delta t$ and $\mu$.
Initially, we focus on controlled datasets denoted as A, B, and C. These
datasets adhere to specific criteria: both $\Delta t$ and $\mu$ are uniformly
sampled, with $\mu$ following a uniform distribution in the range $[1/3,3]$ in
all sets, and the minimum allowed $\Delta t$ is 5, 4, and 3 days respectively
for set A, B, and C. The outcomes for these controlled sets are provided
below.

  * •

We have observed 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 the 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 $\sim 90\%$ with a FPR $\sim 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 $\sim 4-7\%$. Consequently, confident detection
of unresolved lensed system 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, reliable detection of lensed cases is
feasible with FPR $\lesssim 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.
101010Note that all the results presented above, listed as separate points,
are confined to the controlled sets where we exclude low $\Delta t$ and $\mu$
systems from training and testing; see Table 2.

Finally, we test the models, which are trained on these controlled sets (A, B
and C), on the full LSST sample (which has no restriction on the $\Delta t$
and $\mu$). 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, approximately $60\%$ and
$80-90\%$ of the systems with time delays ranging between 9 and 12 days, and
greater than 15 days, respectively, could be recovered. 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,
characterized by biases and significant uncertainties, even for doubles. This
difficulty stems from the predominance of low $\Delta t$ and $\mu$ 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 Lochner et al. (2016); Alves et al. (2022), demonstrate a
high level of reliability in detecting SNe Ia among other transients based on
photometry.

All in all, 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
the rest those could not be sufficiently monitored, e.g. because of season
gaps. Since these criteria were mostly ad hoc, we have the flexibility to
relax them, allowing for the inclusion of a larger fraction of unresolved
systems. This adjustment enables 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.

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
would not necessarily appear point-like; rather, it could be elongated in the
difference imaging, possibly with its centroid shifting over time as different
image dominates the observed light. Furthermore, for type Ia SNe, estimating
the magnification is feasible as they can serve as standardizable 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 input data, we can utilize 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’s worth noting that unresolved supernovae can be observed up to higher
luminosity distances, presenting an opportunity to study the properties of
supernovae at high redshifts. For instance, they can be utilized to
investigate whether there is any evolution in the light curves of type Ia
supernovae.

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

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

Figure 9: 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 Figure 9 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, their smaller
angular separations require larger lens redshifts compared to resolved
systems.

## Acknowledgement

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. SB acknowledges the insightful
discussions with Georgios Vernardos, Mikhail Denissenya and Wuhyun Sohn on
various aspects of this project. SB 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. SS 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. AS 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.

## Data availability

The catalog of unresolved lensed SNe Ia simulated using OM10 simulation
pipeline will be made public upon acceptance.

## References

  * Aghanim et al. (2020) Aghanim, N. et al. 2020, Astron. Astrophys., 641, A6
  * Alves et al. (2022) Alves, C. S., Peiris, H. V., Lochner, M., et al. 2022, ApJS, 258, 23
  * Arendse et al. (2023) Arendse, N., Dhawan, S., Sagués Carracedo, A., et al. 2023, arXiv e-prints, arXiv:2312.04621
  * Bag et al. (2021) Bag, S., Kim, A. G., Linder, E. V., & Shafieloo, A. 2021, Astrophys. J., 910, 65
  * Bag et al. (2022) Bag, S., Shafieloo, A., Liao, K., & Treu, T. 2022, Astrophys. J., 927, 191
  * Bag et al. (2023) Bag, S., Sohn, W., Shafieloo, A., & Liao, K. 2023, ApJ, 950, 158
  * Bellm et al. (2019) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002
  * Birrer et al. (2020) Birrer, S., Shajib, A. J., Galan, A., et al. 2020, A&A, 643, A165
  * Birrer & Treu (2021) Birrer, S. & Treu, T. 2021, A&A, 649, A61
  * Biswas et al. (2020) Biswas, R., Daniel, S. F., Hložek, R., et al. 2020, ApJS, 247, 60
  * Bonvin et al. (2017) Bonvin, V., Courbin, F., Suyu, S. H., et al. 2017, MNRAS, 465, 4914
  * Chan et al. (2021) Chan, J. H. H., Rojas, K., Millon, M., et al. 2021, A&A, 647, A115
  * Chen et al. (2022) Chen, W., Kelly, P. L., Oguri, M., et al. 2022, Nature, 611, 256
  * Delgado & Reuter (2016) Delgado, F. & Reuter, M. A. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9910, Observatory Operations: Strategies, Processes, and Systems VI, ed. A. B. Peck, R. L. Seaman, & C. R. Benn, 991013
  * Delgado et al. (2014) Delgado, F., Saha, A., Chandrasekharan, S., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9150, Modeling, Systems Engineering, and Project Management for Astronomy VI, ed. G. Z. Angeli & P. Dierickx, 915015
  * Denissenya et al. (2021) Denissenya, M., Bag, S., Kim, A. G., Linder, E. V., & Shafieloo, A. 2021, arXiv e-prints
  * Denissenya & Linder (2022) Denissenya, M. & Linder, E. V. 2022, MNRAS, 515, 977
  * Dimitriadis et al. (2019) Dimitriadis, G., Rojas-Bravo, C., Kilpatrick, C. D., et al. 2019, ApJ, 870, L14
  * Ding et al. (2021) Ding, X., Liao, K., Birrer, S., et al. 2021, Monthly Notices of the Royal Astronomical Society, 504, 5621
  * Foxley-Marrable et al. (2018) Foxley-Marrable, M., Collett, T. E., Vernardos, G., Goldstein, D. A., & Bacon, D. 2018, MNRAS, 478, 5081
  * Frye et al. (2023) Frye, B., Pascale, M., Cohen, S., et al. 2023, Transient Name Server AstroNote, 96, 1
  * Goldstein & Nugent (2017) Goldstein, D. A. & Nugent, P. E. 2017, ApJ, 834, L5
  * Goldstein et al. (2018) Goldstein, D. A., Nugent, P. E., Kasen, D. N., & Collett, T. E. 2018, ApJ, 855, 22
  * Goobar et al. (2023) Goobar, A., Johansson, J., Schulze, S., et al. 2023, Nature Astronomy, 7, 1098
  * Goobar et al. (2017) Goobar, A. et al. 2017, Science, 356, 291
  * Grillo et al. (2024) Grillo, C., Pagano, L., Rosati, P., & Suyu, S. H. 2024, arXiv e-prints
  * Huber & Suyu (2024) Huber, S. & Suyu, S. H. 2024, arXiv e-prints, arXiv:2403.08029
  * Huber et al. (2022) Huber, S., Suyu, S. H., Ghoshdastidar, D., et al. 2022, A&A, 658, A157
  * Huber et al. (2019) Huber, S., Suyu, S. H., Noebauer, U. M., et al. 2019, A&A, 631, A161
  * Huber et al. (2021) Huber, S., Suyu, S. H., Noebauer, U. M., et al. 2021, A&A, 646, A110
  * Kelly et al. (2015) Kelly, P. L. et al. 2015, Science, 347, 1123
  * Kelly et al. (2023) Kelly, P. L. et al. 2023, Science, 380, abh1322
  * Kodi Ramanah et al. (2022) Kodi Ramanah, D., Arendse, N., & Wojtak, R. 2022, MNRAS, 512, 5404
  * Kromer & Sim (2009) Kromer, M. & Sim. 2009, MNRAS, 398, 1809
  * Lochner et al. (2016) Lochner, M., McEwen, J. D., Peiris, H. V., Lahav, O., & Winter, M. K. 2016, ApJS, 225, 31
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, arXiv e-prints, arXiv:0912.0201
  * LSST Science Collaboration et al. (2017) LSST Science Collaboration, Marshall, P., Anguita, T., et al. 2017, arXiv e-prints, arXiv:1708.04058
  * More et al. (2017) More, A., Suyu, S. H., Oguri, M., More, S., & Lee, C.-H. 2017, ApJ, 835, L25
  * Naghib et al. (2019) Naghib, E., Yoachim, P., Vanderbei, R. J., Connolly, A. J., & Jones, R. L. 2019, AJ, 157, 151
  * Nomoto et al. (1984) Nomoto, K., Thielemann, F. K., & Yokoi, K. 1984, ApJ, 286, 644
  * Oguri (2007) Oguri, M. 2007, The Astrophysical Journal, 660, 1
  * Oguri (2018) Oguri, M. 2018, MNRAS, 480, 3842
  * Oguri & Marshall (2010) Oguri, M. & Marshall, P. J. 2010, MNRAS, 405, 2579
  * Pakmor et al. (2011) Pakmor, R., Hachinger, S., Röpke, F. K., & Hillebrandt, W. 2011, A&A, 528, A117
  * Pakmor et al. (2012) Pakmor, R., Kromer, M., Taubenberger, S., et al. 2012, ApJ, 747, L10
  * Pascale et al. (2024) Pascale, M., Frye, B. L., Pierel, J. D. R., et al. 2024, arXiv e-prints, arXiv:2403.18902
  * Pierel et al. (2024) Pierel, J. D. R., Newman, A. B., Dhawan, S., et al. 2024, arXiv e-prints, arXiv:2404.02139
  * Quimby et al. (2014) Quimby, R. M., Oguri, M., More, A., et al. 2014, Science, 344, 396
  * Refsdal (1964) Refsdal, S. 1964, Monthly Notices of the Royal Astronomical Society, 128, 307
  * Refsdal & Bondi (1964) Refsdal, S. & Bondi, H. 1964, Monthly Notices of the Royal Astronomical Society, 128, 295
  * Riess et al. (2022) Riess, A. G. et al. 2022, Astrophys. J. Lett., 934, L7
  * Rodney et al. (2021) Rodney, S. A., Brammer, G. B., Pierel, J. D. R., et al. 2021, Nature Astronomy
  * Röpke et al. (2012) Röpke, F. K., Kromer, M., Seitenzahl, I. R., et al. 2012, ApJ, 750, L19
  * Saha et al. (2006) Saha, P., Coles, J., Macciò, A. V., & Williams, L. L. R. 2006, ApJ, 650, L17
  * Seitenzahl et al. (2013) Seitenzahl, I. R., Ciaraldi-Schoolmann, F., Röpke, F. K., et al. 2013, MNRAS, 429, 1156
  * Shafieloo et al. (2006) Shafieloo, A., Alam, U., Sahni, V., & Starobinsky, A. A. 2006, Mon. Not. Roy. Astron. Soc., 366, 1081
  * Shu et al. (2021) Shu, Y., Belokurov, V., & Evans, N. W. 2021, MNRAS, 502, 2912
  * Sim et al. (2010) Sim, S. A., Röpke, F. K., Hillebrandt, W., et al. 2010, ApJ, 714, L52
  * Suyu et al. (2024) Suyu, S. H., Goobar, A., Collett, T., More, A., & Vernardos, G. 2024, Space Sci. Rev., 220, 13
  * Suyu et al. (2020) Suyu, S. H., Huber, S., Cañameras, R., et al. 2020, A&A, 644, A162
  * The LSST Dark Energy Science Collaboration et al. (2018) The LSST Dark Energy Science Collaboration, Mandelbaum, R., Eifler, T., et al. 2018, arXiv e-prints, arXiv:1809.01669
  * Treu et al. (2022) Treu, T., Suyu, S. H., & Marshall, P. J. 2022, Astron. Astrophys. Rev., 30, 8
  * Vernardos & Fluke (2014) Vernardos, G. & Fluke, C. J. 2014, Astron. Comput., 6, 1
  * Vernardos et al. (2014) Vernardos, G., Fluke, C. J., Bate, N. F., & Croton, D. 2014, ApJS, 211, 16
  * Vernardos et al. (2015) Vernardos, G., Fluke, C. J., Bate, N. F., Croton, D., & Vohl, D. 2015, ApJS, 217, 23
  * Weisenbach et al. (2021) Weisenbach, L., Schechter, P. L., & Pontula, S. 2021, ApJ, 922, 70
  * Wojtak et al. (2019) Wojtak, R., Hjorth, J., & Gall, C. 2019, MNRAS, 487, 3342
  * Wong et al. (2020) Wong, K. C. et al. 2020, Mon. Not. Roy. Astron. Soc., 498, 1420