© The Authors 2025

1. Introduction

Type Ia supernovae (SNe Ia) span a range of peak absolute magnitudes that can be standardised using properties of their light curves around the peak (e.g. Phillips 1993; Phillips et al. 1999). However, besides normal SNe Ia, there are also events that are spectroscopically and/or photometrically different, creating their own subclasses and that cannot be standardised in the normal way. One such subclass are those that are thought to be interacting with circumstellar material (CSM), called SN Ia-CSM events. The first SN Ia-CSM reported was SN 2002ic which had Hα and Hβ emission in its spectrum around the peak, which was suggested to originate from an interaction with CSM (Hamuy et al. 2003a,b). The optical light curve of SN 2002ic was found to be declining much more slowly than was expected for normal SNe Ia (Wood-Vasey et al. 2004), and it was found to be consistent with about 1.3 M_⊙ of CSM being present around the progenitor system (Nomoto et al. 2005).

It has been suggested that there is a link between the Ia-CSM and 91T subclasses (bright SNe Ia with a distinct pre-peak spectroscopic evolution, showing a blue pseudo-continuum with two strong Fe III absorption multiplets instead of intermediate mass element lines, and predominantly occurring in younger stellar populations; Filippenko et al. 1992; Ruiz-Lapuente et al. 1992) due to similarities in their peak brightness and spectroscopy before the start of the interaction (Leloudas et al. 2015). In some cases, a 91T-like SN could start interacting with CSM hundreds of years after the explosion, as has been suggested for Kepler’s SN (Patnaude et al. 2012; Katsuda et al. 2015).

Currently the Ia-CSM subclass contains several dozen members (Aldering et al. 2006; Silverman et al. 2013; Sharma et al. 2023). In all cases, the CSM interaction started within about two months of explosion, significantly altering the light curve and leading to the (re-)classification as a Ia-CSM after spectroscopic confirmation. The amount and duration of the interaction varies quite significantly from one object to another. While some mainly feature an unusually slow decline rate (e.g. SN 2018gkx and SN 2020xtq; Sharma et al. 2023), other objects have a plateau for several 100 d before starting to fade away slowly (e.g. SN 2020aekp; Sharma et al. 2023).

SN 2011km (PTF11kx) is a well-observed Ia-CSM event, and showed signs of a complex CSM consisting of multiple shells with which the SN ejecta interact. Dilday et al. (2012) showed that its photometry is similar to 91T-like objects before the CSM interaction begins. They explained the geometry of the system using a symbiotic nova progenitor system. Another interesting Ia-CSM is SN 2020eyj, which Kool et al. (2023) present as the first detection of a Ia-CSM interacting with He-rich material. Up until then all members of the subclass had shown strong Hα emission and weak He signatures. SN 2020eyj, however, showed little to no H present in the CSM. This was also the first time a SN Ia was detected in the radio. Non-detections in normal SNe Ia suggest a clean environment for the ejecta to expand in, while in this case there is a lot of material present. Considering this, Kool et al. (2023) suggest the progenitor system to have been a He star and white dwarf (WD) binary.

The progenitor systems and mechanisms responsible for triggering the explosions of SNe Ia are still not clear. The progenitor system could be single degenerate, with a main sequence or evolved non-degenerate star accompanying a WD (Whelan & Iben 1973; Nomoto 1982). Mass is transferred from the companion to the WD until the explosion is triggered. In the double degenerate scenario, both stars are WDs and merge or interact (Iben & Tutukov 1984; Webbink 1984). The different progenitor scenarios can be separated into two categories of explosion models. In classical models the WD comes close to or reaches the Chandrasekhar mass (M_Ch ∼ 1.4 M_⊙, Chandrasekhar 1931) by accreting matter from the secondary star before exploding (Whelan & Iben 1973), possibly through a delayed detonation (Khokhlov 1991; Mazzali et al. 2007). The second category assumes a situation where the explosion is triggered in a lighter WD, resulting in a sub-Chandrasekhar mass explosion. This can be done in, for example, the double detonation models where an accreted layer of surface material ignites and explodes, compressing the WD and igniting the core causing a second explosion (Taam 1980; Livne & Arnett 1995; Shen & Bildsten 2009; Fink et al. 2010), or in mergers with the core of an evolved star in the core-degenerate scenario (Kashi & Soker 2011). Other models involve the collision or violent merger of two WDs (Rosswog et al. 2009; Pakmor et al. 2010, 2012).

The CSM identified to be present in a SN Ia-CSM could be created by different mechanisms depending on the progenitor system, and its composition depends on the type of donor star present in the system. In the single degenerate scenario, CSM rich in hydrogen can be created by a WD generating a fast wind, blowing away a part of the material it received from the mass transfer (Nomoto et al. 2005). In the double-degenerate scenario, part of the tidally disrupted secondary WD becomes unbound from the system, creating (H-poor) CSM and may be able to produce detectable signatures depending on the time between the tidal disruption and the SN Ia explosion (Raskin & Kasen 2013). When the SN ejecta reach and start to sweep up the CSM surrounding the system, the interaction is revealed as an additional source of light that brightens and alters the light curve of the SN. Depending on the distance of the CSM to the explosion site, there may be a significant delay between the explosion and the start of the interaction.

In an effort to systematically search for late-time CSM interaction, Graham et al. (2019b) looked at old (≥1 year) SNe using the Hubble Space Telescope (HST). They focussed their search on subclasses that are associated with CSM interaction, such as 91T-like SNe Ia (Leloudas et al. 2015). Out of 72 targets, only ASASSN-15og and SN 2015cp were found to show late-time CSM interaction. ASASSN-15og is a Type IIn SN with detected CSM interaction around its peak, and was used as a control object. SN 2015cp had been classified as a 91T-like SN Ia without signs of CSM interaction around its peak. This showed that CSM interaction may start much later after the explosion, and may be systematically missed due to SNe Ia not being actively followed at these phases. From a progenitor point of view this means that material can be ejected from the system prior to the explosion, potentially giving it time to travel further before being caught up by the SN ejecta.

Dubay et al. (2022) use archival UV-band data from the Galactic Evolution Explorer (GALEX) to look for late-time CSM interaction in SNe Ia. Out of a sample of 1080 SNe Ia, 4 were detected in the UV near peak, but none showed signs of late-time CSM interaction. They show that this type of CSM interaction is rare, occurring between 500 to 1000 d after the initial discovery of the SN in < 5% of the SNe Ia at a strength similar to SN 2015cp, and a decreasing percentage as the interaction gets stronger.

With today’s large sky surveys such as the Zwicky Transient Facility (ZTF; Bellm et al. 2019; Graham et al. 2019a; Masci et al. 2019; Dekany et al. 2020) and Asteroid Terrestrial-impact Last Alert System (ATLAS; Tonry et al. 2018), transient events are discovered and followed automatically until they fade below the detection limit. This strategy is extremely efficient in finding and cataloguing transients and is a reliable method to find rare subclasses and interesting objects, assuming their defining features can be identified before they fade away. However, depending on the phase of a Ia-CSM at which the interaction begins (potentially more than one year after the peak), it may be systematically missed as the SN is no longer actively followed at these phases, or the detections could be close to the detection threshold and thus not necessarily be associated with the original SN.

ZTF has covered the entire northern sky above −30° declination every 2–3 nights in three broad optical bands (gri) up to limiting magnitudes of ∼20.5 mag since early 2018. A survey this deep and extensive both in space and time coverage may have detected late-time CSM interaction in a SN Ia that has gone unnoticed due to being close to the detection limit. We attempt to push this limit as much as possible by binning the post-SN observations together for each confirmed SN Ia observed with ZTF before December 2020. While binning data will reduce our temporal resolution, it is traded for deeper detections and upper limits.

In Sect. 2, we introduce the sample we use in our search for optical signals for late-time CSM interaction in the ZTF data stream. In Sect. 3, we present our custom pipeline for identification of late-time flux excesses and set up a simulation to test it and estimate the detection efficiency of our pipeline. Sect. 4 shows the result of running our sample through our pipeline, and provides further investigation on some interesting objects. These results are discussed in Sect. 5, and we conclude in Sect. 6. A flat ΛCDM cosmology for H₀ = 67.7 km s⁻¹ Mpc⁻¹ and Ω_m = 0.310 (Planck Collaboration VI 2020) is assumed where required.

2. Data

Our aim is to look for late-time (> 100 d after peak brightness) signatures of CSM interaction in the largest sample of SNe Ia to date. This has been obtained by the ZTF. We are particularly interested in events that appear to be normal SNe Ia from their spectra and light curves around the peak but may display signs of late-time interaction, as seen in SN 2015cp (Graham et al. 2019b). Our starting sample is 3628 events that were discovered by ZTF from March 2018 to October 2020 (hereafter the ZTF data release 2, ZTF DR2). Each event is spectroscopically classified as a SN Ia or one of its sub-classes. An overview of the ZTF DR2 is presented in Rigault et al. (2025), including the sample definition, properties, and use for cosmology. In this study, since we are searching for likely rare signatures of interaction in the ZTF light curves, we are as inclusive as possible in our sample definition and include all SNe Ia in the DR2 covering a redshift out to z = 0.288.

2.1. ZTF light curve data

ZTF observes in three optical bands gri on a 2–3 day cadence. Reference images, mainly made using observations at the start of ZTF, are subtracted from the science images using the ZOGY image subtraction algorithm (Zackay et al. 2016) to produce difference images. We use forced photometry at the transient location on the difference images using ZTFFPS (Reusch 2020) to get a measure of the observed flux at each epoch. This includes non-detections before each SN was first detected and after each SN has faded below detection limits.

Light curve quality cuts on specific light curve points are applied as in Smith et al. (2025). We do not correct the light curves for Milky Way extinction in our initial analysis but do consider it when focussing on specific objects of interest in Sect. 4.

Another approach for extracting photometry at the transient location is by using Scene Modelling Photometry (SMP; Holtzman et al. 2008). We extract SMP for a few selected objects of interest in Sect. 4 to test if the identified late-time detections are independent of our approach. When using SMP, one has to define an ‘off’ time and ‘on’ time. The observations taken during the off time are used to create a model of the region, or scene, where the SN occurs. This is then used as a template during the on time to calculate and remove host contributions to the photometry, leaving just the transient itself (Lacroix et al. in prep.). The advantage of this method is the significantly lower uncertainty in the model compared to the difference imaging technique, allowing us to find fainter detections. Since we assume that a signal from late-time interaction could occur at any point after the SN, we define the off time as everything up to shortly before the SN explodes, and the on time as everything after this moment.

2.2. Baseline correction

Issues in the construction of the reference images such inaccurate flat-fielding or artefacts in the images that are subsequently co-added can result in a systematic offset in the forced photometry light curve made using different images. The technique of baseline correction is used to correct for this (Yao et al. 2019; Miller et al. 2020).

To estimate the necessary baseline correction, we calculated the weighted mean of the flux of all data points up to 40 days before the estimated SN peak (which is assumed to be the highest flux detection deemed real), and subtracted it from the light curve. Baseline corrections are done separately for each combination of band (gri), field (telescope pointing), and rcid (part of the camera, which is arranged in 4 × 4 charge-coupled devices (CCDs) with four readout channels each, giving a total of 64 readout channels) as each of these combinations uses different, unique reference images. To be able to apply a baseline correction, at least two observations are needed. If this is not possible all observations with that band, field, rcid combination are removed.

Since we are interested in post-SN detections, our baseline correction method using only pre-SN data differs from the one used in Smith et al. (2025) with both pre- and post-SN data. A comparison between the methods found that for most objects our corrections agree with the ones used in Smith et al. (2025) within the uncertainties. This is as expected as objects with late-time flux excesses are expected to be rare, meaning that the two baseline correction methods should give the same result for most objects.

If there is insufficient data to perform a baseline correction, the relevant data (based on field, filter and rcid) is removed from the light curve. If this includes data around the peak position, the peak position in the light curve may change and therefore, the position of the peak is recalculated.

3. Analysis

To systematically search for objects with late-time flux excesses, a custom pipeline was developed. In Sect. 3.1, we describe how the late-time photometry for each object is binned to reach deeper magnitude limits, as well as the scheme used to select objects with robust, significant detections. In Sect. 3.2, we identify and remove bright nearby SNe Ia whose late-time detections are due to the SN radioactive decay tail. In Sect. 3.3, we describe our method of visually inspecting images of potentially interesting sources, and in Sect. 3.4, we detail our use of SIMSURVEY to simulate SNe Ia with late-time interaction signatures. Table 1 shows our sample size after each step of analysis that is discussed in the subsequent sections.

3.1. Binning and filtering programme

