© The Authors 2024

1 Introduction

The light emitted by astrophysical sources detected from Earth is subject to interaction with matter that lies along the line of sight between the source and the observer. Studying the effects of this line-of-sight material on the background light source gives information on the nature of that material, whether it is close to the light source, for example circumstellar material (CSM), within galaxies, that is, the interstellar medium (ISM), or between galaxies, that is to say the intergalactic medium (IGM). It is important to understand these processes: CSM for stellar mass loss, and ISM and IGM for galaxy enrichment and evolution. At the same time, this intervening material is also extremely important for understanding the light sources themselves as it absorbs the light, especially dimming it at optical wavelengths and re-emitting it at longer wavelengths. It thus impacts our understanding of the brightness and properties of the emitters.

This material in the line of sight may be composed of dust and gas. Dust grains are a product of stellar evolution and are responsible for the effect of absorption and reddening of optical starlight (e.g. Draine 2011). Molecular gas is known to exist, for example in the form of hydrogen (H₂) and carbon monoxide (CO) (Dame et al. 2001), but also in the form of diffuse interstellar bands (DIBs) which trace unidentified molecules as of yet (e.g. Lan et al. 2015). On the other hand, the atomic gas phase component of the intervening material is mostly composed of hydrogen. Nevertheless, a fraction of metals is also present and provides valuable information on the chemistry, ionization state, and gas temperature of the line-of-sight matter. These different species absorb the light at characteristic wavelengths, generating detectable absorption lines in the spectroscopic observations.

Optical interstellar absorption lines in the form of sodium, calcium, potassium, and DIBs have long been observed in spectroscopic binaries (e.g. Heger 1919a,b; Young 1922), early-type stars (Struve 1928; Wilson & Merrill 1937; Merrill 1937; Sanford 1937; Evans 1941), novae (Merrill 1935), galaxies (e.g. Heckman & Lehnert 2000), active galactic nuclei (e.g. Baron et al. 2016), and supernovae (SNe; e.g. Penston & Blades 1980). In contrast to material moving at a high velocity, as is the case for SN ejecta, these intervening lines appear much narrower than the characteristic P-Cygni profiles of fast-expanding material. Today, narrow interstellar lines are routinely detected in the optical spectra of a great number of astrophysical transients (e.g. Galazutdinov 2005; Megier et al. 2009; Chen et al. 2010; Park et al. 2015). The relation of the strength of these lines with a colour excess and extinction has been studied in depth (e.g. Merrill & Wilson 1938; Spitzer 1948; Buscombe & Kennedy 1968; Hobbs 1974), providing empirical relations between the amount of dust extinction and the strength of the gas lines (e.g. Richmond et al. 1994; Munari & Zwitter 1997; Turatto et al. 2003; Poznanski et al. 2012; Murga et al. 2015), with caveats for the sodium line since it saturates (Munari & Zwitter 1997; Poznanski et al. 2011). Alternative lines such as potassium (e.g. Hobbs 1974; Munari & Zwitter 1997; Galazutdinov 2005; Phillips et al. 2013) and DIBs (Kos & Zwitter 2013; Phillips et al. 2013; Krełowski et al. 2019) have been claimed to be superior extinction tracers.

Supernovae offer a particularly good laboratory to investigate the material in the line of sight as they are bright and are found in many different environments. On the other hand, understanding this material is also key to understanding SN progenitor systems and their explosions, as it can give us information on stellar evolution and mass loss. The study of narrow line-of-sight lines in transient objects has indeed increased during the last few years, partly due to the detection of time-variable absorption features in high-resolution spectra. While intermediate-luminosity red transients embedded in dusty cocoons show clear signs of sodium line strength evolution (Byrne et al. 2023), the only SN coming from the core collapse of a massive star showing this evolution is the broad-lined type Ic, SN 2012ap (Milisavljevic et al. 2013), with changes in the strength of the DIBs but no variations in the sodium lines. The authors conclude that the material producing these variations is nearby, and the SN interacts with it. More surprisingly, a select number of type Ia SNe – the thermonuclear runaway of white dwarfs in binary systems – have shown evolving sodium lines in high-resolution spectra: SN 2006X (Patat et al. 2007), SN 1999cl (Blondin et al. 2009), SN 2007le (Simon et al. 2009), SN 2013gh (Ferretti et al. 2016), and SN 2014J with time-varying potassium instead of sodium (Graham et al. 2015). This has significant implications for the debate surrounding the SN progenitor system: it has been suggested that such variable lines could originate from the CSM ejected by the progenitor system, in principle favouring the single degenerate scenario of a binary system composed of a white dwarf and a main sequence or red giant (e.g. Moore & Bildsten 2012). However, double degenerate models composed of two white dwarfs could potentially also exhibit CSM characteristics (Raskin & Kasen 2013; Shen et al. 2013; Levanon & Soker 2019). Such CSM claims are further emphasized with the observation of a statistical excess of blueshifted (vs. redshifted) sodium absorption lines in high-resolution SN Ia spectra (Sternberg et al. 2014; Maguire et al. 2013; Phillips et al. 2013; Clark et al. 2021), which may also be related to other SN properties such as ejecta velocity and SN colour (Foley et al. 2012).

On the other hand, measurements of narrow lines from intervening material in SN Ia spectra of a lower resolution have also shown intriguing correlations between line strength and intrinsic SN properties such as the nebular velocity shifts (Förster et al. 2012) and the late colour decline rate (Förster et al. 2013), arguing again for a possible interaction of nearby material with a sub-group of SNe Ia. This has been recently confirmed by Wang et al. (2019) also with low- and mid-resolution spectra and they additionally find possible line evolution in the interacting SNe Ia.

A significant caveat of the above studies is that lower resolution and signal-to-noise ratio (S/N) spectra may affect or even prevent proper measurements. Can low and mid-resolution spectra be used to infer narrow line equivalent widths (EWs) and velocities? Can weaker lines such as DIBs and potassium reliably be detected in such spectra? Can the presence of the high-velocity, broad features from the SN ejecta interfere with the narrow line measurements? Are there any biases that could simulate an apparent evolution or blueshift? Given the large amount of SN spectra gathered in the literature and the scarcity of methods to investigate the material in the light of sight of the SN, the capability to robustly use these spectra for narrow line studies could provide powerful, statistically meaningful results on the presence and distribution of these species in the intervening matter.

In order to investigate narrow lines, we used simulations and a large sample of observed SNe from the literature that consists of more than 1000 objects of different types – including SNe Ia, SNe II (classical plateau and linearly declining hydrogen-rich SNe), stripped-envelope SNe (SESNe) with SNe Ic, SNe Ib, and SNe IIb, and interacting SNe including SNe IIn, SNe Ibn, and SNe Icn. The spectra and simulations include various spectral resolutions that allowed us to test several different possible biases in the measurements of the EWs and velocities of narrow line-of-sight lines. After defining the limitations of our technique, we also looked into the likelihood that these lines in SN spectra evolve over time. In forthcoming papers, we will investigate the differences in these lines for various SN types, SN properties, and environmental characteristics.

The paper is organized as follows. Section 2 introduces the lines from intervening matter analysed in this study. We present the data sample in Sect. 3 and the automated line measurement in Sect. 4. The systematic biases are presented in Sect. 5, while the evolution analysis and discussion are in Sect. 6. We conclude in Sect. 7.

2 Line-of-sight absorption lines

The properties of the ISM of a SN host galaxy and, more specifically, the environment where the SN exploded – or of the CSM ejected by the progenitor prior to explosion – can be studied with spectroscopy. The spectral resolution, usually given by R = λ/∆λ reveals the capacity of the instrument to resolve narrower lines; we define here high-resolution spectra as those with R > 10000 (or Δv < 30 km s⁻¹; Appenzeller 1986) and mid-resolution as 5000 < R < 1000 (or 30 km s⁻¹ < Δv < 60 km s⁻¹). In nearby objects exploding in reddened environments, it is possible to detect narrow lines close to the host galaxy’s rest-frame line wavelength, but also due to the MW close to the observer-frame line wavelength, at a variety of spectral resolutions and signal-to-noise ratios, as will be shown later (Sects. 5.2 and 5.3).

We outline here the main optical absorption lines that can be measured and utilized to investigate the ISM and the extinction towards SNe.

• Sodium (Na I D) lines are very prominent and the typical tracers of the ISM. They were discovered in the solar spectrum by Fraunhofer in 1814. They are the doublet from the fine structure splitting of the excitation states of neutral sodium at λ = 5891.58 Å for D₂ and λ = 5897.56 Å for D₁. Neutral sodium has an ionization potential of 5.14 eV (Morton 2003).

• Calcium (Ca II H&K) lines were also seen by Fraunhofer in sunlight and, similarly, they correspond to the doublet of the fine structure splitting of the excitation states of the singly ionized calcium at λ = 3934.8 Å for K and λ = 3969.6 Å for H, with an ionization potential of 11.87 eV.

• Potassium (K I) lines were first detected by Dunham (1937). This resonant doublet from the fine structure splitting of the neutral potassium excitation states occurs at λ = 7664.90 Å for K₁ and λ = 7698.97 Å for K₂. Neutral potassium has an ionization potential of 4.34 eV.

• The diffuse interstellar bands (DIBs) were first discovered in stellar spectra by Heger (1922). Today, hundreds of DIBs from unidentified molecules in the ISM have been discovered (see Herbig 1995). Proposed carriers include polycyclic aromatic hydrocarbons (PAHs), fullerenes and other hydrocarbons ions and molecules (Salama et al. 1999; Kroto et al. 1985; Motylewski et al. 2000). We choose here three of the strongest and more isolated DIBs at λ = 4428.2, 5780.5, 6283.8 Å.

Table 1 summarizes these narrow lines and their characteristics. We also include the ‘integration window’, a velocity interval around the central wavelength in which the line will be measured (see Sect. 4.1). It is related to the strength of the line and therefore the wavelength range it may cover.

SNe in low-reddening environments may only display the strongest line-of-sight lines, which normally are the Na I D absorption features. If the strength of the absorption line is very low, it indicates negligible material in the line ofsightand very low reddening. Figure 1 shows the location of the narrow Na I D features for some low-redshift (ɀ < 0.1) SNe observed with a low (top panels) and high (bottom panels) resolution instrument. In SNe II (left panels), the narrow Na I D lines are on the top of the He I λ5876 broad line from the ejecta at early phases, but after 30–40 days, they are around the peak of the broad Na I D emission line (from the SN). In SESNe (middle panels), they are between the He I/Na I D and the Si II/He I or Hα (in SNe IIb), while in SNe Ia (right panels), the narrow Na I D features are between the Si II lines. The effect these broad lines have in the narrow lines will be further discussed in Sect. 5.

