© The Authors 2025

1. Introduction

Strong gravitational lensing is a consequence of general relativity, where light from a distant source is deflected and magnified due to the mass of an intervening object. The intervening massive object is the lens, and it acts like a gravitational telescope to observe the high-redshift universe. Over the last few decades, there have been many observations of strongly lensed galaxies and quasars, but a rapidly developing area of research is the observation of lensed transients. Type Ia supernovae (SNe Ia), the explosions of white dwarf stars in binary star systems, are bright, standardisable candles. Therefore, strongly lensed type Ia supernovae (glSNe Ia) are extremely valuable astrophysical tools.

In terms of supernova physics, glSNe Ia allow us to probe distant SNe Ia and study their redshift evolution (Petrushevska et al. 2017; Johansson et al. 2021). This is important for verifying that high-redshift SNe Ia are standardisable and can be included in the Hubble diagram. In terms of lensing physics, glSNe Ia can be used to determine the matter distribution in lensing systems. Furthermore, from the recent discovery of SN Zwicky (Goobar et al. 2023), the authors suggested that glSNe Ia could be vital for the study of previously unknown compact lens systems. In terms of cosmology, glSNe Ia can be utilised in the method of time-delay cosmography to measure the Hubble constant (H₀), which was first proposed by Refsdal (1964). Strong lensing produces multiple images of the supernova explosion, which travel different paths with different gravitational potentials to reach the observer, leading to time delays between the images. From measurements of the time delays and a model of the lens galaxy’s gravitational potential, an absolute distance (known as the ‘time-delay distance’) for the system can be determined. This distance provides constraints on H₀ and, to a lesser degree, the mass density parameter, Ω_M, and the dark energy equation of state, w(z) (Goobar et al. 2002).

The method of time-delay cosmography has previously been employed with quasar sources (e.g. Wong et al. 2020; Shajib et al. 2020; Millon et al. 2020). One of the benefits of using time-delay cosmography with glSNe Ia is that SNe Ia are standardisable candles, with a known peak absolute magnitude that provides absolute magnification constraints. This allows us to constrain the mass-sheet degeneracy (Falco et al. 1985); a systematic effect whereby adding a constant surface mass density to the lens system results in a degeneracy in which all of the lensing observables remain unchanged, except for the time delays. This degeneracy can be broken by knowledge of the unlensed apparent brightness (Birrer et al. 2022). Another benefit is that SNe Ia have predictable light curve shapes with a duration of ∼30 days, whereas the timescale required to monitor time delays from lensed quasars is several years due to their stochastic variability. Additionally, the host and lens galaxy can be studied without any contamination from the supernova once it has faded, which is not always possible for lensed quasars.

At the time of publication, only three glSNe Ia have been detected with ground-based telescopes: PS1-10afx with the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chornock et al. 2013; Quimby et al. 2013), iPTF16geu with the Intermediate Palomar Transient Factory (iPTF; Goobar et al. 2017), and SN Zwicky with the Zwicky Transient Facility (ZTF; Goobar et al. 2023; Pierel et al. 2023). All of these events had a high magnification (with a total magnification factor, μ, of greater than 20), allowing them to be followed up and classified by surveys that obtain spectra of nearby transients. This implies that a larger sample of less magnified glSNe Ia could exist, which may have already been observed by surveys such as ZTF, but remain undiscovered because they would not be targeted by spectroscopic follow-up programmes. Additionally, with upcoming large-scale optical surveys such as the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST; Ivezić et al. 2019), it is necessary for us to develop a systematic approach for detecting all glSNe Ia, instead of discovering only the extremely luminous objects. Most transients found by LSST will be too faint for spectroscopy, even with 8 m-class telescopes. Thus, with the lack of available spectroscopic follow-up resources, it would be valuable to devise a method to identify candidate glSNe Ia before a spectrum of the transient is obtained.

In this work, we present a systematic search of the ZTF archive for glSNe Ia. We developed a pipeline that is capable of filtering alerts and distinguishing candidate glSNe Ia from other variable sources. Initial cuts were applied to the ZTF alert photometry (e.g. constraints on peak absolute magnitude) before forced photometry was obtained for the remaining candidates. Additional cuts were applied to the remaining candidates from their light curve colours, lens galaxy colours, and parameters from light curve fits. We also cross-matched the candidates to the Dark Energy Spectroscopic Instrument (DESI) survey catalogue.

The structure of this article is as follows. Section 2 discusses the methods previously employed to search for lensed supernovae, how these methods were used as guidelines when implementing the search pipeline, and predictions for what we should expect to find in the ZTF archive. The data from ZTF used in this analysis is summarised in Sect. 3.1. Section 3.2 outlines the real-time alert processing platform AMPEL, which allowed us to analyse the ZTF data and structure our pipeline. In Sect. 4, we describe our method to identify probable glSNe Ia, starting with the initial alert pipeline described in Sect. 4.1. The candidates remaining after the initial selection method were analysed in parallel in Sect. 4.2, in which we cross-match to the DESI spectroscopic catalogue, and Sect. 4.3, in which we continue to filter the candidates based on more restrictive cuts. Section 4.4 summarises the sample of candidates that were selected by the two methods. In Sect. 5, we report on the properties of this sample by discussing the possible contaminants (Sect. 5.1) and by comparing to the Bright Transient Survey (Sect. 5.2) and to previously observed glSNe Ia (Sect. 5.3). Our seven most likely candidates (labelled the ‘gold sample’) are analysed in more detail in Section 6. The two objects we present as likely glSNe Ia, ZTF19abpjicm and ZTF22aahmovu, are examined in Sects. 6.4.1 and 6.4.2, respectively. Finally, in Sect. 7, we elucidate the findings of our study and what further work is required to confirm that our candidates are truly glSNe Ia.

2. Strategy for detecting lensed supernovae

The strategy for finding lensed supernovae depends on whether the instrument used is capable of resolving the individual images. Lensed supernovae discovered from spatially resolved data include SN Refsdal (type II, Kelly et al. 2015) and SN Requiem (likely a type Ia, but unconfirmed, Rodney et al. 2021), which were recovered from Hubble Space Telescope (HST) data that distinguished the separate images of each SN. More recently, the James Webb Space Telescope (JWST) discovered the triply imaged SN Ia H0pe (Frye et al. 2024) and the lensed SN Ia Encore (found in the same host galaxy as SN Requiem; Pierel et al. 2024).

To detect unresolved lensed supernovae, a variety of strategies can be employed. Goldstein & Nugent (2017) put forward the idea of searching for magnified sources (with a peak B-band absolute magnitude greater than the typical −19.4 mag for SNe Ia) close to elliptical galaxies. Their reasoning is that normal SNe Ia are the brightest objects in quiescent galaxies, so there is a high likelihood that the magnified object is a lensed transient from a background galaxy. Wojtak et al. (2019) calculated glSNe discovery rates as a function of survey depth and found that shallower, pre-LSST surveys such as ZTF would not be able to detect multiple images from glSNe Ia and the only viable method would be to find magnified sources.

Other proposed methods include detection from light curve information alone (Bag et al. 2021; Denissenya et al. 2022), monitoring known galaxy-galaxy strong-lens systems (Shu et al. 2018; Craig et al. 2024), or known cluster systems (Saini et al. 2000; Sullivan et al. 2000). However, these methods are not utilised in our following analysis. Firstly, we would struggle to characterise glSNe Ia from light curves alone due to occasional gaps in the light curve and limited filters (usually the g- and r-bands are available, whereas the i-band is limited in ZTF), so instead we must combine the parameters we extract from the light curve with additional information (e.g. from a catalogue-matched lens or host galaxy). Secondly, the three glSNe Ia discovered by ground-based telescopes had lens galaxies that were previously unknown, as they were compact (with an Einstein radius, θ_E, of less than 0.5 arcseconds) and could not be detected from ground-based galaxy surveys. Therefore, monitoring known galaxy-galaxy strong-lens systems might not be the optimal search method as it would bias the sample to massive lenses. Thirdly, monitoring cluster systems with shallow ground-based surveys has so far only been successful with the discovery of weakly lensed supernovae (e.g. Goobar et al. 2009; Patel et al. 2014; Nordin et al. 2014; Rodney et al. 2015; Petrushevska et al. 2016), which are not the focus of this work.

Quimby et al. (2014) were the first to suggest using light curve colours as a metric to distinguish lensed supernovae from unlensed. Specifically, they proposed that the observed r − i versus i-band distribution of unlensed SNe during the pre-peak phase differed significantly from the colours of PS1-10afx (the glSNe Ia detected by Pan-STARRS). We expect glSNe to be observed at higher redshifts (due to the probability of lensing and the number of SNe increasing with redshift) and, as a result, they will be redder than nearby unlensed SNe at similar apparent magnitudes (Quimby et al. 2014). This concept is also demonstrated in simulations by Sagués Carracedo et al. (2024), which we shall hereafter refer to as SC24. By generating realistic light curves from SNe Ia templates, the authors of SC24 show how the distributions of g − r, g − i, and r − i at different light curve epochs differ for lensed and unlensed SNe Ia.

In SC24, the authors examine the sample of glSNe Ia that would have been detected by ZTF specifically, by incorporating the observing logs for ZTF, so the results can be directly applied to our systematic search. The findings of SC24 formed a part of our pipeline; in particular, the light curve colours discussed above, but also the simulated distributions of host and lens redshifts. From the distributions of the host and lens galaxy redshifts in SC24, we do not expect to find any glSNe Ia in host galaxies in the range z < 0.1. Additionally, less than 3% of identifiable glSNe have lens galaxies in the range z < 0.1.

SC24 also performed SALT2 light curve fits to simulated glSNe Ia light curves. The SALT2 model is a light curve template for SNe Ia based on two observable parameters: the stretch, x₁, and the colour, c. This model is based on empirical observations, due to the current lack of understanding about SN Ia progenitor explosions (Guy et al. 2007). For typical SNe Ia (e.g., the ones that are used in cosmological analyses), we expect x₁ < |3| and c < |0.3|. As is shown in SC24, simulated glSNe Ia display a wider distribution of stretch values (approximately −3 ⪅ x₁ ⪅ 10) and are redder in colour (c ⪆ 0). They demonstrate that a cut of c > 0 removes approximately half of the unlensed SN Ia sample and most of the lensed core-collapse supernova contaminants, while also retaining 90% of the glSN Ia sample. SC24 also determined the peak B-band absolute magnitude distribution from the SALT2 fits, assuming the redshift of the lens galaxy is the supernova redshift (the motivation behind this assumption will be explained in Sect. 4.1). They show that we can exclude almost all unlensed SNe Ia with a cut of M_B > −20 mag, although this removes approximately 20% of the glSN Ia sample.

In addition to the simulations of SC24, our pipeline utilised a relationship between x₁ and peak brightness for glSNe Ia. There are two effects that contributed to the width and peak brightness of the observed light curves in our study. The first effect is due to the superposition of the unresolved images of the lensed transient. We expect that objects with shorter time delays will have narrower unresolved light curves with greater peak brightness because the peak values of the individual images will be observed at approximately the same observation epoch. Conversely, we expect that objects with longer time delays will have wider unresolved light curves with lower peak brightness¹. The second effect is due to lensing in the system causing an increase in flux, which is known as magnification. In a strong lensing system, the time delays and magnifications are negatively correlated. The combination of these two effects results in a trend in which unresolved glSN Ia light curves with a smaller x₁ value will have a larger peak absolute brightness. We implemented a cut for this relationship in our methodology to be sensitive to objects with short and long time delays.

The most challenging part of an archival study is estimating the true brightness of the candidates in the sample. Our analysis pipeline allowed us to process and filter the entire ZTF alert stream, which was necessary to consider faint, unclassified transients that have not been studied previously. Because we did not have spectra for these unclassified transients, we utilised the redshifts from catalogue-matched nearby galaxies (which, in the case of an actual lensed supernova, would likely be the lensing galaxy). Current spectroscopic surveys are limited to brighter, lower redshift galaxies, meaning that we rely on photometric redshift estimates for distant galaxies, which have larger uncertainties but are more complete.