After pre-processing, the late-time observations are binned in phase to push the detection limit beyond the limit of the individual observations. We define late-time observations as being at least 100 days after the estimated date of SN peak brightness in the observer frame. We remove all SNe Ia that have no data in any band beyond this phase. The exact choice of 100 days is arbitrary but was chosen as a balance between minimising spurious detections due to the light curves being dominated by SN light at earlier times and maximising the phase range over which the interaction can be searched for.

Binning of the light curves is performed in each band separately. Larger bins are better for pushing the magnitude limit as deep as possible, but they sacrifice temporal sensitivity. To balance the time sensitivity and magnitude limit, we use bins with widths of 100, 75, 50, and 25 days. To make sure that the placement of the bin edges does not affect our results, we repeat the binning four times for each bin size while adjusting the starting epoch of the bins by decreasing the size of the first bin by 25% in each iteration. This results in a total of 16 trials for each band.

Each bin starts at the position of a data point to avoid empty bins. A gap in the data that is larger than the bins being used can cause the bins after the gap to always be placed in the same location despite the size modifications of the first bin. To avoid this we re-apply this modification of adjusting the size of the first bin after the data gap to trial different start positions. Lastly, if a bin would only contain one or two points, and the phases of these points occur no later than 10% of bin size of the previous bin (e.g. if the bin size is 25 d, they would occur within 2.5 d of the end of the bin), the previous bin is increased in size to include these points. An example of the bin placement is shown in Fig. 1.

Fig. 1. First 250 days of the gri-band light curves of SN 2019hbb in flux scaled to the peak flux (top panel) and magnitude (bottom panel) space. Binning starts 100 days after the estimated peak date (vertical black dashed and solid lines), using 25 d bins. The g-, r-, and i-bands are shown green, orange, and red, respectively. Before 100 days, we show the unbinned detections with their uncertainties (coloured circles) and non-detections (inverted triangles). After 100 days we show the bins as horizontal lines to show their size, a circle to show their mean value, and the shaded region showing the 1σ uncertainty (dashed regions are non-detections). A bin is deemed a non-detection if the flux f < 5σf. The 5σ magnitude limit is calculated and shown as a downward arrow. In both the g- and r-bands, the first bin is a detection and there are multiple adjacent bins with detections, triggering the tail-fitting procedure (see Sect. 3.2). The resulting tail fits are shown in the green and red lines, respectively, with their 1σ uncertainties as hashed regions. The half-life times are t1/2, g = 70 ± 6 d (χred2 = 0.6) and t1/2, r = 27 ± 4 d (χred2 = 1.4). This tail is therefore deemed to be a normal SN Ia tail.

For each bin the weighted mean and uncertainty of the observed flux are calculated. When binning flux measurements taken from difference images, the uncertainty of the reference images used in the difference imaging procedure has to be considered, as this will limit the depth of the binned observations (Strotjohann et al. 2021) and this is added in quadrature to the weighted uncertainty of the binned flux.

After the binning procedure each light curve has undergone 16 trials per band across four bin sizes and four bin placements. An attempt in a specific trial is considered significant if it has two or more adjacent bins with ≥5σ detections. A late-time detection is considered ‘robust’ if at least four out of 16 attempts have significant detections suggesting that the detections are insensitive to bin placement and/or size. We make this choice of robust detection in at least four attempts to ensure that we are not dominated by spurious detections but that we can still detect a long but faint signal that can only be picked up in the four trials involving the largest (100 day) bins.

3.2. Removal of SN Ia radioactive tail detections

For some nearby SNe Ia, the normal light curve tail powered by the radioactive decay of ⁵⁶Co → ⁵⁶Fe is still visible at the phases investigated here (> 100 d after the peak), possibly triggering false positives in our pipeline. To test if the fading tail is the reason for detections after 100 d, we checked if the detections follow a declining power law in flux space consistent with that of a normal SN Ia tail. We take all bins, normalise to the brightest data point, and fit a declining power law. To ensure that the tail matches the earlier data points, we include the unbinned observations between 60 and 100 d after the peak, making sure that if there are N bins only the latest N/2 unbinned detections are used to ensure that the fit focusses on the bins and not the unbinned points.

A successful fit of a declining normal SN Ia tail has a reduced chi-square of χ_red, fit² < 5 and fitted half-life of t_1/2 with uncertainty σ_t1/2, satisfying t_1/2 − 5σ_t1/2 ≤ 50 d. We chose a threshold of 50 d as Dimitriadis et al. (2017) showed that this is the approximate decay time scale for a normal SN Ia at these phases. A fit with a high χ_red, fit² value could have failed due to bad or uncertain data, or due to the late-time detections not following a power law decay. Fits with a t_1/2 significantly larger than that of a normal SN Ia tail suggest an additional luminosity source contributing to the light curve at these phases. Figure 1 shows an example where this tail fitting procedure determines the late-time detections in an object to be a normal declining SN Ia tail. 432 SNe Ia that were flagged as having late-time detection are discounted from further discussion because their light curves can be explained by a normal fading SN Ia tail.

3.3. SuperNova Animation Programme (SNAP)

After performing the binning and filtering, and removing SNe Ia with contamination from the radioactive tail, we are left with 134 SNe Ia with robust late-time detections (see Table 1). Since the binning and filtering programme is designed to handle a large quantity of light curves and cannot be tailored specifically to suit the peculiarities of a single object, it is possible that there are objects remaining with issues in the data or data processing (e.g. cosmic rays, bad subtractions), resulting in false positive detections. Therefore, we manually check the difference imaging to search for potential issues. To do this efficiently, we use SNAP¹.

Using ZTFQUERY (Rigault 2018), SNAP takes all difference images of the requested sky position during the requested time period(s) in the requested band(s) and shows them in chronological order in an animation. At the start of each animation the reference images in all bands are shown. The programme can show the image in grey-scale, a three-dimensional wire-frame representation of the intensities measured per pixel, the averaged values along both axes of the image, the observation date and duration, the peak and mean pixel values of the shown region, the last spectrum taken before the currently shown image, and highlight the resulting forced photometry point in the light curve corresponding to the plotted images.

Using SNAP, issues in the difference images can be identified, including SN ghosts (the SN is visible in the reference image, leaving a negative imprint in difference images after it has faded), cosmic rays, and bad pixels (NaN, or a large negative number). Variability of a separate source can also be seen, which, when close by, can contaminate the forced photometry at the SN location, for example, an active galactic nucleus (AGN).

3.4. Simulated interaction recovery fractions

To make sure the binning programme works as expected and estimate its detection efficiency in finding late-time signals, we simulated an observing campaign using SIMSURVEY (Feindt et al. 2019a,b), a python package designed to simulate large scale time domain surveys such as ZTF. To successfully simulate an observing campaign, the programme needs to be told what, when, where, and how something is observed, and under what conditions. For this, we need a model of the SN Ia-CSM that is going to be observed, an explosion rate as a function of redshift, and a time range for these explosions to occur. We also require an observing log specifying which part of the sky is being observed at a specific time, the length of the observations, and the weather conditions during the observations. Lastly, we require details of the camera that is used to carry out the observations. The SN model needs to be in a similar format to the SNCOSMO models (a Python package made for supernova cosmology, Barbary et al. 2021). This means we need a set of spectra over the entire phase range, all having the same wavelength spacing and range. Since no such model exists, we built our own model as described in the next section.

3.4.1. The interacting SN model

We chose SN 2011fe as the template of a normal SN Ia as it is well observed and close by (in M101 at a distance of 6.4 Mpc, Shappee & Stanek 2011). Optical spectra of SN 2011fe were obtained between phases of −18 to 1017 d relative to the peak. We made a custom model using SNCOSMO and spectra found on WISeREP (Yaron & Gal-Yam 2012)², which are listed in Table A.1. The spectra were flux calibrated to match the observed coeval broadband magnitudes. Spectra up to 45 d after the peak were flux calibrated using a SALT2 (Guy et al. 2007) fit of the PTF48g and PTF48R-band photometry (Law et al. 2009; Rau et al. 2009), which is used to estimate the flux in the g- and r- bands at these phases. Spectra between 45 and 400 d are flux calibrated using a cubic spline interpolation of photometry in the PTF48g- and PTF48R-bands. The interpolation extends up to 600 d after the peak in the PTF48R-band, there is no PTF48g photometry used between 400 and 600 d after the peak. Three interpolated photometry points from the SALT2 fits were used as anchor points to connect these two parts of the calibration. Dimitriadis et al. (2017) show that there is a slight kink in the light curve tail around 600 d after the peak. This is replicated in the model by calibrating all spectra more than 600 d after the peak using a cubic spline interpolation of photometry from the Large Binocular Telescope (LBT; Hill et al. 2006) in the Bessel R-band (Shappee et al. 2017).

After flux calibration, the spectra were dereddened to remove dust extinction effects, using the Cardelli et al. (1989) extinction law with A_V = 0.04 mag (Patat et al. 2013). The spectra were rebinned, and any wavelength region that was not covered in all spectra was removed as required by SIMSURVEY. Lastly, the model is corrected for distance, redshift and time dilation. The resulting model is that of a normal SN Ia, which exploded at a distance of 10 pc without any dust between the source and observer.

The best late-time detection of CSM in a normal/91T-like SN Ia was for SN 2015cp, where Hα emission was identified in its spectra at 664 d after light-curve peak (Graham et al. 2019b). To model potential CSM interaction signals similar to that of SN 2015cp, we add a narrow Hα line with a Gaussian profile to the SN 2011fe model. Due to the rareness of interaction in otherwise normal SNe Ia, we do not have good constraints on the diversity of interaction signatures and simulate a broad parameter space. The interaction was chosen to start at 100, 200, 300, or 500 days after the peak, last for 100, 300, or 500 days, and have a similar brightness to the observed signal in SN 2015cp (Graham et al. 2019b), as well as 10 times weaker or 10 times stronger than it. All possible combinations of these values are used, and a simulation without any interaction is also used as a control test, giving a total of 37 simulations.

Figure 2 shows an example of our model spectra at 300 d and SN 2015cp at 694 d in the rest frame (Graham et al. 2019b). It also shows the model redshifted to z = 0.07, where the Hα line is partly shifted into the i-band. Figure 3 shows the absolute magnitude ri-band light curves of the SN 2011fe model in the rest frame, as well as light curves with different strengths of Hα emission. Our interaction model will only generate a late-time interaction signal in the r-band (or i-band at z > 0.06) as we only add a Hα emission line. This is enough to test the binning programme but is likely too simple to reflect the actual late-time signal seen in SN 2015cp (which also showed O I and Ca II in the restframe i-band) or potential other events.

Fig. 2. Model spectrum at 300 days (SN 2011fe with the added Hα line) is shown in magenta overlaid on a rest-frame spectrum of SN 2015cp at 694 days in grey. The model flux has been scaled to the distance of SN 2015cp for comparison. The green, orange, and red shaded regions are the bandwidths of the g-, r-, and i-bands, respectively. The transmission profiles are plotted in the same colours for each band. The model is also shown shifted to z = 0.07 (and offset up in flux), where the Hα line has just started to be in the i-band.

Fig. 3. r (orange) and i (red) absolute magnitude light curves of the SN 2011fe model used in the simulations in the ZTF bands as a function of phase from rest-frame B-band peak (Mazzali et al. 2014). The bumpiness in the models is because the underlying SNCOSMO model class interpolates in flux space but fails to find an exponential decay. The added rest-frame CSM interaction model based on Hα emission (starting at a phase of 100 d) is shown with dotted lines for the r-band. Once the interaction becomes the dominant source, it smooths out the bumps from the underlying tail. The black line shows a radioactive decay with t1/2 = 50 d, typical of a declining normal SN Ia tail.

3.4.2. Simulating the observing campaign

By specifying the object to observe, as well as the telescope details and observation schedule, an observation campaign can be performed using SIMSURVEY resulting in a collection of observed light curves. For a deeper explanation of SIMSURVEY we refer the reader to Feindt et al. (2019b). The parameters used as input are listed in Appendix B. In each SIMSURVEY run, 10⁵ SNe Ia are simulated to produce observed light curves and meta-data such as redshift, observed peak date, etc. To ensure that the SNe Ia are similar around the peak to those recovered, we require that the SN Ia light curves must have at least three detections of ≥5σ and are brighter that 19 magnitude at the peak. This reduces the sample to ∼40 000 objects per simulation. These are sent through the binning and filtering programme (Sect. 3.1) as if they were real observed light curves to determine the recovery efficiency.

