A Systematic Study of Ia-CSM Supernovae from the ZTF Bright Transient Survey

Among the supernovae (SNe) that show strong interaction with a circumstellar medium (CSM), there is a rare subclass of Type Ia supernovae, SNe Ia-CSM, which show strong narrow hydrogen emission lines much like SNe IIn but on top of a diluted Type Ia spectrum. The only previous systematic study of this class identified 16 SNe Ia-CSM, eight historic and eight from the Palomar Transient Factory (PTF). Now using the successor survey to PTF, the Zwicky Transient Facility (ZTF), we have classified 12 additional SNe Ia-CSM through the systematic Bright Transient Survey (BTS). Consistent with previous studies, we find these SNe to have slowly evolving optical light curves with peak absolute magnitudes between −19.1 and −21, spectra having weak Hβ and large Balmer decrements of ∼7. Out of the 10 SNe from our sample observed by NEOWISE, nine have 3σ detections, with some SNe showing a reduction in the red wing of Hα, indicative of newly formed dust. We do not find our SN Ia-CSM sample to have a significantly different distribution of equivalent widths of He i λ5876 than SNe IIn as observed in Silverman et al. The hosts tend to be late-type galaxies with recent star formation. We derive a rate estimate of 29−21+27 Gpc−3 yr−1 for SNe Ia-CSM, which is ∼0.02%–0.2% of the SN Ia rate. We also identify six ambiguous SNe IIn/Ia-CSM in the BTS sample and including them gives an upper limit rate of 0.07%–0.8%. This work nearly doubles the sample of well-studied Ia-CSM objects in Silverman et al., increasing the total number to 28.


INTRODUCTION
When it comes to supernovae (SNe) interacting with circumstellar material (CSM), a number of sub-types of corecollapse SNe (CCSNe) show signs of strong interaction, like SNe IIn (Schlegel 1990;Filippenko 1997), SNe Ibn (Pastorello et al. 2008;Foley et al. 2007;Chugai 2009;Hosseinzadeh et al. 2017) and most recently SNe Icn (Gal-Yam et al. 2021, 2022;Perley et al. 2022).SN IIn progenitors are generally thought to be massive stars (like Luminous Blue Variables, LBVs) that lose their hydrogen envelopes to winddriven mass loss and outbursts (Gal-Yam et al. 2007;Gal-Yam & Leonard 2009;Kiewe et al. 2012;Taddia et al. 2013;Smith 2014).Helium-rich but hydrogen-deficient CSM in the case of SNe Ibn (Pastorello et al. 2008;Foley et al. 2007;Chugai 2009) and both hydrogen and helium deficient CSM in SNe Icn (Gal-Yam et al. 2022;Perley et al. 2022;Pellegrino et al. 2022) are thought to arise from high-velocity wind mass loss or stripping of the envelope in binary configurations of massive Wolf-Rayet (WR) like stars.For SNe IIn in most cases, the mass-loss rate derived from the CSM velocity is consistent with estimates from LBV-like eruptive mass loss.
However, there exists a rare sub-type of thermonuclear supernovae (SNe Ia) which also interacts strongly with CSM i.e.SNe Ia-CSM.This class poses a challenge to the progenitor debate of SNe Ia.There is some consensus on there being at least two major progenitor channels for SNe Ia; the doubledegenerate (DD) channel (Webbink 1984;Iben & Tutukov 1984) which is the merging of two C/O white dwarfs and the single-degenerate (SD) channel (Whelan & Iben 1973) where the white dwarf accretes enough material from a nondegenerate companion to explode.Although there are more arguments for the DD scenario from observations of nearby SNe Ia (Nugent et al. 2011;Li et al. 2011;Brown et al. 2012;Bloom et al. 2011), the strongest observational evidence for the SD scenario are SNe Ia with CSM.
Indications of CSM around SNe Ia ranges from detection of time varying narrow Na ID absorption lines (Patat et al. 2007;Blondin et al. 2009;Simon et al. 2009) in highresolution spectra (found in at least 20% of SNe Ia in spiral hosts, Sternberg et al. 2011;Maguire et al. 2013;Clark et al. 2021), to strong intermediate and narrow Balmer emission features in the spectra and large deviations of the light curves from the standard shape.The latter phenomena have been named SNe Ia-CSM (Silverman et al. 2013), but were earlier referred to as "SNe IIna" or "SNe Ian" due to the strong similarity between their spectra and those of SNe IIn.The first two examples of this class studied in detail were SNe 2002ic (Hamuy et al. 2003;Deng et al. 2004;Wang et al. 2004;Wood-Vasey et al. 2004;Kotak & Meikle 2005;Chugai et al. 2004) and 2005gj (Aldering et al. 2006;Prieto et al. 2007), but for a long time there was ambiguity regarding their thermonuclear nature (Benetti et al. 2006).These SNe were dominated by interaction from the first spectrum and were quite over-luminous compared to normal SNe Ia.The first clear example of a thermonuclear SN Ia-CSM was PTF11kx (Dilday et al. 2012;Silverman et al. 2013).It looked like a luminous SN Ia (99aa-like) at early phases but started showing interaction at ∼ 60 days from explosion and thereafter strongly resembled SNe 2002ic and 2005gj at late times.Higher resolution spectra taken at early times indicated multiple shells of CSM with some evacuated regions in between.Dilday et al. (2012) suggested a symbiotic nova progenitor involving a WD and a red giant (similar to RS Ophiuchi) could produce such CSM distribution, however later studies argued that the massive CSM of PTF11kx was inconsistent with the mass-loss rates from symbiotic nova systems (Silverman et al. 2013;Soker et al. 2013).
Ever since, a handful of SNe of this class have been studied in detail to investigate their progenitors and to distinguish them from their spectroscopic cousins, the Type IIn SNe.Both SN Ia-CSM and SN IIn spectra share a blue quasicontinuum, a strong Hα feature with an intermediate and a narrow component, and often a broad Ca NIR triplet feature, but they differ with regards to the line strength of Hβ, strength/presence of helium and presence of emission lines from intermediate mass elements often found in CCSNe.There are some individual SNe with unclear type often referred to as SN Ia-CSM/IIn, like SN 2012ca for which some papers argue for core-collapse (Inserra et al. 2014(Inserra et al. , 2016) ) and others for a thermonuclear origin (Fox et al. 2015).This ambiguity becomes more dominant as the underlying SN flux gets smaller compared to the interaction power (Leloudas et al. 2015).Silverman et al. (2013, hereafter S13) is the only study to analyze a sample of SNe Ia-CSM, 16 objects in total including 6 previously known, 3 re-discovered (reclassified SNe IIn) and 7 new from the Palomar Transient Factory (PTF).Their paper presents the common properties of optical light curves, spectra and host galaxies and contrast them against SN IIn properties.In this paper, we present 12 new SNe Ia-CSM discovered as part of the Zwicky Transient Facility's (ZTF; Bellm et al. 2019;Graham et al. 2019;Dekany et al. 2020) Bright Transient Survey (BTS;Fremling et al. 2020;Perley et al. 2020) and analyze their optical light curves, spectra, hosts and rates.Throughout this paper, we have compared the results derived from our sample to the ones in S13.
This paper is organised as follows; we first discuss the sample selection criteria, the photometric and spectroscopic data collection in §2, then the analysis of light-and color-curves and the bolometric luminosities is done in §3.1.The analysis of early and late-time spectra and emission line identification is presented in §3.2, and analysis of the host galaxies is provided in §3.3.The rates are estimated from the BTS survey in §3.4.We end with a discussion about the nature of SN Ia-CSM progenitors and a summary in §4 and §5.

OBSERVATIONS AND DATA REDUCTION
In this section, we outline our selection criteria, and present the optical photometry and spectroscopic observations of the 12 SNe Ia-CSM in our sample.