In a previous work by Magee et al. (2023), the authors performed a similar archival search for lensed supernovae from the sample of transients that were observed by ZTF and reported to the Transient Name Server (TNS). However, this public sample was limited to brighter objects that would have been reported by the ZTF brokers, and the analysis of the transients did not utilise ZTF forced photometry to study objects at deeper magnitudes. The methods of Magee et al. (2023) included cross-matching with known lens systems, applying colour cuts to identify red objects as in Quimby et al. (2014), selecting light curves that are inconsistent with SN Ia near elliptical galaxies as is suggested in Goldstein & Nugent (2017), and identifying intrinsically luminous objects according to the host redshift. In this work, we focus on two of these methods: applying colour cuts and selecting for intrinsically luminous objects. Magee et al. (2023) identified 132 unique transients across their selection methods, but all of their candidates are consistent with non-lensed transients (and more than half are already spectroscopically classified). We find no overlap between their 132 candidates and our shortlisted candidates, which is likely due to their limited sample selected from public brokers.

Several studies have forecasted the rates of gravitationally lensed supernovae that ZTF would find per year. For example, in a study by Goldstein et al. (2019), the authors predicted a discovery rate of 1.23 glSNe Ia per year and at least 7.37 other types of supernovae per year. Therefore, across the 3.25 years in this archival search, we expect there to be approximately 28 glSNe in our sample, with four of those being SNe Ia. However, Goldstein et al. (2019) assumed that observations taken in the same filter in a single night would be stacked, while our analysis uses forced photometry data instead. Additionally, Wojtak et al. (2019) predicted that ZTF would observe approximately two glSNe Ia per year (that is, there should be six or seven total objects in our sample). More conservatively, Sainz de Murieta et al. (2023) forecasted that ZTF should have observed approximately 0.2 glSNe Ia per year, which is consistent with SN Zwicky being the only discovery. However, when incorporating the ZTF observing logs, SC24 predicted that ZTF would detect glSNe Ia with at least five 5-sigma detections around the peak at a rate of approximately four per year (Goobar et al. 2024). After requiring that the object is measured as over-luminous (assuming the object is at the lens redshift), the simulations of SC24 predict that approximately 1.4 glSNe Ia per year will remain. As a result of this study, we can expect to find up to four archival candidates in our analysis.

3. Analysing ZTF data using AMPEL

3.1. The ZTF archive

The Zwicky Transient Facility is a wide-field, optical survey for time-domain astronomy (Bellm et al. 2019b; Graham et al. 2019). The survey utilises a camera with a field of view of 47 degrees squared on the Palomar 48-inch (P48) Telescope that observes in three filters (g-, r- and i-bands) (Dekany et al. 2020). The observation plan consists of a public survey in the g- and r-bands that scans the northern sky every two days and a partnership survey with a higher cadence and additional i-band photometry over a smaller region of the sky (Bellm et al. 2019a). In this study, we had access to both public and partnership data. On average, ZTF produces between 600 000 and 1.2 million alerts on a complete night of observation (Masci et al. 2019; Patterson et al. 2019).

Owing to the rapid scanning and extensive sky coverage of the survey, ZTF is an excellent precursor to LSST and is capable of detecting rare transients such as lensed supernovae. The Bright Transient Survey (BTS, Fremling et al. 2020; Perley et al. 2020) is a spectroscopic survey within ZTF that aims to acquire spectra of all transients brighter than 18.5 magnitudes, which constitute the vast majority of classified objects in ZTF. Previously observed glSNe Ia may be unclassified in the ZTF archive if they were significantly less magnified than discoveries such as SN Zwicky. Therefore, a systematic search of the ZTF archive is a promising avenue to identify less magnified glSNe Ia.

In this analysis, we have demonstrated a method to systematically identify probable unresolved glSNe Ia, by reducing a total of 98 567 163 unique alerts in ZTF between June 1 2019 and September 1 2022 to a sample of seven glSN Ia candidates. This approach will also be utilised in the final months of ZTF for a live search for lensed transients.

3.2. The AMPEL platform: Alert processing and catalogue-matching

To systematically process all archival ZTF data, we require a platform that can ingest and filter large datasets. We also require the ability to cross-match our large sample of alerts to various photometric and spectroscopic redshift catalogues, so that we can estimate the redshift of each transient’s lens or host galaxy. To achieve these objectives, we used AMPEL², which is a publicly available real-time alert processing platform developed by Nordin et al. (2019). As well as being a broker for ZTF, it will also be one of the brokers for LSST, and is designed to be flexible enough to process a variety of astronomical data streams.

Within AMPEL, it is possible to filter alerts before they are ingested into a database, perform matches to numerous catalogues, and apply analysis methods (such as fitting to the SALT2 template) in a single pipeline. It is also possible to apply this pipeline to both alert and forced photometry from ZTF. We developed our glSNe Ia identification pipeline using the AMPEL infrastructure, specifically the astronomy analysis units from the AMPEL-HU-astro respository³. An example Python notebook that matches the structure of the pipeline is publicly available⁴.

Obtaining accurate lens (and host) redshifts for our candidates was the most crucial part of the study because the redshifts are necessary to estimate the absolute luminosity (and therefore the magnification). Several catalogues are available within AMPEL and the complete list is publicly available via the API⁵. The specific catalogues that we queried are as follows: Sloan Digital Sky Survey (SDSS) DR10 (Brescia et al. 2015); NASA/IPAC Extragalactic Database (NED⁶), accessed through the catsHTM tool (Soumagnac & Ofek 2018); Galaxy List for the Advanced Detector Era (GLADE) v2.3 (Dálya et al. 2018); WISExSCOS Photometric Redshift Catalogue (WISExSCOSPZ) (Bilicki et al. 2016); 2MASS Photometric Redshift catalogue (2MPZ) (Bilicki et al. 2014); Legacy Survey (LS) DR8 (Duncan 2022); Pan-STARRS1 Source Types and Redshifts with Machine Learning (PS1-STRM) (Beck et al. 2021).

From the LS catalogue, we obtained the galaxy colours in the g-, z-, and W₁-bands as well as the photometric redshift for each object. In addition to the catalogues offered by AMPEL, we cross-matched our more promising transients with the DESI spectroscopic survey. This is currently the largest spectroscopic survey operating, with the ability to observe high-redshift galaxies. Our catalogue matching to this survey is discussed in more detail in Sect. 4.2.

4. Candidate glSNe Ia from the ZTF archive

4.1. Alert pipeline

Here, we describe the functional steps of all modules of the workflow (as was mentioned above, this pipeline can be reproduced using the publicly available AMPEL interface). A flowchart to summarise the entire pipeline is shown in Fig. 1. We note that the term ‘alert’ refers to a transient event, whereas ‘object’ refers to a persistent source that is associated with the alert.

Fig. 1. Flowchart to illustrate the stages of the method, as is described in Sect. 4. The red numbers to the bottom right of the boxes indicate the number of transients that remained after that stage of the analysis. The asterisk represents a stage that has not been completed during our study and is left for future work.

Alert photometry from June 1 2019 to September 1 2022 was ingested based on a query with basic initial criteria. These criteria were: a point spread function (PSF) apparent magnitude of greater than 18; greater than six previous detections but fewer than 30; and a positive flux value.

We applied a filter that allowed for SN-like alerts. The filtering criteria included a real bogus (RB) value of greater than 0.3, loose cuts on image quality (such as the full width at half maximum, elongation, and number of bad pixels), distance to known solar system objects, and an alert time history of less than 60 days. Furthermore, probable stars were rejected using a match to the Pan-STARRS PS1 and Gaia DR2 catalogues. The details of this filter are given in the DecentFilter unit of AMPEL⁷. The number of alerts remaining after this stage was 31930.

All remaining alerts were cross-matched to the catalogues listed in Sect. 3.2 to establish whether they are associated with a galaxy. From the catalogues, we determined properties of the associated galaxy, such as the redshift and the angular separation to the alert location. If there was more than one catalogue match for a particular alert, the catalogue with the better ranking (based on accuracy⁸) was chosen. If the catalogues are ranked similarly, an average of the redshifts and distances to the galaxy were used. We note that the catalogue-matched redshift might belong to either the host or lensing galaxy. However, it will most likely belong to the lensing galaxy since this is closer and more luminous (it is unlikely the host galaxy will be sufficiently magnified). This was assumed in the following analysis.

Using the Python library SNCosmo (Barbary et al. 2022)⁹, we performed a preliminary SALT2 fit to the alert light curves (using version 2.4, trained by Betoule et al. 2014). From this fit, we obtained a lower limit for the peak B-band absolute magnitude. This value is a lower limit because we used the catalogue-matched redshift, which we assume belongs to the potential lens.

We applied ‘base cuts’ to find transients that are most likely to be lensed. Initially, all alerts already classified as something other than a glSN by the ZTF Bright Transient Survey (BTS) were removed. We applied the following cuts on the remaining candidates:

1. Peak B-band absolute magnitude brighter than −19.5 mag. This cut selects transients brighter than normal SNe Ia. The simulations of SC24 show that most glSNe Ia are in this magnitude region (regardless of whether the catalogue-matched redshift belongs to the host or the lens galaxy).

2. A redshift greater than 0.1 (also based on simulations of lens and galaxy distributions from SC24). Additionally, the probability of lensing increases at higher redshifts, due to the larger number of possible lensing galaxies, and thus a higher chance of strong lensing alignment.

3. Less than 3″ radius to a catalogue-matched galaxy, which ensures that the candidate is close to the lens galaxy core, where strong lensing is likely to occur. This is also supported by the simulated catalogues of Goldstein et al. (2019), which we utilise in Fig. 2 to show the distribution of the median angular image separation. From this simulation study, we calculate that 98% of the glSNe Ia have a median angular image separation of less than 3 arcseconds. For a lensing system in perfect alignment, the image separation is double the separation of the image from the lens. Therefore, we can assume that a 3″ radius is a reasonable upper limit for the galaxy separation.

   Fig. 2. Histogram illustrating the distribution of median angular image separation of simulated glSNe Ia, using data from the catalogues of Goldstein et al. (2019). The dashed black line indicates a 3″ cut that selects 98.4% of the distribution.

4. At least five detections around the peak (within a range of 10 days before and 20 days after the peak). This is required to distinguish the nature of the transient from the light curve.

The number of glSN Ia candidates remaining at this stage was 7075. Next, we obtained forced photometry (FP) to recover data points at the lower detection limit of ZTF. This was done using PSF fitting (Masci et al. 2023), with the forced photometry pipeline fpbot developed by Reusch (2023)¹⁰. Forced photometry allows us to differentiate between supernovae and long-lasting transients such as active galactic nuclei (AGN) and follow the tail of said supernovae. We removed the baseline of the forced photometry data using the AMPEL unit ZTFFPbotForcedPhotometryAlertSupplier¹¹. Objects where it was not possible to remove the baseline, for example, due to long-lasting stochastic variability, were removed from the sample. This accounted for approximately 13% of the sample at this stage. The data were then truncated using the unit T2PhaseLimit¹², which selects the portion of the light curve with a significant supernova-like peak. The forced photometry data was utilised in the following two analysis methods of Sects. 4.2 and 4.3 to further narrow down our candidates.

We note that the following analysis is based on point source forced photometry, whereas our aim is to detect multiple point sources (e.g. two or four lensed images) that appear unresolved. For glSNe Ia with larger angular separations between the images, this can result in biased light curves. Additionally, it is likely that glSNe Ia with larger image separations (and therefore, displaying blended or multiple images) would be rejected by the ZTF real bogus classifier at the stage where we filter the alerts. Therefore, a different search pipeline (likely image-based) would be required to detect glSNe Ia with larger image separations, which is not the focus of this work.

Executing the archive search pipeline produced 7075 candidate events. This set was analysed further based on whether a DESI spectroscopic redshift was available (Sect. 4.2) or not (Sect. 4.3).