The volumetric rate used as input in the simulations favours more distant SNe, which results in very few SNe at extremely low redshift values and hence larger uncertainties. To mitigate this, we split 0 ≤ z ≤ 0.015 into bins of size 0.001 and simulate an additional 100 SNe in each bin using the same parameters as in the original simulations. Introducing these additional events does not impact the recovery efficiencies because we are comparing the number of recovered events relative to the input number in each redshift bin and therefore, are insensitive to the input rate of events.

3.4.3. Simulated interaction recovery

Our aim is to determine from our simulations how many SNe Ia with signatures of late-time interaction similar to that of SN 2015cp would have been detected by our pipeline. For each of the simulations, we binned the SNe based on their redshift and looked at the fraction of SNe that were reported by the pipeline to show late-time excesses. Figure 4 shows the recovery fractions as a function of redshift for an example simulation when the interaction starts at 500 d and lasts for 500 d for simulations of no CSM interaction, late-time interaction with the same strength as SN 2015cp, and interaction 10 times as strong as SN 2015cp (strong interaction). As expected the recovery fraction drops off with increasing redshift for both the SN 2015cp equivalent strength and the interaction that is 10 times stronger, with the strong interaction recoverable out to a higher redshift.

Fig. 4. Fraction of SNe Ia for one of our simulations (interaction occurring between 500 and 1000 d after the peak) where the CSM interaction was recovered per redshift bin of size 0.002. The simulations are shown for interaction strengths of zero (grey), similar to SN 2015cp (red), and 10 times stronger than SN 2015cp (blue). In the simulation without CSM interaction, the recovery fraction should be interpreted as the fraction of false positives. The simulations with normal ZTF quality reference images are shown with dots and fitted sigmoid functions with solid lines. Simulations where one magnitude deeper reference images were assumed are shown in triangles, with their fitted sigmoid functions in dashed lines.

The recovery fraction of the simulations with CSM interaction does not reach 100% in the lowest redshift bins. This is because the radioactive tails of these SNe Ia tend to be bright out to hundreds of days after explosion. Therefore, depending on the cadence and uncertainties of the simulated photometry, the SN light can dominate over the CSM interaction and the CSM interaction signal does not alter the shape of the SN decay tail enough to be flagged as CSM interaction.

In the simulation without CSM interaction, the recovery fraction is non-zero at small redshifts, meaning that some objects are falsely identified as having late-time excess. For these very bright and high signal-to-noise SN Ia light curves, our decaying tail model for normal SNe Ia proves to be too simple. Our analysis pipeline detects real deviations of the SN light curve evolution from our simple decay tail model. This only occurs at the lowest redshifts and nearby SNe Ia are rare, with only 0.6% of our observed SN Ia sample at z ≤ 0.01. This means that contamination of our sample due to normal SNe Ia tails that cannot be fit by our simplified tail fit model is very low.

We fitted a sigmoid function to the recovery fractions of each simulation, using the total amount of objects in each bin as its weight (Fig. 4). A sigmoid function is an oversimplification (the recovery fraction is underestimated at the low redshifts) but it allows us to easily estimate the redshift limit where CSM interaction can be recovered. We define our redshift limit where CSM interaction can be recovered as z₅₀, the redshift where 50% of the SN interactions are recovered. These values are listed for all simulations in Table A.2.

As discussed in Sect. 3.4.1, we simulated 36 models with interaction signatures starting at 100, 200, 300 and 500 d post peak, lasting for 100, 300 and 500 d, and with strengths the same as SN 2015cp, 10 times weaker and 10 times stronger. We also simulated a model without any late-time CSM interaction. For the models with an interaction strength similar to SN 2015cp (Graham et al. 2019b), when the interaction is short and early (starting 100 days after the peak and lasting for 100 days), the interaction cannot be distinguished from a normal SN Ia decaying light curve and the recovery fraction is as low as the no interaction simulation. When the interaction is longer, or if it starts later, the light curve flattens enough to be identified as deviating from a normally declining SN Ia tail. This pushes the redshift boundary where 50% of the interaction would be recovered by ZTF to z₅₀ = 0.0105 ± 0.0003 for the longest and latest interaction (500–1000 days after peak).

If the CSM interaction is 10 times weaker than that of SN 2015cp, the decaying SN Ia tail generally dominates over the interaction signature and the light curve shows little deviation from a normal decaying SN Ia tail. Even in the best case scenario of the longest and latest CSM interaction simulation, the 50% recovery thresholds lies at z₅₀ = 0.0050 ± 0.0005. In the simulations where the interaction is 10 times stronger compared than that of SN 2015cp, the shortest and earliest interaction (lasting from 100 to 200 days after peak) has z₅₀ = 0.0091 ± 0.0014. For the longest and latest interaction, the 50% recovery rate is at z₅₀ = 0.0323 ± 0.0004.

3.4.4. Impact of reference image depth

The mean limiting magnitude of the ZTF reference images is ∼21.8 mag and as discussed in Strotjohann et al. (2021), this is the limiting factor for recovering faint signals from binned light curve data. To test the improvement of deeper reference images, the assumed limiting magnitude was changed to be 0.5 and 1 mag deeper. The recovery fraction for one magnitude deeper is shown for comparison in Fig. 4. As expected, deeper reference images allows the interaction signatures to be detected to higher redshift, although the increases in z₅₀ values are modest (see Table A.2). For example, for the latest onset and longest interaction duration interaction, z₅₀ increases from 0.0323 ± 0.0004 to 0.0407 ± 0.0009.

4. Results

We run our custom detection pipeline on the ZTF DR2 light curves in the same way as it was performed on the simulated light curves in Sect. 3.4. In 1932 light curves, nothing is detected in any of the 16 trials discussed in Sect. 3.1, in 432 light curves the late-time detections are attributed to declining SN Ia tails, and in 1020 light curves the detections were not considered robust (< 4 successful trials). These are the three largest cuts in our sample, as can be seen in Table 1, and leave us with 134 light curves that pass the pipeline. In Sect. 4.1, we describe the light curves of these events and discuss how some light curves fit into known classes of events (e.g. known Ia-CSM, late-time SN Ia tail detections). In Sect. 4.2, we describe the additional tests that were performed on the remaining promising 10 events to determine if their late-time excesses are due to CSM interaction or other scenarios.

4.1. Initial summary of detected events

As can be seen in Table 1, the result of the pipeline is a list of 134 objects that require visual inspection after passing the detection cuts of positive > 5σ detections in adjacent light curve bins in at least four of the 16 bin size and placement combinations. In 47 of these cases, by visual inspection we identify that an incorrect baseline caused false positives. In five cases, the peak date estimation failed and estimated the peak to be over 100 days before the actual SN explosion. Because of this the SN itself was detected as a late-time signal. Furthermore, in 29 cases there was evidence of the host galaxy being improperly subtracted or showing signs of activity, which interfered with the forced photometry at the SN location. Finally, in 20 cases the tail fit test was unable to show the nature of the tails due to various reasons (e.g. the fits did not converge properly or there was a gap in the observations while the tail was visible). After this step, 33 objects were remaining in the sample.

Further details of these remaining 33 SNe Ia are shown in Table 2 and can be split up into four main groups: (i) known Ia-CSM events, (ii) transient siblings, where a second transient event occurs near the identified SN Ia causing its light to (partially) be picked up during forced photometry at the first SN location, (iii) nearby objects whose tail could not be fitted by our simple model, and (iv) objects that do not fall in the first three groups. In the following sections, we describe the first three of these groups that are not due to potential CSM interaction at late times.

4.1.1. Known Ia-CSM

The first group are the 13 known Ia-CSM, defined as those objects that already had a Ia-CSM classification. These objects started interacting relatively soon after the explosion and remained active long enough to be picked up by our pipeline beyond the 100 day threshold. Figure 5 shows the light curves of the recovered SNe Ia-CSM in absolute magnitude space (uncorrected for extinction). Even if the peak identified by our code is not the real peak due to it not being observed (e.g. for SN 2018evt, SN 2019agi, and SN 2019ibk), the CSM interaction persists for long enough for it to be picked up by our pipeline.

Fig. 5. Binned late-time observations of the recovered known SNe Ia-CSM. All objects are shown in absolute magnitude and over the same time range for easy comparison. All objects are detected beyond 100 days after the peak without using the binning technique. We do not show these individual data points to increase readability. The tail fits are shown as solid lines with the hashed region denoting their 1σ uncertainties. For comparison, SN 2020ssf (ZTF20abyptpc) in the top left corner is a normal SN Ia with a normally declining tail with t1/2, g = 53 ± 1 days and t1/2, r = 26 ± 1 days. The fitted tails for the SNe Ia-CSM are significantly shallower.

Ten of the known Ia-CSM SNe are presented in Sharma et al. (2023), who search for SNe Ia-CSM discovered in the ZTF Bright Transient Survey from May 2018 to May 2021 (BTS; Fremling et al. 2020; Perley et al. 2020). They find two objects that are not in our sample: SN 2020abfe (ZTF20acqikeh) and SN 2020aekp (ZTF21aaabwzx). SN 2020abfe is in the DR2 sample, but due to a combination of a gap in the observed light curve and the interaction not altering the declining tail sufficiently, our tail fit procedure is unable to distinguish it from a normal declining SN Ia tail. SN 2020aekp was first detected after the final date for objects to be included into our sample. Two of the events in our sample (SN 2020eyj and SN 2020kre) are not presented in Sharma et al. (2023). SN 2020eyj was excluded as Sharma et al. (2023) focussed on interaction with H-rich material and this object showed He emission lines suggesting interaction with He-rich material (Kool et al. 2023). SN 2020kre is not in the BTS sample and therefore, not included in Sharma et al. (2023). However, it was confirmed with spectroscopy to have Hα emission in its peak spectra.

Out of the 13 events SN 2020aeuh is an outlier, due to its distinct light curve. While the other 12 known Ia-CSM events detected in our sample have decline tails whose slopes are significantly shallower than for a normal SN Ia or steepen over time, SN 2020aeuh brightens significantly, having a double peaked nature with the second peak at around 100 d after the first. Even though the light curve suggests the SN to be interacting, no H emission (Hα or other lines) are observed. Kool et al. (in prep.) present a detailed analysis of this object.

4.1.2. Siblings

Siblings are transients that occur in the same host galaxy as each other and can be useful for understanding differences in local environments (e.g. Biswas et al. 2022; Graham et al. 2022). In some cases, the siblings occur in (almost) the same place on the sky, only differing in explosion time. This can be either due to the two transients being physically close together, or a projection effect due to the inclination of the host. However, the result is the same: forced photometry at the location of one sibling will result in a (partial) recovery of the other. Assuming that the first transient is a SN Ia in our sample and the second transient is fainter, our pipeline will flag the late-time rebrightening as a late-time excess in one of our objects.

Careful examination of the images using SNAP and cross-referencing using Fritz (an alert broker, van der Walt et al. 2019; Duev et al. 2019; Kasliwal et al. 2019; Coughlin et al. 2023) and the Transient Name Server³ (TNS) showed that there are five objects in our shortlist whose late-time detections are due to a sibling. Figure 6 shows the binned light curves of these objects in flux space. In each light curve there is a sudden significant spike in the detected fluxes in all observed bands, which falls back down again after a short period of time. Table 3 lists the name and type of each sibling if known, as well as their sky separation. In some cases the siblings are close enough together that they are not automatically recognised as separate events, resulting in them having the same name. In the case of SN 2019tzj the sibling (ZTF18aanhpii) exploded close to the nucleus, which had some spurious detections in 2018. This caused the sibling to have a 2018 ZTF name, although it exploded in 2022.

Fig. 6. Binned late-time observations in flux space of the five events with a detected sibling, with the flux normalised to the found peak flux. All objects are plotted on the same flux scale for easy comparison except for SN 2018big, as its late-time detections are much weaker compared to the original SN peak magnitude due to the larger distance offset between the siblings.

Figure 7 shows the detection of a sibling (SN 2019nvm) in the late-time light curve of SN 2018big in magnitude space. SN 2019nvm is slightly offset (∼4″) from the location of the original SN. The photometry pipeline forces the point spread function (PSF) fit at the position of SN 2018big. As the position of SN 2019nvm is slightly offset, only some of the total flux of SN 2019vnm is captured in the fit. Besides these five siblings, we identified two other pairs of siblings with SNAP: SN 2019gcm and SN 2021fnj, and SN 2020jgs and SN 2021och. These siblings were too far apart to be picked up with the forced photometry (9.6 and 10.8″, respectively) but were found while inspecting using SNAP. For a complete list and study on the siblings found in the ZTF DR2, we refer the reader to Dhawan et al. (2024).

Fig. 7. Light curves of SN 2018big and its sibling SN 2019nvm in magnitude space using bins of 25 d. The g (green) and r (orange) bins follow the tail of SN 2018big until it disappears in the noise. About 450 d after the peak of SN 2018big, new detections are identified in the binned photometry. The individual observations remain upper limits, although their shape hint to the true nature of these late-time detections.