Selection Criteria
To carefully curate our sample of SNe Ia-CSM, we used the BTS sample and its publicly available BTS Sample Explorer1 website to obtain the list of all classified Type Ia subtypes during the period 2018-05-01 to 2021-05-01.We then filter out oddly behaving Type Ia SNe based on their lightcurve properties.We used two criteria; the primary being rest-frame duration considering flux above 20% of peak flux, and the second being change in magnitude after 30 days from peak (∆m 30 ).We calculated these two properties from either g or r-band light curves (whichever had maximum number of detections) grouped into 3-day bins and used Gaussian Process Regression2 to interpolate the light curves where coverage was missing.For the first filtering, we calculated the mean (µ ≈ 35 days) and standard deviation (σ ≈ 16 days) of the duration distribution and selected everything that had a duration greater than µ + 3σ.Given the large sample size (N = 3486), the standard error on the mean is ∼ 0.5 days, hence our duration cut of 3σ is suitable.This filtering selected 41 out of 3486 BTS SNe Ia.Then from these 41 SNe, we calculated the mean and standard deviation of the ∆m 30 distribution and removed SNe that were more than 1σ away from the mean on the higher side to reject the relatively steeply declining long SNe, which resulted in 35 SNe being kept.Again, the mean and standard deviation of ∆m 30 distribution of these 41 long duration SNe are 0.48 mag and 0.27 mag respectively and the standard error on mean is ∼ 0.04, making our 1σ cut suitable.Finally, we manually inspected the 35 selected SNe Ia to confirm their classification.20 out of the 35 SNe that passed the above filtering criteria were just normal SNe Ia either caught late or missing some postpeak coverage in ZTF or had spurious detections that resulted in long duration estimates, 2 had incorrect duration estimate due to an interpolation error and were recalculated and 1 (AT2020caa; Soraisam et al. 2021) had some detections before the SN explosion which could be connected to a different SN (i.e. a sibling; Graham et al. 2022).
The remaining 12 long-duration SNe Ia all turned out to be spectroscopically classified SNe Ia-CSM in BTS, and none of the classified BTS SNe Ia-CSM were missed in this filtering.No other SNe apart from these stood out in particular, indicating the classification reliability of the BTS sample.During the same period, 9 SNe Ia-CSM were reported to the Transient Name Server (TNS), out of which 7 are already in our sample, 1 was detected by ZTF but did not meet the BTS criteria, and 1 was not detected in ZTF as the transient location fell too close to the field edges and was masked by the automated image subtraction pipeline.Yao et al. (2019) presented early photometric observations of one SN Ia-CSM in our sample, SN 2018crl.Table 1 summarizes the coordinates, redshifts, peak absolute magnitudes, durations, host galaxy information and Milky Way extinction for the 12 SNe Ia-CSM in our sample.
Furthermore, we re-checked the classifications of 142 SNe IIn classified in BTS during the same period as above, in case any SN Ia-CSM was masquerading among them and found 6 to have ambiguous classifications.These are discussed further in Appendix A.

Discovery
All SNe Ia-CSM were detected by the ZTF (Bellm et al. 2019;Graham et al. 2019;Dekany et al. 2020) and passed the criteria for the BTS (Fremling et al. 2020;Perley et al. 2020) automatic filtering, i.e. extra-galactic real transients with peak magnitudes brighter than 19 mag.These were saved and classified as part of BTS which aims to classify all transients brighter than 18.5 magnitude, and reported to the Transient Name Server3 (TNS) during the period 2018-05-01 to 2021-05-01.Out of the 12 SNe, 6 were first reported to TNS (i.e.discovered) by ZTF (AMPEL, Nordin et al. 2019;Soumagnac & Ofek 2018 and BTS), 3 were first reported by GaiaAlerts (Hodgkin et al. 2021), 2 by ATLAS (Smith et al. 2020) and 1 by ASAS-SN (Shappee et al. 2014).For classification, 9 were classified by the ZTF group, 1 by ePESSTO (Smartt et al. 2015;Stein et al. 2018a), 1 by SCAT (Tucker et al. 2018;Payne et al. 2019) and 1 by the Trinity College Dublin (TCD) group (Prentice et al. 2020).The follow-up spectral series for these SNe were obtained as part of the BTS classification campaign as many were difficult to classify with the ultra-low resolution spectrograph P60/SEDM (Blagorodnova et al. 2018) and hence were followed up with intermediate resolution spectrographs.The SEDM spectra were helpful in determining an initial redshift but the template matches were unclear (matched to SN IIn as well as SN Ia-CSM and SN Ia-pec templates, some matched poorly to SN Ia/Ic at early times).SNe 2019agi (classification and spectrum taken from TNS), 2019rvb, 2020onv, 2020qxz and 2020uem were classified as Ia-CSM ∼ 1 − 2 month after discovery using spectra at phases of 42, 26, 38, 45 and 51 days respectively.SNe 2018crl, 2018gkx and 2019ibk were classified ∼ 2 − 3 months after discovery using spectra at phases of 92, 75 and 103 days respectively.SNe 2018evt, 2020abfe and 2020aekp were classified ∼ 4 − 5 months after discovery using the spectra at phases of 144, 146 and 132 days respectively.SN 2020xtg immediately went behind the sun after its first detection in ZTF therefore its first spectrum (using SEDM) was taken at 91 days since explosion which was dominated by strong Hα emission, and thus SN 2020xtg was initially classified as a Type II.As this SN was exhibiting a long lasting light curve, an intermediate resolution spectrum was taken at 340 days which matched very well to SN Ia-CSM and therefore its classification was updated.SNe 2020uem and 2020aekp showed peculiar features and were followed up for more optical spectroscopy for single object studies (to be presented in future papers).

Optical photometry
To assemble our sample light curves, we obtained forced PSF photometry via the ZTF forced-photometry service (Masci et al. 2019;IRSA 2022) in g, r and i bands and also added data from ATLAS (Tonry et al. 2018;Smith et al. 2020) forced-photometry service in c and o bands.The high cadence ZTF partnership survey in i band contributed some photometry to SNe 2018crl, 2018gkx, 2019agi, 2019ibk and 2019rvb.The ZTF and ATLAS data were supplemented with data from the Rainbow camera (RC, Ben-Ami et al. 2012) on the robotic Palomar 60-inch telescope (P60, Cenko et al. 2006) and the Optical wide field camera (IO:O) on the Liverpool telescope (LT, Steele et al. 2004).The P60 data was processed with the automatic image subtraction pipeline FPipe (Fremling et al. 2016) using reference images from SDSS when available, and otherwise from Pan-STARRS1.The IO:O data was initially reduced with their standard pipeline4 then image subtraction was carried out using the method outlined in Taggart (2020).For SN 2018evt, some early time data available from ASAS-SN (Shappee et al. 2014;Kochanek et al. 2017) in the V band was obtained through their Sky Patrol5 interface.
We corrected all photometry for Milky Way extinction with the Python package extinction (Barbary 2016) using the dust extinction function from Fitzpatrick (1999), the Schlafly & Finkbeiner (2011) dust map, and an R V of 3.1.Then we converted all measurements into flux units for analysis and considered anything less than a 3σ detection an upper limit.There is moderate to good coverage in g, r, c and o bands for all SNe in our sample.Figure 1 shows a multipaneled figure of the light curves of the objects in our sample.The transients were observed during the ongoing NEO-WISE all-sky mid-IR survey in the W 1 (3.4 µm) and W 2 (4.5 µm) bands (Wright et al. 2010a;Mainzer et al. 2014).We retrieved time-resolved coadded images of the field created as part of the unWISE project (Lang 2014a;Meisner et al. 2018).To remove contamination from the host galaxies, we used a custom code (De et al. 2020) based on the ZOGY algorithm (Zackay et al. 2016) to perform image subtraction on the NEOWISE images using the full-depth coadds of the WISE and NEOWISE mission (obtained during 2010-2014) as reference images.Photometric measurements were obtained by performing forced PSF photometry at the transient position on the subtracted WISE images until the epoch of the last NEOWISE data release (data acquired until December 2021).Further analysis of the mid-IR photometry is presented in §3.1.4