4.2. Selection method 1: Spectroscopic sample with DESI cross-matching (S1)

The Dark Energy Spectroscopic Instrument (DESI; Levi et al. 2013; DESI Collaboration 2016a,b, 2022, 2024b,a; Silber et al. 2023; Miller et al. 2024; Guy et al. 2023; Hahn et al. 2023; Bailey et al. in prep.) is a spectroscopic instrument designed to measure the impact of dark energy on the expansion of the Universe. The 5-year DESI survey is currently ongoing and is collecting redshifts for tens of millions of galaxies and quasars (QSOs). The DESI Y1 Key Project publications are DESI Collaboration (2024c,d,e). As part of our study, collaborators working with the DESI survey contributed redshifts (including data not yet public) for galaxies and quasars that were less than 10″ away from our archival transients, and with redshift values greater than 0.15.

We submitted the 7075 candidates that passed the ‘base cuts’ in the pipeline described in Sect. 4.1, for which DESI had lens or host redshifts for 1269 of the candidates. Of these, 519 were labelled as galaxies (and not QSOs) and within 3″ of the ZTF archival transient. We used the DESI redshifts where available to redo the SALT2 fits, as is described in Sect. 4.1. A total of 461 objects had a converging SALT2 fit and were utilised in the next stage of the analysis. The objects that did not pass this stage were likely not supernova-like, with long-duration light curves and stochastic variability. This long-term variability may not have been present during the alert stage, but is visible in the forced photometry data. After reapplying the base peak B-band absolute magnitude cut of M_B < −19.5 mag, 257 transients remained out of the 461. From this sample, we can assert that approximately 56% of the photometrically over-luminous candidates in our sample (with a converging SALT2 fit) were also spectroscopically over-luminous (in this case, over-luminous means peak M_B < −19.5 mag). This is a useful statistic for when we consider the sample with only photometric redshifts in Sect. 4.3.

Next, a stricter peak B-band absolute magnitude cut of M_B < −20 mag was applied. This was done primarily to select objects that are significantly magnified, and therefore are more likely to be glSNe Ia (as opposed to belonging to over-luminous SNe Ia subclasses such as type Ia-91T, Ia-CSM, or Ia-03fg). From the simulations of SC24, we can see that this cut removes almost all unlensed SNe Ia and approximately 20% of real glSNe Ia events. We applied this cut after obtaining forced photometry because the additional data points close to the magnitude limit should allow the SALT2 fit to better characterise the light curve shape (including the peak). In other words, we were more confident in our measurement of M_B with the forced photometry data and could apply a stricter cut.

The remaining 117 candidate light curves were inspected by at least two scanners. Initially, objects were excluded if they did not look supernova-like (i.e. with long-term AGN-like variability). From this inspection, 50 candidates were selected as possible glSNe Ia. At this point, we requested the DESI galaxy spectra for these candidates, to confirm that the redshifts were consistent with the spectral emission lines. The parameters for this sample of 50 are included in Table A.2 (available at the CDS).

From our spectroscopic sample of 50, we narrowed down our candidates by visual inspection once again. This time, objects were excluded if they did not look Ia-like (e.g. if they looked similar to superluminous supernovae, tidal disruption events, or AGN; refer to Sect. 5.1) or if they were previously classified on TNS. Additionally, redshifts were updated as DESI data reprocessing and re-observations occurred throughout this programme, which meant that seven objects were removed from the sample. Consequently, we had 27 final candidates from the spectroscopic sample, which are discussed in Sect. 4.4 and shown in Table 3.

4.3. Selection method 2: Photometric sample from selection with additional cuts (S2)

The majority of the objects that passed the base cuts described in Sect. 4.1 did not have a spectroscopic redshift from DESI. It was not possible to manually inspect these remaining objects, so we applied further cuts motivated by the simulations of SC24 and observational data. For example, Fig. 3, illustrates the LS g − z versus g − W₁ lens galaxy colours for the sample of candidates remaining after the base cuts within the redshift bin 0.2 < z < 0.3. By comparing to catalogues of known lens galaxies¹³ and QSOs (from the WISE catalogue of AGN; Assef et al. 2018), it is evident that a large fraction of the candidates’ associated galaxies do not have the expected colours of a lensing galaxy.

Fig. 3. Plot of the LS g − z versus g − W1 galaxy colours for the candidates remaining after the base cuts (orange), compared to contours from confirmed and candidate lens galaxies (blue), and QSOs detected by WISE (green). The data was acquired by cross-matching the AGN and lens catalogues with LS (Duncan 2022). The dashed lines indicate the colour cuts which were devised to select likely elliptical and lens galaxies, shown in Table 2. A single redshift bin of 0.2 < z < 0.3 is shown for clarity. The two star markers show the lens galaxy colours for iPTF16geu (grey) and SN Zwicky (black), which had lens redshifts of 0.216 and 0.226, respectively. For SN Zwicky, the g − z versus g − W1 colours from Pan-STARRS are shown instead, because the lens galaxy was not present in the LS catalogue).

By implementing stricter cuts, we aimed to increase the likelihood that a transient was truly a lensed supernova, even though photometric redshifts are less reliable. We required that the transient more closely resembled the expected features of a glSN Ia. These characteristics include the light curve colour evolution, the lens galaxy colours, and the SALT2 fit parameters. We note that the additional cuts were applied to all remaining candidates at the end of Sect. 4.1, not just the objects without a DESI galaxy match in Sect. 4.2.

First, we refitted the light curves with the SALT2 template to recalculate the peak absolute magnitude values from the new forced photometry data. If a SALT2 fit could not converge, or the new peak absolute magnitude was dimmer than −20 mag, the candidate was removed from the sample. The justification for this stricter absolute magnitude cut is the same as in Sect. 4.2.

The remaining candidates were selected from a set of additional cuts:

1. A selection of g − r, g − i, and r − i colour cuts at light curve phases of t₀ − 7, t₀, and t₀ + 7 days (where t₀ is the time at peak B magnitude) were determined based on simulations from SC24. The values for these cuts are summarised in Table 1. The motivation for these cuts is that they differentiate normal SN Ia from lensed ones, given that the latter are usually significantly redder (when observed at the same peak apparent magnitude). In SC24, the simulations show that this is the colour space where we expect to find 90% lensed and 10% unlensed SNe Ia. Candidates were labelled as passing this cut if they passed at least one of the colour cuts at a certain epoch (i.e. they only needed to pass one out of nine of the cuts). The leniency of this cut is due to some light curves missing bands – usually the i-band – and/or certain epochs. Additionally, we used light curve information from a window of ±3 days around the phases of t₀ − 7, t₀, and t₀ + 7 days, so that we did not penalise objects that were observed with a cadence of two or three days.

   Table 1.

   Light curve colour cuts for the g-, r-, and i-bands applied at a phase of t₀ − 7, t₀, and t₀ + 7 days.

2. Using galaxy colours from LS, we implemented colour cuts that aimed to exclude AGN and blue star-forming galaxies whilst retaining lens-like galaxies. For this cut, we assumed that the catalogue-matched galaxy is the lens or a blend of the lens and the SN host galaxy (the blended systems, similar to iPTF16geu, will be more in the upper right corner, and so will also be preserved with the cuts). The values are chosen from the empirical observations in Fig. 3 and are given in Table 2. We note that the population of lenses that host glSNe Ia in surveys such as ZTF may be different that than the population of lenses found in galaxy surveys (as is discussed in Goobar et al. 2023). However, this is not possible to quantify with such a small sample of known glSNe Ia. Additionally, a minimum galaxy brightness of greater than 21.07 in the LS z-band (corrected for Milky Way extinction) was required. This was based on simulations of lens galaxies from Wojtak et al. (2019) and Arendse et al. (2024) with K-corrections from Lenspop (Collett 2015)¹⁴. The simulations show that we expect 90% of lens galaxies to be brighter than this magnitude limit.

3. Based on the simulated distributions shown in SC24, we predict that 90% of glSNe Ia light curves have a SALT2 colour parameter, c, of greater than 0 (which is synonymous with being redder). Thus, a cut of c > 0 was applied to the sample.

4. As was discussed in Sect. 2, we expect that glSNe Ia with narrower light curves will have a greater total magnification (and conversely, wider light curves will have a smaller total magnification). This means that there are likely two populations of glSNe Ia: less magnified objects with larger x₁ values that may have poor SALT2 fits, and more magnified objects with smaller x₁ values that more closely resemble normal SNe Ia. The transients observed with greater magnification will have a better line-of-sight alignment with the core of the lensing galaxy, meaning that the distance to the catalogue-matched galaxy should be less (the actual value depends on the Einstein radius of the lensing system, but for typical compact systems we would expect a separation of θ < 1″). This information suggests that we could study the two populations separately so that we could remove contamination by normal SNe Ia more easily. To select the highly magnified glSNe Ia with smaller time delays, we applied a SALT2 stretch cut of −3 < x₁ < 3, which would select the objects that look like normal SNe Ia. Additionally, we required that the object was less than 1″ away from the catalogue counterpart, to increase the likelihood that the object was a highly magnified glSN Ia. Conversely, to select the less magnified glSNe Ia, we applied a SALT2 stretch cut of 3 < x₁ < 20, which would exclude the normal SNe Ia. This cut is motivated by the distribution of x₁ parameters for simulated lensed and unlensed SNe Ia in SC24. We note that a small angular separation to the lensing galaxy was used as a proxy for magnification, instead of the peak absolute magnitude, to minimise the impact of inaccurate photometric redshifts.

The candidates were inspected by at least two scanners. Objects were excluded if they did not look Ia-like (i.e. if they looked similar to other transients, such as the contaminants discussed in Sect. 5.1), if they were classified as something else, or if the data was too limited or noisy to determine what they were.

After this manual inspection stage, 27 photometric candidates were remaining. Seven of the candidates also had a DESI cross-match and were identified in the spectroscopic sample from Sect. 4.2. The 27 photometric candidates are discussed in detail in Sect. 4.4 and shown in Table 4.

A discussion of the accuracy of the various photometric catalogues utilised in this study is beyond the scope of this paper. However, from our analysis in Sect. 4.2, we can assume that at least half of the photometrically over-luminous candidates identified in this selection method would also be spectroscopically over-luminous if a spectrum for the lens or host was obtained.

The methodology described in Sect. 4 is summarised as a flowchart in Fig. 1, which includes the number of transients that remained after each stage in the analysis (indicated in red).

4.4. Summary of candidates

Table 3 summarises the 27 candidates from S1 described in Sect. 4.2, with spectroscopic redshifts and angular separations from a catalogue match with DESI. This selection of spectroscopically over-luminous candidates is interesting since we expect that the majority of the events are bright enough to be glSNe Ia. However, if they do not pass the cuts specified in Sect. 4.3, it is more likely that they are superluminous supernova candidates. Despite this, we cannot exclude the possibility that some of them are glSNe Ia (or other lensed transients) and the reason(s) that they did not pass the cuts are:

1. The fit to the SALT2 template fails to model the light curve well enough, which may be the case if the lensed object is not a typical SN Ia or it is a SN II. This could mean that it does not pass the x₁, c, and light curve colour cuts, or that a SALT2 fit does not converge at all.

2. The catalogue-matched galaxy information is not present in the LS catalogue, the LS catalogue galaxy is dimmer, or the colours are different than what we would expect from a lens galaxy. This would mean that it does not pass the galaxy cuts.

To eliminate these possibilities, careful inspection of each of the candidates is necessary. To increase the likelihood that a candidate is a glSN Ia, we could fit the light curves to a combined SALT2 model with multiple images, and check that the resulting fit parameters are reasonable compared to the results from simulation studies.