4.1.3. Kinked tails

This group consists of five objects where the simple tail model, based on the typical decline rate of SNe Ia based on SN 2011fe (Dimitriadis et al. 2017), failed to fit the observations at later times. The reason for this failure in four of the events is that there is a slow-down in the r- and i-band decline rates at ∼200–250 d after peak, which the model does not take into account. The fifth event, SN 2020sck (Dutta et al. 2022), is a known SN Iax, a subclass known for having lower ejecta velocities and luminosities, suggesting that the explosion did not necessarily fully disrupt the star (Jordan 2012; Kromer et al. 2013). This event was flagged because of a slow-down in its decline rate roughly 80 days after the identified peak. The change in slope is visible in all bands and significantly longer than the assumed t_1//2 = 50 d of normal SNe Ia. The presence of a bound remnant has been suggested to be the cause of similar late-time signatures seen in other SNe Iax (Kawabata et al. 2018; McCully et al. 2022; Camacho-Neves et al. 2023).

The four SNe objects that deviate from the simple tail model are very nearby (z ≤ 0.009) compared to the majority of ZTF DR2 sample and when using the binned observations they are bright enough to be detected up to (nearly) a year after their first detection. As Rigault et al. (2025) are using these nearby events to calibrate their H₀ measurement, their peak magnitudes are currently blinded. Figure 8 shows the g-, r-, and i-band light curves of these objects in magnitude space. The R-band light curve of SN 2011fe (Zhang et al. 2016), which had a similar change in decline slope, is also shown for comparison. As discussed in Dimitriadis et al. (2017), the radioactive decays produce γ-rays and positrons, as well as X-rays and electrons that can be thermalised, depositing their energy in the expanding SN ejecta. As the ejecta expand over time they become more transparent, shortening the delay between thermalisation of the deposited energy and the emission of optical radiation. This results in the SN tail slope changing as the opacity changes.

Fig. 8. Five objects with kinks in their tails that start to deviate from the assumed decline rate at ∼200–250 d post peak are shown in magnitude space as a function of days since peak. As most of these objects have their peak magnitude blinded, no scaling is shown. A normal radioactive decay model was fitted to these tails, shown as solid straight lines with their 1σ uncertainty as dashed regions. But as the ejecta opacity changes over time so does the half-life time of the tail, causing a kink seen in the bins which is not reproduced by the model. The arbitrarily normalised R-band light curve of SN 2011fe (known not to have CSM interaction from detailed spectral studies) (Zhang et al. 2016) is shown in blue, showing the same shift in decline slope at a similar phase.

To be able to observe a change in decline slope such as this, a SN has to be both bright and well observed during the time it is visible. The other SNe Ia in our sample at similarly low redshifts have gaps in their observations or the change in slope is not strong enough for the tail fits to fall outside the allowed range of reduced chi-squared values (χ_red² > 5) and therefore, are not flagged by the pipeline due to this.

4.2. Additional tests of promising events

After performing the tests discussed in the previous sections on each event, there are ten objects remaining with an unexplained late-time excess. Light echoes produced by SN light scattering off of nearby dust clouds were considered, but were ruled out as these are typically ≥10 mag fainter than the SN at its peak (Patat 2005; Graur et al. 2016). We performed several additional tests on these ten events to try to determine the origin of their late-time excess. The first is using SMP as described in Sect. 2.1. The second is testing for coincidence with an AGN and third is more detailed comparisons with known transient classes. These tests are discussed below and a summary is shown in Table 4.

4.2.1. Scene modelling photometry

To test for issues in data processing, we ran the ten events through the scene modelling pipeline of Lacroix et al. (in prep.) as discussed in Sect. 2.1. Two SNe Ia (SN 2020awr and SN 2020kzd) were found to have issues with their difference imaging and forced photometry light curves, causing false detections. The detections in both these events are close to the detection limit and therefore, are impacted by even small errors in the baseline placement. The difference imaging (left panel) and scene-modelling (right panel) light curves of these two events are shown in Fig. 9.

Fig. 9. Two objects whose late-time detections were revealed to be caused by the photometry extraction. The colours are as in Fig. 8. The left side shows the binned forced photometry light curve, and the right side shows the binned SMP light curve. Bins with 5σ detections are shaded solid, while the non-detections are hashed.

For SN 2020awr, ∼300 d after the SN peak, the i-band observations jump up to detections at ∼20.3 mag. In contrast to other objects where baseline issues were found, here the offset occurs only for part of the light curve, while the pre-SN baseline has no visible issues. We could not identify a clear reason for this when inspecting the images with SNAP, and the SN is too far away from the host nucleus for it to be host activity. For the scene-modelling version of the photometry, the jump in the i-band observations has disappeared completely, showing that there was indeed an unidentified issue with the i-band data for this object. Since the late-time detections are determined to be spurious, this object is ruled out from having late-time detections.

In the case of SN 2020kzd, late-time detections are present in the g- and r-bands for hundreds of days (see Fig. 9). The SN is in a complex environment with three galaxies close to its sky position, which likely complicates the image subtraction and baseline correction. When SMP is performed on the event the detections disappear and average flux at late times is consistent with zero, showing that the binned forced photometry light curve likely suffered from a wrongly determined baseline correction.

4.2.2. AGN contamination

If a SN explosion site is coincident with a host galaxy that has an AGN, host activity is a likely cause of the late-time detections. Using data from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), Hviding et al. (2022) present a criterion to test if a galaxy hosts an AGN based on the WISE W1–W2 and W2–W3 colours, which we apply to our events. One object (SN 2019vzf) is 4.86 kpc of its host centre and its host is a known AGN, with WISE colours of W1–W2 = 0.55 ± 0.03 and W2–W3 = 2.90 ± 0.04 mag. In SNAP, the late-time signal appears to cover both the SN and AGN locations. The AGN contamination is too strong to put any meaningful constraints on the late-time flux at the SN location. We kept the object in our sample until now to test if it is possible to use scene modelling to reduce the AGN contamination. However, this is not possible. We attribute the late-time signal to host activity and disregard it in future discussion. The light curves of SN 2019vzf are shown in Fig. 10.

Fig. 10. Objects whose late-time detections are explained in the additional tests. The top four rows show the light curves of the objects where a previously undetected sibling transient as an explanation for the late-time observations could not be ruled out, with forced photometry light curves in the left-hand panels and the scene modelling photometry light curves in the right-hand panels. The 5σ detections are shown as bins with solid uncertainty regions and bins with hashed uncertainty regions are non-detections. The object whose late-time detections are caused by the host galaxy AGN, SN 2019vzf, is shown in the bottom row. The colours are as in Fig. 8 with g band in green, r band in orange, and i band in red.

4.2.3. Presence of a sibling close to the SN location

To test if a previously unidentified sibling transient is causing the late-time detections of the remaining seven objects, we compared their late-time light curves to known classes of transients, including a SN Ia, core-collapse SNe (Type Ib, Type Ic, Type IIP) and two tidal disruption events (TDE). Firstly, we have estimated the amount of potential host extinction from the main SN peak by assuming it was a normal SN Ia with a typical g-band peak of −18.8 to −19.3 mag after correcting for the distance to the SN and for Milky Way extinction. These estimated host extinction values are given in Table 4, assuming R_V = 3.1. The bright end of the absolute peak magnitude gives an upper limit for the host extinction, and the faint end gives a lower limit. After correcting for this range of additional host galaxy extinction, the late-time excesses have mean absolute r-band magnitudes of −15.4 to −17.5 mag (see Table 4).

After estimating the allowed extinction for each primary SN Ia, we initially assumed that if the late-time excess is due to another transient then it will have the same extinction along the line-of-sight. For these other transients, we used examples of a SN Ia, Ib, Ic, two IIPs, and two TDEs to compare against, with details of the comparison objects described in Table 5. The transients chosen to represent their category are all in the typical magnitude range for their type. Two TDEs and SN IIPs were chosen to represent the upper and lower end of the range of peak magnitudes expected for these transients. We use our model for SN 2011fe to represent SNe Ia and test the lower and upper edge of the range of normal SN Ia peak r-band magnitude.

These comparison objects were chosen as they have well-sampled light curves in the ZTF filters, and literature values for the host extinction. We correct the light curves of the comparison objects using their literature redshifts and their extinction values before correcting for the redshift and extinction of each SN Ia in our sample with a potential late-time excess. We then compare this transient light curve to the found late-time detections to see how well they match. A good match will have a similar magnitude, colour, and duration.

It could be the case that the suspected sibling was in the same line-of-sight direction, but had a different amount of extinction due to, for instance, exploding behind a cloud that adds additional extinction. We check this by adding enough extinction to match the r-band detections between the different comparison events and the observed late-time detection and again check if the colour and duration match up, as the observed colour is affected by the extinction. For the TDE comparisons, we allow a host galaxy E(B – V) of up to one magnitude, as was estimated for the ZTF TDE sample of Hammerstein et al. (2023).

For four events (SN 2019mse, SN 2019rqn, SN 2020alm, and SN 2020pkj), the late-time detections are consistent with at least one of the comparison classes, as detailed in Table 4 and shown in Fig. 10. We describe them individually in the following sections. The three remaining events (SN 2018grt, SN 2019ldf, and SN 2020tfc) cannot be explained by the presence of a sibling transient and are discussed further in Sect. 4.2.4.

SN 2019mse. SN 2019mse has late-time detections starting at 450 d after the peak and lasting ∼250 d in the gri-bands with an absolute r-band magnitude of −17.5 mag during the excess. Careful re-examination of the difference images show that the late detections are slightly offset (about one pixel) from the SN location, and appear to be on the host nucleus location instead. However, with its WISE colours being W1 – W2 = 0.25  ±  0.04 and W2 – W3 = 2.44  ±  0.09, the host is determined to not contain an AGN. The scene modelling version of the light curve shows a similar late-time excess, showing that this is not an artefact from the chosen photometric analysis method.

The late-time signal is detected in all three bands, and its behaviour is very similar in all of them (see Fig. 10). Its rise and decline time scales, absolute magnitude, and colours agrees well with the ranges seen for TDE. Together with the observation that the late-time detection are at the host nucleus location, this suggests that a nuclear transient explains the late-time signal adequately.

SN 2019rqn. In the case of SN 2019rqn, there is a short period of detections at 950–1050 d in the r-band after a gap in the observations, declining and fading below the detection significance within 100 d of the first detection. Nothing is detected at a ≥5σ level in the g or i-band observations. The i-band SMP light curve had its host contribution not completely subtracted, causing a flux offset in the data points. We therefore do not consider the SMP i-band further. The data points in the g- and r-band light curve have slightly larger uncertainties in the SMP version, causing the main SN light curve to fall below the 5σ threshold at an earlier epoch resulting in fewer individual points visible (Fig. 10). Similarly, larger uncertainties for the SMP prevents the detection of a late-time r-band signal.

Assuming a sibling exploded during the gap in the observations between 870 and 950 d post peak, SNe Ib and Ic with E(B – V) ≤0.3 mag extinction in the r-band can fit their tail to match the observed detections in the r-band without being excluded by the g-band and i-band non-detections. Therefore, we cannot rule out a sibling as the source of the detected late-time signal.

SN 2020alm. The late-time signal in SN 2020alm is seen in all three bands, beginning at ∼750 d after the peak and lasting for at least 300 d. There is a gap of 80 d in the observations immediately before the period of activity. The detections slowly rise to a plateau. Our initial analysis only included data for SN 2020alm up to ∼1000 d after the peak. However, when this object was identified as having late-time detections, the light curve pipeline was rerun and it was found to be still bright at later times, with significant detections in the r- and i-bands, but not above 5σ detections in the g-band. The i-band SMP light curve contained a significant flux offset, as the host was not fully subtracted. We therefore do not further consider this band. The binned SMP r-band light curve fails to reproduce the late-time detections found in the forced photometry light curve. However, the g-band detections are recovered in the SMP. The SN is close to the host nucleus at 0.66″ (0.85 kpc) offset at the redshift of the SN, but the WISE colours of W1–W2 = 0.06 ± 0.04 and W2–W3 = 2.42 ± 0.12 place it far outside the AGN region.

As SN 2020alm was still active while our analysis was on-going, we obtained two spectra using the Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy (OSIRIS) instrument on the Gran Telescopio CANARIAS (GTC) at Roque de los Muchachos in La Palma on 26 July 2023 using the R1000R grism. As the spectrum is heavily dominated by the host galaxy, we subtracted a rebinned spectrum of the host taken by the Sloan Digital Sky Survey (SDSS; York et al. 2000; Adelman-McCarthy et al. 2006; Gunn et al. 2006; Smee et al. 2013) in 2003, well before the SN occurred. We confirmed successful host subtraction by checking for residual NaID and Mg Iλ5175 absorption lines and found that no residual features were present. A more detailed explanation is given in Appendix C.