Optical spectroscopy
The main instruments used for taking spectra and the software used to reduce the data are summarized in Table 2. Additionally, the spectrum Reguitti (2020) obtained using the Asiago Faint Object Spectrograph and Camera (AFOSC) on the 1.8 m telescope at Cima Ekar, and the spectrum Stein et al. (2018b) obtained using the ESO Faint Object Spectrograph and Camera version 2 (EFOSC2) on ESO New Technology Telescope (NTT) were taken from TNS.
The details for all optical spectra (61 for the sample in total) presented in this paper are provided in Table 3.Furthermore, all spectra were corrected for Milky Way extinction using extinction and the same procedure as for the photometry.The SN redshifts were derived using narrow host lines for the objects which did not already have a host redshift available in the NASA/IPAC Extragalactic Database6 (NED).Photometric calibration was done for all spectra i.e. they were scaled such that the synthetic photometry from the spectrum matched the contemporaneous host-subtracted ZTF r-band data.For SN 2018crl, a host galaxy spectrum taken using P200/DBSP was available, which was subtracted from the P200/DBSP SN spectrum taken at +92 days.For SN 2020aekp, more spectra beyond ∼ 350 days were obtained but will be presented in a future study of the object (34 additional spectra up to ∼600 day).
These processed spectra were used for the rest of the analysis as detailed in §3.2 and will be available on WISeREP7 (Yaron & Gal-Yam 2012).For the purpose of this paper, the 'explosion time' simply refers to the time when optical flux rises above the zero-point baseline (i.e.first light).We used pre-peak g, r, i-band ZTF photometry and c, o-band ATLAS photometry (binned in 1day bins), when available, for our analysis.For each SN, the light curve was interpolated using Gaussian process regression to obtain the peak flux epoch, then a power-law (PL) model was fit using epochs from baseline to 60% of peak brightness in each band following Miller et al. (2020).The PL fits converged in at least one band for 6 out of 12 BTS SNe Ia-CSM.For the rest, we simply took the middle point between the first 5σ detection and the last upper limit before this detection as the explosion epoch with half of the separation between these two points as the uncertainty.
The explosion time estimates, light curve bands used for the PL fits and the 1σ uncertainties on explosion times are listed in Table 4.The unfilled 'PL fit filters' column in the table are the SNe for which the PL fit did not converge and averages were used.For the PL fits this typically constrains the time of explosion to within a fraction of a day.Given the high cadence of the ZTF survey, even in the cases where we use only the last non-detection the uncertainty range is typically less than 3 days.Only for SN 2020uem is the date of explosion virtually unconstrained (±57 days) as it was behind the sun at the time of explosion.
Although for SN 2019ibk the explosion time is formally constrained with a ±3 day uncertainty, this estimate was derived using only ATLAS o-band data right after the SN emerges from behind the sun.There is not a clear rise observed over a few epochs but two non-detections before a 5σ detection.It is possible that the actual peak of this SN occurred earlier while it was behind the sun and the rising o-band points after it emerged are due to a second peak or bump (similar to SN 2018evt, in that case the actual rise was caught before the SN went behind the sun in ASAS-SN data).If the former explosion epoch estimate from o-band is to be believed then SN 2019ibk would be the most sub-luminous among the SNe Ia-CSM, peaking at −17.5.

Duration and absolute magnitudes
Figure 2 shows the SNe Ia-CSM (colored squares) in our sample in the duration-luminosity and duration-∆m 30 phase space.In the top panel, the x-axis is duration above half-max and the y-axis is the peak absolute magnitude (see Table 1) when we have photometric coverage both pre-peak and postpeak.For SNe missing the pre-peak coverage, their discovery magnitude is taken to be the upper limit to peak absolute magnitude and the duration from discovery the lower limit Table 4. Explosion time epoch estimates derived from pre-peak multi-band light curves.For 6 out of 12 SNe Ia-CSM, we were able to fit a power-law model to multi-band data following Miller et al. (2020).For the remaining 6 SNe, the explosion epoch was estimated by taking the mean of the first 5σ detection and last upper-limit before the first detection.to duration above half-max (marked by arrows in Figure 2).The BTS SN Ia sample is shown in gray points, and we also show the SNe Ia-CSM presented in S13 with empty triangles for comparison in the top panel.In the bottom panel, the x-axis is duration above 20% of peak flux (∆t 20 ) and the y-axis is ∆m 30 , the two parameters used in the selection criteria.Most of the SNe Ia-CSM lie on the longer duration and brighter luminosity side, and are even more distinctly separated in the ∆t 20 -∆m 30 phase space.This makes the SN initial decline rate and duration useful tools for identifying thermonuclear SNe potentially interacting with CSM, if they have not revealed themselves already in their early time spectra.The gray points lying in the same phase space as SNe Ia-CSM are the false positive cases described in §2.1.Also worth noting is that the duration calculated by taking the flux above half of peak flux value does not capture the true duration of the light curve when the plateau phase falls below half-max as is the case for SN 2020aekp (> 500 days light curve) but ∆t 20 and ∆m 30 do.

Light and color curves
We have good pre-peak coverage in ZTF data for 8 of the 12 SNe in our sample 8 .SN 2018evt was discovered by ASAS-SN on JD 2458341.91 (Nicholls & Dong 2018) and classified by ePESSTO the next day (Stein et al. 2018a), around 115 days before the first detection in ZTF when the SN came back from behind the sun.Hence we have only one Change in magnitude 30 days after peak (∆m30) vs. rest-frame duration above 20% of peak-flux for BTS SNe Ia and SNe Ia-CSM.These criteria were used to filter out potential SNe Ia-CSM from all SNe Ia and demonstrate that SNe Ia-CSM occupy a distinct portion in this phase space.However some gray points (not SN Ia-CSM) remain on the longer duration side and are the false positive cases described in §2.1.
epoch of pre-peak photometry and one early spectrum for SN 2018evt.
Our mixed bag of SNe Ia-CSM show post-maximum decline rates ranging from 0.5 to 2.0 mag 100d −1 in the r band from peak to ∼ 100 days post peak.The median decline rate is 1.07 mag 100d −1 , which is much slower than the decline rates of normal SNe Ia.We see a variety of changes in decline rates after around 100 days from peak.Two SNe (2020onv and 2020abfe) show no change and have a constant slow decline throughout.Four SNe (2018gkx, 2019agi, 2019ibk and 2019rvb) evolve to a shallower slope going from ∼ 0.6-1 mag 100d −1 to ∼ 0.2-0.5 mag 100d −1 .Three SNe (2018crl, 2020qxz and 2020aekp) show a major change in decline rate with the light curves becoming almost flat, and SN 2020aekp shifts back to a slow decline from this plateau after ∼ 200 days.In three cases, the decline rate actually becomes steeper, SN 2018evt goes from 0.52 mag 100d −1 to 1.4 mag 100d −1 , SN 2020uem goes from 0.52 mag 100d −1 to 1.25 mag 100d −1 and SN 2020xtg seems to go from 0.61 mag 100d −1 to 1.35 mag 100d −1 (even though there is only one epoch at late times to measure this change).The 3 SNe with fastest initial decline rates ( 1.5 mag 100d −1 in the r band) are similar to SN 2002ic (initial decline of 1.66 mag 100d −1 in V ) and PTF11kx (initial decline of 3.3 mag 100d −1 in R) and coincidentally are also the ones that evolve into a plateau.The rest of the sample have initial decline rates comparable to SN 1997cy (0.75 mag 100d −1 ) and SN 2005gj (0.88 mag 100d −1 ) (Inserra et al. 2016).From these observations, we can conclude that SNe Ia-CSM exhibit a range of slow evolution indicating that there exists a continuum of phases at which strong CSM interaction begins to dominate the powering of the light curves for these SNe.It is, however, difficult to pinpoint the exact phase when interaction starts from the light curve without modeling.CSM interaction could be affecting the peak brightness significantly even in cases where interaction only appears to dominate after a few weeks (SNe 2018crl, 2020qxz 2020aekp).Considering the average peak phase to be ∼ 20 days past explosion from the light curves and assuming an ejecta velocity of ∼ 20000 km s −1 , the CSM is located at ∼ 3.5 × 10 15 cm.This estimate can be refined by considering the phase of the earliest spectrum that shows interaction signatures (see §3.2).At late times, all the decline rates are slower than that expected from Cobalt decay (0.98 mag 100d −1 ), confirming that the power from CSM interaction dominates the light curve behaviour for a long time.
Figure 3 shows the g − r color evolution of our sample SNe as a function of phase (rest-frame days from rband maximum), comparing them with some famous SNe Ia-CSM (SNe 2005gj, 1997cy, 1999E), and SNe 2012ca (Ia-CSM/IIn), 2010jl (IIn) and 1991T (over-luminous Type Ia).The color evolution of normal SNe Ia from ZTF (Dhawan et al. 2022) is shown in grey lines.We use g − r colors when available, otherwise we estimate the g − r color by fitting Planck functions to estimate the black-body temperatures from the V − R colors.Our SNe Ia-CSM show similar color evolution as the older Type Ia-CSM/IIn interacting SNe, i.e. the g − r color increases gradually for about 100 days and then settles onto a plateau or slowly declines, and one object (SN 2019ibk) becomes redder at late times similar to SN 2012ca.The interacting SNe are redder at late times compared to the normal SNe Ia.