Table 4 summarises the 27 candidates from S2 described in Sect. 4.3, with the photometric redshifts and angular separations from catalogue matches within AMPEL. The uncertainty on the peak absolute magnitude is dominated by the uncertainty on the photometric redshifts. As a result of the large uncertainties, the majority of the candidates from S2 are consistent with a normal SN Ia within 3σ, so it is hard to assign a confidence level to the claim that a candidate is truly a lensed supernova. Therefore, the photometric sample in Table 4 will need to be targeted in future spectroscopic follow-up missions to get the lens and host redshifts before we can make any conclusions about them being glSNe Ia. At the time of publication, this was not yet possible, but it is something that we aim to do in future work.

5. Sample properties

5.1. Sources of contamination

Five main sources of contamination may have passed the cuts in our pipeline: AGN, tidal disruption events (TDEs), superluminous supernovae (SLSNe), lensed core-collapse SNe, and unlensed SNe Ia (either due to their peculiar over-luminosity or an incorrect photometric redshift).

Firstly, AGN are luminous enough to look similar to lensed events and constitute a large fraction of the 7075 transients that pass the base cuts with alert photometry. However, they were likely removed by the cuts applied in Sects. 4.2 and 4.3 after we obtained forced photometry (because the FP should contain any long-term variations that were slightly below the detection threshold for the alert photometry). Even if there were AGN that passed these additional cuts, the long-term variability in the light curves would be clear enough to remove during the visual inspection stages. An example of this in our study is ZTF20abgdrea, shown in Fig. 4. This object was present in both selection methods, but was filtered out during the visual inspection stage. The dotted black lines indicate the section of the light curve that was selected by the pipeline for the SALT2 fit, which explains why it passed our stringent cuts from selection method 2 (S2). Furthermore, ZTF20abgdrea is 0.7″ away from a candidate lens system in the SuGOHI VI catalogue (Sonnenfeld et al. 2020), which means that it could also be a lensed AGN or QSO. However, the DESI spectrum shows no evidence of AGN activity, so this is not possible to confirm without further study of the lens system.

Fig. 4. Light curve of ZTF20abgdrea, a candidate that passed S1 and S2, but was removed during visual inspection due to the long-term variability that suggests it is a likely AGN or QSO. The dotted black lines indicate the section of the light curve that was selected by the pipeline for the SALT2 fit. The baseline was not removed from this forced photometry, to more clearly illustrate the variability.

Secondly, TDEs reside close to galaxy cores and can also be very luminous, so they satisfy most of our criteria. However, they are quite rare events (although not as rare as glSNe Ia) and their light curves have distinct features, such as an approximate power-law decline, and are sometimes accompanied by a dust echo (for example, detected in the WISE infrared bands) (Reusch et al. 2022). This distinctive light curve shape allowed us to identify and remove likely TDEs from our sample. An example of this is ZTF20acyxxfo, shown in Fig. 5. This object (with a DESI redshift of z = 0.2 which implies a peak M_B < −20  mag) was present in both selection methods but was filtered out during the visual inspection stage, as it has a good fit to a TDE template light curve (Reusch et al., in prep.).

Fig. 5. Light curve of ZTF20acyxxfo, a candidate that passed S1 and S2, but was removed during visual inspection because it was identified as a likely TDE due to the approximate power-law decline.

Thirdly, SLSNe are a large contaminant because they are very luminous (with a peak absolute magnitude of approximately −21) and some of their light curves closely resemble the light curves of glSNe Ia with longer time delays. However, we can exclude some SLSNe with longer durations that would not be feasible for glSNe Ia (i.e. a rest-frame duration of greater than 100 days). Additionally, the galaxy brightness cut that we describe in Sect. 4.3 should exclude the hosts of hydrogen-poor superluminous supernovae (SLSNe-I), which are typically bluer and fainter (Lunnan et al. 2014; Leloudas et al. 2015; Schulze et al. 2018). To fully exclude SLSNe as contaminants, a careful examination of each candidate and all the data available is necessary; we perform this examination in Sect. 6 for our most likely candidates. An example of a possible SLSN that we removed during the visual inspection stage of S2, ZTF20acudrme, is shown in Fig. 6. This object had a photometric redshift of z = 0.25, which allowed us to estimate a peak M_B < −20 mag. However, the long duration and large separation from its galaxy counterpart (at 2″) suggest that ZTF20acudrme is a SLSN candidate.

Fig. 6. Light curve of ZTF20acudrme, a candidate that passed S2, but was removed during visual inspection because the long duration indicates it is a likely SLSN.

Additionally, some lensed core-collapse supernovae may be present in this sample. Typically, core-collapse SNe have a peak B absolute magnitude of approximately −17 (Perley et al. 2020), which means that they would have to be extremely magnified to pass our final peak cut of M < −20 mag. Furthermore, the simulations of SC24 showed that applying a c > 0 cut should remove almost all lensed core-collapse SNe. This significantly reduces the probability of their detection by our algorithm, although their presence is not ruled out. If we did find something magnified but looked more similar to a core-collapse SN than a SN Ia, we would not exclude it from our analysis (as this would also be an interesting discovery). However, none of the candidates in our sample have the obvious characteristics of a lensed core-collapse supernova.

Transients that interact with the circumstellar material (CSM) are possible contaminants. Supernova classes like type Ia-CSM and type IIn can be more luminous than normal SNe Ia due to the supernova ejecta interacting with the material that surrounds the progenitor, causing an increase in the brightness we observe (Sharma et al. 2023). However, typically the light curve will display signatures of CSM interaction, such as multiple bumps or plateauing phases. Fig. 7 illustrates an example candidate that was removed during the visual inspection of S1, due to the multiple peaks of varying widths, which could not be caused by a multiply imaged SN.

Fig. 7. Light curve of ZTF20acqhhyr, a candidate that passed S1, but was removed during visual inspection because of the presence of multiple bumps in the light curve of not-equal width, which suggests that CSM interaction is causing the high luminosity.

Other peculiar subclasses of SNe Ia can be over-luminous, such as 03fg-like (or super-Chandrasekhar) and 91T-like SNe Ia. SNe Ia-03fg are believed to be the result of white dwarf progenitors that exceed the Chandrasekhar limit before they explode (Howell et al. 2006). These objects can surpass a peak M_B of −20 mag; however, they are typically found in low-mass galaxies with a lower luminosity (for example, Howell et al. 2006; Childress et al. 2011; Hsiao et al. 2020; Lu et al. 2021). As a result, it is unlikely that they would pass our galaxy cuts in Sect. 4.3. Examples of 03fg-like SNe found in more massive galaxies (e.g. Taubenberger et al. 2011; Chakradhari et al. 2014) are further from the centre of the host and so would not pass our requirement for less than 3″ from a catalogue-matched galaxy. Additionally, 03fg-like SNe are typically less red than what we expect from glSNe Ia, so only a small fraction would pass our light curve colour cuts. Similarly, SNe Ia-91T are believed to arise from non-typical progenitor systems that cause the SN explosion to be more luminous. However, SNe Ia-91T are typically 0.2 mag brighter than normal SNe Ia (Yang et al. 2022), so it is unlikely that they would pass our peak magnitude cut.

Finally, we note that some objects with erroneous photometric redshifts will be present in the photometric sample in Table 4. These are likely to be normal SNe Ia at redshifts of z ∼ 0.1.

5.2. Comparison with the Bright Transient Survey (BTS)

BTS is a spectroscopic supernova survey within ZTF that aims to acquire spectra of all transients brighter than 18.5 magnitudes (excluding galactic sources, AGN, or moving objects), using the Spectral Energy Distribution Machine (SEDM) that operates on the Palomar 60-inch telescope (Blagorodnova et al. 2018). The survey aims to create the largest unbiased, brightness-limited sample of supernovae. Additionally, they target objects brighter than 19 magnitudes when the spectroscopic resources are available. The first spectrum of SN Zwicky was obtained as a result of BTS scanners submitting it to the SEDM for spectroscopic classification (Goobar et al. 2023).

To quantify the amount of contamination from unlensed supernova-like sources in this study, a comparison to the BTS was performed. We obtained forced photometry for 4127 transients in the BTS catalogue¹⁵ with a redshift (including the classes SLSN, SN Ia, SN Ib/c, TDE, and their subclasses), and processed the data as is described in Sect. 4.1. A SALT2 template fit was performed using the redshift from the BTS catalogue. Figure 8 shows a histogram with the derived peak M_B from the SALT2 fits for all the objects in the sample.

Fig. 8. Histogram illustrating the peak B absolute magnitude distribution for the 5 main classes of transients (SN Ia, SN II, SN Ib/c, SLSNe, and TDEs) recorded by the BTS. The dashed black line illustrates the objects with MB < −20 mag.

First, we applied the base cuts as is described in Sect. 4.1 (M_B < −19.5 mag, z > 0.1, less than 3″ distance from a catalogue-matched galaxy, at least 5 detections around the peak). In the sample remaining after base cuts, we retained 43% of the original SLSNe, 29% of the TDEs, 2% of the SNe II (of which the majority were SN IIn), 2% of the SNe Ia, and less than 1% of the SNe Ib/c. After applying the stricter absolute magnitude and redshift cuts (M_B < 20 mag and z > 0.15) described in Sects. 4.2 and 4.3, the only contamination in the sample was due to SLSNe (24% remaining), TDEs (14% remaining), and a small contribution from SNe II (less than 1% remaining). After applying all the cuts described in Sect. 4.3, only one object remained (ZTF21aahfjrr, classified as a SN IIn, but likely a SLSN-IIn). This study illustrates that we should expect non-negligible contamination from SLSNe and TDEs from the selection method of Sect. 4.2, but very little contamination from the method of Sect. 4.3. The most effective cuts from the method of Sect. 4.3 for removing objects were the galaxy colour cuts, as they removed over half of the remaining sample.

However, it is important to note that the brightness range of transients in the BTS does not represent the range expected in this study (as we are considering objects up to the detection limit of ZTF). As a result, we expect more contamination from intrinsically bright objects at larger redshifts.

5.3. Comparison with known glSNe Ia

It is valuable to compare the statistics of our candidates with the glSNe Ia already observed by ZTF and iPTF. Figure 9 shows the distribution of the peak B-band absolute magnitude versus the redshift for the candidates in the sample, assuming that the object is at the lens redshift. This peak absolute magnitude would be a lower limit on the true peak if the object were truly lensed. The blue contour displays the transients that passed the base cuts described in Sect. 4.1 and the data points illustrate the candidates passing the two selection methods (S1 and S2) described in Sects. 4.2 and 4.3, respectively. The corresponding values for SN Zwicky and iPTF16geu are also displayed. From Fig. 9, it is apparent that SN Zwicky and iPTF16geu match the brightest magnification tail. This implies that, if our candidates are truly glSNe Ia, they would belong to a less magnified population. Excluding a single high redshift candidate from the S1 sample (ZTF20acirhoc, z = 1.4), the majority of the candidates are found to have lens galaxies within the range 0.2 < z < 0.4.

Fig. 9. Distribution of the peak B-band absolute magnitude at the lens redshift versus the lens redshift for the transients that passed the base cuts described in Sect. 4.1 (blue contours), the candidates passing the two selection methods described in Sects. 4.2 and 4.3 (green and orange crosses, respectively), and for the already discovered glSNe Ia SN Zwicky and iPTF16geu (black and grey stars, respectively). The upper limit of the x axis has been adjusted for clarity, omitting one higher redshift object from S1 (ZTF20acirhoc).

Figure 10 shows the distributions for peak B-band apparent magnitude (m_B), peak B-band absolute magnitude (M_B), lens redshift (z), SALT2 x₁, SALT2 c, and angular separation values for the transients that passed each stage of the cuts (the base cuts, S1 and S2). For comparison, the dotted line with a star marker indicates the corresponding values for SN Zwicky for each distribution parameter, using the values given in Goobar et al. (2023) (with the caveat that the SN Zwicky parameters were fitted using a SALT2 time-delay model with four separate images, while our pipeline is only fitting for one resulting SALT2 light curve).