The resulting spectrum shows an excess that is stronger towards longer wavelengths (Fig. C.1). This is consistent with the broadband photometry finding an brighter excess in the redder bands, while the g-band remains within the noise after binning the observations. There is some excess in the narrow [N II] λλ6548, 6583, Hα and [S II] λλ6716, 6730 emission lines, but we lack the resolution to check if this is significant. There is no visible additional Hα component in the spectrum, which would be indicative of CSM interaction. Integrating the spectrum over the r- and i-band efficiencies gives an r − i colour of 0.6 mag, which is within 3σ of the value found in the latest photometry bin.

Hammerstein et al. (2023) show that TDEs generally have a g–r colour of zero and can have featureless spectra. We approximate a TDE by a flat line in order to estimate the amount of extinction needed to generate a red excess similar to the spectrum. We find that the general shape of the spectrum can be approximated with 0.6 < E(B − V)_host < 1 mag. This would mean an absolute r-band magnitude of −18.8 > M_r > −19.8 mag. Hammerstein et al. (2023) show that both this amount of host extinction and late-time brightness are possible for TDEs. They also show that it is possible for a TDE to rise and fall back down within the 80 d gap in observations, although this is seen in fainter TDEs than corresponding to our estimated absolute magnitude range.

The TDE sample of Hammerstein et al. (2023) did not contain a single object that matches our late-time detections in duration and luminosity in SN 2020alm. However, it is possible to combine parts of different TDEs together to make a TDE that peaked and decayed within the 80 d gap and levelled out by the time it became observable again. Based on this, a TDE is a plausible explanation for the late-time signal detected in this object.

SN 2020pkj. In the case of SN 20120pkj, the first r-band bin with a detection is the end of the normally declining tail, but after the first bin the detections rise slightly in the next r-band bin (Fig. 10). None of the binned photometry for the g- and i-bands give significant detections. Unfortunately, there is a gap in the observations immediately after the r-band detections preventing us from following its evolution closely at these phases. When it became observable again at > 200 d, no significant detections were found in any band. The duration of the transient is at least 75 d and could be up to 150 d. In the SMP light curve of this object, a similar rise at the same epochs is recovered, though it is found in the g-band instead of the r-band. This is likely due to the low significance of these detections at just 5.8-σ in the r-band in the forced photometry and 5.2-σ in the g-band in the SMP, showing that they are just on the detection limit of our binning technique.

With a redshift of z = 0.02456, this object is nearby enough for the end of the tail to be visible in the bins at 100 days. This could explain why the first bin is a detection, but the slight increase is still unexpected for a normal SN Ia decline tail. While the second bin is marginally consistent with our tail fit, a slight but steady increase can be seen in the unbinned flux values as well, suggesting that the brightening is real.

Assuming that the detected part of the late-time signal is the brightest part of a sibling transient, the tested transients need significantly more extinction of ∼3 to 6 mag in the r-band (depending on the comparison transient) than was found for the main SN Ia peak. However, the ∼80 d gap is long enough for a Type I sibling SN to peak higher during the gap and dim again before observations resumed. The other transient could have just started when it stopped being observable, and declined below the detection limit when the location became observable again. Since there is no constraint on the magnitude, this could work with any amount of extinction. Therefore, we cannot rule out conclusively that the late-time detections are due to another transient.

4.2.4. Late-time interaction candidates

Finally, we present the three objects (SN 2018grt, SN 2019ldf, and SN 2020tfc) whose late-time detections could not be explained by any of the explanations discussed above. Their light curves are shown in Fig. 11, and the colours of the late-time detections in Fig. 12 and the SMP light curves for SNe 2019ldf and 2020tfc in Fig. D.1. A scene modelling analysis could not be performed for SN 2018grt due to it exploding early in the ZTF survey. The forced photometry pipeline uses images obtained during the ZTF commissioning phase as templates but the SMP does not. We present these three objects as late-time interaction candidates as there are no spectra of these objects at these late times to confirm the presence of features consistent with CSM interaction.

Fig. 11. Three candidate objects, shown in magnitude and flux space. All three have significant detections (≥5σ) after a period of observations consistent with zero flux. From the alternative explanations the best fitting alternate transients are shown in dotted lines. For SN 2018grt this is the Type IIP SN 2017gmr, for SN 2019ldf and SN 2020tfc this is the TDE AT 2018hco.

Fig. 12. Colour curves of the three candidate objects, together with the colours of the best fitting alternate transients. The top and bottom panels show g–r and r–i, respectively. The first bin for each object starts at zero days, but the bins can be shifted horizontally in an attempt to better fit the colour curve of the transient compared against (given that this is allowed by the rest of the light curve). Bins whose mean observation dates are closest to each other are used to calculate the colour, provided that these bins overlap in time. If there is a detection in only one band used to calculate the colour while the other is a non-detection, the result is a lower or upper limit.

By using the estimated mean absolute r-band magnitude after removing all extinction effects (see Table 4) we can estimate the required Hα flux assuming that it is the source of all the detected flux in the r-band, and that the emission line has the same width as was used during the simulations. The identified strength of the signal is compared to the Hα signal detected in SN 2015cp for each of the events.

SN 2018grt. The late-time signal of SN 2018grt is only detected in the r-band, the g-band stays around zero flux, and there are no i-band observations at these phases. The first detection in the r-band begins at 1350 d post peak with a magnitude of 21.4  ±  0.2 (absolute magnitude of −16.4 mag). It varies little over the ∼100 d period where it is detected, after which it returns to zero flux within 50 d. The SN is close to the host nucleus with an offset of 0.35″ (0.32 kpc) but its host colours place it well outside the expected AGN region. Checking the difference images with SNAP shows that the host nucleus and SN location differ by ∼1 pixel.

The late-time detections are ∼2.5 mag below the main SN peak. If the late-time signal is due to another SN Ia or a TDE, it would require a significantly higher extinction value than was found for SN 2018grt itself. In addition, this object shows a sudden drop in the r-band, which is not normal behaviour for a TDE. Most of the other transient types cannot reproduce the plateau followed by the sharp decline only detectable in the r-band, apart from a SN IIP, where the plateau of SN 2017gmr has nearly the same time span. However, to fit the observed magnitude with a IIP SN, it would require a host E(B − V) that is three times higher than what was found for SN 2018grt, making this scenario less likely to be the case. Therefore, we conclude that late-time CSM interaction is a plausible scenario for the late-time signal in this object.

If we assume 0.21 ≤ E(B − V)_host ≤ 0.36 mag and that the r-band signal is produced only by an Hα emission line with a similar width to the one used in our simulations in Sect. 3.4, we estimate the strength of the emission to be much stronger than SN 2015cp, at 60 to 100 times its emission strength. However, there are examples of SNe Ia-CSM with interaction strengths this strong, for example, SN 2020eum was within this range (Sharma et al. 2023).

SN 2019ldf. SN 2019ldf has late-time detections in the r-band beginning at 1050 d after the peak and lasting for about 100 d, with an additional increase in brightness towards the end. There is a single 5σ detection in the i-band but there are a number of lower significance detections coeval with the r-band detections. These detections are directly after a long period without observations due to the object being behind the Sun. Nothing is detected in the g-band during the time of the rise in the r band. The binned SMP light curve recovers these late-time r-band detections, and a single i-band detection, showing that these detections are not specific to the photometry method.

We compared the properties of the late-time detections to those of our comparison transient objects. Even if we assume that a SN exploded during the gap in observations in order to avoid needing a significantly larger E(B − V)_host value, SNe evolve too much over a period of 100 days to explain the detections. In addition, detections would also be expected in the g-band, which have not been found.

A TDE could fit the detections if it was intrinsically bright but heavily extincted, as this could explain the red colour and absence of signal in the g-band. However, SN 2019ldf is offset from the host nucleus by 0.65″ (0.78 kpc) and inspection of the difference images using SNAP shows that the late-time signal is more consistent with the SN location than the host nucleus location, disfavouring the TDE explanation. Therefore, we conclude that the late-time signal could be due to late-time CSM interaction.

The late-time detections persist until the end of the observation window. To determine if there were still signs of interaction once it was visible again, we obtained g- and r-band photometry with the EFOSC2 imaging spectrograph (Buzzoni et al. 1984) on the ESO New Technology Telescope (NTT) in La Silla, Chile on 2023 May 19 as part of the extended Public ESO Spectroscopic Survey of Transient Objects+ (ePESSTO+; Smartt et al. 2015).

To examine whether the SN is still detected in our images from May 2023, we used image subtraction techniques. Due to the lack of reference images in the g- and r-band filters from EFOSC2, we used images from the DESI Legacy Imaging Surveys Data Release 9 (Dey et al. 2019). After aligning the images, we subtracted them from each other with the High Order Transform of Psf ANd Template Subtraction code version 5.11 (HOTPANTS; Becker 2015). We measured the brightness in the difference images using aperture photometry. The photometry was calibrated against stars from DESI Legacy Imaging Surveys. The EFOSC2 and DESI Legacy Imaging Surveys filters are not identical which might add an unknown systematic to the reported photometry. The 5σ upper limits are m_g = 24.7 and m_r = 24.3 mag at 1397 days after the peak, with no detection in either band. This means that the signal has disappeared at this time, and thus could have lasted, at most, for about 500 days.

Assuming that the r-band signal is entirely due to the Hα emission, we estimate it to be ∼60 times as strong as the late-time interaction found in SN 2015cp. However, this assumption is very simplistic, as it completely disregards the rise in the i-band and therefore it is only a first order estimate.

SN 2020tfc. This object has late-time gri-band detections, beginning at 550 d after the peak and lasting for at least 250 days. While the r- and i-bands are at more or less the same magnitude, the g-band detections are about 1.3 mag fainter. This immediately poses an issue for any alternate transient considered, as either there is a low amount of extinction and a weak g-band signal, or there is a high amount of extinction and the intrinsic signal is even brighter in the i-band. This, combined with the fact that the signal lasts for several hundreds of days with little variation, disfavours a SN as an alternate transient explanation. The intrinsic colour also heavily disfavours a TDE as these objects tend to have a similar intrinsic brightness in the g-, r-, and i-bands. In the binned SMP light curve the late-time detections were confirmed in the g-band. The SMP i-band data points have a large scatter, most likely due to uncertain background removal because of a low number of available images for the SMP template. Therefore, we do not consider them further in our analysis. In the r-band SMP light curve, the data points are below our 5σ cut-off for detections, although one is very close to our limit at 4.6σ. However, the g-band detections seen in the forced photometry are confirmed by SMP suggesting a real signal is present at late-times in at least this band, and at lower significance in the r band.

Similar to SN 2020alm, the late-time signal was on-going during our analysis but unfortunately there is no archival host galaxy spectrum available to compare to. As the host dominates the late-time signal (the SN is at a distance of 0.28″ (0.21 kpc) from the host nucleus) and cannot be removed using difference imaging or SMP as for the photometry, this prevented us from taking a spectrum of the late-time signal. With all alternate explanations ruled out or severely challenged by observations, the late-time CSM interaction remains as a plausible explanation for these late-time detections.

If we assume that the r-band signal is due to Hα emission only, the interaction is estimated to be 110 to 150 times as strong as the interaction found in SN 2015cp. This is by far the strongest of the three, but again this simple assumption is unrealistic as it completely ignores the measured g- and i-band signal that suggests a contribution from a continuum or other spectral lines to the late-time signal.

4.3. Late-time CSM interaction rates based on our candidate objects

We identified three objects (SNe 2018grt, 2019ldf, and 2020tfc) with potential late-time CSM interaction signatures. For SN 2018grt and SN 2019ldf, these detections were in the r-band (with no i-band data available for SN 2018grt and low significance i-band detections for SN 2019dlf). For SN 2020tfc, late-time detections were found in all three bands (gri). In our initial simulations to determine our recovery efficiency (Sect. 3.4), we made the assumption that any CSM interaction would be dominated by Hα emission that would be present in only the r-band up to z ∼ 0.07 and in the i-band beyond that. This was based on the dominant interaction signatures seen in SNe Ia-CSM and also the one event with late onset interaction, SN 2015cp (Graham et al. 2019b).