Mid-IR brightness comparison
Out of 12 SNe in our sample, only one observed (SN 2020abfe) did not have 3σ detections post explosion in the unWISE difference photometry light curves and two (SNe 2019rvb and 2020qxz) did not have coverage post explosion.The unWISE light curves for the rest of the SNe Ia-CSM having > 3σ detections in W1 (3.3 µm) and W2 (4.6 µm) bands are shown in Figure 4 (black and red stars) along with Spitzer IRAC survey data of SN 2008cg (indigo and magenta empty triangles), SN 2008J (indigo and magenta empty squares) (both Ia-CSM) and some SNe IIn (blue and orange crosses) taken from Fox et al. (2011).The most nearby SN in our sample, SN 2018evt, is among the brightest (∼ 17 AB mag) in MIR at least until ∼1000 days after explosion and has a bumpy light curve.SNe 2019ibk and 2018crl however are the most luminous with an absolute magnitude of −18.7 mag in the W1 band.The brightness of the BTS SNe Ia-CSM is comparable with other interacting SNe and span a similar range (−16 to −19).However, SNe IIn have been detected until even later epochs (up to 1600 days) than SNe Ia-CSM, probably due to the larger number of SNe IIn at closer distances.SN 2020abfe has upper limits around ∼ −18 in W1 band and ∼ −18.5 in W2 band up to ∼300 days post explosion shown with upside down filled triangles.As the mid-IR luminosity can be fainter than these limits for SNe Ia-CSM (as can be seen for other nearby SNe in this sample) and SN 2020abfe is at a redshift of 0.093, it might just be out of reach for WISE.
This brightness of SNe Ia-CSM in mid-IR can be indicative of existing or newly formed dust.A clear signature of new dust is reduced flux in the red wing of the Hα emission line at late phases as the new dust formed in the cold dense shell behind forward shock absorbs the far-side (redshifted) intermediate and narrow line emission (see bottom panel of Fig. 7).For our sample, this reduction in Hα red wing is the most pronounced for SN 2018evt.

Bolometric luminosity
As the SN Ia-CSM luminosity is dominated by CSM interaction, their spectra comprise of a pseudo-continuum on the blue side and strong Hα emission on the red side, hence a blackbody fit to multi-band photometric data is not appropriate to estimate the bolometric luminosity.Instead we calculate a pseudo-bolometric luminosity from the available multiband optical data by linearly interpolating the flux between the bands and integrating over the optical wavelength range spanned by the ATLAS and ZTF bands.The individual band light curves are first interpolated using Gaussian process regression to fill in the missing epochs.This estimate places a strict lower limit on the bolometric luminosity.
All BTS SNe Ia-CSM show a slow evolution in bolometric luminosity, inconsistent with the decay of 56 Co to 56 Fe.The sample's overall luminosity decline rates are comparable to those of SNe 1997cy and 2013dn, as shown in Figure 5.Only SNe 2018crl and 2020aekp seem to show early decline in their pseudo-bolometric light curves similar to SN 1991T for about 40 days after peak like SN 2002ic and PTF11kx.Another BTS interacting SN Ia, ZTF20aatxryt (Kool et al. 2022), was found to follow the PTF11kx light-curve evolution very closely and as its light curve fell into a plateau the SN started showing signs of interaction with a heliumrich CSM and evolved into a helium-rich SN Ia-CSM.We have excluded ZTF20aatxryt from the sample as we focus on typical SNe Ia-CSM interacting with hydrogen-rich CSM in this study.At late phases (∼ 300 days), the SNe Ia-CSM are approximately 100 times brighter than normal SNe Ia at the same epoch.Therefore, at these late phases, the luminosity and spectral features of SNe Ia-CSM are entirely dominated by CSM-interaction with little emergent SN flux.From the pseudo-bolometric light curves, we place a lower limit on the total radiated energy for SNe Ia-CSM to be 0.1-1.5 ×10 50 erg.This is well below the thermonuclear budget (E kin ∼ 10 51 erg), but as this is a lower limit and some SNe in the sample have unconstrained peaks, the true total radia-tive energy might come close to the thermonuclear budget, requiring high conversion efficiency to achieve their luminosity.1991T, 1997cy, 1999E, 2002ic, 2005gj, 2013dn, and PTF11kx from literature.The light curves in each filter having more than 10 epochs were interpolated using Gaussian process regression to fill in the missing epochs, and at each epoch the fluxes between the bands were linearly interpolated and integrated over the optical wavelength range spanned by ZTF and ATLAS filters to get the pseudobolometric luminosity.For BTS SNe, the phases are with respect to the estimated explosion epochs, while for comparison SNe the phases are with respect to discovery.