Fig. 10. Set of histograms displaying the distributions for peak mB, peak MB, lens z, SALT2 x1, SALT2 c, and angular separation values for the transients that passed the base cuts described in Sect. 4.1 (orange bars), and the candidates passing the two selection methods described in Sects. 4.2 and 4.3 (solid black and red lines, respectively). The dotted line with a star marker indicates the corresponding values for SN Zwicky for each distribution parameter.

Additionally, we applied our pipeline to the alert photometry for SN Zwicky and found that it passed all of the base cuts, as well as all of the additional cuts described in Sect. 4.3, except for the galaxy colour cuts (as they rely on Legacy Survey data, which is not available for this area of sky). However, as is shown in Fig. 3, a catalogue counterpart for SN Zwicky observed by Pan-STARRS would have passed our galaxy colour cuts.

6. Gold sample of glSN Ia candidates

There is a subsample of candidates that pass all the cuts described in Sect. 4.3 and also have a cross-match with a galaxy counterpart observed by DESI. Table 5 summarises the candidates that were found by both selection methods. Hereafter we shall refer to this sample of candidates as the ‘gold sample’. The light curves for the seven candidates in the gold sample are displayed in Fig. 11, which includes their DESI galaxy redshift and peak absolute magnitude (assuming the object is located at the catalogue redshift). Additionally, the Legacy Survey (LS) field cutouts for the seven candidates are shown in Fig. 12, where the location of the ZTF transient is indicated by a green circle. This sample represents our most likely glSN Ia candidates because we are confident that they are more luminous than normal SNe Ia and they have the expected light curve and lens galaxy characteristics of glSNe Ia.

Fig. 11. Light curve plots for the seven candidates in the gold sample that had a spectroscopic redshift from the DESI catalogue and passed all the cuts described in Sect. 4.3. Each plot shows the light curve, the lens redshift, and peak absolute magnitude at the lens redshift (from the DESI cross-match).

Fig. 12. Image cutouts of the LS field at the location of the transient (indicated by a green circle) for the seven gold sample candidates (the ZTF name for each transient is shown in the top left of the cutout). The galaxies we see are likely the lensing galaxies, and the host is more distant and faint.

In the following subsections, we examined the seven gold sample candidates further by fitting to a simple two-image SN Ia model (Sect. 6.1), calculating the rise, decline, and duration of the light curves (Sect. 6.2), and by obtaining the galaxy photometry (Sect. 6.3). Each object is also discussed in depth, starting with the most likely candidates (Sect. 6.4).

6.1. Two-image combined SALT2 fit

We fitted the light curves of the gold sample to a combined SALT2 model with two images. While simulations from Goldstein et al. (2019) predict that most events detected by ZTF would have four images (at a percentage of 62%, versus 22% for two images), constraining the parameters for four images simultaneously is challenging without additional information about the system (e.g. the flux ratios of the images to constrain each magnification) and also not possible with the limited number of detections we have for some candidates. By including a second component to the SALT2 light curve, we can show that the transient differs from a typical unlensed SN, which strengthens the argument that the object is multiply imaged. The combined model, in flux space, is given as

$$\begin{array}{c}\hfill f_{\mathrm{model}}(t)=F_1\times f_{\mathrm{SALT}2}(t)+F_2\times f_{\mathrm{SALT}2}(t-\mathrm{\Delta }t),\end{array}$$(1)

where F₁, F₂ are the relative image magnifications, Δt is the second image time delay (in the observer-frame), t is the phase of the light curve, and f_SALT2 is the flux model for a single SALT2 template. The SALT2 model also depends on the following parameters: the time of the light curve peak (t₀), the amplitude (x₀), colour (c), and stretch (x₁) of the light curve, and the redshift of the transient (z). We note that these image magnifications are relative, and not absolute values. If we assume that the lenses can be modelled as a singular isothermal sphere or halo (SIS), then the absolute total image magnification (which is the sum of the individual image magnifications) for the two-image system is given by

$$\begin{array}{c}\hfill {\mu }_{\mathrm{tot}}=\frac{2(r+1)}{r-1},\end{array}$$(2)

where we define the flux ratio as r = F₁/F₂¹⁶.

We allowed the redshift parameter to vary between a lower limit of the DESI z, and two times this value. This is because we expect the DESI z to belong to the lens, and the actual supernova to belong to a more distant background galaxy. To avoid degeneracy with the image magnifications, the intrinsic amplitude of the SNe Ia light curve (x₀) is fixed to an arbitrary value. We implemented bounds of [ − 2, 2] on x₁ and [ − 0.5, 0.5] on c, to match the parameters of a normal SN Ia. We limited the time delay for the second image to values between zero and 60 days (this value is motivated by the simulations of SC24 and Goldstein et al. 2019, and also implied by the absence of a distinct secondary peak).

Of the seven candidates in the gold sample, only three of them converged to reasonable two-image combined SALT2 fits: ZTF19abpjicm, ZTF20abjyrxf, and ZTF22aahmovu. Figure 13 illustrates the best fit two-image model for each candidate and Table 6 displays the fit parameters (as well as the derived parameters μ_tot, μ₁, and μ₂) for each model. The best fit time delays are larger than what we observed for SN Zwicky (where the time delays for all the images were less than a day) but are not unexpected, according to the simulations of SC24. They find that the median maximum time delay for glSNe Ia with a peak apparent magnitude of m > 19 is 8.9 days in the rest-frame (which would correspond to 12.5 days in the observer-frame, assuming a glSN at z = 0.4). Additionally, they predict that approximately 20% of glSNe Ia in this magnitude regime would have maximum time delays of greater than 25 days in the rest-frame.

Fig. 13. Two-image combined SALT2 template fits for the three candidates in the gold sample that had a reasonable fit; ZTF19abpjicm, ZTF20abjyrxf, and ZTF22aahmovu. Each plot shows the light curve in g-, r-, and i-bands, the best fit combined SALT2 model (solid grey line), and the model for each lensed image (dashed grey lines). The fit parameters are given in Table 6.

The remaining four objects in the gold sample have poor combined two-image SALT2 fits because they are too wide to only display a single peak (ZTF21aablrfe, ZTF22aabifrp, and ZTF22aadeqlh) or they are too noisy to provide a convincing fit (ZTF21abcwuhh). The light curves and corresponding fits are shown in Fig. B.1 of the Appendix. The two-image fits for the wide objects produced larger time delays, for which we would have expected two distinct light curve peaks. However, this is not evident from any of the light curves. We note that a model with three or four images could explain the wider light curves and the absence of bumps, but this is not something that we explicitly modelled (due to the difficulty in constraining the extra model parameters).

Additionally, we calculated the expected image separation (Δθ) from our best fit two-image SALT2 parameters. Once again assuming the lens to be an SIS, we have the following equation for image separation, (derived from e.g. Mörtsell et al. 2020):

$$\begin{array}{c}\hfill \mathrm{\Delta }\theta =0.27″\sqrt{\left(\frac{\mathrm{\Delta }t}{\mathrm{days}}\right)\left(\frac{500\mathrm{Mpc}}{2D_l}\right)\left(\frac{D_s}{D_{\mathit{ls}}}\right)\frac{(r+1)}{(r-1)}\frac{1}{(1+z_l)}},\end{array}$$(4)

where z_l is the lens redshift, and D_l, D_s, D_ls are the angular diameter distances to the lens, source, and between the lens and the source, respectively. Furthermore, the image separation scales with the square of the velocity dispersion, v, of the lensing galaxy, according to the following equation:

$$\begin{array}{c}\hfill \mathrm{\Delta }\theta =1.15″{\left(\frac{v}{200\mathrm{km}{\mathrm{s}}^{-1}}\right)}^2.\end{array}$$(5)

We combined Equation 4 and 5 to calculate v in terms of our fitted parameters:

$$\begin{array}{c}\hfill v=47\mathrm{km}{\mathrm{s}}^{-1}{\left[\left(\frac{\mathrm{\Delta }t}{\mathrm{days}}\right)\left(\frac{500\mathrm{Mpc}}{2D_l}\right)\left(\frac{D_s}{D_{\mathit{ls}}}\right)\frac{(r+1)}{(r-1)}\frac{1}{(1+z_l)}\right]}^{1/4}.\end{array}$$(6)

The calculated values for Δθ and v are displayed in Table 7 for the three objects with reasonable two-image combined SALT2 fits, derived from the parameters in Table 6.

6.1.1. Fitting a two-image SALT2 model to superluminous supernovae and tidal disruption events

To illustrate the difference between the light curves of our candidates and those of our typical contaminants, we also performed fits of the two-image combined SALT2 model to the forced photometry light curves of SLSNe and TDEs from BTS.

For a sample of 55 SLSNe (34 SLSNe-I and 21 SLSNe-II), we found that only two (3.6%) objects could be fit with a reduced χ² of less than five. The ZTF names of the SLSNe are ZTF19aavouyw and ZTF21aaarmti, with reduced χ² values of 3.6 and 1.4, respectively. Both of these objects were SLSNe-I, and we observed a trend that the SLSNe-II fits had higher reduced χ² values and poorer fits by eye. The light curves and corresponding fits are shown in Fig. C.1 of the Appendix. It should be noted that these fits do not look convincing, and the pull values are larger than we would expect for lower reduced χ² values. This is because there are large uncertainties in the resulting SALT2 model (particularly in the g-band) that are included in the χ² calculation. As a result, we conclude that both SLSN-I and SLSN-II are not well modelled by the two-image SALT2 model.

For a sample of 18 TDEs, we found that only one (5.6%) object could be fit with a reduced χ² of less than five. The ZTF name of this TDE is ZTF22aagvrlq, with a reduced χ² value of 2.8. The light curve and corresponding fit are shown in Fig. C.2 of the Appendix. Once again, it is notable that the fit is very poor by eye. This is not surprising, given the characteristic light curve shape of TDEs. Therefore, we conclude that it is unlikely that a typical TDE would be well modelled by the two-image SALT2 model.

This adds further evidence in favour of a glSN Ia interpretation for our candidates ZTF19abpjicm, ZTF20abjyrxf, and ZTF22aahmovu.

6.1.2. Comparison with SLSN and TDE model fit

To strengthen the argument that our gold candidates are glSNe Ia, we also performed fits to SLSN and TDE models as a comparison for the three objects that are well fit by the two-image SALT2 template. To do this, we used the software package redback¹⁷ (Sarin et al. 2024), which relies on the Bayesian parameter inference tool Bilby¹⁸ (Ashton et al. 2019). Within redback, it is possible to fit a wide range of transient models with user-defined priors and use Bayesian inference to determine the best fit parameters for the model.

Although the exact powering mechanism behind SLSNe is unknown, several models exist to explain the high luminosities that we observe. One possible power source is strong circumstellar interaction (CSI) between the ejecta and the material surrounding the transient, which was likely ejected by the progenitor before the explosion. It is possible that the high luminosities of hydrodogen-rich SLSN-II are caused by CSI, which is usually evidenced by strong hydrogen emission lines similar to the narrow Hα lines shown in the spectra of SN IIn (e.g. Smith et al. 2007). Another possible powering mechanism is a central engine, such as the spin-down of a magnetar. There is evidence that the spin-down magnetar model (Kasen & Bildsten 2010; Woosley 2010), where energy is transferred from a strongly magnetised neutron star to the ejecta, is the cause of the high lumninosities of hydrogen-poor SLSN-I. However, since SLSN-II are the only SLSN subclass that have been found in luminous, Milky Way-like galaxies (Gal-Yam 2012), this subtype is the most likely contaminant in our sample. Conversely, SLSN-I are found in lower mass, lower luminosity, metal-poor hosts (Lunnan et al. 2014; Perley et al. 2016; Chen et al. 2017), which would not pass our galaxy cuts described in Sect. 4.3.

Therefore, we performed a redback fit to an SLSN model that is consistent with the characteristics of a hydrogen-rich subtype. Our model combines CSI and nickel-cobalt decay (the csmni model described in MOSFiT; Guillochon et al. 2018). We used the nested sampler nessai¹⁹ (Williams 2021; Williams et al. 2021, 2023) and implemented broad uniform priors. Figure D.1 shows the resulting models for the three promising candidates.