In Graham et al. (2019b), the Ca II NIR triplet emission was also identified in its spectra, along with emission consistent with Mg Iλ5175. Harris et al. (2018) speculated that although the Ca II NIR region was very noisy, the Ca II emission may be a similar strength to the Hα emission. For SN 2019ldf at a redshift of 0.057, the Ca II NIR triplet, if present, would be partially shifted out of the i-band, which could potentially explain its low-significance i-band detections. SN 2020tfc is at lower redshift (0.031) so the Ca II NIR triplet would fall completely in the i-band and could potentially be the source of the i-band detections. The strength of the potential Mg Iλ5175 line in SN 2015cp was weak and not well constrained (Graham et al. 2019b) but could potentially result in a weak signature in the g-band. However, without spectral confirmation this is only speculation.

The other difference between our initial simulations for determining recovery efficiency and the observed detections in our three candidate late-time CSM objects is that the late-time CSM interaction signature is much stronger in our observed events. In our three candidates, the late-time signal is between 2 to 3 magnitudes below the SN peak magnitude, while in the simulations with strongest interaction the signal was about 4.4 magnitudes below the peak.

It is clear from our candidate events, that there is significant diversity (detected bands, timescales and strengths) in their potential CSM signatures and without knowledge of their underlying spectra (e.g. emission lines that are present), it is difficult to develop spectral models and run simulations covering this full diversity. Therefore, although very simplistic and ignoring other potential emission lines occurring in the g- and i-bands as were seen in SN 2015cp (Graham et al. 2019b), we focussed on Hα emission dominated models and the corresponding r-band detection efficiency to estimate a rate of late-time Hα emission dominated CSM interaction based on our three candidate objects. We use the same underlying spectral model (SN 2011fe combined with Hα emission) as detailed Sect. 3.4 in further SIMSURVEY simulations but constrain the strength and timescale of the CSM signature based on our three events. In Table 6, we show the range of start epoch, duration and strength (dependent on assumed host galaxy extinction) simulated for each event. For each of the three objects, we take the worst case scenario (shortest, weakest, latest start epoch) and the best case scenario (longest, strongest, earliest start epoch) allowed by the constraints from the data to estimate a worst and best case scenario recovery efficiency for each object.

Since the available ZTF survey plan for our simulations only spans 1004 d, we cannot set the interaction signal to start at phases after this as the interaction would never be observed. Therefore, we choose to start the interaction in our simulations no more than 750 d after the peak. This is done to ensure there is a decent chance to observe an event until after the interaction has started, while still being able to apply a good baseline correction. As the interaction is by far the dominant source of light at these epochs, changing the start date of the interaction in this way has little effect on the total brightness of the SN. However, it does affect the number of objects observed by allowing more objects to have observations at later times, when the interaction occurs. This is a redshift independent effect, but will impact the resulting recovery fraction function. To account for this, we only consider objects with observations after the start of the interaction, and remove those where the simulation limits prevent any possible late-time detection. While changing the interaction start date has increased the number of objects that satisfy this condition, there are still some objects that do not due to e.g. being in a sparsely observed field.

Not all objects in the DR2 have detections up to the time after the peak where the interaction was detected in the three final candidates. Similar to the simulation, in order to get a proper estimate of the rate we can only consider the objects that could have been found interacting (N_sample), and remove those without observations this late after the peak. This gives us an effective sample size that is smaller than the full DR2, and is listed in Table 6.

The expected number of observed objects with late-time interaction in a sample N_recov(N_sample, η(z)) is a function of the efficiency and the intrinsic rate, which in turn is a fraction of the total SN Ia rate. If we assume the SN Ia rate from Frohmaier et al. (2019) as in our initial simulations and that the fraction of these showing late-time interaction, f, is redshift independent, then the only z dependency is in the efficiency of our pipeline. In a similar approach to Prajs et al. (2017) and Frohmaier et al. (2018), we run a Markov chain Monte Carlo (MCMC) simulation with 10⁷ realisations for each scenario to find the fraction of interacting SNe that best explains our findings.

In each realisation of the MCMC, we draw N_sample objects, and assign a redshift to each object according to the distribution from Frohmaier et al. (2019). For each object we also draw two random numbers from a flat distribution between zero and one. The first random number is used to decide if the object was interacting, and the second to decide if any interaction would have been recovered. Interaction is true if the random number is below the chosen value for f for that realisation. Recovery is true if the second random number is below η(z) at the z of the object (e.g. if the recovery efficiency is 0.7 at the SN redshift, a drawn number below 0.7 will result in a detection if the object was interacting). The number of recovered objects with interaction can be found by counting the objects for which both interaction and recovery is true, while the actual number of interacting objects can be found by counting the objects where interaction is true.

Figure 13 shows per bin the fraction of realisations that resulted in one object being recovered as a function of late-time CSM interaction fraction with the best and worst scenario parameters for the three discovered events. We approximate them by fitting a skewed normal distribution and use the fit to estimate the peak of the distribution and the 68% confidence interval on either side. These values are quoted as a fraction of the total SN Ia rate and as the late-time CSM interaction rate in the last two columns of Table 6.

Fig. 13. Fraction of MCMC realisations per bin that resulted in one object being recovered in a sample with the same size as the effective DR2 sample size as a function of the late-time CSM interaction fraction. The best and worst case scenario for each object is shown in blue and red, respectively. A skewed normal distribution fit with a 1σ uncertainty band is shown for each scenario, and the dashed and dotted lines give maximum and 68% confidence interval of these distributions, respectively. The distributions continue on the right side of each plot.

5. Discussion

In this section, we first discuss the overall rate of CSM interaction signatures found in our analysis (Sect. 5.1) and some of the properties of the three SNe Ia displaying these unexplained late-time detections (Sect. 5.2).

5.1. Late-time interaction is rare in SNe Ia

In our study of the ∼3500 SNe Ia in the ZTF DR2 with photometry at > 100 d, we identified three objects (SN 2018grt, SN 2019ldf, and SN 2020tfc) for which late-time CSM interaction is the best explanation for their detected late-time flux excesses. However, as each of these has significantly different parameters for the interaction, they have to be treated separately when trying to estimate the rate of late-time CSM interaction in SNe Ia (Sect. 4.3). Except for the worst case scenario of SN 2018grt, the identified rate for all cases agree with one another within the uncertainty, with the fraction of SNe Ia showing late-time Hα-dominated CSM interaction between $0.{03}_{-0.02}^{+0.08}$ and $0.{23}_{-0.11}^{+0.38}$%. This is similar to what was found by Sharma et al. (2023) in their Ia-CSM sample, and agrees with the upper limit set by Graham et al. (2019b) with their discovery of the late-time interaction in SN 2015cp. Some objects were disregarded as we could not definitively determine that they were not due to a sibling transient. Therefore, there could be a handful of additional objects displaying strong late-time Hα-dominated CSM interaction that were not included in our rate estimates but this would double the rate at most. Therefore, we estimate that the intrinsic rate of late-time Hα-dominated CSM interaction in SNe Ia is < 0.5%.

The identified late-time signals that are consistent with CSM interaction occurred between 1.5 to over 3.5 years after the original SN peaked. If we assume that the SN ejecta have a velocity of the order of 10⁴ km s⁻¹, the distance at which the CSM shell resides is of the order of 10¹⁷ cm. As CSM moves slowly, this means it must have been ejected from the system a long time ago. The further out the shell is from the progenitor system, the more mass is contained in even a thin shell. This is especially true when one considers that in order to get a detectable signal from the interaction, the CSM density cannot be too low. At these distances, light travel time effects also become significant, and a short interaction with a thin shell will be smeared out over a long time. Unfortunately, since we have only partial constraints on the interaction timescales and CSM mass estimation from an interaction signature is far from straightforward, we do not attempt to provide a CSM mass estimate.

An interesting aspect of the late-time detections in our three candidates is their inferred strength. While our model is likely an oversimplification by trying to explain the entire r-band signal using only a narrow Hα line, it is clear that the signal is much stronger compared to the one found in SN 2015cp, especially when there are other band(s) at a similar magnitude. Interaction this strong is not unprecedented though, SN 2020aekp for instance, starts to plateau around 50 days after peak at an absolute magnitude M ∼ −18.5 in all three ZTF bands (Sharma et al. 2023) and holds this plateau for several hundreds of days. Most of the known SNe Ia-CSM presented in Sharma et al. (2023) have detections in multiple bands, and the spectra presented show more than just an Hα line coming from the CSM interaction. CSM around a SN Ia also does not have to be spherically symmetric, especially in a scenario where material has been stripped and partially lost from the donor star in the progenitor system. Smith et al. (2015) show for SNe IIn that if the CSM mostly resides in a disk, only the ejecta travelling in the direction of the disk will be slowed down and interact while ejecta at higher inclinations continues to expand normally.

The colour curves for the late-time detections in our three objects (Fig. 12), show that g − r∼1.7 mag for SN 2020tfc, with limits of g − r > 0.5 and 1.5 mag for SN 2018grt and SN 19ldf, respectively. For the two objects with i-band observations, r − i < 0.5 mag. The SN Ia-CSM sample of Sharma et al. (2023) had late-time (> 300 d) g − r colours of 0 − 0.5 mag, which are similar to SN 2018grt but bluer than SN 2019ldf and SN 2020tfc. Kool et al. (2023) find g − r ∼ −1 mag and r − i ∼ 0.5 mag at late times for the He-interacting Ia-CSM SN 2020eyj. When taking into account the host galaxy extinction that they identify, it becomes even bluer and is much bluer than our three objects.

5.2. Properties of SNe Ia with late-time interaction

The three SNe Ia with late-time flux excesses that cannot be explained by other scenarios have SALT2 light curve fitter x₁ and c values that are generally typical of normal SNe Ia (Rigault et al. 2025). However, the c value of SN 2018grt of 0.61 ± 0.03 is at the high end and would be excluded from cosmological analyses. We also identified significantly more host extinction for this object than for the other two based on their luminosities at the peak.

All three candidates are found at small projected distances from their host nuclei. While the small sample size might be used to explain part of this observation, the mean projected distance of all objects in the ZTF DR2 is 6.3 kpc (Rigault et al. 2025). The chance for our three objects to be at most 0.8 kpc from the host is < 1%, suggesting that objects with late-time interaction have a preference for small host separations, the detections are caused by nuclear variability, or they are caused by bad host subtraction. The location of the late-time signal is too close to the host nucleus to distinguish which location the late-time signal is more consistent with by using SNAP. The host galaxy colours are inconsistent with being AGN and the properties of the late-time detections cannot be easily explained by known nuclear transient classes such as TDEs. It is possible that by going to deeper limiting magnitudes with our light curve binning that we have identified previously unstudied nuclear variability but this cannot be confirmed. Imperfect template subtraction can occur to galaxy centres, resulting in a dipole artefact. However, these are easily recognisable with SNAP, and this possibility was ruled out for our three candidates. Therefore, we conclude that that a likely interpretation is that there may be something intrinsic to SNe Ia exploding in these environments that produces late-time CSM interaction but further samples are required to confirm.

The morphologies of the host galaxies of the candidate events are broadly elliptical or spheroidal in nature with no signs of spiral arms. The galaxies have stellar masses of 4 × 10⁹, 6 × 10¹⁰, and 2 × 10¹⁰ M_⊙ for SN 2018grt, SN 2019ldf, and SN 2020tfc, respectively. These host masses are all within the bulk of the masses of the ZTF DR2 SN Ia sample (Rigault et al. 2025). The hosts of the known H-rich Ia-CSM sample have stellar masses of ∼3 × 10⁸ to 3 × 10¹⁰ M_⊙, consistent with our objects (Sharma et al. 2023). SN 2020eyj, the He-interacting and radio detected SN Ia-CSM (Kool et al. 2023) also had a very small offset from its host of 0.57 ± 0.02″ (0.36 ± 0.01 kpc). Its host only has a WISE W3 upper limit so it cannot be excluded from having an AGN present. It is a compact star forming galaxy with a mass of ∼6 × 10⁷ M_⊙, which is lower than the typical Ia-CSM range (Sharma et al. 2023).

The hosts of the three candidate events all have WISE W1 – W2 ≈ 0 mag. SN 2020tfc has W2 – W3 of ∼0.8 mag, while the other two have only limits with W2 – W3 > 1.2 mag. Based on Fig. 2 of Irani et al. (2022), this places them in the overlap region between elliptical and spiral host properties from the galaxy sample of Lintott et al. (2011). The vast majority of Ia-CSM sample of Sharma et al. (2023) has W2 – W3 colours of > 1 mag, again broadly consistent with our sample. However, the morphologies of our candidate events do not show evidence for spiral arms and their overlap between ellipticals and spirals in the W2 – W3 parameter space suggest they may come from different, older stellar populations than the known Ia-CSM events (Kool et al. 2023; Sharma et al. 2023). However, to prove that these are different populations requires a larger sample.

5.3. Limitations of the analysis