Spectroscopy
Figure 6 displays the spectral series obtained for the BTS SNe Ia-CSM.Most of the early time spectra were taken with the SEDM, the BTS workhorse instrument (R ∼100), which is not able to resolve the narrow CSM lines.Therefore, these SNe were followed up with higher resolution instruments to get more secure classifications.For each spectrum in Figure 6, the phase is provided with respect to the explosion epoch estimate given in Table 4.We have spectra ranging from a few to around 470 days from explosion.Considering the well constrained explosion times of SN 2018evt, presence of narrow Hα in its first spectrum at 8 days since explosion and assuming a typical ejecta velocity of ∼20000 km s −1 , this implies that the CSM interaction start as close as ∼1.4×10 15 cm.
Figure 7 shows the early time (left) and late time (right) spectral behaviour of the BTS SNe Ia-CSM together with a few historical SNe for comparison, namely SNe Ia-CSM SN 2011jb (Silverman et al. 2013)  and [Fe II/III] line regions, and vertical dashed lines mark the Balmer emission lines.The sample spectra have been multiplied by a constant factor to magnify relevant spectral features.In the following paragraphs, we compare the observations of some of the spectral features with previous analysis of this class (Silverman et al. 2013;Fox et al. 2015;Inserra et al. 2016).A few of our early time SNe Ia-CSM show underlying SN Ia absorption features like PTF11kx and SN 2002ic (most are, however, quite diluted and also affected by the low resolution and signal-to-noise ratio (SNR) of the SEDM spectra), the most notable being SNe 2018evt, 2020qxz and 2020aekp.SNe 2020qxz and 2020aekp also have among the fastest initial post-peak decline rates in the sample, similar to PTF11kx, while coverage around peak is not available for SN 2018evt.On the other hand, SNe with slower decline rates similar to SN 1997cy and SN 2005gj have more SN IIn-like early time spectra dominated by blue pseudo-continuum and Balmer emission.The faster decline rate suggests we are still seeing some of the emission from the ejecta at those phases.To unveil the nature of the progenitor of interacting SNe, it is therefore necessary to obtain some spectroscopic follow-up before peak light.Spectroscopic data at the phase of transition to interaction-dominated luminosity would also help in deducing the extent and density structure of the optically thick CSM.
Late time spectra of SNe Ia-CSM look very similar to those of SNe IIn, heavily dominated by Hα emission.The CSM interaction masks the underlying SN signature and we instead see late-time spectra riddled with photoionized CSM lines.In some cases, the photosphere might lie in an optically thick cold dense shell (CDS) formed between the forward and reverse shocks which obscures the ejecta completely (Smith et al. 2008;Chugai et al. 2004).The continuum is also enshrouded under a blue quasi-continuum from a forest of irongroup element lines (S13) as identified and analyzed for SNe 2012ca and 2013dn by Fox et al. (2015).
The blue quasi-continuum blend of iron lines ([Fe III] lines around ∼4700 Å and [Fe II] around ∼5200 Å) in the spectra of the BTS SN Ia-CSM sample (see Figure 7 top right panel) is the dominant feature blue-ward of 5500 Å but the ratio of [Fe III]/[Fe II] is much weaker compared to for SNe Ia (like SN 1991T).This feature is more apparent in the SNe Ia-CSM like PTF11kx and SN 2002ic that became interactiondominated later than for other SNe Ia-CSM such as SNe 1997cy, 1999E and SN 2012ca (SN Ia-CSM/IIn, for which a clear type has not been established).Inserra et al. (2014) argues for a core-collapse origin for SN 2012ca given this low amount of [Fe III] along with the detection of blueshifted Carbon and Oxygen lines (which however, were later argued to be [Fe II] lines by Fox et al. 2015).S13 instead argues in favor of a thermonuclear origin given the presence of this blue quasi-continuum, despite [Fe III] being weaker.Fox et al. (2015) points out that a similarly suppressed ratio of [Fe III]/[Fe II] is observed in some SNe Ia, particularly the super-Chandra candidate SN 2009dc, for which the explanation was suggested to be a low ionization nebular phase owing to high central ejecta density and low expansion velocities (Taubenberger et al. 2013).Fox et al. (2015) argue that in the case of SNe Ia-CSM, a lower ionization state could arise owing to the deceleration of ejecta by the dense CSM explaining the Fe line ratio suppression.Since Ca has lower first and second ionization potentials than Fe, the detection of [Ca II] λλ7291, 7324 would be consistent with this low ionization, which Fox et al. (2015) confirms for SNe 2012ca and 2013dn.Indeed, we find clear evidence of [Ca II] emission for 8 out of 12 SNe in our sample and moderate to weak signal for the remaining 4.Although this does favor the argument for a thermonuclear origin, a similar blue quasicontinuum is also observed in other interacting SN types like SNe Ibn (SN 2006jc, Foley et al. 2007) and SNe IIn (SNe 2005ip and 2009ip), making Fe an incomplete indicator of the progenitor nature (see detailed discussion in Fox et al. 2015).
We do not find strong evidence of O I λ7774 or [O I] λλ6300, 6364 emission in our sample, although they might be present at very weak levels in some SNe (e.g.SN 2020uem).SN 2020uem has strong emission lines at 6248, 7155 and 7720 Å which are consistent with being iron lines and were also observed in SNe 2012ca, 2013dn and 2008J.S13 note that the very broad emission around 7400 Å can be due to a blend of [Ca II] λλ7291, 7324 and [O II] λλ7319, 7330, however we note that this broad emission is likely to be from calcium as O II is harder to excite than O I which is either very weak or absent in our spectra.The broad Ca NIR triplet feature resulting from electron scattering is the next strongest feature after the Balmer emission and is present in all mid to late-time spectra of the SNe in our sample where the wavelength coverage is available.We observe it increasing in relative strength with phase, at least for a year, after which we no longer have spectral coverage.
The bottom panel of Figure 7 shows the line profile of Hα, with the blue side reflected over the red side at the maximum flux after continuum removal.We do see evidence of diminished flux in the red wing of Hα at late phases in some SNe (most notable in SNe 2018evt and 2020uem), which can indicate formation of new dust in the post-shock CSM.S13 claim to observe this for all non-PTF SNe Ia-CSM in their sample starting at ∼75-100 days, while for the PTF SNe Ia-CSM they do not have spectra available post that phase range.For some BTS SNe Ia-CSM, we also do not have spectra available post 100 days which limits any analysis of this phenomenon for a large enough sample.and SNe 2018crl, 2018gkx, 2018evt, 2019agi, 2019ink, 2019rvb, 2020onv, 2020qxz in left panel.
The spectra were reduced and processed as outlined in §2.5 for the emission line analysis, the results of which are described in the next section.We used only good SNR SEDM spectra and intermediate resolution spectra for line identification and analysis.

Hα, Hβ and He I emission lines
To analyze the Hα line emission, we first fit the continuum level using the fit continuum function of the specutils Python package, where the continuum is estimated by a cubic function fitted on regions on each side of the line.We remove this continuum level and then fit the Hα line with a broad and a narrow component Gaussian function using the fit lines function of specutils which returns the best fit Gaussian model and the 1σ uncertainty on the model parameters.We generate 1000 sample models within 1σ uncertainties of the parameters centered around the best-fit values and calculate the intensity, flux and velocity (FWHM) of the broad and narrow components for each model.Then we take the median and standard deviation of the intensity, flux and velocity FWHM distributions to get their final best value and 1σ uncertainty.The equivalent width was also calculated for the Hα line using the model fit as well as directly from the data, and the difference between the values derived from model and data is reported as the error on the EW.All values are reported in Table 5.For 3 SNe in our sample, we have a series of intermediate resolution spectra through which we can trace the evolution of the Hα line with phase.Figure 8 shows this trend of the Hα line parameters (integrated flux in the top panel and equivalent width in the bottom panel) versus phase for all SNe in our sample.The un-filled markers represent the narrow emission while the filled markers represent the broad emission.For SNe where this analysis could be done on multiple spectra, we see that the Hα equivalent width generally increase over time, with some SNe showing fluctuations up to 100 days possibly due to interaction of ejecta with multiple CSM shells of varying density.For SN 2018evt, Yang et al. (2022) analyzed Hα line properties from a comprehensive spectral series data, which are plotted in Figure 8 in gray circles and seem to agree well with our analysis at comparable epochs.
From the Gaussian profile line fitting analysis of the Hα emission line, we found that the broader component has velocities ranging from ∼1000 to ∼4000 km s −1 (intermediate width) and the narrow component has velocities of about ∼200 km s −1 to ∼1000 km s −1 (see Figure 9).The narrow component could only be resolved down to ∼300 km s −1 limited by the mediocre resolution of the spectrographs used (KeckI/LRIS R∼800, P200/DBSP R∼1000, NOT/ALFOSC has R∼360).While we know that the narrow lines originate in the unshocked ionized CSM, the exact origin of the intermediate components is uncertain.They could arise from the post-shock gas behind the forward shock or from the shocked dense clumps in the CSM (Chugai & Danziger 1994).
The luminosities of the Hα line measured from the BTS SNe Ia-CSM lie in the range 2.5-37×10 40 erg s −1 which are comparable to the values from S13 who reported most of their SNe in the 1-10×10 40 erg s −1 range except one object that had a luminosity of 39×10 40 erg s −1 .From the broad Hα luminosity, we did a simple estimate of the mass-loss rate assuming spherically symmetric CSM deposited by a stationary wind ρ ∝ r −2 having velocity v w (Chugai 1991;Salamanca et al. 1998).The mass-loss rate Ṁ can be related to the broad Hα luminosity L Broad Hα as (Salamanca et al. 1998, their Eq.2) where v s is the shock velocity (obtained from the broad component velocity of the Hα line).We used a value of 100 km s −1 considering previous high resolution spectral studies of SNe Ia-CSM (Kotak & Meikle 2005;Aldering et al. 2006;Dilday et al. 2012) for v w as we cannot fully resolve the narrow component and a maximum value of 0.1 for the efficiency factor Hα (Salamanca et al. 1998).The mass-loss rates were estimated from the available spectra and are shown in Figure 10 as a function of years before explosion (t w = vst vw , where t is the phase of the spectra).For most SNe in the sample, the mass-loss rates lie between 0.001-0.02M yr −1 , except for SN 2019rvb which has ∼0.07 M yr −1 lost within 2 years prior the explosion.These rates are much higher than what could be attained from a red giant superwind (∼ 3 × 10 −4 M yr −1 ) but are comparable to previous estimates (calculated through multiple methods) for SNe Ia-CSM and require some unusual mechanism to reach such persistently higher mass-loss rates in the decades prior to explosion.Also to consider is that the simplistic assumption of spherical symmetry likely does not apply for SNe Ia-CSM.Evidence of multiple thin shells and asymmetric CSM was observed for PTF11kx (Dilday et al. 2012) and light curve modeling of SNe 1997cy and 2002ic suggested a better fit to a flat density profile rather than stationary wind (Chugai & Yungelson 2004).An asymmetric or clumpy CSM might be the norm for SNe Ia-CSM (and some SNe IIn) rather than the exception.
The same analysis as for the Hα line was also carried out for Hβ and He I λ5876 with a one component Gaussian fit.
For cases where a Gaussian model could not fit the data, we integrate the flux value in a 100 Å region centered at 5876 Å for He I.The Na ID absorption lines are also prevalent in some spectra and blend with the He I line, resulting in positive EWs for some SNe.The cumulative distributions of Hβ and He I equivalent widths are shown in the top and bottom panels of Figure 11 respectively.The Hβ median EW measured from the BTS SN Ia-CSM sample is 7.1 Å , close to the S13 value of ∼6 Å and quite weak compared to what S13 measured for SNe IIn (∼13 Å ).The overall cumulative distribution of Hβ EW is also comparable to the S13 SNe Ia-CSM rather than to the S13 SNe IIn.For the He I λ5876 line, the median EW measured for our BTS SN Ia-CSM sample, considering only significant emission features, is 2.4 Å .This is close to the value of ∼2 Å reported in S13, and again significantly different from their SN IIn value of ∼6 Å (∼4 Å with upper limits), however the overall distribution seems to be closer to the S13 SNe IIn (but still weaker) rather than to the S13 SNe Ia-CSM.This indicates that perhaps He I is not as good a discriminant between the populations compared to Hβ.Among the most Herich SNe in our sample are SNe 2019ibk, 2020uem, 2020xtg, 2020aekp and 2018evt, and these SNe also have the higher Hα equivalent widths in the sample.
Figure 12 plots the cumulative distribution of the Balmer decrements ( F Hα F Hβ ) measured for our sample SNe.The higher Balmer decrement values (>15) have large errors associated to them because of low SNR of the spectra from which they were derived, particularly near the Hβ line.Consistent with the results of S13, the SNe Ia-CSM from this sample also have a high median Balmer decrement value of ∼7 (∼5 in S13), indicating that the emission line mechanism is probably collisional excitation or self-absorption rather than recombination, from which the expected Balmer decrement value is ∼3.In the case of SNe Ia-CSM, if the CSM distribution consists of multiple shells as suggested for PTF11kx, moderately high densities could be created when fast moving ejecta overtake slowly moving thin dense CSM shells creating large enough optical depth in the Hα line which results in the Hβ transition decaying as Paα + Hα (Xu et al. 1992).For some individual SNe where multiple spectra are available, the Balmer decrement is observed to first increase and later on decrease with phase.