For the TDE case, we performed a redback fit to the tde_fallback model (which is based on the tde model described in MOSFiT; Guillochon & Ramirez-Ruiz 2013; Guillochon et al. 2018; Mockler et al. 2019). We used the nested sampler nestle²⁰ (Barbary 2021) and implemented broad uniform priors. Figure D.2 shows the resulting models for the three promising candidates.

To compare models, we calculated the Bayesian information criterion (BIC; Schwarz 1978) for the light curve fits to the individual SALT2 template, the two-image SALT2 template, SLSN-II model, and the TDE model, which is shown in Table 8.

For the case of ZTF19abpjicm, we see that a single SALT2 fit is the preferred model, closely followed by a two-image SALT2 template. This is likely due to the BIC penalising the two-image model for the additional fit parameters. Additionally, the brightness of ZTF19abpjicm excludes a typical SN Ia explanation, as we shall discuss in Sect 6.4.1. The larger BIC values for the SLSNe and TDE models support our hypothesis that this object is more likely to be a glSN Ia. Additionally, we can see in Figs. D.1 and D.2 that the SLSNe and TDE models fail to capture the behaviour of the light curve.

For ZTF20abjyrxf, the BIC values suggest that a TDE or SLSN interpretation is preferred. The posterior resulting from the TDE model predicts a blackhole mass (M_BH) of $9.{75}_{-0.99}^{+0.89}\times {10}^6{\mathrm{M}}_{\odot }$ and a disrupted star mass (M_star) of $12.{62}_{-2.50}^{+1.57}{\mathrm{M}}_{\odot }$, which are on the upper limits of the mass ranges of known TDEs (Wevers et al. 2017; Ryu et al. 2020; Chang et al. 2024) but are still consistent with expectations. However, without any additional evidence supporting a TDE interpretation (such as the presence of a dust echo), this conclusion is less likely. Additionally, we can see in Fig. D.2 that the model fails to capture the behaviour of the light curve at the pre-peak epoch. Considering the SLSN interpretation, the posterior predicts a high ejecta mass (M_ej) of $5.{42}_{-0.78}^{+1.13}{\mathrm{M}}_{\odot }$ and a high nickel fraction (f_Ni) of $0.{81}_{-0.19}^{+0.13}$. This would correspond to a nickel mass of greater than a solar mass, which would suggest a pair instability supernova (Gal-Yam et al. 2009). This is also a very rare transient event, with a small number of possible candidates (e.g. Schulze et al. 2024), and is therefore very unlikely. As a result, we expect that the two-image SALT2 model for ZTF20abjyrxf is the most likely in this comparison.

Finally, the BIC values for ZTF22aahmovu also suggest that a TDE or SLSN interpretation is preferred. Once again, the posterior for the TDE model predicts a high blackhole mass (M_BH) of $18.{99}_{-1.16}^{+0.68}\times {10}^6{\mathrm{M}}_{\odot }$ and a disrupted star mass (M_star) of $11.{29}_{-1.64}^{+1.50}{\mathrm{M}}_{\odot }$ which are physically reasonable. However, the TDE model once again fails to capture the behaviour of the light curve at the pre-peak epoch in Fig. D.2. The posterior parameters of the SLSN model produce unrealistic results; the model predicts a very low M_ej of $3.{53}_{-0.06}^{+1.43}\times {10}^{-4}{\mathrm{M}}_{\odot }$ which would suggest a failed supernova or a different type of stellar event. As this is not compatible with the observed luminosity, this SLSN-II model does not accurately describe ZTF22aahmovu, and we propose that the two-image SALT2 model is the most likely in this comparison.

6.2. Rise, decline, and duration calculations

We calculated the time of rise, decline, and duration in the r-band for each of the seven gold sample candidates, shown in Table 9. This is a useful characteristic for comparison with contaminants like SLSNe, which typically have a longer rise time and duration. However, we note that glSNe Ia with a brighter second bump would display a longer rise and a shorter decline (such as in the case of ZTF20abjyrxf in Fig. 13).

We defined the rise time as the time it takes the transient to rise from 50% of the peak flux value to the peak (t_peak − t_{50%, r}), and conversely for the decline time (t_{50%, d} − t_peak). Hence, the duration was defined as the sum of the rise and decline times. The r-band was chosen as this band usually had the most detections, due to the regular cadence and redness of the objects. For simplicity, and to avoid introducing any biases from models, we assumed that the largest flux value was the peak of the light curve. The values t_{50%, r} and t_{50%, d} were determined by averaging the time of detection before and after the flux values crossed the 50% threshold. If there was no detection before (or after) the 50% peak flux threshold, then the first (or last) detection date was chosen. From our crude estimations, we expect uncertainties on the rise or decline times of approximately ±1 day (half the typical survey cadence).

A comparison of the duration of the glSN Ia candidates in the gold sample with the duration of SLSNe-I and SLSNe-II from the Bright Transient Survey (Perley et al. 2020) is shown in Fig. 14. We do not see a clear distinction between our glSN candidates and the SLSNe, but we note that the selection of candidates that had a reasonable two-image SALT2 fit displayed a shorter duration.

Fig. 14. Plot comparing the rest-frame duration of the seven gold candidates (red stars for candidates with a reasonable two-image SALT2 fit and gold stars for the remaining objects) with SLSNe-I and SLSNe-II from the ZTF Bright Transient Survey. In this plot, we assume that the catalogue-matched redshift for our candidates belongs to the host (which would be the case if we treated the candidates like potential SLSNe).

6.3. Candidate lens galaxy properties

Next, we estimated the stellar masses of the gold sample candidate’s galaxies to determine whether they are consistent with SLSN host galaxies or lens galaxies. For each candidate in the gold sample, we retrieved science-ready co-added images from the Galaxy Evolution Explorer general release 6/7 (GALEX; Martin et al. 2005), the Sloan Digital Sky Survey DR 9 (SDSS; Ahn et al. 2012), the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS, PS1) DR1 (Chambers et al. 2016), DESI Legacy Imaging Surveys (Legacy Surveys, LS) Dey et al. (2019), the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), and images from the Wide-Field Infrared Survey Explorer (WISE; Wright et al. 2010) processed by Lang (2014). We measured the brightness of the potential lensing galaxy using LAMBDAR²¹ (Lambda Adaptive Multi-Band Deblending Algorithm in R; Wright et al. 2016) and the methods described in Schulze et al. (2021). Tables E.1 to E.7 in the Appendix display the measurements in the different bands.

We modelled the observed spectral energy distribution with the software package Prospector version 1.1 (Johnson et al. 2021)²². We assume a Chabrier IMF (Chabrier 2003) and approximate the star formation history (SFH) by a linearly increasing SFH at early times followed by an exponential decline at late times [functional form t × exp(−t/τ), where t is the age of the SFH episode and τ is the e-folding timescale]. The model is attenuated with the Calzetti et al. (2000) model. The priors of the model parameters are set identically to those used by Schulze et al. (2021). Table 10 summarises the properties of the galaxy.

To put these measurements in context of hosts of other SN classes, we show in Fig. 15 the mass versus star formation rate of the candidate glSNe Ia lens galaxies and the hosts of hydrogen-poor superluminous supernovae (SLSNe-I), interaction-powered supernovae (SNe IIn) and superluminous supernovae (SLSNe-IIn) from the PTF survey (Schulze et al. 2021). The lens galaxies of the known glSNe Ia, SN Zwicky and iPTF16geu, are also plotted for comparison. The potential lens galaxies of the glSN Ia candidates are significantly more massive than the hosts of H-poor SLSNe, but they are consistent with the properties of SLSNe-IIn and SNe IIn hosts. However, they are also consistent with previously observed galaxy lenses, so we cannot exclude them as potential lensing galaxies.

Fig. 15. Star formation rate and stellar mass of the glSN Ia candidate lens galaxies (black), known glSNe Ia lens galaxies (red), and SLSNe-I (green), SLSNe-II (blue) and SNe IIn (grey) host galaxies. The lens galaxy of the glSN Ia candidates are similar to hosts of SLSNe-IIn and SNe IIn, suggesting that some glSN Ia candidates could be SLSNe-IIn. However, the masses are also consistent with previously observed galaxy lenses.

6.4. Individual candidate discussion

6.4.1. ZTF19abpjicm

ZTF19abpjicm has a good two-image combined SALT2 fit (with a reduced χ² of 0.86), and a relatively short duration of 31.4 days. From our BIC comparison to SLSN and TDE models in Sect. 6.1.2, the two-image combined SALT2 model is the preferred interpretation. Therefore, we state that the light curve is not typical for contaminants such as SLSNe or TDEs.

Fitting the light curve with a single SALT2 template also produces a reasonable fit, with a larger x₁ of 2.40 ± 0.62 and a reduced χ² of 1.04, which indicates the object could be an unlensed SN Ia. The photometric and spectroscopic redshift values agree, which suggests that we can trust that this object is over-luminous. This means that it cannot be a normal SN Ia. The only SN Ia subclasses known to surpass a peak M_B of −20 mag are Ia-CSM and 03fg-like. As was discussed in Sect. 5.1, 03fg-like SNe are commonly found in low-mass, star-forming galaxies, which is not consistent with the galaxy properties for ZTF19abpjicm given in Table 10. Also, we do not see obvious signs of CSM interaction in the light curve (there are no plateau phases, late-time bumps, or rebrightening in the forced photometry) for a SN Ia-CSM. However, SNe Ia-CSM are a diverse subtype, and it is possible that we do not see CSM signatures such as bumps and plateaus due to the SN being close to the detection limit.

From the relative magnification of the two fitted images alone, we calculated a total magnification of μ_tot = 3.6 ± 1.3 in Sect. 6.1. We can also estimate the total magnification by comparing the peak apparent magnitude to the apparent magnitude we would expect from a normal SN Ia (with M_B = −19.4 mag) at a redshift of 0.44 (as predicted by the combined SALT2 model). Although this requires the assumption of a cosmological model (here we assume a flat ΛCDM model with H₀ = 67.4 km s⁻¹Mpc⁻¹ and Ω_M = 0.315; Planck Collaboration VI 2020), and relies on the fitted redshift, it is independent of the previous method that relies solely on the flux ratios of the images. From this method, we calculate a total magnification of μ_tot = 4.7. This is consistent with the μ_tot estimated from the parameters of the combined SALT2 model, where the deviation is within one standard deviation of the predicted value.

Also in Sect. 6.1, we calculated an image separation (Δθ) of 1.1″ for the case of ZTF19abpjicm, from the estimate of the source redshift and the time delays. This is within the range of values we would expect for typical lensing systems and agrees with our observation that the separate images were unresolved by ZTF. We also calculated a lens galaxy velocity dispersion (v) of 190 km s⁻¹, which is typical for strong lens galaxies with the stellar mass quoted in Table 10.

As a result, we have justified that our two-image model fit parameters are physical and consistent with our photometric measurements of the light curve and the suspected lens.

6.4.2. ZTF22aahmovu

ZTF22aahmovu has a good two-image combined SALT2 fit with a reduced χ² of 1.38. From our comparison to SLSN and TDE models in Sect. 6.1.2, the TDE model failed to accurately model the light curve in the pre-peak epoch and the SLSN model predicted an unphysical ejecta mass. As a result, we claim that the light curve is not typical for contaminants such as SLSNe or TDEs. Additionally, the light curve does not resemble that of a normal SN Ia, so we can also exclude over-luminous Ia subclasses as being a contaminant.

From the relative magnification of the two fitted images alone, we calculated a total magnification of μ_tot = 4.3 ± 0.6 in Sect. 6.1. As with ZTF19abpjicm, we also estimated the total magnification by comparing the peak apparent magnitude to the apparent magnitude that we would expect from a normal SN Ia at a redshift of 0.35 (assuming the same flat ΛCDM cosmology). From this method, we calculate a total magnification of μ_tot = 2.8. These values are consistent within 3σ, which is reasonable when we consider that there are additional uncertainties we have not accounted for (such as microlensing by intervening stars, which could increase or decrease the magnification from the expected value).