Since there is no good model of a Type Ia SN interacting with CSM at late times that could be used for our SIMSURVEY simulations, we had to make our own based on several assumptions. Our main assumptions were that the SN Ia looks normal until the moment the interaction starts, and we assumed that this interaction was with hydrogen-rich material showing itself primarily as an Hα emission line. This is motivated by the fact that Hα is found in most SNe Ia-CSM and the late-time interaction in SN 2015cp was confirmed through the observation of an Hα emission line. However, Graham et al. (2019b) also found other emission lines that are associated with the SN, such as the Ca II triplet near 8500 Å and a tentative detection of Mg Iλ5175, suggesting that our model is too simplistic. Similarly, for SN 2020tfc, we identified late-time detections in the gri-bands, which cannot be explained by Hα emission alone.

This method of binning the late-time light curves allows to push the detection limit down to a limiting magnitude m = 21.44 mag. Since the references used in ZTF have a mean limiting magnitude of ${\overline{m}}_{\lim }\sim 21.8$ mag, this becomes the leading uncertainty preventing our binning technique to go deeper. We have shown through our simulations that deeper references will allow for the binning technique to go deeper and recover fainter interaction signatures in a larger redshift volume. We have also used SMP to generate light curves that did not rely on forced photometry to test for issues coming from the photometry measurement. In most cases the late-time detections were recovered with both methods, showing their robustness. However, the SMP has been found to have issues with identifying a baseline flux when the number of images used in scene model template creation is small. This highlights the requirement for a long off time to model the underlying galaxy light sufficiently.

The pipeline focusses more on the historical light curve of the SNe than their current state. Since we require two adjacent bins to be detections in order to avoid false positives and the smallest bins we use are 25 d, the interaction needs to be active for over a month before it is picked up. Only after this, can it be followed up with spectroscopy if identified fast enough and the interaction signature is still occurring. Unfortunately, in the cases of our three candidates this could not be done due to a combination of recovering the interaction late, the objects being blocked by the Sun for a section of the year, and not having a suitable reference spectrum of the host galaxy to subtract and isolate the interaction spectrum. Adapting the pipeline to run in real-time would help with shortening the timescale to detection but if the host galaxy contribution is strong, then an archival reference spectrum is required.

5.4. By-product pipeline detections

The conservative, catch-all approach outlined in our work poses the challenge of identifying exactly why each object is flagged by the pipeline. Our approach allows for different kinds of objects that become interesting at late times to be caught. We recovered nearly all known Ia-CSM SNe through their non-normal SN Ia tails, as well as a number of known and unknown siblings, along with evidence of a change in the declining tail slope of close-by, normal Type Ia SNe around 200 days after the peak. In our sample, 4% of the objects needed to be visually inspected, upon which the cause of the (false) positive became clear quickly for most objects. Some of the checks, such as using the WISE colours to identify an AGN host close to the SN location, have been automated, reducing the workload in future attempts to use this method for searching for late-time signals, especially when using large datasets such as the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST; Ivezić et al. 2019).

The pipeline, as presented in this paper, is specifically tailored for finding SNe Ia that deviate from normal behaviour at > 100 d after the peak. However, the method of binning late-time observations can be used for any type of transient, but the checks to ensure that the resulting detections are not due to expected behaviour are SN Ia specific. When other types of transients (e.g. SN Ib, Ic, IIP) are used as input, the result will likely vary. For instance, a long lasting plateau in an SN IIP will be unexpected by the pipeline as it does not follow the decline tail of a normal SN Ia. However, modifying the pipeline to work for different kinds of transients mostly requires revision of these checks.

6. Conclusions

We have presented a search of the ZTF DR2 SN Ia light curve sample for signatures of late-time CSM interaction. This is the first systematic search for signatures of late-time interaction in SNe Ia in a large optical survey. We made a custom pipeline to calibrate and bin the light curve data at more than 100 d after the peak, using bins with sizes between 25 and 100 d. Our analysis was based on searching for Hα emission (as was seen in SN 2015cp and Ia-CSM) in the r-band at lower redshift and i-band above a redshift of z = 0.07. We performed simulations with SIMSURVEY to determine the efficiency of our search, as well as the intrinsic rate of potential Hα-dominated CSM interaction in the sample. Our main conclusions are:

1. Our pipeline returned 134 SNe that were potentially interesting based on their late-time light curves. Visual inspection of these objects was performed using our visualisation programme, SNAP, to inspect the difference images, removing 101 objects as false positives.

2. Of the remaining 33 objects, we identified 13 out of the 14 known Ia-CSM objects in the DR2, five siblings close to the position of the original SN Ia, four very nearby events whose late-time light curves were not captured by our simple radioactive tail model and one Iax with a late-time excess in all three bands, consistent with the presence of a bound remnant.

3. Out of our final shortlist of ten candidate events, we identified three SNe Ia (SN 2018grt, SN 2019ldf, and SN 2020tfc) that displayed late-time detections beginning 550–1350 d after peak and lasting at least 100–250 d, which could not be explained by data issues, AGN activity, or other transient events exploding at a similar location.

4. For SN 2018grt, these late-time detections were only in the r-band (no coeval i-band data was available). For SN 2019ldf, detections were made in the ri-bands and for SN 2020tfc in all three bands suggesting potential contributions with Ca II NIR emission or other Mg I as was identified in SN 2015cp (Graham et al. 2019b).

5. The r-band magnitudes of the late-time interaction are −16.5, −16.4, and −16.8 mag for SN 2018grt, SN 2019ldf, and SN 2020tfc, respectively. At their respective redshifts, this corresponds to Hα interaction strengths of 60–150 times that of SN 2015cp (depending on the extinction correction used). The strong nature of this signal could suggest we might only have found the high end of the late-time interaction strength distribution.

6. Using SIMSURVEY simulations of the ZTF survey, we estimated the intrinsic rate of strong Hα-dominated late-time (> 100 days after the SN peak) interaction to be occurring in < 0.5% of SNe Ia. This translates to absolute rates of $8_{-4}^{+20}$ to ${54}_{-26}^{+91}$ Gpc⁻³ yr⁻¹, assuming a constant SN Ia rate of 2.4 × 10⁻⁵ Mpc⁻³ yr⁻¹ for z ≤ 0.1.

The rarity of late-time interaction (occurring in < 0.5% of SNe Ia) highlights the importance of a large dataset of objects that have been observed for multiple years. The late-time detections occurred at different epochs for each object (from 550–1350 d post peak), showing that the phase at which SNe Ia will start to show signs of CSM interaction is highly variable. Therefore, the only viable strategy is to keep observing SNe even after they have faded beyond the detection limit and binning the late-time light curves. The interaction strength that we are sensitive to is heavily dependent on both the science and reference image depth. Future improvements to the analysis would include the use of deeper reference images to detect fainter signatures of late-time interaction, as well as running the pipeline in real time so that additional photometry and spectroscopy can be obtained to further characterise the late-time excesses. The deeper magnitude limits of LSST would be ideal for this study but cannot be immediately performed when the survey starts because of the requirements of deep reference images, as well as the need to wait up to more than three years after the SN Ia peak for the interaction to occur.

Data availability

The late-time spectrum of SN 2020alm 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/A11 The ZTF DR2 photometry is available on GitHub. The binning programme is available at https://github.com/JTerwel/late-timelcbinner. SNAP is available at https://github.com/JTerwel/SuperNova_Animation_Program. The ePESSTO+ photometry is available on the ESO archive.

Acknowledgments

JHT, KM and MD acknowledge support from EU H2020 ERC grant no. 758638. L.G. acknowledges financial support from the Spanish Ministerio de Ciencia e Innovación (MCIN), the Agencia Estatal de Investigación (AEI) 10.13039/501100011033, and the European Social Fund (ESF) “Investing in your future” under the 2019 Ramón y Cajal program RYC2019-027683-I and the PID2020-115253GA-I00 HOSTFLOWS project, from Centro Superior de Investigaciones Científicas (CSIC) under the PIE project 20215AT016, and the program Unidad de Excelencia María de Maeztu CEX2020-001058-M. This work was funded by ANID, Millennium Science Initiative, ICN12_009. T.E.M.B. acknowledges financial support from the Spanish Ministerio de Ciencia e Innovación (MCIN), the Agencia Estatal de Investigación (AEI) 10.13039/501100011033, and the European Union Next Generation EU/PRTR funds under the 2021 Juan de la Cierva program FJC2021-047124-I and the PID2020-115253GA-I00 HOSTFLOWS project, from Centro Superior de Investigaciones Científicas (CSIC) under the PIE project 20215AT016, and the program Unidad de Excelencia María de Maeztu CEX2020-001058-M. MN is supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 948381) and by UK Space Agency Grant No. ST/Y000692/1. Y.-L.K. has received funding from the Science and Technology Facilities Council [grant number ST/V000713/1]. TWC acknowledges the Yushan Young Fellow Program by the Ministry of Education, Taiwan for the financial support. 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 was supported by the GROWTH project (Kasliwal et al. 2019) funded by the National Science Foundation under Grant No 1545949. The Gordon and Betty Moore Foundation, through both the Data-Driven Investigator Program and a dedicated grant, provided critical funding for SkyPortal. The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID #2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. Pipeline processing and analyses of the data were supported by NOIRLab and the Lawrence Berkeley National Laboratory (LBNL). The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS Web Site is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington. Based on observations made with the Gran Telescopio Canarias (GTC), installed at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, on the island of La Palma. This work is (partly) based on data obtained with the instrument OSIRIS, built by a Consortium led by the Instituto de Astrofísica de Canarias in collaboration with the Instituto de Astronomía of the Universidad Autónoma de México. OSIRIS was funded by GRANTECAN and the National Plan of Astronomy and Astrophysics of the Spanish Government. The ztfquery code was funded by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement n°759194 – USNAC, PI: Rigault). This work made use of Astropy (http://www.astropy.org) a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022). Many people have contributed to lmfit. The attribution of credit in a project such as this is difficult to get perfect, and there are no doubt important contributions that are missing or under-represented here. Please consider this file as part of the code and documentation that may have bugs that need fixing. Some of the largest and most important contributions (in approximate order of size of the contribution to the existing code) are from: Matthew Newville wrote the original version and maintains the project. Renee Otten wrote the brute force method, implemented the basin-hopping and AMPGO global solvers, implemented uncertainty calculations for scalar minimizers and has greatly improved the code, testing, and documentation and overall project. Till Stensitzki wrote the improved estimates of confidence intervals, and contributed many tests, bug fixes, and documentation. A. R. J. Nelson added differential_evolution, emcee, and greatly improved the code, docstrings, and overall project. Antonino Ingargiola wrote much of the high level Model code and has provided many bug fixes and improvements. Daniel B. Allan wrote much of the original version of the high level Model code, and many improvements to the testing and documentation. Austen Fox fixed many of the built-in model functions and improved the testing and documentation of these. Michal Rawlik added plotting capabilities for Models. The method used for placing bounds on parameters was derived from the clear description in the MINUIT documentation, and adapted from J. J. Helmus’s Python implementation in leastsqbounds.py. E. O. Le Bigot wrote the uncertainties package, a version of which was used by lmfit for many years, and is now an external dependency. The original AMPGO code came from Andrea Gavana and was adopted for lmfit. The propagation of parameter uncertainties to uncertainties in a Model was adapted from the excellent description at https://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html#confidence-and-prediction-intervals, which references the original work of: J. Wolberg, Data Analysis Using the Method of Least Squares, 2006, Springer. Additional patches, bug fixes, and suggestions have come from Faustin Carter, Christoph Deil, Francois Boulogne, Thomas Caswell, Colin Brosseau, nmearl, Gustavo Pasquevich, Clemens Prescher, LiCode, Ben Gamari, Yoav Roam, Alexander Stark, Alexandre Beelen, Andrey Aristov, Nicholas Zobrist, Ethan Welty, Julius Zimmermann, Mark Dean, Arun Persaud, Ray Osborn, @lneuhaus, Marcel Stimberg, Yoshiera Huang, Leon Foks, Sebastian Weigand, Florian LB, Michael Hudson-Doyle, Ruben Verweij, @jedzill4, @spalato, Jens Hedegaard Nielsen, Martin Majli, Kristian Meyer, @azelcer, Ivan Usov, and many others. The lmfit code obviously depends on, and owes a very large debt to the code in scipy.optimize. Several discussions on the SciPy-user and lmfit mailing lists have also led to improvements in this code. Other software: numpy (Harris et al. 2020), pandas (The pandas development team 2020; McKinney et al. 2010), matplotlib (Hunter 2007), sncosmo (Barbary et al. 2021), simsurvey (Feindt et al. 2019a).

References

Appendix A: Tables