Host galaxies
We retrieved science-ready co-added images from the Galaxy Evolution Explorer (GALEX) general release 6/7 (Martin et al. 2005) Equivalent Width (Å) Line Velocity (km s 1 ages (Wright et al. 2010b) from the unWISE archive (Lang 2014b) 9 .We used the software package LAMBDAR (Lambda Adaptive Multi-Band Deblending Algorithm in R) (Wright et al. 2016) and tools presented in Schulze et al. (2021), to measure the brightness of the host galaxy.The spectral energy distribution (SED) was modelled with the software  2017) model for the ionized gas contribution.The priors were set as described in Schulze et al. (2021).
Figure 13 shows the log of star formation rate (SFR) as a function of stellar mass for hosts of BTS SNe Ia-CSM.We also use a Galaxy-zoo (Lintott et al. 2011) sample of elliptical and spiral galaxies (randomly sampled in the redshift range z = 0.015−0.05),and BTS SN Ia hosts as comparison samples collected by and used for comparison in Irani et al. (2022).We find the SN Ia-CSM host galaxy population to be consistent with late-type spirals and irregulars with recent star formation history. 4 out of 12 SNe have clearly spiral hosts, 3 have edge-on host galaxies, 4 seem to have irregulars as hosts and 1 has an unclear host type.Host galaxies of 10 out of 12 SNe have w2 − w3 measurements available which are all > 1 mag, putting them in late-type category (Irani et al. 2022), 1 (SN 2019rvb) does not have W3 measurement but has N U V − P S1 r ∼ 1 mag again putting it towards late-type and 1 (SN 2020abfe) does not have any of the above information available except the P S1 r band magnitude of 20.766, which is the faintest host galaxy (absolute SDSS r-band magnitude of −17.4) in our BTS SN Ia-CSM sample.As noted in S13, the SN Ia-CSM hosts of their sample had generally low luminosities (−19.1 < M r < −17.6) except MW like spiral hosts.Our BTS SN Ia-CSM host luminosities lie in the range of −21.8 < M r < −17.4 covering low to MW like luminosities.SN Ia-CSM rate:

Rates
where T is the duration of the survey, N is the number of transients that pass the quality cut, D max,i is the distance out to which the i th transient with peak absolute magnitude M i can be detected above the survey magnitude limit m lim (=19 mag for BTS SNe Ia-CSM) at peak light without any extinction, f sky is the average active survey coverage as a fraction of full sky, f ext is average reduction in effective survey volume due to Galactic extinction, f rec is the average recovery efficiency for a detectable transient within the survey coverage area, and f cl,i is the classification efficiency dependent on apparent magnitude.
The duration of the survey in which these 12 SNe Ia-CSM were detected is from 2018-05-01 to 2021-05-01, i.e.T = 3 years.We calculate f sky during this time period by averaging the sky area coverage of the public MSIP survey considering 3 day cadence for ZTF Phase I (2018-05-01 to 2020-10-31) and 2 day cadence for ZTF Phase II (since 2020-11-01), which turns out to be 12505 deg 2 for Phase I and 14831 deg 2 for Phase II, giving a mean f sky = 0.32.We use the same value of 0.82 for f ext as calculated in Perley et al. (2020) given there has not been any change in the number and positions of ZTF fields.
To estimate f rec , we consider SNe Ia-CSM brighter than −18.5 peak absolute magnitude and brighter than 18 apparent magnitude (total 5) of which 4 pass the quality cut, giving an f rec of 0.8.We take classification completeness of 0.75 at 19 mag, 0.9 at 18.5 mag and 1 at 17.2 mag and linearly interpolate in between these values to get f cl,i .
Then using H 0 = 70 km s −1 Mpc −1 , ignoring cosmological effects11 as in Perley et al. (2020) and applying a uniform K-correction (K = 2.5×log 10 (1 + z)), we get a rate of 29.35 +27.53  −21.37 Gpc −3 yr −1 for SNe Ia-CSM.We also calculate a SN Ia rate of 2.88 +0.28  −0.25 × 10 4 Gpc −3 yr −1 from SNe Ia observed in the same period following the same method, which is close to the value of 2.35×10 4 Gpc −3 yr −1 calculated in Perley et al. (2020).This puts SNe Ia-CSM to be 0.02-0.2% of SNe Ia.However this rate estimate should be considered a lower limit given various caveats in the correct identification of SNe Ia-CSM (see discussion §4.3).If the ambiguous classification cases outlined in Appendix A are considered to be SN Ia-CSM and included in the rate calculation, we obtain a rate upper limit of 97.7 +135.8  −77.3 Gpc −3 yr −1 , which is 0.07-0.8% of SNe Ia.