In high-resolution spectra, the narrow Na I D lines can be resolved. It is possible to detect four lines in the spectra (see, for example, bottom panel of Fig. 1), two from the MW (λ5890 and λ5896) and two from the host galaxy (λ5890 × (1 + ɀ) and λ5896 × (1 + ɀ), where ɀ is the SN redshift), as long as a single cloud component is present. For low-resolution spectra, only two lines can be detected (one that corresponds to the MW (λ5893, on average) and another related to the host galaxy (λ5893 × (1 + ɀ)). When the SN is very nearby (ɀ < 0.04; e.g. SN 1987A), the MW and the host galaxy lines can be blended, and just a single narrow absorption line is visible. These cases are not included in the current analysis (see the selection criteria in Sects. 3 and 5).

The narrow line-of-sight lines play an important role in estimating the reddening in the line of sight of SNe. Theoretically, this is somewhat expected as the line EW of a given gas species X is given by: $${\text{EW}}_{\text{X}}={\displaystyle \int \left(1-e^{-{\tau }_X(\lambda )}\right)}\text{d}\lambda ,$$(1)

where τ_X is the optical depth, which is a function of the column density of the species, N_X, multiplied by the absorption coefficient, κ₀. In turn, the gas column density depends on the gas-to-dust ratio, ƒ_gd, the mass fraction of the gas, ƒ_X, and the ionization fraction, ƒ_ion. For low optical depths, EW_X ∝ τ_X ∝ N_X, so that: $${\text{EW}}_{\text{X}}\propto {\tau }_X={\kappa }_0{N}_X={\kappa }_0\left(1-f_{\text{ion\hspace{0.17em}}}\right)f_X{f}_{\text{gd}}{N}_{\text{dust\hspace{0.17em}}},$$(2)

where the dust column density, N_dust, is directly proportional to the extinction A_V: $${\text{EW}}_{\text{X}}\propto {\kappa }_0\left(1-f_{\text{ion\hspace{0.17em}}}\right)f_X{f}_{\text{gd}}{A}_V.$$(3)

As the optical depth increases, the line profile becomes saturated, and the equivalent width does not increase linearly with column density. Indeed, the strength at the central wavelength stalls whereas the Doppler broadening increases, so that ${\text{EW}}_X\propto \sqrt{\ln N_X}$. At even larger optical depths, the collisional broadening dominates and ${\text{EW}}_X\propto \sqrt{N_X}$.

The linear relation between EW and A_V at low optical depths (Eq. (3)) depends strongly on other factors that are quite uncertain and vary substantially in our MW. For example, the low ionization potentials of the gaseous metals (below the Hydrogen ionization potential of 13.6 eV, see Table 1) mean that they are easily photo-ionized, changing the ionization fraction. Multiple gas clouds in the line of sight can also create several spectral components whose varying velocities cannot be disentangled in mid- and low-resolution spectra and thus contribute to the saturation of the line.

Despite these uncertainties in the theoretical relation, it has been shown empirically, as mentioned above, that the strength of the narrow Na I D intervening lines correlates with dust extinction, based either on nearby stars (e.g. Munari & Zwitter 1997), SNe (Turatto et al. 2003; Sollerman et al. 2005) or large samples of extragalactic sources whose light is absorbed in the MW gas and dust (Poznanski et al. 2012; Murga et al. 2015). The trends are limited by the large dispersion and a distinct saturation for EW values larger than 0.5 Å in low-resolution data (Munari & Zwitter 1997; Poznanski et al. 2011), and 1.8 Å in higher resolution spectra (1 Å for Na I D₂ and 0.8 Å for Na I D₁; Poznanski et al. 2012) respectively. Poznanski et al. (2011) suggested that the EW estimated from low-resolution spectra is a bad estimator of E(B − V), while Phillips et al. (2013) argued that the Na I D gives an overabundance of E(B − V) for several SNe Ia and the DIB at 5780 Å is a better tracer of reddening than Na I D. Indeed, Baron et al. (2015) demonstrate that eight different DIBs correlate well with extinction for a large sample of extragalactic objects from the Sloan Digital Sky Survey.

In the context of SNe, the presence of the high-velocity, wide features from the SN makes the measurement of the EW difficult (see Sect. 5.4). Additionally, time-variable Na I D absorption lines detected in a few SNe Ia (Patat et al. 2007; Blondin et al. 2009; Simon et al. 2009) also complicate their use as reddening tracers. However, considering their easy identification and the scarcity of other methods to constrain dust extinction along the line of sight, their use is still widespread. That is why a robust methodology and a careful study of systematic biases, as the one presented here, are needed.

3 Sample

The data analysed in this paper consists of a large heterogeneous sample from the literature obtained through large web compilations such as SUSPECT (Richardson et al. 2001), WIS-EREP¹ (Yaron & Gal-Yam 2012), the Open Supernova Catalog (Guillochon et al. 2017), and from individual journal publications. It contains SNe Ia from the thermonuclear runaway of white dwarfs in binary systems, including the peculiar 91bg-like, 91T-like, 02es-like and Iax; it contains hydrogen-rich SNe (SN II) from the core-collapse of massive stars that ave retained their hydrogen envelopes, including both old light-curve classifications, SNe IIP (plateau) and SNe IIL (fast-decline); it contains stripped-envelope (SESNe) from the core-collapse of massive stars that have lost their outer layers: IIb with hydrogen shown only initially in their spectra but not later on, Ib without hydrogen throughout their evolution, Ic without hydrogen nor helium in their spectra and Ic-BL with broad lines from a central engine; and it contains various transients related to stellar outbursts and death of massive stars with clear signs of interaction with circumstellar material like SNe IIn, SNe Ibn, SNe Icn, and SN impostors, collectively called here “SNe-int”.

This dataset consists of 12 257 spectra of 1788 SNe observed between 1937 and 2023. 2538 of these spectra correspond to 327 SNe II, 2053 spectra to 328 SESNe, 1269 spectra to 116 SNe-int and 6397 spectra to 1017 SNe Ia. The sample also includes 174 spectra with mid- and high-resolution. Most of them are from SNe Ia. Only 19 are from SESNe, and 60 are from SNe II. All spectra are corrected for MW extinction with the maps of Schlafly & Finkbeiner (2011) and brought to the rest frame. Due to the heterogeneous nature of the sources, this sample is not complete, neither in volume nor in magnitude. The sample details can be found in Table D.1.

Figure 2 (left panel) shows the redshift distribution of our sample. The mean redshift value is 0.026, while the median is 0.019. The nearest object is SN 1987A with a redshift of 0.00001, while the farthest object is iPTF-13ajg with a redshift of 0.74. In the right panel of Fig. 2, we present the distribution of spectra per object. 720 SNe have at least two spectra, 873 SNe have at least three spectra, and 937 SNe have between 3 and 10 spectra. The objects with the most data are the type Ia SN 2005cf (115), SN 2006X (85) and SN 2001V (82). On average, we have 6.7 spectra per SN. For the sample of SNe II, the average is 7.6, 6.4 for SESNe, and 6.4 for SNe Ia.

To give a better idea of the spectral quality of this large and diverse sample, we show in the top panel of Fig. 3 the average resolution of all spectra, where we define this resolution as the average sampling around 10 000 km s⁻¹ of the narrow line (in this case Na I D): Δv_res = ∆λ_pc/λ_rest and we also show the corresponding pseudo-resolution, R_p = λ_rest/∆λ_P, in the upper axis. ∆λ_P is the wavelength sampling provided in the spectral files instead of the full-width half-maximum of unresolved lines, which characterizes an instrument resolution. This is why we call R_p a pseudo-resolution, which is only a proxy for the real spectral resolution of the instrument. In particular, some spectra have higher wavelength spacing without meaning the spectral resolution is better. This is especially true for many old SN spectra from years before 1995 that have a narrower spacing data format but are not of particularly high resolution. These represent most of the first peak shown in the distribution, although some real high-resolution spectra are also present. Other peaks in the distribution are mostly but not exclusively populated by the spectra of large surveys with its instruments and their characteristic wavelength spacing: the supernova group at the Harvard-Smithsonian Center for Astrophysics² (CfA; e.g. Blondin et al. 2012; Modjaz et al. 2014; Hicken et al. 2017), the UC Berkeley SN group (Berkeley; e.g. Silverman et al. 2012; Shivvers et al. 2019; Stahl et al. 2020) and Carnegie Supernova Project³ (CSP; e.g. Folatelli et al. 2013; Gutiérrez et al. 2017; Stritzinger et al. 2023) are the most characteristic ones.

On the bottom panel of Fig. 3, we also show the average S/N in the same wavelength regime calculated with the inverse of the root mean square (RMS) of the difference between the original spectrum and a smoothed spectrum with a cosine kernel. The distribution is exponentially falling with a mode around 30 and a median of 150.

With he aim of doing a systematic analysis and decreasing the bias in our sample, we will define several criteria for our data selection (see Table 2) that will be discussed in the next sections.

4 Automated line measurement

In this section, we outline the automated methodology we have developed to measure EWs and velocities in a diverse set of spectra.

4.1 Equivalent width measurement

To measure the equivalent width of the various line-of-sight narrow lines from the host galaxy and/or the MW, we use a fully automated code inspired by Förster et al. (2012). First, a continuum is found by defining several nodes around the centre wavelength of the narrow line with separations between 50 and 2000 km s⁻¹. The nodes are evenly spaced (in velocity space), and the flux at each node is obtained by smoothing locally the spectrum with a cosine kernel.

Once the continuum is defined, we measure the EW: $$\text{EW}={\displaystyle \sum _i\frac{c_i-f_i}{c_i}}\Delta \lambda ={\displaystyle \sum _{i}1}-\frac{f_i}{c_i}\Delta \lambda ,$$(4)

where the flux at each pixel, f_i, is compared to the continuum at that pixel c_i, and ∆λ is the wavelength separation between pixels in Å.

We sum the area under/over the continuum in a fixed velocity range around the central wavelength, similar to Förster et al. (2012). The velocity interval changes for each line that we consider and is given in Table 1. These intervals were estimated by a visual inspection of the spectra with the most pronounced narrow features, e.g. SN 2Ol3fc (Kangas et al. 2016) or SN 2OO3cg (Elias-Rosa et al. 2006); it is largest for the stronger Na I D(+1000 km s⁻¹) which includes two lines. We note that the EW measured in this way sums all contributions in the given range, meaning that, in principle, more than one feature could be present. Additionally, it is important to remember that emission lines are taken into account giving possible negative values of EWs. The net sum should be close to zero when there is only noise.

The error associated with the EW of a single spectrum has four components: (i) the flux error, σ_f, which is calculated from the RMS between the flux and the continuum outside the range of the line (but within +10 000 km s⁻¹), and an average of this outside error is assumed within the range of the line (see solid blue lines in e.g. Figs. 4 and C.2); (ii) the error in the continuum, σ_c, which comes from the dispersion obtained when assuming a difference of +25% in the node separation that determines the continuum (see red dotted lines in Figs. 4 and C.2); (iii) the error in the integration window for which we re-calculate the EW at ±25% of the range given in Table 1 for each line; and (iv) a systematic error, σ_sys, from the dispersion obtained when measuring the line at eight different positions where no line is expected (see vertical blue dotted lines in Figs. 4 and C.2) and the EW should be zero. The first two sources of error are propagated in Eq. (4) to obtain the error on the EW, and together with the latter, they are added in quadrature.

The node separation for the continuum plays a crucial role in this technique (see Figs. 4 and C.2), so the optimal separation is found by maximizing the absolute value of the signal-to-noise of the EW (see Fig. 5). This is done by calculating the EW for various node separations in bins of 50 km s⁻¹ and through a secondary loop in bins of 10 km s⁻¹ centred at the highest |EW|/σ_EW found in the first iteration. The absolute value is used to make sure that we do not necessarily force positive values on the EW. We also require that the profile of the narrow line intersects with the continuum within the velocity window considered.

To verify the accuracy of the automated method, we compare the results obtained by measuring the EW of Na I D manually using IRAF (Tody 1986, 1993) and with our automatic software. Figure 6 shows the comparison for a sub-sample of 1145 low-resolution spectra and 23 high-resolution spectra. The methods give consistent results with a median absolute deviation of 0.01±0.22 Å.

4.2 Velocity measurement

The velocity measurement is trickier compared to the EW: for weak or absent lines, it cannot be done, while for stronger lines, there might be multiple components that are only resolvable with high-resolution spectroscopy. Obtaining an unbiased estimate of the velocity similar to the EW is perhaps impossible, but we attempt two approaches. The first is model independent and consists simply of a line-weighted average wavelength of the area under/over the continuum measured above (grey shaded areas in e.g. Fig. 4) compared to the rest wavelength of the line: $$\text{v}\left({\text{km\hspace{0.17em}s}}^{-1}\right)=\frac{c}{{\lambda }_{\text{rest\hspace{0.17em}}}}\left[\frac{{\sum }_i\left(1-f_i/c_i\right){\lambda }_i}{{\sum }_i\left(1-f_i/c_i\right)}-{\lambda }_{\text{rest\hspace{0.17em}}}\right],$$(5)

where c is the speed of light in vacuum and λ_rest the restfame wavelength of the line.

The first source of uncertainty stems from the wavelength calibration, which may vary substantially across our heterogeneous sample and is impractical to check individually given the lack of, e.g. sky spectra. Instead, we measure the emission lines from the gas phase of the environment, sometimes present in SN spectra, and calculate the average velocity dispersion between those lines. We obtain a characteristic dispersion of 35 km s⁻¹ (see Appendix A), which we take as an uncertainty floor from the wavelength calibration. Additional sources of uncertainty come from re-measuring the velocities when (i) the flux error is added/subtracted to the original flux, (ii) when the continuum is changed by ±25%, and (iii) when the integration window is changed by ±25%.

This simple technique works well when the narrow line is relatively symmetric but otherwise suffers from several biases. For example, in the case of the Na I D, if one of the two doublet components is stronger, the flux-weighted wavelength will be affected, as seen for SN 2002bo (see Fig. 7). The bias may be even stronger in the case of multiple components with different velocities.

Another approach, commonly used in the literature, consists of a parametric profile fitting with, e.g. one or multiple Gaussians. For this, one has to assume the number of components (i.e. number of Gaussians), which one may not be able to distinguish in the case of low-resolution spectra. In Fig. 7, we show a fit with two Gaussians with a fixed wavelength separation of 6 Å between the two sodium lines and a common Gaussian width. We directly fit the ratio of the flux-to-continuum using the continuum obtained in the previous section. For this particular case, the line-weighted average estimates a slightly bluer velocity because of the stronger line of the doublet. A comparison of the automated EW and flux-weighted velocity measurements of the sodium doublet with the values obtained from a double Gaussian fit is shown in Fig. 8. The match of the EWs is remarkable, while there is a larger dispersion for the velocity – still consistent within the large uncertainties. This dispersion is mostly driven by cases where the EW is very weak, and the measurement of the velocity becomes very uncertain for both methods. Spectra with stronger lines (see green points) agree better and present lower uncertainties. The two approaches are overall quite consistent, and we use our automated technique as our primary method but compare, when necessary, with the Gaussian approach.

4.3 SN average of multiple spectra

In the case where multiple spectra of a SN are available, and no evolution is expected or found (see Sect. 6), one can get a better estimate of the EW and velocity of the intervening line by averaging the individual values found in the previous sections (Sects. 4.1 and 4.2). Nonetheless, given the heterogeneous set of spectra taken with different instruments and with different underlying SN ejecta throughout the SN evolution, we find that a better approach consists of taking the median of the individual flux-to-continuum ratios to generate a stacked ratio with higher signal-to-noise. This approach is similar to the work done by Poznanski et al. (2012) for NaiD measured in galaxy spectra. We show an example of 12 stacked flux-to-continuum spectra of SN 2003cg around the DIB-5780 Å line in Fig. 9. After the stacking, the EW and velocity are calculated in the same way as for individual spectra. To measure the corresponding error of this stacked value, we do a bootstrap with 100 iterations in which the same number of flux-to-continuum spectra are randomly taken from the original sample (some spectra may be repeated or missing in each iteration), and the EW and velocity are measured each time. The final asymmetric uncertainties are given by 16 and 84% of the distribution. The blue lines in Fig. 9 contain the 1σ bootstrap realizations for illustration purposes. As for individual spectra, we also include the systematic errors from the continuum, the integration window and the dispersion of measuring the line at other wavelengths. A comparison of the EW measured from stacked spectra with the median or weighted averages of individual EW measurements is shown in Fig. C.1.

5 Spurious temporal evolution of narrow lines and systematic biases

Although the methodology laid out in the previous section correctly measures the EW of narrow lines of spectra of various resolutions, there are some general caveats inherent to the use of low-resolution spectra that we investigate further. Some of these biases become evident when studying the behaviour of the narrow lines over time: we see a conspicuous evolution of the EW for some SNe. In particular, in Fig. 10, we show the evolution of the Na I D EW for the MW line for different SN types measured with the automated technique. All objects show variations in their EWs; however, in SN l999em and SN 2007gr, a constant decrease is striking. Narrow intervening lines from the MW are expected to be generally constant as a function of time since they come from material in our own Galaxy that does not interact with the SN; therefore, variations in their EWs suggest the interference of external factors. A similar, albeit more extreme, behaviour was observed for the narrow lines from the host galaxy (see Sect. 5.4). This is worrying since it can lead to false physical interpretations like interaction with circumstellar material. These results motivate us to explore the possible biases that affect the strength of the narrow line-of-sight lines. They are described in detail in the following subsections.

5.1 Redshift

In order to ensure that the narrow Na I D absorption line from the MW and the host galaxy are separately resolved (i.e. to ensure they are not blended), we set a lower limit cut at 0.004 on the redshift of our sample (independently of the spectral resolution). This eliminates a small number of very nearby SNe (see Table 2). As will be shown later in Sect. 5.4, the redshift also plays an important role in the location of the broad lines of SNe with respect to the narrow intervening lines affecting their measurement.

5.2 Spectral resolution

In this section, we explore the impact of varying spectral resolutions in measuring the EW and velocity of the lines from intervening matter. We pick high-resolution spectra (∆υ_res < 50 km s⁻¹) and smooth them with a cosine kernel (see Sect. 4) with increasing velocity windows to simulate lower resolutions (Fig. C.3). We try velocity windows spanning 50–500 km s⁻¹, which conservatively cover our sample’s average velocity resolutions (see Fig. 3). As long as the continuum is relatively flat (Sect. 5.4), we find excellent agreement in the measured EW and velocity of the line-of-sight line for all resolution windows with the original high-resolution values. This is shown in Fig. 11 for several SNe of different Na I D line strengths. We highlight that the uncertainties in the EW become larger, as expected, for worse resolutions. Therefore, no cut on the resolution is applied to our sample.

5.3 Signal-to-noise

We now evaluate the performance of the automated technique when subject to varying S/N ratios in the spectral range of the narrow line. The S/N, as explained in Sect. 3, is estimated by calculating the RMS of the observed spectrum with respect to a smoothed version. We perform a simulation by taking spectra with high signal-to-noise and randomly perturbing them according to various decreasing S/N ratios (see Fig. C.4). For each given S/N, we do 10 simulations and take the mean and standard deviation to show the effect in Fig. 12. The measured EWs are remarkably consistent with decreasing S/N, and the uncertainties measured by the automated code (dashed error bars) are always larger than the standard deviation (solid error bars) of the simulations. Nevertheless, we do see that for weak lines (e.g. those in the spectra of SN 2002an and SN 2007oc), the average EW starts to increase at very low S/N, albeit always within the errors. Because of this, we add a conservative additional minimal S/N cut of 15 that rejects around 1500 spectra (see Table 2).

The recovered velocities are also consistent, but we note that, as expected, there is a large dispersion for weak lines as there is not much line flux to obtain the flux-weighted average wavelength. However, the uncertainties on the velocities are always much larger than the standard deviations of the simulated spectra (measured uncertainties are larger than 300 km s⁻¹ and are not shown in the figure).

5.4 Broad line interference

When measuring the narrow line-of-sight lines in individual spectra of SNe and investigating their evolution, it becomes evident that several objects have a seemingly temporal change in the strength of Na I D. Particularly for SNe II, there is often a strong apparent decrease in the EW as a function of time, but variations can also be seen for SNe Ia (see, e.g. Fig. 13). Upon closer examination of the SN spectra, these changes appear to be strongly linked to the appearance and evolution of the emission and absorption components in the P-Cygni profile of the broad Na I D (for SNe II) and Si II 6355 Å (for SNe Ia). If the variations are connected to the development of broad lines (see Gutiérrez et al. 2017 for the appearance and evolution of P-Cygni lines in SNe II), they do not constitute real intrinsic variations. In that case, one will infer a wrong evolution and, more generally, EW measurements that depend on the underlying profile and are thus not a reliable reflection of the column density of the intervening material.

To investigate such bias in more detail, we create a simple simulation in which we input fake narrow lines into various differing P-Cygni profiles and calculate the recovered EW. For these simulations, we take both observed optical spectra of seven SNe II, three SNe Ia, and five SESNe, together with two SN radiative transfer models: the delayed-detonation model of SNe Ia (Blondin et al. 2015), and the ml5mlt3 for SNe II (Dessart & Hillier 2011; Dessart et al. 2014). Most of the selected observed SNe do not show evident narrow Na I D lines; however, when the lines are clearly detected, they are removed from the spectra by interpolating the continuum to the line wavelengths so that we can then add simulated lines of our choice.

Once the spectra are clean from the Na I D lines, we smooth them using a moving average box by convolution with a cosine kernel. Although most spectra have a good S/N, the smoothing process provides easier identification of the continuum and the narrow lines. Following the smoothing, the spectra are normalized at the average position of the narrow lines (λ5893). In parallel, we generate narrow Na I D absorption features with similar characteristics to those observed in SN spectra. These features are created for two cases: (i) for low-resolution spectra that assume the Na I D doublet is blended and just one single feature is detected, and (ii) for high-resolution spectra that produce the Na I D doublet (two resolved lines). For the two cases, we test different strengths of the narrow line.

To see the effect of the narrow lines in the SN spectra, we add the narrow Na I D feature at different locations offset from the central wavelength in all SN types to simulate different ejecta velocities and redshifts. Figure 14 shows the type II SN 2003bn spectra at 100 days from explosion centred on the broad Na I D P-Cygni profile with the added narrow Na I D feature at different locations. Here, the narrow Na I D features were added for low (single line) resolution spectra. As one can see, the strength of the narrow Na I D lines at low-resolution visually changes depending on the location in the spectrum; for example, at l₂ and l₃, the line looks weaker. This same procedure of adding the narrow lines at different positions was also done for other observed SNe and for the aforementioned SN Ia and SN II models.

When we then calculate the EW of these simulated narrow lines, we find significant changes and a non-intrinsic evolution, as can be seen for observed and modelled spectra of SNe II in Fig. 15. It is evident that as the P-Cygni profile of Na i D develops with time during the SN photospheric evolution, the narrow line measurement is deeply affected: the EW we measure is lower and drops as a function of the steepness of the underlying continuum. Not only does this affect the EW measured for the narrow line, but it may also affect the photospheric or nebular line from the SN, e.g. producing an apparent redshift in the broad Na I D emission line at late times in CC-SNe (>100 days from explosion; see Fig. C.5). This effect depends on the SN cosmological redshift and ejecta expansion velocity, which determine the location of the narrow line with respect to the P-Cygni profile. At different locations, the inferred behaviour changes accordingly. However, if we do the same simulations for high-resolution narrow lines, we do perfectly recover the EW, no matter the underlying P-Cygni profile.

Similarly, we show in Fig. 16 the results for a simulation of SESN models and observed spectra. We also find that the Na I D P-Cygni profile of the SN ejecta may swallow portions of the narrow line at later epochs in low-resolution spectra and consequently induce smaller measurements of the EW. Regarding SNe-int, there are no available models in the literature, and simulating narrow lines on top of a few objects will misrepresent a very heterogeneous class. Nonetheless, we expect – and visually confirm – the effect to be smaller, since strongly interacting SNe generally do not show prominent P-Cygni profiles.

Although less pronounced than SNe II, SNe Ia can also exhibit an apparent evolution from contamination by the broad ejecta features. This is seen for observed and simulated spectra (see Fig. 17). In this case, the P-Cygni profile of Si II 6355 Å at early times, has very high expansion velocities, and part of it may reach the rest wavelength of the narrow Na I D; later as the ejecta slows down, the narrow line is no longer engulfed by the absorption of the P-Cygni profile, and the EW is well recovered. The exact behaviour naturally depends strongly on the expansion velocity that moves different parts of the broad line into the narrow Na I D, and the effect will be stronger for higher-velocity objects (bottom panel of Fig. 13). Again, this is not seen for high-resolution spectra.

It is thus no surprise that using low-dispersion spectra, Wang et al. (2019) find the evolution of Na I D lines in high-velocity SNe Ia when not properly considering this effect. Their evolution is strikingly similar to our simulated case shown in Fig. 17. Therefore, variations in the strength of Na I D narrow lines in SNe Ia at early times are not always indications of circumstellar interaction. In most cases, when low-resolution spectra are used, variations are produced by the interference of SN broad lines; in this case, by the Si II line.

Having seen the clear important bias of the interference of the broad line in the measurement of the narrow lines, is there a way to identify these dubious cases? Since the bias stems basically from a steep continuum from the P-Cygni profile, we can measure for a given spectrum the average absolute slope of the inferred continuum, m_c, from the automatic method (see Sect. 4.1) by doing: $$m_c=\frac{1}{\langle f_c\rangle }\frac{1}{N}{\displaystyle \sum _i^N\frac{\left|f_c^{i+1}-f_c^i\right|}{{\lambda }_c^{i+1}-{\lambda }_c^i}},$$(6)

where the sum is over all N nodes of the continuum within ±50Å of the line, and 〈ƒ_c〉 is a normalization average of the continuum flux in the same range. The normalized slope continuum m_c[Å⁻¹ ] is a measure of the steepness of the line in the immediate vicinity of the narrow line. Such slopes clearly reflect the presence of underlying broad profiles, as can be seen in Figs. 15 and 17. The discrepancy with the real simulated EW increases for larger continuum slopes, as can also be seen in Fig. 18. We find that a cut of continuum slopes of m_c < 0.002 Å⁻¹ ensures deviations of less than 30% in the EW measurement for all simulated cases. We strongly encourage such a cut in m_c for the EW measurement of narrow lines in mid- and low-resolution spectra.

5.5 Slit size and orientation

The stability of the strength of narrow lines could be affected by different external instrumental factors (besides the resolution and S/N), such as the size and orientation of the slit used for spectroscopy or the amount of host-galaxy light entering into the slit. These factors could introduce changes in the EWs that could be easily confused with intrinsic variations, given the diversity of instruments and observing conditions. To estimate the effect of some of these potential biases in the strength of the narrow lines, we measure the EW of the Na I D lines in different environments and with different simulated slits. For this aim, we use data obtained with MUSE at the Very Large Telescope (VLT) as a part of the All-weather MUse Supernova Integral field Nearby Galaxies (AMUSING, Galbany et al. 2016) collaboration.

To evaluate the effect of the amount of host-galaxy light entering the slit, the impact of the slit orientation and its size, we define six different positions, with four distinct slit angles for each location (0, 45, 90 and 135 degrees with respect to the horizontal west-east axis) and two widths for the slit, one of 1″ and another one of 1.5″. The length of the slit is taken to be 15". Figure 19 shows an example of the procedure for NGC 988, the host galaxy of the type II SN 2017gmr. As one can see, the chosen positions are distributed along the galaxy, from outer regions to the centre. We note that position 1 actually includes the SN 2017gmr in the slit, making this a particularly good study case. We then integrate all spectra within a rectangular aperture of the slit dimension. The resulting profiles centred on the narrow Na I D lines are presented in Fig. 20 for all different slit orientations and positions. We show the profiles obtained using the slit of 1″ width on the left panels and the slit of 1.5″ on the right panels. The average and standard deviations of the Na I DEWs for every position, displayed in each panel, show small differences (the largest deviation for the slit size is ~3% and for orientation is ~13% at position 3). Although positions 2–6 do not enclose a supernova, we expect the effect of just integrated host light to be an upper limit on the EW since the light contribution of the SN will dominate. In fact, position 1 with the SN does present only minor changes of less than 3% with slit size and orientation. Similar results were also found using several other galaxies targeted by AMUSING, particularly the largest error associated with these components is about 0.20 Å, which is generally always lower than the uncertainty estimate of our automated technique and thus not a major concern henceforth.

5.6 Minimum number of spectra

In this section, we evaluate how the automated measurement of the EW changes with the number of spectra used in the stacking procedure of Sect. 4.3. We take non-evolving SNe with more than 20 spectra, passing the previous cuts (in particular, the continuum slope). We then explore the effect of using 1–15 spectra for calculating the EW. For each number of spectra used, we draw 50 random combinations of spectra and calculate the median and standard deviation, as shown in Fig. 21. Besides the dispersion of all draws, we also show the typical median uncertainty obtained from our method with dashed error bars. The bottom plot shows that the uncertainties are conservative when using more than two or three spectra. The standard deviation indicates a high dispersion in the measured values for one and sometimes two spectra, even beyond the measured uncertainties. The bottom panel of Fig. 21 shows that this underestimation can reach a factor of 2– 10 when using only one spectrum. We thus prefer to discard SNe with only one spectrum because of the large dispersion of 1–2 Å of some cases.

Our final set of spectra and SNe is considerably reduced after the four cuts (see Table 2): redshift, signal-to-noise in the spectra, continuum slope, signal-to-noise in the EW and minimum number of spectra. However, we believe these cuts provide a clean sample suitable for the automated measurement of narrow lines.

6 Na I D EW evolution analysis and discussion

The evolution of narrow lines from intervening material has been largely studied in the literature, and variations in the EW have been reported for various transients (e.g. Byrne et al. 2023; de laeger et al. 2015). Such observations have important implications for the presence and possible interaction with circumstellar matter or nearby patchy interstellar clouds. In high-resolution spectra of SNe Ia, SN 2006X (Patat et al. 2007), SN 2007le (Simon et al. 2009), SN 20l3gh (Ferretti et al. 2016), and SN 20141 (Graham et al. 2015) showed changes in the EWs and profiles of the narrow intervening features. In the previous section, we showed with simulations that the EW measurement of narrow lines in high-resolution spectra is not affected by systematic biases; therefore, changes in profiles and intensities presented in these works probably represent true changes in the absorption of material in the line of sight of those SNe. On the other hand, we have shown that variations in low-resolution spectra should be analysed with care. In particular, the effect of the varying morphology of the underlying ejecta P-Cygni profile must be considered. Wang et al. (2019) found variable Na I D absorption lines in low-resolution spectra of high-velocity SNe Ia at early phases. We found a similar behaviour at similar phases in our low-resolution simulation: a decrease in the EW, followed by an increase, finally settling onto a constant evolution. As shown in Fig. 17, the changes in the EW are closely related to the slope of the continuum. When a continuum cut is considered, the variability in most objects disappears. This does not mean that the variability of the narrow absorption lines in low-resolution spectra cannot be real, but we suggest that careful cuts must be considered before performing such analysis and drawing conclusions.

After including the cuts described in previous sections (see summary in Table 2), we evaluate the behaviour of the strength of the narrow Na I D absorption lines as a function of time. Since EW measurements from an individual spectrum can be unreliable (Sect. 5.6), we require at least two spectra per epoch range for a given SN; we then stack them and recalculate the EW (see Sect. 4.3). For this analysis, we explore various possibilities, using at least two or at least three spectra per SN in intervals of 10 to 20 days. While stacking SN spectra in time intervals will wash out shorter timescale variations, it ensures that we can study the time evolution of the EW with our methodology and low-resolution spectra. Once the spectra are stacked, we look at the evolution of the narrow lines both statistically and individually. Particularly, we focus on the EW of the line because it is more robust than the velocity (see Sect. 4.2).

6.1 Statistical evolution of EW

We study the main SN types, SNe Ia, SNe II, SESNe and SNe-int, for statistical signs of evolution. Finding population-wide signs of evolution could give us important insight into the progenitor systems and nearby material of a given population. For instance, a large number of SNe II are known to be subject to very early interaction with CSM (e.g. Bostroem et al. 2023; Martinez et al. 2024) that could perhaps be revealed in the narrow absorption line evolution. SNe-int are by definition subject to interaction with nearby material and thus good candidates to observe evolution in the narrow absorption lines, as is the case for SN 2011A (de Jaeger et al. 2015); however, the line is often absent despite strong CSM interaction(e.g. Ofek et al. 2010; Pastorello et al. 2018). To investigate this, we require a common reference epoch for each group. We use the maximum B-band dates obtained through multi-wavelength light-curve fits with spectrophotomet-ric templates of each SN type (see Appendix B). The epochs are normalized to this reference epoch and put in rest-frame. All spectra of a given SN in time intervals of, e.g. 10 days, are stacked, and the EW is measured. We investigate the evolution of the σ-deviations per epoch interval i with respect to the all-epoch stacked EW measurement of each SN, i.e. $\overline{\text{EW}}$, so that: $$\delta {\sigma }_i=\frac{{\text{EW}}_i-\overline{\text{EW}}}{\sqrt{{\sigma }_{{\text{EW}}_i}^2+{\sigma }_{\overline{\text{EW}}}^2}}.$$(7)

Figure 22 shows the evolution of the distribution of these deviations for the Na I D lines of the four SN types considered. In general, we can see that there is no large evolution of these groups within the uncertainties. For SESNe, there seems to be a decrease after 50 days post maximum, although this is likely due to low-number statistics. Similarly, SNe Ia appear to increase within 30 and 50 days post maximum, but with a small number of SNe at these epochs, this is mainly driven by one SN (SN 2012dn), as we will see in the next section. SNe-int also have very low statistics making any conclusion rather difficult. When investigating the evolution with different epoch intervals between 10–15 days and also requiring a minimum of three spectra per SN per interval instead of two, we obtain similar results.

Additionally, we also investigate narrow lines from other species, such as Ca II and K I finding no evidence of evolution. However, it is possible that smaller sub-groups might show signs of changes in time. Given the small statistics after all our cuts, we leave this for a future investigation. In the next section, we analyse the evolution of individual well-sampled SNe.

6.2 Individual evolution of EW

To analyse the evolution of the EW of narrow absorption lines for individual objects, we take the stacked EW of a minimum of two spectra in intervals of 10 days, as in the previous section. Besides a visual inspection, we apply two methods to identify possible candidates for evolution automatically:

• Linear slope: we fitted the EW of all intervals of each SN to a simple line with a Bayesian fit from which the posterior distribution gives us the probability for the slope to be different from zero (Kelly 2007). A significantly non-zero slope (p > 95%) indicates that there is a general monotonic trend in the evolution of the SN.

• Cumulative sum: to try to capture more complicated evolution, we also do a cumulative sum analysis (CUSUM, Page 1954) that monitors deviations from a target, in this case, the EW value obtained from the stack of all the spectra (at all epochs) of a given SN, $\overline{\text{EW}}$. Then the algorithm measures accumulated deviations, chosen here as a drift of half the uncertainty on the stacked EW, $0.5{\sigma }_{\overline{\text{EW}}}$, to both positive and negative directions, so that: $$\begin{array}{l}{H}_0=0\\ L_0=0\\ H_i=\max \left[0,H_{i-1}+\left({\text{EW}}_i-{\sigma }_{{\text{EW}}_i}\right)-\left(\overline{\text{EW}}+0.5{\sigma }_{\overline{\text{EW}}}\right)\right]\\ L_i=\min \left[0,L_{i-1}+\left({\text{EW}}_i+{\sigma }_{{\text{EW}}_i}\right)-\left(\overline{\text{EW}}-0.5{\sigma }_{\overline{\text{EW}}}\right)\right]\end{array}$$

  where H_i and L_i indicate the positive and negative cumulative deviations updated at each time step i. We also take into account the uncertainty in each individual EW_i. This measures the cumulative deviations starting from the earliest epoch interval of the SN and increasing in time. In order to not give higher weight to the initial time bins, we also start the algorithm from the latest epoch interval going in decreasing time steps. The final cumulative H and L are the average of both sequential cumulative deviations. SNe with more than two cumulative deviations above the error on the stacked EW are considered viable candidates.

Figure 23 shows these two methods applied to the measurements of SN 2006X (SN Ia) and SN 2014G (SN II) after the cuts (filled black points). Doing a Bayesian linear fit to the stacked EWs per epoch interval (red points) results in a slope consistent with zero for SN 2006X and no significant points in the cumulative deviation analysis. SN 2014G, on the other hand, has a slope with a probability of 93.7% to be non-zero and one significant cumulative deviation point. Both of these SNe are considered to be not evolving in time according to our diagnostics. However, it is important to observe that if we were to do our analysis on all data points without cuts (empty black points), these SNe would have non-zero slope probabilities of 88.4 and 100%, with 4 and 12 significant cumulative deviation points, respectively. Both SNe would be considered to have an evolution. This shows the importance of our cuts, particularly the continuum slope that prevents strong blending with the broad P-Cygni profile (see Sect. 5.4).

When analysing various narrow lines for all SNe after cuts and stacking at least two spectra per 10 or 15-day intervals, we find that very few SNe are tagged as possibly evolving according to our two diagnostics. For every case, we ensure that no equivalent evolution is seen for the same lines of the MW and inspect them visually. The few SNe tagged (e.g. SN 2001N, SN 2012dn) have few points in the nebular phase with stronger EW. However, these are a very small number and, in some cases, still noisy after stacking per interval. In the case of SN 2012dn, a 03fg-like “super-Chandra” SN Ia (Chakradhari et al. 2014; Taubenberger et al. 2019), the narrow Na I D at ~40 days after the maximum gets mixed with features due to iron-group elements and it is difficult to disentangle the two.

SNe known to show variation in their lines with high-resolution spectra, e.g. SN 2006X, SN 2007le, SN 2013gh or SN 2014J, are not picked up by our culls, but we have predominantly low-resolution spectra for which such small variations remain within the uncertainties. We have no access to high-resolution spectra for SN 2006X, SN 2013gh and only one for SN 2007le. For SN 2014J, if we restrict our analysis to the 10 high-resolution spectra we have and do not make any cuts on the continuum slope, we find no evolution of the EW for Na I D and no evolution of the three spectra that cover the wavelengths of K I. Graham et al. (2015) found evolution in K I lines in other spectra of even higher resolution (R > 100 000). In the case of the type IIn SN 2011A (de Jaeger et al. 2015), we do see the increase in EW with time that they report. Nonetheless, with only 5 spectra after cuts (and only two stacked points per 15-day interval), the linear slope is not significantly positive, and the stacked EW of all spectra is too influenced by the late spectra with increasing EW, such that no significant deviation points are found.

We conclude that the current dataset, after all our cuts and with our methodology, is consistent with no variation in the EW of narrow lines over time within the uncertainties. Subtler variations are naturally possible but necessitate better sampling and, foremost, higher spectral resolution.

7 Conclusions

In this paper, we have presented a statistical analysis of the evolution of narrow line-of-sight lines from intervening material (Na I D, Ca II H&K, K I, and DIBs) visible in nearby SNe. We have analysed over 11 000 spectra of more than 1600 SNe with various resolutions. To quantify the properties of these narrow lines, we developed a robust automated methodology allowing for the estimation of their EWs and velocities. When searching for possible EW evolution, we found systematic variations in the Na I D line from the MW and the host galaxies. Changes in the Na I D depth from the MW narrow interstellar lines are in principle not expected because the material in our own Galaxy is at sufficiently large distances from the SN. The strength of the narrow absorption lines in low- and mid-resolution spectra may thus be affected by different external factors that we explore: the signal-to-noise of the spectra, the spectral resolution, the slit size and orientation, the number of spectra used, and, most importantly, the interference of the broad P-Cygni lines from the SN. We found that these biases, particularly those associated with the presence of the broad lines, induce a conspicuous apparent evolution (variability) over time in many SNe. For instance, in SNe II, the narrow sodium line overlaps with the broad sodium from the SN ejecta, and in SNe Ia the narrow sodium overlaps with the P-Cygni profile of Si II. We do indeed confirm this with simulations in which we input fake narrow lines into different locations of the broad SN P-Cygni lines and calculated the recovered EW. In our simulations of low-resolution spectra, we found large significant changes in the EW depending on the underlying profile. However, in our high-resolution simulations, the strength of the narrow lines was not affected by the interference of the broad SN P-Cygni profiles.

Based on these findings, we presented an easy way to detect and exclude cases of high contamination of the SN profile in the narrow line: when the normalized slope of the continuum is too high, that is, m_c > 0.002 Å⁻¹, there is a broad component underneath blending with the narrow line. We show that excluding these cases with a steep continuum is essential in the study of narrow lines from intervening material in the line of sight. Not taking this into account may lead to important biases in the frequently used estimations of dust extinction through the EW ∝ A_V relation (Eq. (3)) in low-resolution spectra, and it may lead to wrong inferences of the evolution of the lines and the presence of nearby material. Previous claims of the evolution of narrow lines in low-resolution spectra should be revisited in light of this study.

We also strongly encourage the use of more than one low-resolution spectrum in the estimation of the EW. Fluctuations of the S/N, the spectral resolution, and the morphology of the underlying P-Cygni profile over time can generate changes in these lines that are washed out with a larger number of spectra stacked to obtain a more accurate EW. We show that at least two or three spectra are needed to measure robust EWs with conservative uncertainties.

Finally, after considering all of these effects, we analysed the temporal evolution of line-of-sight features in a large sample of nearby SNe to detect any possible statistical variation in their EWs. We found no overall evolution for three principal SN groups: SNe II, SNe Ia, and SESNe, in all narrow lines considered. Significant individual SN variations were also ruled out based on our methodology and low-resolution spectra.

In this work, we have paved the way for the robust measurement of narrow absorption lines in a variety of spectra from different instruments and resolutions. Armed with this methodology, we can now use a large sample to study the variations of the SN-averaged narrow lines according to the SN type, SN properties, and host galaxy characteristics. A wealth of information on the pervading intervening material is contained in those lines, providing essential clues as to the immediate and farther vicinity of SNe and their progenitors, with far-reaching consequences for our understanding of the ISM cycle of galaxies and the persistent question of extinction in supernova cosmology.

Acknowledgements

We sincerely thank the referee for very useful comments that improved the quality of the paper. We heartily thank all the research groups and observers for making their data public. This research has made use of the Berkeley SNDB database as well as the CfA Supernova Archive, which is funded in part by the National Science Foundation through grant AST 0907903. S.G.G. acknowledges the financial support from the visitor and mobility program of the Finnish Centre for Astronomy with ESO (FINCA), funded by the Academy of Finland grant nr 306531. S.G.G. and A.M.M. acknowledge support by FCT under Project CRISP PTDC/FIS-AST-31546/2017 and Project No. UIDB/00099/2020. C.P.G. acknowledges financial support from the Secretary of Universities and Research (Government of Catalonia) and by the Horizon 2020 Research and Innovation Programme of the European Union under the Marie Skłodowska-Curie and the Beatriu de Pinós 2021 BP 00168 programme. C.P.G. and L.G. recognize the support from the Spanish Ministerio de Ciencia e Innovación (MCIN) and the Agencia Estatal de Investigación (AEI) 10.13039/501100011033 under 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. L.G. also thanks the financial support from the European Social Fund (ESF) “Investing in your future” under the 2019 Ramón y Cajal program RYC2019-027683-I. J.P.A. acknowledges support by ANID, Millennium Science Initiative, ICN12_009. A.M.G. acknowledges financial support by the European Union under the 2014-2020 ERDF Operational Programme and by the Department of Economic Transformation, Industry, Knowledge, and Universities of the Regional Government of Andalusia through the FEDER-UCA18-107404 grant. S.M. acknowledges support from the Research Council of Finland project 350458.

Appendix A Wavelength calibration uncertainty

The heterogeneous dataset used in our analysis has been wavelength-calibrated by different groups for different instruments and observing conditions. Here we estimate an average uncertainty of the wavelength calibration for the entire sample. For this, we randomly take 1556 spectra of 180 core-collapse SNe from our sample and use our automated technique (see §4.2) around three narrow emission lines characteristic of gas-rich environments: Hα (6563Å), Hβ (4681Å) and [OIII] (5007Å). We use the following integration windows: ± 750 km/s, ± 750 km/s and ± 900 km/s, respectively. An example measurement of Hα for a spectrum of SN 2009bz is shown in Figure A.1. The dispersion of the three line velocities for this spectrum gives 33.2 km s⁻¹. Of the 1556 spectra, we select 696 that clearly show at least two of the three lines, i.e. requiring that EW/σ_EW > 1.3 and EW > 0.5Å; and we calculate the median dispersion for the sample finding 35.0 km/s. Doing the same analysis using Gaussian fits instead, we obtain proper fits with EW> 0.5 in at least two lines for 767 spectra for which the average of the velocity dispersion gives 36.6 km/s. This is in remarkable agreement with the automated value confirming that this dispersion is not dependent on the discussed methodologies. This value could reveal calibration uncertainties or real physical differences in the elements of the gas clouds, or a combination of both. We conservatively take it here as an intrinsic calibration uncertainty that is added in quadrature to other uncertainties of the method.

Appendix B B-band maximum dates

To obtain the B-band maximum reference epochs used in section 6, we perform multi-band light-curve fits to every SN. We use a light-curve fitter commonly used for SNe Ia (Conley et al. 2008) but with different spectrophotometric templates for each type: the template by Hsiao et al. (2007) for SNe Ia, the template of Nugent et al. (2002) for SESNe and a model of SNe II used in the PLAsTiCC challenge (Kessler et al. 2019) based on empirical spectral energy distributions that are a linear combination of three ‘eigenvector’ components. For SNe-int, we try both, the template for SNe II and for SESNe, choosing visually the best obtained fit. Example fits are shown in Figure B.1. All fits are visually inspected, and SNe with poor data or bad fits are discarded. Although core-collapse SNe are not as homogeneous as SNe Ia, and this fitter was devised for SNe Ia, this simple approach gets rough estimates that allow us to have epochs in the same reference system. We also attempted to use the explosion epochs as a reference for SNe II (Gutiérrez et al. 2017), confirming no general evolution of EW with time in all three SN types.

Appendix C Figures

We present here additional figures: Figure C.1 shows the comparison of the EW measured from the stacked flux-to-continuum ratios of several spectra of a given SN (see §4.3) versus the median and the weighted average of spectra for a given SN (the uncertainties come from the median absolute deviation and the weighted error, respectively); Figure C.2 displays the impact of the continuum node separation on the EW measurement (see §4.1); Figures C.3 and C.4 shows how degrading the spectral resolution (§5.2) and the signal-to-noise (§5.3) affect the narrow line and its EW measurement; and Figure C.5 demonstrates that the narrow line at nebular times can emulate a blueshift in the line from the SN ejecta (see §5.4).

Appendix D Tables

In Table D.1 we present all SNe used in this analysis with their type, redshift, number of spectra and the references for their spectra.

Additional references

Ailawadhi, B., Dastidar, R., Misra, K., et al. 2023, MNRAS, 519, 248

Arbour, R., Li, W., Mostardi, R., et al. 2007, IAU Circ., 8829, 2

Arcavi, I., Gal-Yam, A., Yaron, O., et al. 2011, ApJ, 742, L18

Aretxaga, I., Benetti, S., Terlevich, R. J., et al. 1999, MNRAS, 309, 343

Barbon, R., Buondí, V., Cappellaro, E., & Turatto, M. 1999, A&AS, 139, 531

Benetti, S. 1991, in European Southern Observatory Conference and Workshop Proceedings, Vol. 37, European Southern Observatory Conference and Workshop Proceedings, 339

Benetti, S., Cappellaro, E., Danziger, I. J., et al. 1998, MNRAS, 294, 448

Benetti, S., Cappellaro, E., Turatto, M., et al. 1994, A&A, 285, 147

Benetti, S., Chugai, N. N., Utrobin, V. P., et al. 2016, MNRAS, 456, 3296

Benetti, S., Turatto, M., Balberg, S., et al. 2001, MNRAS, 322, 361

Benetti, S., Turatto, M., Cappellaro, E., Danziger, I. J., & Mazzali, P. A. 1999, MNRAS, 305, 811

Benetti, S., Turatto, M., Valenti, S., et al. 2011, MNRAS, 411, 2726

Black, C. S., Milisavljevic, D., Margutti, R., et al. 2017, ApJ, 848, 5

Blanco, V. M., Gregory, B., Hamuy, M., et al. 1987, ApJ, 320, 589

Blanton, E. L., Schmidt, B. P., Kirshner, R. P., et al. 1995, AJ, 110, 2868

Blaylock, M., Branch, D., Casebeer, D., et al. 2000, PASP, 112, 1439

Blondin, S., Challis, P., & Kirshner, R. 2008, Central Bureau Electronic Telegrams, 1454, 1

Borish, H. J., Huang, C., Chevalier, R. A., et al. 2015, ApJ, 801, 7

Botticella, M. T. et al. 2009, MNRAS, 398, 1041

Branch, D., Falk, S. W., Uomoto, A. K., et al. 1981, ApJ, 244, 780

Branch, D. & Venkatakrishna, K. L. 1986, ApJ, 306, L21

Bufano, F., Pignata, G., Bersten, M., et al. 2014, MNRAS

Bufano, F. et al. 2012, ApJ, 753, 67

Cappellaro, E., Danziger, I. J., della Valle, M., Gouiffes, C., & Turatto, M. 1995, A&A, 293, 723

Casebeer, D., Branch, D., Blaylock, M., et al. 2000, PASP, 112, 1433

Childress, M. J., Hillier, D. J., Seitenzahl, I., et al. 2015, MNRAS, 454, 3816

Childress, M. J., Tucker, B. E., Yuan, F., et al. 2016, PASA, 33, e055

Ciabattari, F., Mazzoni, E., Rossi, M., et al. 2013, Central Bureau Electronic Telegrams, 3576, 1

Clocchiatti, A., Benetti, S., Wheeler, J. C., et al. 1996, AJ, 111, 1286

Corsi, A., Ofek, E. O., Gal-Yam, A., et al. 2012, ApJ, 747, L5

Cristiani, S., Cappellaro, E., Turatto, M., et al. 1992, A&A, 259, 63

Davis, K., Johnson, J., Tinyanont, S., et al. 2022, Transient Name Server Classification Report, 2022-362, 1

de Jaeger, T., Anderson, J. P., Pignata, G., et al. 2015, ApJ, 807, 63

de Jaeger, T., Zheng, W., Stahl, B. E., et al. 2019, MNRAS, 490, 2799

D’Elia, V., Pian, E., Melandri, A., et al. 2015, A&A, 577, A116

Deng, J. S., Qiu, Y. L., Hu, J. Y., Hatano, K., & Branch, D. 2000, ApJ, 540, 452

DerKacy, J. M., Paugh, S., Baron, E., et al. 2023, MNRAS, 522, 3481

Drake, A. J., Djorgovski, S. G., Graham, M. J., et al. 2013, Central Bureau Electronic Telegrams, 3494, 1

Drout, M. R., Milisavljevic, D., Parrent, J., et al. 2016, ApJ, 821, 57

Drout, M. R., Soderberg, A. M., Mazzali, P. A., et al. 2013, ApJ, 774, 58

Dugas, A., Fremling, C., & Sharma, Y. 2019, Transient Name Server Classification Report, 2019-2840, 1

Elias-Rosa, N., Pastorello, A., Maund, J. R., et al. 2013, MNRAS, 436, L109

Elias-Rosa, N., Van Dyk, S. D., Li, W., et al. 2010, ApJ, 714, L254

Ergon, M., Sollerman, J., Fraser, M., et al. 2014, A&A, 562, A17

Faran, T., Poznanski, D., Filippenko, A. V., et al. 2014a, MNRAS, 442, 844

Faran, T., Poznanski, D., Filippenko, A. V., et al. 2014b, MNRAS, 445, 554

Fassia, A., Meikle, W. P. S., Chugai, N., et al. 2001, MNRAS, 325, 907

Filippenko, A. V. 1988, AJ, 96, 1941

Filippenko, A. V. 1997, ARA&A, 35, 309

Filippenko, A. V., Porter, A. C., & Sargent, W. L. W. 1990, AJ, 100, 1575

Filippenko, A. V., Richmond, M. W., Branch, D., et al. 1992a, AJ, 104, 1543

Filippenko, A. V., Richmond, M. W., Matheson, T., et al. 1992b, ApJ, 384, L15

Filippenko, A. V. & Sargent, W. L. W. 1986, AJ, 91, 691

Foley, R. J. 2015, MNRAS, 452, 2463

Foley, R. J., Filippenko, A. V., & Jha, S. W. 2008, ApJ, 686, 117

Foley, R. J., Jha, S. W., Pan, Y.-C., et al. 2016a, MNRAS, 461, 433

Foley, R. J., Pan, Y.-C., Brown, P., et al. 2016b, MNRAS, 461, 1308

Foley, R. J., Smith, N., Ganeshalingam, M., et al. 2007, ApJ, 657, L105

Foley, R. J., Van Dyk, S. D., Jha, S. W., et al. 2015, ApJ, 798, L37

Foley, R. J. et al. 2012, ApJ, 744, 38

Foley, R. J. et al. 2013, ApJ, 767, 57

Fraser, M., Ergon, M., Eldridge, J. J., et al. 2011, MNRAS, 417, 1417

Fraser, M. et al. 2013, MNRAS, 433, 1312

Gagliano, R., Post, R., Weinberg, E., Newton, J., & Puckett, T. 2016, Transient Name Server Discovery Report, 2016-568, 1

Gal-Yam, A., Arcavi, I., Ofek, E. O., et al. 2014, Nature, 509, 471

Gal-Yam, A., Bruch, R., Schulze, S., et al. 2022, Nature, 601, 201

Gal-Yam, A., Bufano, F., Barlow, T. A., et al. 2008, ApJ, 685, L117

Gal-Yam, A., Kasliwal, M. M., Arcavi, I., et al. 2011, ApJ, 736, 159

Gal-Yam, A., Leonard, D. C., Fox, D. B., et al. 2007, ApJ, 656, 372

Galbany, L., Moreno-Raya, M. E., Ruiz-Lapuente, P., et al. 2016, MNRAS, 457, 525

Gerardy, C. L., Meikle, W. P. S., Kotak, R., et al. 2007, ApJ, 661, 995

Gomez, G. & Lopez, R. 1995, AJ, 109, 737

Gómez, G. & López, R. 1998, AJ, 115, 1096

Gómez, G. & López, R. 2000, AJ, 120, 367

Gomez, G., Lopez, R., & Sanchez, F. 1996, AJ, 112, 2094

Gorbikov, E., Gal-Yam, A., Ofek, E. O., et al. 2014, MNRAS, 443, 671

Graham, M. L., Kumar, S., Hosseinzadeh, G., et al. 2017, MNRAS, 472, 3437

Graham, M. L., Sand, D. J., Valenti, S., et al. 2014, ApJ, 787, 163

Green, D. W. E. 2007, Central Bureau Electronic Telegrams, 974, 2

Gutiérrez, C. P., Bersten, M. C., Orellana, M., et al. 2021, MNRAS, 504, 4907

Gutierrez, C. P., Cartier, R., & Yaron, O. 2017, Transient Name Server Classification Report, 2017-468, 1

Gutiérrez, C. P., González-Gaitán, S., Folatelli, G., et al. 2016, A&A, 590, A5

Hamuy, M., Maza, J., Pinto, P. A., et al. 2002, AJ, 124, 417

Han, X., Zheng, W., Stahl, B. E., et al. 2020, ApJ, 892, 142

Harutyunyan, A. H. et al. 2008, A&A, 488, 383

Hendry, M. A., Smartt, S. J., Crockett, R. M., et al. 2006, MNRAS, 369, 1303

Hicken, M., Friedman, A. S., Blondin, S., et al. 2017, ApJS, 233, 6

Hiramatsu, D., Arcavi, I., Burke, J., et al. 2019a, Transient Name Server Classification Report, 2019-1137, 1

Hiramatsu, D., Burke, J., Arcavi, I., et al. 2019b, Transient Name Server Classification Report, 2019-738, 1

Hosseinzadeh, G., Kilpatrick, C. D., Dong, Y., et al. 2022, ApJ, 935, 31

Hosseinzadeh, G., Sand, D. J., Valenti, S., et al. 2017a, ApJ, 845, L11

Hosseinzadeh, G., Valenti, S., Arcavi, I., et al. 2017b, Transient Name Server Classification Report, 2017-1514, 1

Hosseinzadeh, G., Valenti, S., McCully, C., et al. 2018, ApJ, 861, 63

Hosseinzadeh, G. et al. 2017c, ApJ, 836, 158

Humphreys, R. M., Davidson, K., Jones, T. J., et al. 2012, ApJ, 760, 93

Inserra, C., Pastorello, A., Turatto, M., et al. 2013, A&A, 555, A142

Inserra, C., Turatto, M., Pastorello, A., et al. 2012, MNRAS, 422, 1122

Irani, I., Schulze, S., Gal-Yam, A., et al. 2019, ApJ, 887, 127

Jacobson-Galán, W. V., Venkatraman, P., Margutti, R., et al. 2022, ApJ, 932, 58

Jencson, J. E., Prieto, J. L., Kochanek, C. S., et al. 2016, MNRAS, 456, 2622

Jha, S. & Brown, W. 2000, IAU Circ., 7394, 2

Jha, S., Challis, P., Kirshner, R., & Calkins, M. 2000, IAU Circ., 7379, 1

Jha, S., Garnavich, P., Challis, P., et al. 1999, IAU Circ., 7107, 2

Jones, D. O., Dimitriadis, G., Pan, Y. C., & Foley, R. J. 2017, Transient Name Server Classification Report, 2017-1464, 1

Kangas, T. et al. 2016, MNRAS, 456, 323

Karamehmetoglu, E., Fransson, C., Sollerman, J., et al. 2021, A&A, 649, A163

Khazov, D., Yaron, O., Gal-Yam, A., et al. 2016, ApJ, 818, 3

Kiewe, M., Gal-Yam, A., Arcavi, I., et al. 2012, ApJ, 744, 10

Kilpatrick, C. D., Foley, R. J., Abramson, L. E., et al. 2017, MNRAS, 465, 4650

Kosugi, G., Mizumoto, Y., Kawai, N., et al. 2004, PASJ, 56, 61

Kumar, B. et al. 2013, MNRAS, 431, 308

Kumar, S., Hsiao, E. Y., Ashall, C., et al. 2023, ApJ, 945, 27

Leonard, D. C., Filippenko, A. V., Barth, A. J., & Matheson, T. 2000, ApJ, 536, 239

Leonard, D. C., Filippenko, A. V., Ganeshalingam, M., et al. 2006, Nature, 440, 505

Leonard, D. C., Filippenko, A. V., Li, W., et al. 2002, AJ, 124, 2490

Li, L., Wang, X., Zhang, J., et al. 2018, MNRAS, 478, 4575

Lin, H., Xiang, D., Rui, L., et al. 2018, The Astronomer’s Telegram, 11205, 1

Lunnan, R., Kasliwal, M. M., Cao, Y., et al. 2017, ApJ, 836, 60

Maeda, K., Tajitsu, A., Kawabata, K. S., et al. 2016, ApJ, 816, 57

Maguire, K., Jerkstrand, A., Smartt, S. J., et al. 2012, MNRAS, 420, 3451

Maguire, K., Sullivan, M., Pan, Y.-C., et al. 2014, MNRAS, 444, 3258

Maguire, K., Sullivan, M., Patat, F., et al. 2013, MNRAS, 436, 222

Mahabal, A. A., Drake, A. J., Djorgovski, S. G., et al. 2012, Central Bureau Electronic Telegrams, 3297, 1

Margutti, R., Milisavljevic, D., Soderberg, A. M., et al. 2014, ApJ, 780, 21

Marion, G. H., Brown, P. J., Vinkó, J., et al. 2016, ApJ, 820, 92

Matheson, T., Challis, P., Kirshner, R., & Calkins, M. 2003a, IAU Circ., 8203, 1

Matheson, T., Challis, P., Kirshner, R., & Calkins, M. 2003b, IAU Circ., 8228, 1

Matheson, T., Filippenko, A. V., Chornock, R., Leonard, D. C., & Li, W. 2000, AJ, 119, 2303

Matheson, T., Filippenko, A. V., Li, W., Leonard, D. C., & Shields, J. C. 2001, AJ, 121, 1648

Matheson, T., Garnavich, P. M., Stanek, K. Z., et al. 2003c, ApJ, 599, 394

Maund, J. R., Fraser, M., Smartt, S. J., et al. 2013, MNRAS, 431, L102

Mazzali, P. A., Danziger, I. J., & Turatto, M. 1995, A&A, 297, 509

Mazzali, P. A., Sullivan, M., Filippenko, A. V., et al. 2015, MNRAS, 450, 2631

Mazzali, P. A., Sullivan, M., Hachinger, S., et al. 2014, MNRAS, 439, 1959

Meikle, W. P. S., Kotak, R., Farrah, D., et al. 2011, ApJ, 732, 109

Melandri, A., Pian, E., Ferrero, P., et al. 2012, A&A, 547, A82

Milisavljevic, D., Margutti, R., Parrent, J. T., et al. 2015, ApJ, 799, 51

Milisavljevic, D., Soderberg, A. M., Margutti, R., et al. 2013, ApJ, 770, L38

Miller, A. A., Magee, M. R., Polin, A., et al. 2020, ApJ, 898, 56

Miller, A. A. et al. 2010, MNRAS, 404, 305

Modjaz, M., Kirshner, R., Challis, P., Matheson, T., & Nutzman, P. 2005, Central Bureau Electronic Telegrams, 131, 1

Modjaz, M., Kirshner, R. P., Blondin, S., Challis, P., & Matheson, T. 2008, ApJ, 687, L9

Morales-Garoffolo, A., Elias-Rosa, N., Bersten, M., et al. 2015, MNRAS, 454, 95

Nagao, T., Kuncarayakti, H., Maeda, K., et al. 2023, A&A, 673, A27

Neill, J. 2016, Transient Name Server Classification Report, 2016-951, 1

Nordin, J., Brinnel, V., Giomi, M., et al. 2019, Transient Name Server Discovery Report, 2019-1287, 1

Nugent, P. E. et al. 2011, Nature, 480, 344

Nyholm, A., Sollerman, J., Taddia, F., et al. 2017, A&A, 605, A6

Oates, S. R., Bayless, A. J., Stritzinger, M. D., et al. 2012, MNRAS, 424, 1297

Ochner, P., Benetti, S., Cappellaro, E., et al. 2015, The Astronomer’s Telegram, 8205, 1

Ofek, E. O., Cameron, P. B., Kasliwal, M. M., et al. 2007, ApJ, 659, L13

Östman, L., Nordin, J., Goobar, A., et al. 2011, A&A, 526, A28

Pan, Y. C., Foley, R. J., Kromer, M., et al. 2015, MNRAS, 452, 4307

Parrent, J. T., Howell, D. A., Fesen, R. A., et al. 2016, MNRAS, 457, 3702

Parrent, J. T., Howell, D. A., Friesen, B., et al. 2012, ApJ, 752, L26

Pastorello, A., Baron, E., Branch, D., et al. 2005, MNRAS, 360, 950

Pastorello, A., Benetti, S., Brown, P. J., et al. 2015a, MNRAS, 449, 1921

Pastorello, A., Hadjiyska, E., Rabinowitz, D., et al. 2015b, MNRAS, 449, 1954

Pastorello, A., Prieto, J. L., Elias-Rosa, N., et al. 2015c, MNRAS, 453, 3649

Pastorello, A., Pumo, M. L., Navasardyan, H., et al. 2012, A&A, 537, A141

Pastorello, A., Turatto, M., Benetti, S., et al. 2002, MNRAS, 333, 27

Pastorello, A., Valenti, S., Zampieri, L., et al. 2009, MNRAS, 394, 2266

Pastorello, A., Wyrzykowski, Ł., Valenti, S., et al. 2015d, MNRAS, 449, 1941

Pastorello, A., Zampieri, L., Turatto, M., et al. 2004, MNRAS, 347, 74

Pastorello, A. et al. 2008, MNRAS, 389, 113

Pastorello, A. et al. 2015e, MNRAS, 454, 4293

Patat, F., Benetti, S., Cappellaro, E., et al. 1996, MNRAS, 278, 111

Patat, F., Cappellaro, E., Danziger, J., et al. 2001, ApJ, 555, 900

Pellegrino, C., Hiramatsu, D., Arcavi, I., et al. 2023, ApJ, 954, 35

Pellegrino, C., Howell, D. A., Terreran, G., et al. 2022, ApJ, 938, 73

Pereira, R., Thomas, R. C., Aldering, G., et al. 2013, A&A, 554, A27

Pessi, T., Prieto, J. L., Monard, B., et al. 2022, ApJ, 928, 138

Phillips, M. M., Heathcote, S. R., Hamuy, M., & Navarrete, M. 1988, AJ, 95, 1087

Pignata, G. et al. 2011, ApJ, 728, 14

Poidevin, F., Perez-Fournon, I., Angel, C. J., et al. 2020, Transient Name Server Classification Report, 2020-2512, 1

Polshaw, J., Kotak, R., Chambers, K. C., et al. 2015, A&A, 580, L15

Poznanski, D., Kostrzewa-Rutkowska, Z., Wyrzykowski, L., & Blagorodnova, N. 2015, MNRAS, 449, 1753

Pozzo, M., Meikle, W. P. S., Fassia, A., et al. 2004, MNRAS, 352, 457

Prentice, S. J., Ashall, C., James, P. A., et al. 2019a, MNRAS, 485, 1559

Prentice, S. J., Maguire, K., Skillen, K., Magee, M. R., & Clark, P. 2019b, Transient Name Server Classification Report, 2019-1306, 1

Pun, C. S. J., Kirshner, R. P., Sonneborn, G., et al. 1995, ApJS, 99, 223

Quimby, R., Höflich, P., Kannappan, S. J., et al. 2006, ApJ, 636, 400

Quimby, R. M., Wheeler, J. C., Höflich, P., et al. 2007, ApJ, 666, 1093

Ransome, C. 2021, Transient Name Server Classification Report, 2021-2830, 1

Reguitti, A., Pastorello, A., Pignata, G., et al. 2019, MNRAS, 482, 2750

Reguitti, A., Pastorello, A., Pignata, G., et al. 2022, A&A, 662, L10

Reynolds, T. M., Fraser, M., Mattila, S., et al. 2020, MNRAS, 493, 1761

Rich, D., Dennefeld, M., Pizzella, A., et al. 2012, Central Bureau Electronic Telegrams, 3227, 2

Sanders, N. E., Soderberg, A. M., Valenti, S., et al. 2012, ApJ, 756, 184

Schmidt, B. P., Kirshner, R. P., Schild, R., et al. 1993, AJ, 105, 2236

Schulze, S., Yaron, O., Sollerman, J., et al. 2021, ApJS, 255, 29

Shivvers, I., Mazzali, P., Silverman, J. M., et al. 2013, MNRAS, 436, 3614

Shivvers, I., Modjaz, M., Zheng, W., et al. 2017a, PASP, 129, 054201

Shivvers, I. et al. 2016, MNRAS, 461, 3057

Shivvers, I. et al. 2017b, MNRAS, 471, 4381

Silverman, J. M., Ganeshalingam, M., Cenko, S. B., et al. 2012, ApJ, 756, L7

Silverman, J. M., Ganeshalingam, M., & Filippenko, A. V. 2013a, MNRAS, 430, 1030

Silverman, J. M., Ganeshalingam, M., Li, W., et al. 2011, MNRAS, 410, 585

Silverman, J. M., Nugent, P. E., Gal-Yam, A., et al. 2013b, ApJS, 207, 3

Silverman, J. M., Pickett, S., Wheeler, J. C., et al. 2017, MNRAS, 467, 369

Silverman, J. M. et al. 2013c, ApJS, 207, 3

Singh, A., Srivastav, S., Kumar, B., Anupama, G. C., & Sahu, D. K. 2018, MNRAS, 480, 2475

Smith, N., Cenko, S. B., Butler, N., et al. 2012a, MNRAS, 420, 1135

Smith, N., Chornock, R., Silverman, J. M., Filippenko, A. V., & Foley, R. J. 2010, ApJ, 709, 856

Smith, N., Li, W., Silverman, J. M., Ganeshalingam, M., & Filippenko, A. V. 2011, MNRAS, 415,773

Smith, N., Mauerhan, J. C., Silverman, J. M., et al. 2012b, MNRAS, 426, 1905

Smith, N., Silverman, J. M., Filippenko, A. V., et al. 2012c, AJ, 143, 17

Sparre, M., Sollerman, J., Fynbo, J. P. U., et al. 2011, ApJ, 735, L24

Spiro, S., Pastorello, A., Pumo, M. L., et al. 2014, MNRAS, 439, 2873

Srivastav, S., Ninan, J. P., Kumar, B., et al. 2016, MNRAS, 457, 1000

Srivastav, S., Smartt, S. J., Huber, M. E., et al. 2023, ApJ, 943, L20

Stahl, B. E. et al. 2020, MNRAS, 492, 4325

Sternberg, A., Gal-Yam, A., Simon, J. D., et al. 2011, Science, 333, 856

Stritzinger, M., Burns, C. R., Phillips, M. M., et al. 2010, AJ, 140, 2036

Stritzinger, M. & Morrell, N. 2009, Central Bureau Electronic Telegrams, 2052, 1

Stritzinger, M. et al. 2009, ApJ, 696, 713

Stritzinger, M. D., Holmbo, S., Morrell, N., et al. 2023, A&A, 675, A82

Stritzinger, M. D., Valenti, S., Hoeflich, P., et al. 2015, A&A, 573, A2

Szalai, T., Vinkó, J., Könyves-Tóth, R., et al. 2019, ApJ, 876, 19

Szalai, T., Vinkó, J., Nagy, A. P., et al. 2016, MNRAS, 460, 1500

Taddia, F., Sollerman, J., Fremling, C., et al. 2019, A&A, 621, A71

Taddia, F., Sollerman, J., Fremling, C., et al. 2016, A&A, 588, A5

Taddia, F., Stritzinger, M. D., Fransson, C., et al. 2020, A&A, 638, A92

Taddia, F. et al. 2013, A&A, 555, A10

Takáts, K., Pignata, G., Pumo, M. L., et al. 2015, MNRAS, 450, 3137

Takáts, K., Pumo, M. L., Elias-Rosa, N., et al. 2014, MNRAS, 438, 368

Tartaglia, L., Elias-Rosa, N., Pastorello, A., et al. 2016, ApJ, 823, L23

Tartaglia, L., Pastorello, A., Benetti, S., et al. 2015, The Astronomer’s Telegram, 8291, 1

Tartaglia, L., Pastorello, A., Sollerman, J., et al. 2020, A&A, 635, A39

Tartaglia, L., Sand, D. J., Groh, J. H., et al. 2021, ApJ, 907, 52

Taubenberger, S., Valenti, S., Benetti, S., et al. 2009, MNRAS, 397, 677

Teffs, J. J., Prentice, S. J., Mazzali, P. A., & Ashall, C. 2021, MNRAS, 502, 3829

Teja, R. S., Singh, A., Sahu, D. K., et al. 2022, ApJ, 930, 34

Terreran, G., Benetti, S., Cappellaro, E., et al. 2016, The Astronomer’s Telegram, 9403, 1

Tomasella, L. 2013, The Astronomer’s Telegram, 5005, 1

Tomasella, L., Benetti, S., Pastorello, A., & Valenti, S. 2012, Central Bureau Electronic Telegrams, 3032, 2

Tomasella, L., Cappellaro, E., Pumo, M. L., et al. 2018, MNRAS, 475, 1937

Trundle, C. et al. 2009, A&A, 504, 945

Tucker, M. A., Ashall, C., Shappee, B. J., et al. 2021, ApJ, 914, 50

Turatto, M., Benetti, S., Cappellaro, E., et al. 1996, MNRAS, 283, 1

Turatto, M., Cappellaro, E., Benetti, S., & Danziger, I. J. 1993a, MNRAS, 265, 471

Turatto, M., Cappellaro, E., Danziger, I. J., et al. 1993b, MNRAS, 262, 128

Turatto, M., Suzuki, T., Mazzali, P. A., et al. 2000, ApJ, 534, L57

Uno, K., Nagao, T., Maeda, K., et al. 2023, ApJ, 944, 204

Valenti, S., Howell, D. A., Stritzinger, M. D., et al. 2016, MNRAS, 459, 3939

Valenti, S., Sand, D., Stritzinger, M., et al. 2015, MNRAS, 448, 2608

Valenti, S. et al. 2012, ApJ, 749, L28

Valerin, G., Pumo, M. L., Pastorello, A., et al. 2022, MNRAS, 513, 4983

Van Dyk, S. D., Zheng, W., Fox, O. D., et al. 2014, AJ, 147, 37

Vinkó, J., Sárneczky, K., Balog, Z., et al. 2009, ApJ, 695, 619

Wang, L. et al. 2003, ApJ, 591, 1110

Wang, Q., Armstrong, P., Zenati, Y., et al. 2023, ApJ, 943, L15

Xi, G., Wang, X., Li, W., et al. 2022, MNRAS, 517, 4098

Yamanaka, M., Maeda, K., Kawabata, K. S., et al. 2015, ApJ, 806, 191

Yamanaka, M., Maeda, K., Tanaka, M., et al. 2016, PASJ, 68, 68

Yaron, O., Perley, D. A., Gal-Yam, A., et al. 2017, Nature Physics, 13, 510

Zelaya, P., Clocchiatti, A., Baade, D., et al. 2017, ApJ, 836, 88

Zhang, J., Wang, X., József, V., et al. 2020, MNRAS, 498, 84

Zhang, J., Wang, X., Vinkó, J., et al. 2018, ApJ, 863, 109

Zhang, K., Wang, X., Zhang, J., et al. 2016, ApJ, 820, 67

Zhang, T. et al. 2012, AJ, 144, 131

Zhang, Y., Zhang, T., Danzengluobu, et al. 2022, PASP, 134, 074201

Zheng, W., Brink, T., & Filippenko, A. V. 2018, Transient Name Server Classification Report, 2018-1219, 1

Zheng, W., Filippenko, A. V., Mauerhan, J., et al. 2017, ApJ, 841, 64

Zheng, W., Silverman, J. M., Filippenko, A. V., et al. 2013, ApJ, 778, L15

References

---

All Tables

All Figures