The calculated image separation from our best fit two-image SALT2 parameters was 1.7″. We also determined a lens galaxy velocity dispersion of 240 km s⁻¹ for this system. These values are typical for strong lensing galaxies with the stellar mass quoted in Table 10.

Therefore, like with ZTF19abpjicm, we have justified that our two-image model fit parameters are physical and consistent with our photometric measurements of the light curve and the suspected lens.

6.4.3. ZTF20abjyrxf

ZTF20abjyrxf also has a two-image combined SALT2 fit (with a reduced χ² of 2.53). The model displays a slight decrease in the g-band brightness as the fit transitions from being dominated by the first image to the second, which is consistent with what we observe in the light curve. Additionally, the light curve is broad with a slower rise time (approximately 58.4 days). This is attributed to the relatively equal flux contributions from the two images, which causes the peak of the model to occur later (at the peak of the second image). However, this slower rise or plateauing feature could be characteristic of an SLSN light curve. From Fig. 15 in Sect. 6.3, the galaxy associated with ZTF20abjyrxf could be consistent with the host of a SLSN-II(n). However, from our comparison to an SLSN model in Sect. 6.1.2, the best fit model predicts a nickel mass of greater than a solar mass, which would suggest a pair-instability supernova (an extremely rare and therefore unlikely event). The TDE model fit from Sect. 6.1.2 also failed to accurately model the light curve in the pre-peak epoch, and so we assert that this object is not a typical TDE.

In Sect. 6.1, we estimated a total magnification of μ_tot = 20.4 ± 6.2. This is quite large, similar to the magnification we observed from SN Zwicky (with μ_tot = 24.3 ± 2.7). For this to be correct, we would expect the source redshift to be ∼2.8 times the value of the lens redshift (assuming the same flat ΛCDM cosmology as for ZTF19abpjicm), which is possible, although not consistent with the best fit redshift parameter from the SALT2 two-image fit. We would also expect larger image separation for ZTF20abjyrxf, although this is highly dependent on the source redshift (it could vary from 2.5″ for z_s = 2.3 z_l, to 6.7″ for z_s = 0.24 as was determined from the SALT2 fit). Considering that the spatial resolution of ZTF is approximately 1″ (combined with the 2″ seeing at Palomar Observatory), it should be possible to distinguish the separate images in the case of a 6.7″ separation. However, we did not find any evidence of spatial variability in the ZTF alert images.

It is also possible that our simplified lens model does not accurately resemble the system of ZTF20abjyrxf. This has been observed before; in the case of SN Zwicky, the observed relative flux contributions for two out of the four images did not agree with the model predictions. Furthermore, we cannot ignore the possibility that four images could be present in our case, which would change the relative image magnifications and the resulting flux ratios. As a result, we conclude that, if ZTF20abjyrxf is a glSN Ia, it is not well described by a simple two-image system with a SIS lens, and further study of the potential lens galaxy is required.

6.4.4. ZTF21aablrfe

It was not possible to fit ZTF21aablrfe with a reasonable two-image SALT2 model, due to the width of the light curve and the longer rise time (approximately 60 days). It also had a poor fit to the single-image SALT2 model (with a reduced χ² of 4.7). Additionally, it is further away from the DESI galaxy, at a separation of 2.4″. The combination of these characteristics suggests that this object could be a SLSN-II(n) candidate.

Despite this, the light curve is sufficiently red to pass the colour cuts and the galaxy counterpart is not consistent with the hosts of SLSNe-I (as is shown in Fig. 15). A four-image SALT2 model might be able to more accurately match the light curve we observe, although this would require us to know more about the lens system.

6.4.5. ZTF22aabifrp

ZTF22aabifrp has an interesting LS cutout, which shows two galaxies close together (less than three arcseconds apart). The DESI catalogue-matched redshift (z = 0.1886) belongs to the object above the green circle plotted in Fig. 11. The photometric redshift for the lower galaxy, which is closer to ZTF22aabifrp at a separation of 0.1″, is given as z = 0.25 ± 0.04 in the LS DR9 catalogue. Therefore, the z = 0.19 object could be a lens for the host at z = 0.25. However, the redshifts of the two galaxies are consistent within 1.5 standard deviations, and so this claim is not statistically significant.

Similarly to ZTF21aablrfe, it was not possible to fit ZTF22aabifrp with a reasonable two-image SALT2 model, due to the width of the light curve. From our fit, we would expect to see two distinct peaks for the separate images; however, there is no evidence for this in the data (unfortunately, there is a large period of no observations where the valley between the two peaks would be). This is another case where a four-image SALT2 model might be able to more accurately match the light curve we observe. Also like ZTF21aablrfe, the light curve is sufficiently red to pass the colour cuts and the closest galaxy is not consistent with the hosts of SLSNe-I. Despite this, we cannot exclude the possibility that this object is a SLSN-II(n).

6.4.6. ZTF22aadeqlh

The LS cutout for ZTF22aadeqlh in Fig. 12 displays two possible galaxy lenses or hosts (a larger blob to the left and a yellower blob to the right). It is not certain which galaxy is the lens or the host, although it appears that the transient is more aligned with the object on the left. The two galaxies are at approximately the same redshift according to DESI (z = 0.32), so they may be a combined two-galaxy lens.

ZTF22aadeqlh was reported to be a candidate SLSN by Gkini et al. (2022) in an AstroNote. This is possible, because of the broad light curve and lack of a reasonable two-image SALT2 fit. Similarly to ZTF21aablrfe and ZTF22aabifrp, the light curve is sufficiently red to pass the colour cuts and the closest galaxy is not consistent with the hosts of SLSNe-I. However, once again, we cannot exclude the possibility that this object is a SLSN-II(n).

6.4.7. ZTF21abcwuhh

ZTF21abcwuhh has a single-image SALT2 fit with a reduced χ² of 1.0. However, the light curve is too faint (and, therefore, too noisy) to provide a convincing two-image combined SALT2 fit. The photometric and spectroscopic redshift values agree, which suggests that we can trust that this object is over-luminous. ZTF21abcwuhh may be a glSN Ia candidate, but the limited data means that we cannot further deduce the true nature of this faint transient.

6.5. Summary of gold sample candidates

We present two candidates, ZTF19abpjicm and ZTF22aahmovu, as likely glSNe Ia. The observed characteristics of the two objects are not typical for contaminants such as SLSNe, TDEs, or over-luminous SN Ia subclasses. They can also be well modelled by a two-image combined SALT2 fit, and the resulting parameters are consistent with our photometric measurements and adhere to the expected physical characteristics of lensing systems.

An additional candidate, ZTF20abjyrxf, was also well modelled by a two-image combined SALT2 fit. However, the resulting parameters suggest that we should have observed a larger image separation (and therefore should have observed the two images as separate sources in ZTF). As this was not the case, we conclude that ZTF20abjyrxf is not well described by a simple two-image system with an SIS lens. Further study of the potential lens galaxy is required for us to present ZTF20abjyrxf as a likely glSN Ia.

ZTF21aablrfe, ZTF22aabifrp, and ZTF22aadeqlh were not well modelled by a two-image combined SALT2 fit due to the width of their light curves (we only see evidence for a single peak in the data; however, the model predicts that we should see two distinct peaks). Despite this, we see evidence of two galaxies at different redshifts (a potential lens and host) in the LS cutout for ZTF22aabifrp. We also observe two galaxies at approximately the same redshift for ZTF22aadeqlh, which could be a combined two-galaxy lens. As a result, we conclude that these candidates are not well described by a simple two-image system with an SIS lens, but a more complex model (e.g. with four images) could potentially provide a more accurate fit to the observed data. As with ZTF20abjyrxf, follow-up of the potential lens and/or host galaxies is required.

7. Conclusions

In this work, we presented a systematic search of ZTF archival data between June 1 2019 and September 1 2022 for strongly lensed type SNe Ia. We created an analysis pipeline using the AMPEL processing platform, which allowed us to query the ZTF alert photometry, apply filters for SN-like objects, and perform SALT2 template fits. Then, we applied some initial cuts and acquired forced photometry for the remaining candidates. Finally, we narrowed down our most likely candidates through two selection methods: by cross-matching to the DESI spectroscopic redshift catalogues, or applying stricter cuts motivated by simulations and previous observations. Additionally, we performed a similar analysis on the ZTF Bright Transient Survey sample of supernovae, to estimate the amount of contamination we should expect from each class.

The conclusions of our analysis and the plans for future work are as follows:

1. We have shown that it is possible to filter the alert photometry of a large optical survey to find transients that match the expected characteristics of glSNe Ia. We devised an initial filtering algorithm to reduce 31930 SN-like alerts to 7075 candidates that passed some preliminary cuts. We reduced the 7075 initial candidates to 50 final glSN Ia candidates by applying additional cuts and matching to the DESI spectroscopic catalogue.

2. Of these final candidates, 27 were the result of a DESI cross-match (spectroscopic sample) and 27 were the result of applying stricter cuts based on simulations (photometric sample). A careful inspection of each of the spectroscopic candidates is necessary to exclude the possibility that they are SLSNe. Both the spectroscopic and photometric samples will be targeted in future spectroscopic follow-up missions to determine the lens and/or host redshifts.

3. The vast majority of the final candidates were found within the lens redshift range of 0.2 < z < 0.4. Furthermore, they are all dimmer than SN Zwicky and iPTF16geu, which suggests that, if they are truly glSNe Ia, they belong to a less magnified population.

4. Seven candidates were present in both the spectroscopic and photometric samples, and these were discussed in more depth in this work. We present two candidates, ZTF19abpjicm and ZTF22aahmovu, as likely glSNe Ia. They are well modelled by a combined two-image SALT2 template and the model fit parameters are consistent with the physical parameters we would expect for strong lensing systems. From this two-image modelling, we estimate time delays of 22 ± 3 and 34 ± 1 days for the two events, respectively, which suggests that we have uncovered a population of glSNe Ia with longer time delays. The remaining five objects are still promising candidate glSNe Ia, but we would require more information about their lensing systems to model them more accurately.

5. The two most promising candidates, ZTF19abpjicm and ZTF22aahmovu, should be targeted in a follow-up programme to confirm that they occurred in a strong lensing system. To do this, we suggest acquiring high-resolution photometry of the system, to search for evidence of a strongly lensed background host galaxy. Alternatively, if the host galaxy is not lensed, we suggest acquiring deeper spectra of the potential lens galaxy, to search for additional lines from a background host galaxy.

In this paper, we have highlighted the importance of lensed transients (particularly standardisable candles like SNe Ia) for cosmology and presented a pipeline that is both efficient and sensitive enough to parse full alert streams from optical observatories to search for glSNe Ia. The methodology presented in this systematic search is currently applied as a live search in ZTF and could be the foundation for future live searches in LSST and beyond. To fully remove the sources of contamination, spectroscopic follow-up while the transient is live is necessary. Despite this, we have shown that it is still possible to eliminate contaminants based on light curve modelling and photometric observations of the galaxy counterpart, and, as a result, we present two likely glSNe Ia.

Data availability

Table A.2 is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/694/A146. All ZTF light curve data is publicly available at https://www.ztf.caltech.edu/ztf-public-releases.html. The instructions for how to install and run the AMPEL platform can be found at https://github.com/AmpelAstro/Ampel-HU-astro/. The specific AMPEL workflow we used to obtain the alert photometry, cross-match to galaxy catalogues, and perform initial filtering is found at https://github.com/AmpelAstro/Ampel-HU-astro/blob/main/notebooks/run_lensing_query.ipynb. The forced photometry pipeline for ZTF data that we utilised can be found at https://github.com/simeonreusch/fpbot. The BTS sample, which we utilised to estimate the contaminants present in our sample, is publicly accessible at https://sites.astro.caltech.edu/ztf/bts/explorer.php. All external data sources are referenced in the main text. Supplemental data are accessible at https://zenodo.org/records/11105340.