Appendix B: Inputs to SIMSURVEY simulations

The specific inputs to SIMSURVEY used in Section 3.4.2 to determine the detection efficiencies for SN 2015cp-like interaction are listed here.

• Model: SN 2011fe + Hα line (Section 3.4.1).

• Sky distribution: RA ∈ [0°, 360°], Dec. ≥ − 30° (The area covered by ZTF, Bellm et al. 2019).

• Volumetric rate: The SN Ia rate is 2.4 × 10⁻⁵ Mpc⁻³ yr⁻¹ for z ≤ 0.1 (Frohmaier et al. 2019). SIMSURVEY uses this to calculate the amount of SNe to generate at a given redshift interval.

• SN peak time distribution: 58 195≤ modified Julian date (MJD) ≤58 487 (between 18 March 2018 and 4 January 2019).

• Galactic extinction: dust maps from Schlegel et al. (1998).

• Host galaxy extinction: Cardelli et al. (1989) extinction law, with E(B – V) drawn from an exponential distribution with exponent λ = 0.11 (Stanishev et al. 2018), the same way as host extinction was added in the original SIMSURVEY paper (Feindt et al. 2019b).

• Telescope specifications: ZTF P48 camera, 4×4 grid of CCDs with four readout channels each, resulting in 64 separate output channels (Dekany et al. 2020).

• Survey plan: ZTF observation logs between 58 197≤ MJD ≤59 211 (between 20 March 2018 and 28 December 2020), ensuring all simulated SNe are followed for a minimum of about 2 years after the peak.

Appendix C: Late-time spectrum of SN 2020alm

After confirming the late-time detections in SN 2020alm to still be ongoing, we obtained a spectrum using OSIRIS on the GTC on 26 July 2023, 1277 days after the estimated peak date of the SN. As the observed spectrum is heavily dominated by the host galaxy, we subtracted a host galaxy spectrum taken by SDSS in 2003 to remove the host contamination. This was done after re-sampling the host spectrum to have the same wavelength spacing as the new spectrum. This left only the spectrum causing the late-time photometry detections, and is shown in Fig. C.1. The subtraction was confirmed to be successful by checking that the Mg Iλ5175 and NaID absorption lines were reduced to the noise level, as these lines are not expected to be due to the late-time signal but purely from the host galaxy. Some of the host galaxy emission lines were not completely subtracted during this process, most noticeably [N II] λ6583 and [S II] λλ6716, 6730, but our resolution is inadequate to draw any conclusions from this.

The only explanation for the late-time detections that uses a second transient at the same location is a TDE. Hammerstein et al. (2023) shows that the intrinsic spectrum of a TDE is flat in this wavelength range, with sometimes some narrow emission lines. Therefore, we model a TDE spectrum as a line of constant flux density and add Milky Way extinction (using the SFD89 dust maps in the direction of the object; Schlegel et al. 1998) and variable host extinction in an attempt to obtain the general shape of the observed spectrum. We find that 0.6 ≤ E(B − V)_host ≤ 1 mag is adequate to reproduce the general spectral shape and suggests that a TDE with approximately constant colour and moderate host extinction can explain the observed spectral excess for this event.

Fig. C.1. Spectrum of the late-time signal in SN 2020alm in its rest frame. The top panel shows the late-time spectrum obtained on 26 July 2023 using OSIRIS on the GTC, and the SDSS spectrum obtained in 2003. The bottom panel shows the late-time excess, obtained by subtracting the SDSS host galaxy spectrum from the observed late-time spectrum. A smoothed spectrum is shown in blue. The smoothing was done using a rolling kernel of size 5 to average over the values. The red lines are a simple TDE model with Milky Way and some amount of host galaxy extinction applied (the amount is shown in the legend), in order to get the approximate shape of the observed spectrum. Narrow emission and absorption lines that were notable in the unsubtracted spectrum are marked with dashed lines. The grey regions are affected by sky lines, and should be ignored. The late-time spectrum of SN 2020alm is available in electronic form at the CDS or upon request to the author.

Appendix D: Binned SMP light curves of the final candidates

Fig. D.1. Binned SMP light curves of the three candidate objects. As no bins ≥5σ are recovered in SN 2018grt, it is not shown here. Both SN 2019ldf and SN 2020tfc do however still have robust detections in some bands. The best fitting alternate transients, shown in dotted lines, are the same as in Fig. 11. The i-band of SN 2020tfc is not shown as the background was not subtracted completely, resulting in a significant flux offset.

All Tables

All Figures

Fig. 1. First 250 days of the gri-band light curves of SN 2019hbb in flux scaled to the peak flux (top panel) and magnitude (bottom panel) space. Binning starts 100 days after the estimated peak date (vertical black dashed and solid lines), using 25 d bins. The g-, r-, and i-bands are shown green, orange, and red, respectively. Before 100 days, we show the unbinned detections with their uncertainties (coloured circles) and non-detections (inverted triangles). After 100 days we show the bins as horizontal lines to show their size, a circle to show their mean value, and the shaded region showing the 1σ uncertainty (dashed regions are non-detections). A bin is deemed a non-detection if the flux f < 5σf. The 5σ magnitude limit is calculated and shown as a downward arrow. In both the g- and r-bands, the first bin is a detection and there are multiple adjacent bins with detections, triggering the tail-fitting procedure (see Sect. 3.2). The resulting tail fits are shown in the green and red lines, respectively, with their 1σ uncertainties as hashed regions. The half-life times are t1/2, g = 70 ± 6 d (χred2 = 0.6) and t1/2, r = 27 ± 4 d (χred2 = 1.4). This tail is therefore deemed to be a normal SN Ia tail.

In the text

Fig. 2. Model spectrum at 300 days (SN 2011fe with the added Hα line) is shown in magenta overlaid on a rest-frame spectrum of SN 2015cp at 694 days in grey. The model flux has been scaled to the distance of SN 2015cp for comparison. The green, orange, and red shaded regions are the bandwidths of the g-, r-, and i-bands, respectively. The transmission profiles are plotted in the same colours for each band. The model is also shown shifted to z = 0.07 (and offset up in flux), where the Hα line has just started to be in the i-band.

In the text

Fig. 3. r (orange) and i (red) absolute magnitude light curves of the SN 2011fe model used in the simulations in the ZTF bands as a function of phase from rest-frame B-band peak (Mazzali et al. 2014). The bumpiness in the models is because the underlying SNCOSMO model class interpolates in flux space but fails to find an exponential decay. The added rest-frame CSM interaction model based on Hα emission (starting at a phase of 100 d) is shown with dotted lines for the r-band. Once the interaction becomes the dominant source, it smooths out the bumps from the underlying tail. The black line shows a radioactive decay with t1/2 = 50 d, typical of a declining normal SN Ia tail.

In the text

Fig. 4. Fraction of SNe Ia for one of our simulations (interaction occurring between 500 and 1000 d after the peak) where the CSM interaction was recovered per redshift bin of size 0.002. The simulations are shown for interaction strengths of zero (grey), similar to SN 2015cp (red), and 10 times stronger than SN 2015cp (blue). In the simulation without CSM interaction, the recovery fraction should be interpreted as the fraction of false positives. The simulations with normal ZTF quality reference images are shown with dots and fitted sigmoid functions with solid lines. Simulations where one magnitude deeper reference images were assumed are shown in triangles, with their fitted sigmoid functions in dashed lines.

In the text

Fig. 5. Binned late-time observations of the recovered known SNe Ia-CSM. All objects are shown in absolute magnitude and over the same time range for easy comparison. All objects are detected beyond 100 days after the peak without using the binning technique. We do not show these individual data points to increase readability. The tail fits are shown as solid lines with the hashed region denoting their 1σ uncertainties. For comparison, SN 2020ssf (ZTF20abyptpc) in the top left corner is a normal SN Ia with a normally declining tail with t1/2, g = 53 ± 1 days and t1/2, r = 26 ± 1 days. The fitted tails for the SNe Ia-CSM are significantly shallower.

In the text

	Fig. 6. Binned late-time observations in flux space of the five events with a detected sibling, with the flux normalised to the found peak flux. All objects are plotted on the same flux scale for easy comparison except for SN 2018big, as its late-time detections are much weaker compared to the original SN peak magnitude due to the larger distance offset between the siblings.
In the text

Fig. 7. Light curves of SN 2018big and its sibling SN 2019nvm in magnitude space using bins of 25 d. The g (green) and r (orange) bins follow the tail of SN 2018big until it disappears in the noise. About 450 d after the peak of SN 2018big, new detections are identified in the binned photometry. The individual observations remain upper limits, although their shape hint to the true nature of these late-time detections.

In the text

Fig. 8. Five objects with kinks in their tails that start to deviate from the assumed decline rate at ∼200–250 d post peak are shown in magnitude space as a function of days since peak. As most of these objects have their peak magnitude blinded, no scaling is shown. A normal radioactive decay model was fitted to these tails, shown as solid straight lines with their 1σ uncertainty as dashed regions. But as the ejecta opacity changes over time so does the half-life time of the tail, causing a kink seen in the bins which is not reproduced by the model. The arbitrarily normalised R-band light curve of SN 2011fe (known not to have CSM interaction from detailed spectral studies) (Zhang et al. 2016) is shown in blue, showing the same shift in decline slope at a similar phase.

In the text

	Fig. 9. Two objects whose late-time detections were revealed to be caused by the photometry extraction. The colours are as in Fig. 8. The left side shows the binned forced photometry light curve, and the right side shows the binned SMP light curve. Bins with 5σ detections are shaded solid, while the non-detections are hashed.
In the text

Fig. 10. Objects whose late-time detections are explained in the additional tests. The top four rows show the light curves of the objects where a previously undetected sibling transient as an explanation for the late-time observations could not be ruled out, with forced photometry light curves in the left-hand panels and the scene modelling photometry light curves in the right-hand panels. The 5σ detections are shown as bins with solid uncertainty regions and bins with hashed uncertainty regions are non-detections. The object whose late-time detections are caused by the host galaxy AGN, SN 2019vzf, is shown in the bottom row. The colours are as in Fig. 8 with g band in green, r band in orange, and i band in red.

In the text

	Fig. 11. Three candidate objects, shown in magnitude and flux space. All three have significant detections (≥5σ) after a period of observations consistent with zero flux. From the alternative explanations the best fitting alternate transients are shown in dotted lines. For SN 2018grt this is the Type IIP SN 2017gmr, for SN 2019ldf and SN 2020tfc this is the TDE AT 2018hco.
In the text

Fig. 12. Colour curves of the three candidate objects, together with the colours of the best fitting alternate transients. The top and bottom panels show g–r and r–i, respectively. The first bin for each object starts at zero days, but the bins can be shifted horizontally in an attempt to better fit the colour curve of the transient compared against (given that this is allowed by the rest of the light curve). Bins whose mean observation dates are closest to each other are used to calculate the colour, provided that these bins overlap in time. If there is a detection in only one band used to calculate the colour while the other is a non-detection, the result is a lower or upper limit.

In the text

Fig. 13. Fraction of MCMC realisations per bin that resulted in one object being recovered in a sample with the same size as the effective DR2 sample size as a function of the late-time CSM interaction fraction. The best and worst case scenario for each object is shown in blue and red, respectively. A skewed normal distribution fit with a 1σ uncertainty band is shown for each scenario, and the dashed and dotted lines give maximum and 68% confidence interval of these distributions, respectively. The distributions continue on the right side of each plot.

In the text

Fig. C.1. Spectrum of the late-time signal in SN 2020alm in its rest frame. The top panel shows the late-time spectrum obtained on 26 July 2023 using OSIRIS on the GTC, and the SDSS spectrum obtained in 2003. The bottom panel shows the late-time excess, obtained by subtracting the SDSS host galaxy spectrum from the observed late-time spectrum. A smoothed spectrum is shown in blue. The smoothing was done using a rolling kernel of size 5 to average over the values. The red lines are a simple TDE model with Milky Way and some amount of host galaxy extinction applied (the amount is shown in the legend), in order to get the approximate shape of the observed spectrum. Narrow emission and absorption lines that were notable in the unsubtracted spectrum are marked with dashed lines. The grey regions are affected by sky lines, and should be ignored. The late-time spectrum of SN 2020alm is available in electronic form at the CDS or upon request to the author.

In the text

Fig. D.1. Binned SMP light curves of the three candidate objects. As no bins ≥5σ are recovered in SN 2018grt, it is not shown here. Both SN 2019ldf and SN 2020tfc do however still have robust detections in some bands. The best fitting alternate transients, shown in dotted lines, are the same as in Fig. 11. The i-band of SN 2020tfc is not shown as the background was not subtracted completely, resulting in a significant flux offset.

In the text