Precursor rates
The ZTF precursor rates were calculated following the method in Strotjohann et al. (2021) which studied the frequency of precursors in interacting SNe found in ZTF.Strotjohann et al. (2021) included 6 of the SNe Ia-CSM presented in this paper in addition to 4 other SNe Ia-CSM not in this paper (see Appendix A for details) for their search but did not find any robust 5σ precursor detections.This non-detection was concluded to be due to the small sample size of SNe Ia-CSM (or that they are more distant) compared to the SN IIn sample, so even if the precursors were as bright or frequent as for SNe IIn, it would be difficult to detect them.
The same search was here carried out for our larger sample by taking the ZTF forced photometry multi-band (g, r, i) light curves generated by the pipeline outlined in Masci et al. (2019) and stacking them in 1, 3 and 7-day long bins to search for faint outbursts.There were 7389 total available pre-explosion epochs for BTS SNe Ia-CSM, the earliest epoch being 1012 days prior to the explosion and the median phase 340 days prior.Hence the results are valid for typical SN Ia-CSM progenitors at about ∼1 year before the SN.We did not find any robust 5σ precursor detections.The upper limits for the precursor rates in different bands are shown in Figure 14, where the solid lines indicate up to what fraction of the time a precursor of a given brightness could have been detected while being consistent with the ZTF non-detections.A precursor of −15 magnitude could occur as frequently as ∼10% of the time given the ZTF non-detections.A continuous search for the precursors as more SNe Ia-CSM are found and classified and their sample size increases could yield a detection if the precursors are as frequent and bright as for SNe IIn.The dense and massive CSM around these objects is close enough to have been deposited within decades prior to the SN but the lack of precursors within 1 year indicates that there is likely no violent event that ejects a lot of mass in that period.Probing for precursors could potentially constrain the progenitor in at least some cases.For example, Soker et al. (2013) predicts for their core degenerate (CD) model for PTF11kx-like SNe release of significant energy (∼10 49 erg) before explosion over timescale of several years, implying a precursor 3-7 magnitudes fainter than the SN explosion spread over several years, peaking in the near-IR.The fastest declining SNe in our sample (SNe 2018crl, 2020qxz and 2020aekp) are also the ones that develop a plateau and show relatively stronger SN Ia-like absorption features in their early spectra.They seem to have a delayed start for the interaction like PTF11kx but not as fast a decline, and thus bridge the gap between PTF11kx and the rest of the strongly interacting SNe Ia-CSM.It remains to be seen how many SNe Ia are weakly interacting where the CSM interaction starts in earnest at timescales of ∼year or more after explosion, this requires searching for faint detections in care- fully calibrated forced photometry light curves (stacked to go fainter), a study currently undertaken by Terwel et al. (in prep).From the current sample, it appears that in addition to SNe Ia-CSM being intrinsically rare, delayed interaction SNe Ia-CSM are even rarer and only constitute about a quarter of all SNe Ia-CSM.This delayed interaction behaviour could also be an effect of asymmetric or clumpy CSM wherein part of the SN ejecta shine through depending on the viewing angle.Observational campaigns that capture the inner boundary of the CSM and the geometry robustly could shed light on the distribution of the inner CSM radius and reveal if it is a continuous distribution or if there are multiple progenitor scenarios within the SN Ia-CSM class.

Implications for progenitor based on observed mass loss
From Figure 10, the estimated mass-loss rates from a simple spherical treatment of the CSM and a stationary wind lie between ∼ 10 −3 to 10 −1 M yr −1 over a period of less than ∼ 60 years before explosion.That gives a total mass loss of ∼ 0.1 to ∼ 1 M .Dilday et al. (2012) estimated ∼ 5 M of CSM around PTF11kx while Graham et al. (2017) revised it to be ∼ 0.06 M .Light curve modeling of SN 1997cy and SN 2002ic by Chugai & Yungelson (2004) resulted in ∼ 5 M estimates for both SNe.Inserra et al. (2016) also fit analytical models to some SNe Ia-CSM and found the CSM mass to lie between 0.4 and 4.4 M .Since from Figure 5, the pseudo-bolometric luminosities of our SNe Ia-CSM lie somewhere between PTF11kx and SNe 1997cy, 2002ic and 2005gj, with SN 1999E somewhere in the middle, we can say that the total CSM mass in our sample of SN Ia-CSM should also be several solar masses.A WD+AGB star system has typically been suggested for historical SNe Ia-CSM to explain this massive CSM.The WD could either gain mass through Roche Lobe overflow (RLOF) from the companion that drives an optically thick wind (OTW) or merge with the core of the AGB star that then explodes in or soon after the common envelope phase.Meng & Podsiadlowski (2019) model WD+MS systems for their common envelope wind (CEW) model and find ∼ 1 M CSM around SNe Ia-CSM.Thus, given the large observed CSM mass range, the nature of the companion cannot be solely determined from total mass lost.High resolution spectroscopy that can resolve the narrow unshocked CSM wind velocity is also needed to determine the compactness of the companion.

Implications for progenitor based observed volumetric rate
Robust observed rate estimates for SNe Ia-CSM have been few and far between.Dilday et al. (2010) found 1 interacting SN Ia (SN 2005gj) in a sample of 79 SNe Ia at z < 0.15 in the SDSS-II SN survey, giving a rate of ∼1%.After the PTF11kx discovery in the Palomar Transient Factory (PTF) survey, the SN Ia-CSM rate was estimated to be ∼0.1% (1 in 1000 classified SNe Ia; Dilday et al. 2012) but without spectroscopic completeness determination.S13 identified 7 more SNe Ia-CSM from the PTF SN IIn sample, bumping up the estimate to ∼0.8%.With this sample we have improved the rate estimate, providing a robust value (along with an uncertainty estimate on that value) from an unbiased survey with high spectroscopic completeness up to 18.5 magnitude.However this rate quite possibly still underestimates the true value for two reasons.The first being possible thermonuclear SNe that are enshrouded so completely by CSM interaction that they are misclassified as SNe IIn in the absence of good early time data.In the BTS SN IIn sample, we found 6 SNe IIn to have ambiguous classifications which could possibly be SNe Ia-CSM and these are described in Appendix A. Including these ambiguous cases in rate estimation results in a rate upper limit of 0.07-0.8%for strongly interacting thermonuclear SNe, while excluding them gives an underestimated rate of 0.02-0.2%.
The second issue with the rates is if there is indeed a continuum of delayed interaction SNe Ia-CSM like PTF11kx, interaction in SNe Ia may present itself hundreds of days later at magnitudes fainter than ZTF's limit (∼20.5)resulting in those SNe not being counted when they may be sharing the same progenitor as the rest of the interacting SNe Ia-CSM.Lastly in some rare cases, the SN might appear normal in its light curve shape and duration (and thus would be missed by the selection criteria used in this paper) but seem to have peculiar narrow Hα in its spectrum or bright mid-IR flux (like in the case of SN 2020aaym; Thévenot et al. 2021).Han & Podsiadlowski (2006) predicted a rate of 0.1-1% for 02ic-like events for their delayed dynamical instability SD model but could not naturally explain the delayed interaction and multiple CSM shells in PTF11kx (which is relevant for some SNe in our sample).A symbiotic nova-like progenitor was suggested by Dilday et al. (2012) for PTF11kx and they quoted the theoretical rates for the same to lie between 1-30%, however the model could not explain the massive CSM.Soker et al. (2013) suggested a core degenerate (CD) scenario in which the explosion is set by the violent prompt merger of the core of the giant companion on to the WD and could naturally explain the massive CSM of PTF11kx (Livio & Riess 2003).Soker et al. (2013) estimated the occurrence of such SNe (M core + M W D 2 M and M env 4 M ) through population synthesis and found it to be 0.002 per 1000 M stars formed.Assuming ∼1-2 SNe Ia occur per 1000 M stars formed (Maoz et al. 2012), this corresponds to 0.1-0.2%,which compares well with our observed rate estimate.
The CEW model by Meng & Podsiadlowski (2019) predicts that the SNe Ia-CSM like objects could arise in the SD CEE scenario when CONe White Dwarfs (WD) steadily accrete material at the base of the CE without quickly spiraling in due to the driving of a CEW wind (10-100 km s −1 ).The WD explodes when it reaches the Chandrasekhar mass (1.38 M ) and could possibly explode within the CE before it is ejected.The CEW model predicts that 25-40% of the SNe Ia from CONe WD in Common envelope evolution with a Main Sequence (MS) companion will show SN Ia-CSM like properties.Meng & Podsiadlowski (2019) also give the ratio of SNe Ia from CONe WDs to normal SNe Ia from CO WDs to be between 1/9 and 1/5 (considering normal SNe Ia only come from CO WD + MS systems).Combining that with the estimate that roughly 10-20% of all SNe Ia may come from the SD scenario (Hayden et al. 2010;Bianco et al. 2011), SNe Ia-CSM from CONe WD according to the CEW model should be 0.28% to 1.6% of all SNe Ia.A spin-down before explosion of the WD (Justham 2011;Di Stefano & Kilic 2012) could also explain the time delay between explosion and interaction.
Soker (2022) estimated the common envelope to explosion delay time distribution (CEEDTD) shortly after the CEE (t CEED < 10 4 yr) from SN in planetary nebula rates and SN Ia-CSM observed rates to be roughly constant rather than having a t −1 dependence, that is the SN explosion could occur very soon after the CEE as well.Our observed rates are on the lower side compared to these theoretical model estimates but compare well within the observational uncertainties, though the CEW model seems to best account for the overall SNe Ia-CSM properties.