Acknowledgments

Based on observations obtained with the Samuel Oschin Telescope 48-inch and the 60-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under Grants No. AST-1440341 and AST-2034437 and a collaboration including current partners Caltech, IPAC, the Weizmann Institute of Science, the Oskar Klein Center at Stockholm University, the University of Maryland, Deutsches Elektronen-Synchrotron and Humboldt University, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National Laboratories, IN2P3, University of Warwick, Ruhr University Bochum, Northwestern University and former partners the University of Washington, Los Alamos National Laboratories, and Lawrence Berkeley National Laboratories. Operations are conducted by COO, IPAC, and UW. This work has been supported by the research project grant “Understanding the Dynamic Universe” funded by the Knut and Alice Wallenberg Foundation under Dnr KAW 2018.0067, Vetenskapsrådet, the Swedish Research Council, project 2020-03444, and the G.R.E.A.T research environment, project number 2016-06012. This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of High-Energy Physics, under Contract No. DE–AC02–05CH11231, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Science and Technology of Mexico (CONACYT); the Ministry of Science and Innovation of Spain (MICINN), and by the DESI Member Institutions: https://www.desi.lbl.gov/collaborating-institutions. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the U. S. National Science Foundation, the U. S. Department of Energy, or any of the listed funding agencies. The authors are honored to be permitted to conduct scientific research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. S. Schulze is partially supported by LBNL Subcontract NO. 7707915. M. W. Coughlin acknowledges support from the National Science Foundation with grant numbers PHY-2308862 and PHY-2117997.

References

Appendix A: Data for candidates that passed selection criteria

Table A.1 shows the photometric redshift errors from the LS and PS1-STRM catalogues for the 27 candidates from selection method S2 (the photometric sample). As is discussed in Sect. 4.1, the redshifts in our analysis are averaged from a selection of photometric galaxy catalogues. All our candidates had a photometeric redshift value from either the LS or the PS1-STRM catalogues, and these were the catalogues with the lowest uncertainties typically. Therefore, we selected the lowest photometric redshift value from either catalogue to propagate to peak absolute magnitude error.

Table A.2 details the measured parameters of the 125 candidates that passed the selection criteria in Sects. 4.2 and 4.3, before the final visual inspection phase. This includes the ZTF ID, right ascension (RA), declination (dec), photometric redshift (phot. z), DESI redshift, angular separation from the catalogue counterpart (ang. sep.) for both the photometric and DESI redshifts, peak absolute B-band magnitude from SALT2 (for both the photometric and DESI redshifts), SALT2 x₁ and associated error, SALT2 c and associated error, whether it passed selection method 1 with DESI cross-matching (detailed in Sect. 4.2) and/or selection method 2 with additional cuts (detailed in Sect. 4.3), and the reason it was rejected from the final candidate sample (if any).

Appendix B: Two-image SALT2 fits to gold sample candidates

As is discussed in Sect. 6.1, we performed fits of the two-image combined SALT2 model to forced photometry light curves of all seven of the gold sample candidates. Three of them converged to reasonable fits: ZTF19abpjicm, ZTF20abjyrxf, and ZTF22aahmovu. The remaining four objects in the gold sample have poor combined two-image SALT2 fits because they are too wide to only display a single peak (ZTF21aablrfe, ZTF22aabifrp, and ZTF22aadeqlh) or they are too noisy to provide a convincing fit (ZTF21abcwuhh). Figure B.1 shows the fits for these four objects.

Fig. B.1. Two-image combined SALT2 model to forced photometry light curves of ZTF21aablrfe (upper left), ZTF21abcwuhh (upper right), ZTF22aabifrp (lower left), and ZTF22aadeqlh (lower right).

Appendix C: Two-image SALT2 fits to SLSNe and TDEs in BTS

As is discussed in Sect. 6.1.1, we performed fits of the two-image combined SALT2 model to forced photometry light curves of SLSNe and TDEs from the BTS. Figure C.1 shows the fits for ZTF19aavouyw and ZTF21aaarmti, the two SLSNe-I with a reduced χ² of less than five. Figure C.2 shows the fit for ZTF22aagvrlq, the only TDE with a reduced χ² of less than five.

Fig. C.1. Two-image combined SALT2 model to forced photometry light curves of the SLSNe-I, ZTF19aavouyw (left) and ZTF21aaarmti (right).

Fig. C.2. Two-image combined SALT2 model to the forced photometry light curve of the TDE, ZTF22aagvrlq.

Appendix D: SLSN and TDE model fit with redback

As is described in Sect. 6.1.2, we performed fits to a SLSN-II model (with a powering mechanism that combines CSI and nickel) and a TDE model for the three candidates that had good fits to the two-image SALT2 model. The model fit light curves are shown in Fig. D.1 and Fig. D.2 for the SLSN-II and TDE models, respectively.

Fig. D.1. CSM and nickel model fits to ZTF19abpjicm (top), ZTF20abjyrxf (middle), and ZTF22aahmovu (bottom) with redback. The grey lines correspond to one hundred random models randomly drawn from the posterior and the maximum likelihood model.

Fig. D.2. TDE model fits to ZTF19abpjicm (top), ZTF20abjyrxf (middle), and ZTF22aahmovu (bottom) with redback. The grey lines correspond to one hundred random models randomly drawn from the posterior and the maximum likelihood model.

Appendix E: Galaxy photometry of gold sample candidates

Tables E.1 to E.7 show the photometry we acquired for the potential lens or host galaxies of the seven candidates in the gold sample (as is described in Sect. 6.3). From this photometry, we calculated the stellar masses and star formation rate for each candidate.

All Tables

All Figures

	Fig. 1. Flowchart to illustrate the stages of the method, as is described in Sect. 4. The red numbers to the bottom right of the boxes indicate the number of transients that remained after that stage of the analysis. The asterisk represents a stage that has not been completed during our study and is left for future work.
In the text

	Fig. 2. Histogram illustrating the distribution of median angular image separation of simulated glSNe Ia, using data from the catalogues of Goldstein et al. (2019). The dashed black line indicates a 3″ cut that selects 98.4% of the distribution.
In the text

Fig. 3. Plot of the LS g − z versus g − W1 galaxy colours for the candidates remaining after the base cuts (orange), compared to contours from confirmed and candidate lens galaxies (blue), and QSOs detected by WISE (green). The data was acquired by cross-matching the AGN and lens catalogues with LS (Duncan 2022). The dashed lines indicate the colour cuts which were devised to select likely elliptical and lens galaxies, shown in Table 2. A single redshift bin of 0.2 < z < 0.3 is shown for clarity. The two star markers show the lens galaxy colours for iPTF16geu (grey) and SN Zwicky (black), which had lens redshifts of 0.216 and 0.226, respectively. For SN Zwicky, the g − z versus g − W1 colours from Pan-STARRS are shown instead, because the lens galaxy was not present in the LS catalogue).

In the text

Fig. 4. Light curve of ZTF20abgdrea, a candidate that passed S1 and S2, but was removed during visual inspection due to the long-term variability that suggests it is a likely AGN or QSO. The dotted black lines indicate the section of the light curve that was selected by the pipeline for the SALT2 fit. The baseline was not removed from this forced photometry, to more clearly illustrate the variability.

In the text

	Fig. 5. Light curve of ZTF20acyxxfo, a candidate that passed S1 and S2, but was removed during visual inspection because it was identified as a likely TDE due to the approximate power-law decline.
In the text

	Fig. 6. Light curve of ZTF20acudrme, a candidate that passed S2, but was removed during visual inspection because the long duration indicates it is a likely SLSN.
In the text

	Fig. 7. Light curve of ZTF20acqhhyr, a candidate that passed S1, but was removed during visual inspection because of the presence of multiple bumps in the light curve of not-equal width, which suggests that CSM interaction is causing the high luminosity.
In the text

	Fig. 8. Histogram illustrating the peak B absolute magnitude distribution for the 5 main classes of transients (SN Ia, SN II, SN Ib/c, SLSNe, and TDEs) recorded by the BTS. The dashed black line illustrates the objects with MB < −20 mag.
In the text

Fig. 9. Distribution of the peak B-band absolute magnitude at the lens redshift versus the lens redshift for the transients that passed the base cuts described in Sect. 4.1 (blue contours), the candidates passing the two selection methods described in Sects. 4.2 and 4.3 (green and orange crosses, respectively), and for the already discovered glSNe Ia SN Zwicky and iPTF16geu (black and grey stars, respectively). The upper limit of the x axis has been adjusted for clarity, omitting one higher redshift object from S1 (ZTF20acirhoc).

In the text

Fig. 10. Set of histograms displaying the distributions for peak mB, peak MB, lens z, SALT2 x1, SALT2 c, and angular separation values for the transients that passed the base cuts described in Sect. 4.1 (orange bars), and the candidates passing the two selection methods described in Sects. 4.2 and 4.3 (solid black and red lines, respectively). The dotted line with a star marker indicates the corresponding values for SN Zwicky for each distribution parameter.

In the text

	Fig. 11. Light curve plots for the seven candidates in the gold sample that had a spectroscopic redshift from the DESI catalogue and passed all the cuts described in Sect. 4.3. Each plot shows the light curve, the lens redshift, and peak absolute magnitude at the lens redshift (from the DESI cross-match).
In the text

	Fig. 12. Image cutouts of the LS field at the location of the transient (indicated by a green circle) for the seven gold sample candidates (the ZTF name for each transient is shown in the top left of the cutout). The galaxies we see are likely the lensing galaxies, and the host is more distant and faint.
In the text

	Fig. 13. Two-image combined SALT2 template fits for the three candidates in the gold sample that had a reasonable fit; ZTF19abpjicm, ZTF20abjyrxf, and ZTF22aahmovu. Each plot shows the light curve in g-, r-, and i-bands, the best fit combined SALT2 model (solid grey line), and the model for each lensed image (dashed grey lines). The fit parameters are given in Table 6.
In the text

Fig. 14. Plot comparing the rest-frame duration of the seven gold candidates (red stars for candidates with a reasonable two-image SALT2 fit and gold stars for the remaining objects) with SLSNe-I and SLSNe-II from the ZTF Bright Transient Survey. In this plot, we assume that the catalogue-matched redshift for our candidates belongs to the host (which would be the case if we treated the candidates like potential SLSNe).

In the text

Fig. 15. Star formation rate and stellar mass of the glSN Ia candidate lens galaxies (black), known glSNe Ia lens galaxies (red), and SLSNe-I (green), SLSNe-II (blue) and SNe IIn (grey) host galaxies. The lens galaxy of the glSN Ia candidates are similar to hosts of SLSNe-IIn and SNe IIn, suggesting that some glSN Ia candidates could be SLSNe-IIn. However, the masses are also consistent with previously observed galaxy lenses.

In the text

	Fig. B.1. Two-image combined SALT2 model to forced photometry light curves of ZTF21aablrfe (upper left), ZTF21abcwuhh (upper right), ZTF22aabifrp (lower left), and ZTF22aadeqlh (lower right).
In the text

	Fig. C.1. Two-image combined SALT2 model to forced photometry light curves of the SLSNe-I, ZTF19aavouyw (left) and ZTF21aaarmti (right).
In the text

	Fig. C.2. Two-image combined SALT2 model to the forced photometry light curve of the TDE, ZTF22aagvrlq.
In the text

	Fig. D.1. CSM and nickel model fits to ZTF19abpjicm (top), ZTF20abjyrxf (middle), and ZTF22aahmovu (bottom) with redback. The grey lines correspond to one hundred random models randomly drawn from the posterior and the maximum likelihood model.
In the text

	Fig. D.2. TDE model fits to ZTF19abpjicm (top), ZTF20abjyrxf (middle), and ZTF22aahmovu (bottom) with redback. The grey lines correspond to one hundred random models randomly drawn from the posterior and the maximum likelihood model.
In the text