SUMMARY
In this paper, we have presented optical and mid-IR photometry, optical spectra and detailed analysis of 12 new SNe Ia-CSM identified in the Zwicky Transient Facility Bright Transient Survey, nearly doubling the total number of such objects discussed previously by Silverman et al. (2013).The properties of the sample extracted in this paper agree very well with similar analysis conducted in S13, particularly the median EW of Hβ is found to be significantly weaker in SNe Ia-CSM compared with SNe IIn and consequently the Balmer decrements are ubiquitously higher in SNe Ia-CSM.The brightness of SNe Ia-CSM in mid-IR is comparable to SNe IIn and observations of reduced flux in the red side of the Hα wing together with the mid-IR brightness points to formation of new dust in the cooling post-shock gas.The host galaxies of SNe Ia-CSM lie towards late-type galaxies with recent star formation.Unlike SNe IIn, no precursors were found within ∼1000 days before explosion for SNe Ia-CSM, which could be an observational bias (less number of SNe Ia-CSM compared to SNe IIn).We provide a robust rate estimate of 0.02-0.2% of all SNe Ia for SNe Ia-CSM on account of the BTS survey being unbiased and spectroscopically highly complete.The simple mass-loss rate estimates from broad Hα luminosity of ∼ 10 −2 M yr −1 are similar to previous estimates from various methods and indicate several solar masses of CSM around these SNe.The observed rate agrees well within the observational uncertainties with the CEW model by Meng & Podsiadlowski (2019) which can also explain the interaction delay and massive CSM.

8Figure 2 .
Figure2.Top: Location of our 12 SNe Ia-CSM in the peak absolute magnitude vs. rest-frame duration above half max phase space.The colored points are the BTS SNe Ia-CSM and the gray points are the rest of the BTS SNe Ia.Also shown with empty triangles are the SNe Ia-CSM from S13.The vertical arrows mark the upper limits to peak absolute magnitudes and horizontal arrows mark the lower limits to durations of SNe not having pre-peak coverage.Bottom: Change in magnitude 30 days after peak (∆m30) vs. rest-frame duration above 20% of peak-flux for BTS SNe Ia and SNe Ia-CSM.These criteria were used to filter out potential SNe Ia-CSM from all SNe Ia and demonstrate that SNe Ia-CSM occupy a distinct portion in this phase space.However some gray points (not SN Ia-CSM) remain on the longer duration side and are the false positive cases described in §2.1.

Figure 6 .
Figure 6.Spectral series of all SNe Ia-CSM presented in this paper.The rest-frame phases are shown alongside the spectra in each subplot and have been calculated using the explosion epoch estimate.The colors depict different instruments used to obtain this data.Major emission lines are marked with vertical dashed lines.

Figure 7 .
Figure 7. Top left: Early-time spectra of BTS SNe Ia-CSM with phases between 0 and 30 days since explosion compared to spectra of SNe 2011jb, 2005gj, 1991T and PTF11kx (phases in days since discovery).Top right: Late-time spectra of BTS SNe Ia-CSM (phases ranging from 40 to 370 days since explosion) compared to spectra of SNe 2011jb, 2005gj, 2010jl and PTF11kx (phases in days since discovery).Bottom left and right: Hα line profiles (post continuum removal) with the blue side reflected across the peak flux, marked by dashed lines.SNe 2020aekp, 2020abfe, 2020xtg and 2020uem in the right panel, and SNe 2018crl, 2018gkx, 2018evt, 2019agi, 2019ink, 2019rvb, 2020onv, 2020qxz in left panel.

Figure 8 .
Figure 8. Integrated fluxes and equivalent widths of Hα emission line with respect to SN phases for the BTS SN Ia-CSM sample.Broad component values are shown with filled markers and narrow component values with un-filled markers.SN 2018evt Hα luminosities and EWs presented in Yang et al. (2022) are also shown in gray circles.

Figure 9 .
Figure 9. Velocity of Hα emission line with respect to SN phases for the BTS SN Ia-CSM sample.Broad component values are shown with filled markers and narrow component values with unfilled markers.

Figure 10 .
Figure10.Mass-loss rates estimated from the luminosity of the broad component of Hα for the BTS SNe Ia-CSM.A value of 100 km s −1 was assumed for the wind velocity.

FollowingFigure 11 .
Figure 11.Cumulative distributions of equivalent width of Hβ and He I λ5876 emission lines calculated from the BTS SNe Ia-CSM (in grey) compared with the respective distributions presented in S13 for SNe Ia-CSM (blue) and SNe IIn (red).Vertical dashed lines mark the median EW of the distributions.

Figure 12 .
Figure 12.Cumulative distribution of Hα/Hβ intensity ratio (Balmer decrement) calculated from intermediate resolution spectra of BTS SN Ia-CSM sample (grey shaded region).The red line is the distribution of Balmer decrement of SNe IIn measured in S13, the blue line is the SN Ia-CSM Balmer decrement distribution from S13.The black circles are a few representative points indicating the high Balmer decrement values and the uncertainties on them.The vertical dashed line is the median Balmer decrement measured from BTS SNe Ia-CSM.

Figure 13 .
Figure13.Host galaxies of BTS SN Ia-CSM (black circles) on SFR vs stellar mass plot with Galaxy-zoo spiral (blue contours) and elliptical (red contours) galaxies for comparison.BTS SN Ia hosts are also shown for comparison in green circles.Equal sSFR lines are marked with grey dashed lines.
Fraction of SNe Ia-CSM with delayed interaction

Figure 14 .
Figure14.Precursor rate as a function of magnitude calculated from BTS SN Ia-CSM pre-explosion ZTF forced photometry stacked in 7-day bins.The different colored shaded regions correspond to different ZTF bands (r-red, g-green, i-grey).The solid lines depict the upper limits on fraction of the time a precursor of the corresponding magnitude would have been detected which is consistent with the ZTF non-detections.
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.The ZTF forced-photometry service was funded under the Heising-Simons Foundation grant #12540303 (PI: Graham).This work was supported by the GROWTH project(Kasliwal et al. 2019) funded by the National Science Foundation under PIRE Grant No 1545949.The Oskar Klein Centre was funded by the Swedish Research Council.Partially based on observations made with the Nordic Optical Telescope, operated by the Nordic Optical Telescope Scientific Association at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofisica de Canarias.Some of the data presented here were obtained with ALFOSC.Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and NASA; the observatory was made possible by the generous financial support of the W. M. Keck Foundation.The SED Machine is based upon work supported by the National Science Foundation under Grant No. 1106171.This work has made use of data from the Asteroid Terrestrial-impact Last Alert System (ATLAS) project.The Asteroid Terrestrial-impact Last Alert System (ATLAS) project is primarily funded to search for near earth asteroids through NASA grants NN12AR55G, 80NSSC18K0284, and 80NSSC18K1575; byproducts of the NEO search include images and catalogs from the survey area.The ATLAS science products have been made possible through the contributions of the University of Hawaii Institute for Astronomy, the Queen's University Belfast, the Space Telescope Science Institute, the South African Astronomical Observatory, and The Millennium Institute of Astrophysics (MAS), Chile.This research has made use of the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology.The Liverpool Telescope is operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias with financial support from the UK Science and Technology Facilities Council.Y. Sharma thanks the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant #1829740, the Brinson Foundation, and the Moore Foundation; her participation in the program has benefited this work.S. Schulze acknowledges support from the G.R.E.A.T research environment, funded by Vetenskapsrådet, the Swedish Research Council, project number 2016-06012.This work has been supported by the research project grant "Understanding the Dynamic Universe" funded by the Knut and Alice Wallenberg Foundation under Dnr KAW 2018.0067,The research of Y. Yang is supported through a Bengier-Winslow-Robertson Fellowship.Fritz (van der Walt et al. 2019; Duev et al. 2019) and GROWTH marshal (Kasliwal et al. 2019) (dynamic collaborative platforms for time-domain astronomy) were used in this work.

Table 1 .
Properties of the 12 BTS SNe Ia-CSM 1 Rest frame duration above 20% of r-band peak flux, uncertainty of ±2 − 3 days from ZTF cadence. 2 Corrected for Galactic extinction.

Table 2 .
Description of spectrographs used for follow-up and the corresponding data reduction pipelines

Table 3 .
Summary of optical spectra Pseudo-bolometric luminosity light curves of BTS SNe Ia-CSM compared with pseudo-bolometric light curves of SNe , SN 2005gj and PTF11kx, the Type Ia SN 1991T and the well-observed Type IIn SN 2010jl.Vertical gray regions mark typical SN Ia absorption features