ZTF SN Ia DR2: High-velocity components in the Si ii $\lambda 6355$

The SN Ia spectra in our sample come from a variety of sources. In this section, we describe the instrument and telescope, as well as data reduction technique for each telescope and instrument setup. Our measurements are sensitive to the wavelength and relative flux calibration across the Si ii $\lambda 6355$ feature, but we do not make use of the absolute values so absolute flux calibration to photometric measurements is not performed. These spectra are publicly released as part of the ZTF DR2 data release, with details on how to access the spectra and metadata provided in Rigault et al. (2024). We perform a number of cuts on our sample to select pre-maxiumum light spectra that reliable phase estimates from their light curves and have sufficient signal to noise in the region of the Si ii $\lambda 6355$ feature (see Section 2.2).

Spectra from the European Southern Observatory’s (ESO) New Technology Telescope (NTT) were obtained with ESO Faint Object Spectrograph and Camera version 2 (EFOSC2; Buzzoni et al., 1984) as part of the ePESSTO/ePESSTO+ collaboration (Smartt et al., 2015) at the La Silla Observatory. These spectra were reduced using a custom built pipeline described in Smartt et al. (2015) to provide wavelength- and flux-calibrated spectra. A number of spectra come from the 2 m Liverpool Telescope (LT; Steele et al., 2004) using the Spectrograph for the Rapid Acquisition of Transients (SPRAT; Piascik et al., 2014) at the Observatorio del Roque de los Muchachos. The spectra were reduced with the pipeline of Barnsley et al. (2012), adapted for SPRAT, along with a custom Python pipeline (Prentice et al., 2018).

Spectra were obtained with the SuperNova Integral Field Spectrograph (Aldering et al., 2002) on the University of Hawai’i 88-inch Telescope (UH88) at the Mauna Kea Observatories. The spectra were reduced using the pipeline of the Spectroscopic Classification of Astronomical Transients (SCAT) Survey (Tucker et al., 2022). This pipeline is insufficient for studies involving absolute spectrophotometric calibration. However, our measurements involve velocities and relative line fluxes of the Si ii $\lambda 6355$ feature and so this is not an issue. We are also not concerned with the region affected by the dichroic crossover ($\sim$5000 - 5200 Å) that require special flat field images that are not applied by this pipeline. Spectra were also obtained at the Las Cumbres Observatory’s FLOYDS Spectrograph on the Faulkes Telescope North (FTN) and on the Faulkes Telescope South (FTS) at Haleakala and Siding Spring, respectively (Brown et al., 2013). These spectra were reduced and calibrated using a custom pipeline¹¹1floyds$\_$pipeline,https://github.com/svalenti/floyds$\_$pipeline described in Valenti et al. (2014).

The KAST Spectrograph (Miller & Stone, 1993) on the Shane 3 m Telescope at the Lick Observatory was used to obtain spectra that were reduced and calibrated using a custom Python pipeline, as detailed in Dimitriadis et al. (2022). Spectra were obtained with the Low Resolution Imaging Spectrometer (LRIS; Oke et al., 1995; McCarthy et al., 1998; Rockosi et al., 2010) on the Keck I telescope and reduced and calibrated using lpipe (Perley, 2019). The spectrum obtained with DEIMOS on the Keck II telescope was reduced using the pypeit software package (Prochaska et al., 2020)²²2https://github.com/pypeit/PypeIt. The Double Beam Spectrograph (DBSP) on the 200-inch Telescope (P200) at Palomar Observatory was used to obtain a number of spectra, these were reduced used a PyRAF-based pipeline, ‘pyraf-dbsp’ (Bellm & Sesar, 2016).

Spectra were obtained with the Asiago Faint Objects Spectrograph and Camera (AFOSC) on the 1.82 m Copernico Telescope at the INAF Osservatorio Astronomico di Padova. They were reduced using standard tasks in iraf, including wavelength calibration with arc lamp spectra and flux calibration using spectrophotometric standard stars. The spectra obtained at Alhambra Faint Object Spectrograph and Camera (ALFOSC) on the Nordic Optical Telescope (NOT) at the Observatorio del Roque de los Muchachos was reduced using a custom pypeit (Prochaska et al., 2020) environment³³3https://gitlab.com/steveschulze/pypeit$\_$alfosc$\_$env. Spectra were obtained using the Dual Imaging Spectrograph (DIS) on the Astrophysical Research Consortium 3.5m Telescope (ARC) at the Apache Point Observatory (APO) and reduced using standard routines in iraf, as in e.g. Sharma et al. (2023). Spectra obtained at the Goodman Spectrograph (Clemens et al., 2004) on the Southern Astrophysical Research Telescope (SOAR) at Cerro Tololo Inter-American Observatory were reduced and calibrated using a custom Python pipeline, as detailed in e.g. Dimitriadis et al. (2022). Spectra were obtained at the Very Large Telescope with the Focal Reducer and Low Dispersion Spectrograph 2 (FORS2) on UT1 of the Very Large Telescope (VLT) at the ESO Paranal Observatory. These spectra were reduced using a custom Python pipeline⁴⁴4forsify, https://github.com/afloers/forsify based on pypeit (Prochaska et al., 2020).

Spectra were observed with the Inamori-Magellan Areal Camera and Spectrograph (IMACS Dressler et al., 2011) mounted on the Magellan-Baade telescope at Las Campanas observatory. The data reduction was performed in iraf⁵⁵5IRAF was distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. following standard reduction procedures. The Low Dispersion Survey Spectrograph 3 (LDSS-3; Stevenson et al., 2016) on the Magellan-Clay Telescope at Las Campanas Observatory was used to obtain one spectrum, which was reduced using standard iraf routines as described in Hamuy et al. (2006) using the same procedure as described for IMACS above.

Spectra were obtained with the Optical System for Imaging and low-Intermediate Resolution Integrated Spectroscopy (OSIRIS; Cepa et al., 2000) on the Gran Telescopio Canarias (GTC) at the Observatorio del Roque de los Muchachos. These data were reduced and calibrated following the method in Piscarrera (in prep.) using custom routines based on pypeit (Prochaska et al., 2020). Some spectra were obtained with the Device Optimized for the Low Resolution (DOLORES) on the Telescopio Nazionale Galileo (TNG) at the Observatorio del Roque de los Muchachos and reduced using pypeit (Prochaska et al., 2020), following Das et al. (2023).

Spectra were obtained with Gemini Multi-Object Spectrographs (GMOS) (Hook et al., 2004; Allington-Smith et al., 2002) on the Gemini North Telescope at the Mauna Kea Observatories. The spectra were reduced using standard iraf/pyraf and Python routines, following Dimitriadis et al. (2022). The Intermediate-Dispersion Spectrograph and Imaging System (ISIS) and the Auxiliary-port Camera (ACAM; Benn et al., 2008) on the William Herschel Telescope (WHT) at the Observatorio del Roque de los Muchachos were used to obtain spectra. These were reduced and calibrated using standard iraf routines. The Robert Stobie Spectrograph (RSS) (Kobulnicky et al., 2003) on the South African Large Telescope (SALT) (Buckley et al., 2006) at the South African Astronomical Observatory (SAAO) was used to obtain spectra. The spectra were reduced using the custom pipeline, PySALT (Crawford et al., 2010) to produce wavelength- and flux-calibrated spectra. Spectra were obtained with the DeVeny spectrograph on the 4.3 m Discovery Channel Telescope (DCT), which was reduced using standard iraf routines, including wavelength and flux calibration (Hung et al., 2017). The Spectral Energy Distribution Machine (SEDm; Blagorodnova et al., 2018; Kim et al., 2022; Lezmy et al., 2022) on the P60 (Cenko et al., 2006) was used to obtain spectra, which were subsequently reduced using pysedm (Rigault et al., 2019).

We do not require absolute flux calibration for our measurements but our measurements are impacted by the relative flux and the wavelength calibration. Not all the spectral reduction techniques used for our sample provide meaningful flux uncertainties. Therefore, as described in Section 3.1.1, we estimate the uncertainty on the flux using the standard deviation of the values in certain continuum regions near to the feature of interest. For SN Ia spectra, where flux uncertainties were available, we cross-checked against the uncertainty estimates from the standard deviation near the feature of interest and found them to be consistent. In Section 2.2, we describe how we performed a cut on the sample based on the signal to noise in the Si ii $\lambda 6355$ region. Each remaining spectrum after this cut was inspected manually while choosing the continuum regions and if any cosmic rays or host galaxy lines were identified they were masked out prior to fitting the Si ii $\lambda 6355$ feature. We include an uncertainty in our fitting of 200 km s^-1 to account for additional velocity offsets associated with motions of the SNe Ia in their host galaxies.

The starting point for defining our sample is the subset of the 3628 confirmed SNe Ia in the DR2 with spectra from the aforementioned facilities for which sufficient data information and reliable reductions could be performed, amounting to 3585 targets with 5028 spectra. 43 SNe Ia that were included in the full DR2 data release are excluded from our sample due to a lack of sufficient information on the observations and data reduction techniques. As our analysis will be restricted to spectra in the pre-peak regime, we require sufficient photometry to produce a reliable estimate of spectral phase. We used the suggested cuts of Rigault et al. (2024) on the SALT2 light-curve fitter (Guy et al., 2007) outputs, of fit probability greater than 10^-5 and uncertainties smaller than the quoted values for the light curve width parameter, $x_{1}$ ($\delta x_{1}<1$), color parameter, $c$ ($\delta c<0.1$), and the time of maximum light, $t_{0}$ ($\delta t_{0}<1$ d). We do not impose any constraints on the measured values of $x_{1}$ or $c$.

To maximise our spectral sample, we investigated all SNe Ia that failed to meet these criteria in order to avoid cutting otherwise acceptable spectra solely due to poor phase estimates from sparse photometry. In many cases this involved fitting supplementary photometry obtained from other surveys. A summary of this investigation and the updated phase estimates can be found in Appendix A. Following this investigation we are left with 4572spectra for which we have reliable phase estimates. As we are specifically interested in the evolution of Si ii $\lambda 6355$ high-velocity components, we limit our search to the pre-peak regime when these features are expected to be present. Cutting any spectra with a phase post-peak, we are left with 2362spectra from 1801SNe Ia.

To assess the spectral quality in the vicinity of the Si ii $\lambda 6355$  we define a local signal-to-noise ratio (SNR) as the ratio of the depth of the Si ii $\lambda 6355$ to the standard deviation in the regions of the continuum to the red and blue (continuum region selection shall be discussed in detail in Sect. 3.1). This measurement is made once the continuum has been removed and is repeated 1000 times with different continuum selection points. The final SNR estimate is taken as the mean of these 1000 values. In making these estimates we identified 432 spectra that either did not cover the wavelength range of the Si ii $\lambda 6355$ feature, or did not have a clearly visible Si ii $\lambda 6355$ feature from which to compute an estimate of local SNR. These spectra were cut from the sample and the measurements of SNR for the remaining objects can be found in Fig. 1, along with the distribution of dispersion (average separation between pixels in Å) and the redshift distribution of the corresponding SNe. The large peak at $\sim$25 Å in the spectral dispersion panel corresponds to the spectra coming from the Spectral Energy Distribution Machine (SEDM; Blagorodnova et al., 2018; Rigault et al., 2019) on the P60 (Cenko et al., 2006) at the Palomar Observatory. The SEDM provides $\sim$60 per cent of the full DR2 spectral sample. The usefulness of the SEDM spectra in the search for HVFs will be determined by the simulations in Sect. 3.2, in which we will place a cut on both dispersion and SNR.

Our aim is to distinguish between those SN Ia spectra that possess a HV component in the Si ii $\lambda 6355$ line and those that do not. The feature fitting with single and double component models shall be performed with an Markov-Chain Monte Carlo (MCMC) framework using the emcee package (Foreman-Mackey et al., 2013), treating the local continuum as linear and leaving the corresponding slope and intercept as free parameters in the fitting.

To investigate the accuracy of the MCMC/BIC method described in Section 3.1 for identifying HV components, as well as its ability to recover the injected parameters, we generated synthetic Si ii $\lambda 6355$ features. The results from these simulations were also used to inform cuts on dispersion and SNR in the observed data. In Section 3.2.1, we describe the construction of our simulated data for testing our method, while in Section 3.2.2, we detail our detection efficiencies and in Section 3.2.3, we discuss the correction for biases in our fitting based on our simulations.

The top panel of Fig. 7 displays the measured velocity evolution for the 329 spectra from our 307 SNe Ia, split into the photospheric components, and high-velocity components where classified. Power-law fits were performed to the velocities as a function of phase for the HV components, as well as three samples of the PV velocities, i) all the PV components, ii) just the PV components with HV counterparts, and iii) just the singular PV components. These three fits are consistent from around $-$10 d onwards, with the PV velocities in objects possessing HV components drifting to higher velocities at early phases. This agrees with the findings of Silverman et al. (2015) in that HVFs tend to be accompanied by higher velocity PVFs. We also plot in grey the velocity evolution from the PTF sample (Maguire et al., 2014) that was used as reference for our Si ii $\lambda 6355$ feature generation in the simulations. These PTF measurements were performed on the overall Si ii $\lambda 6355$ features without distinction between PV and HV components. The PTF velocity evolution therefore likely has some level of contamination from HVFs, and would be expected to lie somewhere between our measured PVF (blue) and HVF (pink) evolution curves, as it does.

In the bottom panel of Fig. 7 we present the measured depths and widths of the singlet lines making up the PV and HV components. As for the velocity evolution we plot the comparison PTF data as grey contours to show how our final measurements align with the distributions chosen to inform the feature generation in the simulations. As expected, the distribution of PVF measurements lie in the same parameter space as the PTF data, but our measurements are slightly more compact in both dimensions which, as with the velocity, is likely due to HVF contamination in the PTF dataset inflating the depths and widths in some cases. For the generation of the HVF components in the simulations, we made the assumption that the widths and depths of the HV components followed the same distribution. However, as clearly visible in Fig. 7, these components tend to cluster more towards shallower depths and narrower widths. While some false positive classifications may be contaminating the HVF cluster, there remains a clear distinction between the HVF and PVF distributions. To evaluate the effect that this discrepancy has upon the calculated detection efficiencies we recalculated the true-positive rates from the simulations, excluding any simulated spectra from the calculations with $a_{\text{HV}}>0.25$ or $c_{\text{HV}}>70$ Å so has to be more in agreement with the measurements from the observations. The residuals between the true positive rates these updated simulations and the original simulations are presented in Figure 8. In most cases, the true-positive rate increases with the cut to only focus on the observed region of the parameter space. This implies that we are less sensitive to finding HV components with similar widths and depths to their PV counterparts - likely due to increased degeneracy in the fit. Therefore, we updated the detection efficiencies with a GP interpolation (following the method of Section 3.2.2) to these new values.

The evolution of the ratio of pEW${}_{\text{HV}}$/pEW${}_{\text{PV}}$ (R${}_{\text{HVF}}$) is seen for our observed sample in the middle panel of Fig. 9. In general, the uncertainties on the measured R${}_{\text{HVF}}$ values are large as many of the two-component fits exhibit high levels of degeneracy and posteriors of the numerator and denominator in the ratio (pEW${}_{\text{HV}}$ and pEW${}_{\text{PV}}$) are inversely correlated.

The SNe with multiple spectra allow us to probe the evolution of this ratio within individual objects, with the overall sample - including single-epoch objects - giving a more global view of how the distribution of R${}_{\text{HVF}}$ changes with phase. For most of our multi-phase objects, we possess only two spectral epochs, separated by less than $\sim$2.5 d, leading to very little evolution in R${}_{\text{HVF}}$ which appears to remain approximately constant on such small timescales. However, for ZTF19aatlmbo (SN 2019ein) and ZTF18abauprj (SN 2018cnw), we have three and four spectra spanning 9.3 d and 10.6 d, respectively. ZTF18abauprj displays a clear decrease in this ratio with time, with the HV component fading away completely somewhere between $-11.4$ and $-5.8$ d, whereas ZTF19aatlmbo exhibits a far flatter evolution; albeit with large uncertainties on the spectrum at $-10.9$ d. The R${}_{\text{HVF}}$ evolution for these multi-spectra objects is consistent with previous studies showing a decline with phase, with HV components starting out strong and fading away over time (Silverman et al., 2015). As clear from these two objects, while we see a decrease in R${}_{\text{HVF}}$ with time in individual SNe, this decay occurs at different phases, and we see many multi-spectra targets still exhibiting HVFs after the HV component in ZTF18abauprj had faded away completely. When combining the multi-spectra and single-spectra objects to get a global view of how R${}_{\text{HVF}}$ varies with phase, we find no clear evolution, with some large and small measurements of this ratio and early and late phases, further supporting the idea that the fading away of HV components occurs at different phases in different objects.

The bottom panel of Fig. 9 presents the evolution of the velocity separation, $\Delta v$, with phase. As for R${}_{\text{HVF}}$, the multi-epoch spectra with small temporal differences do not shed much light upon the phase evolution of $\Delta v$. For ZTF18abauprj, we see an upwards evolution and the converse for ZTF19aatlmbo, implying that there is not a singular strict evolutionary trend for such features. However, when looking at the global evolution of the sample, we observe a dearth in the higher velocity separation spectra as we approach maximum light, with all the separations from HVF spectra after $-4$ d clustering below $\sim$5000 km s^-1. This indicates more separated HV components tend to fade away at earlier phases compared to those that are more entangled with their PV counterpart.

In order to analyse global properties of SNe Ia containing Si ii $\lambda 6355$ HVFs, we are required to understand and mitigate the potential biases present in our sample that make it unrepresentative of the global SN Ia population. To this end we define a ‘control sample’ as the 994 SNe Ia from the DR2 volume-limited (z $\leq$0.06) sample presented in Rigault et al. (2024). At this point we also remove any objects from our sample of 307 SNe that did not pass the initial suggested cuts outlined in Sect. 2.1 (fit probability $>10^{-5}$, $\delta x_{1}<1$, $\delta c<0.1$, $\delta t_{0}<1$ d), reducing our ‘base sample’ to 261 SNe Ia. Through Kolmogorov-Smirnov (KS) tests between this base sample and the control sample, we can identify parameters in which our sample exhibits a significant level of bias and then eliminate these by way of a redshift-based cut. KS p-values below the threshold of 0.05 signify that the difference between the two samples is statistically unlikely to have occurred simply by chance, and points to them coming from separation populations.

The top panel of Fig. 10 presents these KS-test results for a number of measurements between the control sample and the base 261 SNe sample after imposing different redshift cuts, e.g., the data point at $z=0.04$, corresponds to the KS-test between the base sample limited to SNe Ia with $z\leq 0.04$ and the full control sample. The parameters investigated here are the SALT2 stretch and colours parameters, $x_{1}$ and $c$, the phase of first photometric detection, the uncertainty in the date of maximum light ($\delta t_{0}$), the host galaxy mass (measurements from Smith et al. in prep), and the host galaxy redshift (not shown on plot). Our base sample of 261 SNe compared to the control sample has KS p-values that are greater than our bias threshold of 0.05 in all parameters for all redshift cut values, except $c$ and host redshift. As we introduce more restrictive redshift cuts to our sample, we see the KS p-values for $c$ rise, until it crosses our threshold to be considered low bias at a redshift cut of z = 0.06.

The comparison of the redshift distribution between the base and control samples exhibits such small p-values - for all potential redshift cuts - that it falls below the y-axis scaling of the plot in Fig. 10. This difference is not surprising at redshift cuts either far above or far below that of the control sample (z = 0.06). With a redshift cut equal to that of the control sample we see that the difference between the populations is that our sample peaks at lower redshifts, most likely due to the apparent brightness required to obtain a high enough SNR spectrum pre-peak to enter our sample. The redshift distribution describes the geometric distribution of our objects in space, with the intrinsic physics of the objects described by the other parameters. Therefore, the bias in the redshift distribution should not impact our results given that we observe no statistically significant bias in the other parameters. We consequently place a cut at z = 0.06, leaving us with 210 spectra from 190 SNe from which to investigate global HVF properties. The distributions of these parameters in this low-bias cut sample can be seen as step plots in the bottom panels of Fig. 10 along with the control sample as the bars.

In Fig. 11 we present the distribution of velocity separations found in the 64 HVF spectra of the low-bias sample as the solid bar histogram. In general we observe a distribution skewed towards smaller velocity separations with very few spectra showing the extreme separations of 8000 km s^-1. As seen in Fig. 3 our classification method is very sensitive to the larger velocity separations, however its detection efficiency drops off as the features encroach on one another and become increasingly entangled. As an example, we have a true-positive rate of $\sim$25% for features with SNR 8 and dispersion 2 Å/pix that possess a velocity separation of 4000 km s^-1, implying that for every one spectrum that we identify with this parameter set, we miss three others; the true value is the number of identified targets divided by the fractional true-positive rate. Therefore, we can draw efficiency values from our 3D GP interpolation of the true-positive rates to scale the distribution of velocity separations on a spectrum-by-spectrum basis. This is first performed using the measured values of velocity separation, for which a Gaussian KDE of the probability density function (PDF) can be seen in Fig. 11 as the black solid line. In order to probe the potential variation of this distribution we employ a Monte-Carlo method, sampling each of the velocity separations in the sample from their individual distributions - described by their measured medians and upper and lower uncertainties - before then performing the detection efficiency scaling with the resulting $\Delta v$ values. This procedure is repeated 1000 times with each of these PDFs plotted as a faint orange line in Fig. 11. These PDFs are all multiplied by the number of spectra that they represent and therefore, represent spectrum densities instead of probabilities.

In order to probe the potential variation of this distribution we employ a Monte-Carlo method, sampling each of the velocity separations in the sample from their individual distributions - described by their measured medians and upper and lower uncertainties - before then performing the detection efficiency scaling with the resulting $\Delta v$ values. This procedure is repeated 1000 times with each of these PDFs plotted as a faint orange line in Fig. 11. In each iteration we account for the variation due to false positives by randomly reclassifying some HVF spectra with $\Delta v<5000$ km s^-1 as non-HVF in accordance with our adopted conservative false positive rate of 2%. These final PDFs are all multiplied by the number of spectra that they represent and therefore, represent spectrum densities instead of probabilities.

Figure 12 presents the phase distribution of the spectra displaying HVF against our full low-bias sample. The full sample is presented as the faint bars with the HVF spectra as the solid bars. As before, we perform the detection efficiency scaling on a spectrum-by-spectrum basis and plot the resulting density function from the velocity separation values as the solid black curve. The spectrum density function of the full phase distribution is plotted as the dotted black line. As for the velocity separation analysis, we perform Monte-Carlo sampling for these two phase distributions, randomly reclassifying some HVF spectra with $\Delta v<5000$ km s^-1 as non-HVF in accordance with the false positive rate and drawing velocity separations and phases from the corresponding individual distributions to calculate new PDFs. In each of these 1000 iterations we also integrate over three phase bins to calculate the percentage of spectra exhibiting a Si ii $\lambda 6355$ HV component. These three regions are divided up to each hold $\sim$1/3 of the 64 HVF spectra from the low-bias sample, giving bins for $<-$11 d, $-11$ to $-6$ d and $>-$6 d. The percentages of HVF spectra for these three phase ranges are shown in the bottom panel of Fig. 12 along with the shaded regions representing 68% confidence intervals calculated via the Clopper-Pearson method. These calculations indicate that we find Si ii $\lambda 6355$ HV components in 76$\pm^{7}_{9}$% of spectra before $-11$ d, 46$\pm{7}$% of spectra between $-11$ and $-6$ d, and 29$\pm{5}$% of spectra between $-6$ d and maximum light. If we were to use the original detection efficiencies before introducing the cuts to only consider simulations with $a_{\text{HV}}\leq 0.25$ and $c_{\text{HV}}\leq 70$ Å (see Section 4.1), we calculate slightly higher percentages of 83$\pm^{6}_{8}$%, 49$\pm{7}$%, and 32$\pm{5}$% for these three phase intervals. Therefore, Si ii $\lambda 6355$ HV components are close to $\sim$2.5 times more common at early phases, but still appear in approximately one third of spectra in the week before maximum light.

In Fig. 13 we present the distributions of the SALT2 light curve parameter $x_{1}$ for the SNe Ia in our full low-bias sample and only those showing evidence of HVFs as the faint and solid histograms, respectively. This information is then displayed in the form of SN density functions for the full sample (dotted) and those for which we have a spectrum with an identified HVF (solid). As for the $\Delta v$ and phase parameters before, we perform a Monte-Carlo sampling of the individual $x_{1}$ measurements over 1000 iterations to assess the amount of variation in these density curves accounting for the estimated false positive rate and uncertainties in the $x_{1}$ values. For each of these iterations, we integrate under the two density curves over two regions ($x_{1}<0$ and $0\leq$ $x_{1}$) and calculate the percentage of SNe exhibiting HV components in the Si ii $\lambda 6355$. As visible in the bottom panel of Fig. 13, these two percentages (30$\pm{5}$% and 27$\pm^{6}_{5}$%) are consistent with one another and support the idea of HVFs being ubiquitous across the SN Ia population. As before for the phase percentages, these uncertainties represent 68% confidence intervals.

As seen in Fig. 12, Si ii $\lambda 6355$ HVFs are more prevalent at earlier phases and there is a significant decrease in the percentage of SN Ia spectra displaying HVF features as the phase approaches maximum light. By cutting the sample and repeating these percentage calculations with only the earliest phases, we can test the impact of this phase dependence on any potential trend with $x_{1}$. We introduce a phase cut at $-8.7$ d - the median HVF spectral phase - and recalculate the rates of occurrence in HVF in low and high $x_{1}$ SNe Ia at these earlier phases as 43$\pm^{11}_{10}$% ($x_{1}<0$) and 56$\pm{11}$% ($0\leq$ $x_{1}$). While the difference between these two percentages is slightly larger than before (and in the opposite direction), they still overlap in the 1$\sigma$ uncertainty region and therefore support the idea of HVF ubiquity.

We also investigated the relationship between the presence of HVFs and light-curve parameters (peak absolute ZTFg-band magnitude, the g-band decline rate in 15 days post maximum, $\Delta$m_15,ZTFg) measured from GP fitting in Dimitriadis et al. (2024). These are shown in the top panels of Fig. 14. We employ KS tests to quantify the likelihood of the two populations hailing from the same parent population. Resulting p-values that fall below the threshold of 0.05 indicate that the difference between the two samples for a given measurement is unlikely to have occurred simply by chance, and may point to a physical underlying distinction. We observe no significant difference between the HVF and non-HVF populations in terms of the peak ZTFg-band absolute magnitude and $\Delta$m₁₅ parameters with the p values as 0.67 and 0.81 respectively. This again supports the ubiquity of these features across the Ia population.

The host galaxy mass and local g-r colour (Smith et al. in prep.) are shown in the bottom panels of Fig. 14. The p values for the host mass and local colour distributions of the HVF and non-HVF samples are large relative to the 0.05 threshold, with values of 0.65 and 0.59, respectively, again consistent with a lack of a difference in the local environment between SNe Ia with and without a Si ii $\lambda 6355$ HVF.

As before for investigating trends with $x_{1}$, we repeated these statistical tests exclusively with spectra before $-8.7$ d. We once again find large p-values for peak absolute magnitude and $\Delta$m₁₅ in the ZTFg-band (0.91 and 0.52) as well as for the host galaxy mass and local g-r colour (1.0 and 0.94), indicating no difference between the HVF and non-HVF populations in these observables.

Silverman et al. (2015) investigated the properties of SNe Ia both with and without HVF in the context of the velocity-based classification scheme of Wang et al. (2009). The ‘Wang classification scheme’ divides objects into normal velocity (NV${}_{\text{W}}$) and high velocity (HV${}_{\text{W}}$) subclasses depending on whether the photospheric velocity (measured from the Si ii $\lambda 6355$ feature) is less than or greater than 11800 km s^-1 around peak ($-$5 to 5 d), respectively. We chose to add the ‘W’ subscript here to differentiate between the high-velocity Wang subclass (HV${}_{\text{W}}$) and the high-velocity components (HV). While the cutoff of 11800 km s^-1 was employed by Silverman et al. (2015), we chose to adopt a cutoff of 12000 km s^-1 for consistency with other ZTF DR2 studies (Burgaz et al., 2024). Our Wang classifications are drawn from the latest spectrum we have for each object. In the case of the HV${}_{\text{W}}$ class, only spectra later that $-$5 d were used. However, for any object for which the latest pre-peak spectrum has PV $<12000$ km s^-1, regardless of phase, we denote it as NV${}_{\text{W}}$ as the photospheric velocity is not expected to increase and the objects will still have PV$<12000$ km s^-1 around peak. These classifications are then augmented by classifications from (Burgaz et al., 2024) for which they examined a slightly later phase range of $-$5 to 5 d with respect to peak brightness.

Figure 15 presents the pEW ratio, R${}_{\text{HVF}}$, against the photospheric velocity of the Si ii $\lambda 6355$ feature for the low-bias sample of 190 SNe Ia. When investigating the same parameters, Silverman et al. (2015) found a dearth of NV${}_{\text{W}}$ SNe (PV$<12000$ km s^-1) exhibiting HVFs in the Si ii $\lambda 6355$ (R${}_{\text{HVF}}$ $>$0), with these features found more predominantly in HV${}_{\text{W}}$ SNe. Contrary to this, we find HVFs in 26 of the NV${}_{\text{W}}$ classified objects in our low-bias sample (36 in the full sample), largely populating the empty region of the parameter space seen by Silverman et al. (2015). An increase in R${}_{\text{HVF}}$ with increasing photospheric velocity was also observed by Silverman et al. (2015). While this correlation makes intuitive sense for the evolution of individual objects, with R${}_{\text{HVF}}$ and photospheric velocity both decreasing over time, we observe no such global correlation in Fig. 15. As before with the lack of a global downwards trend in the middle panel of Fig. 9, this indicates that the onset of the R${}_{\text{HVF}}$ evolution begins at different phases for different objects.

As described above, Wang classifications (NV${}_{\text{W}}$ vs. HV${}_{\text{W}}$) are generally measured using the Si ii $\lambda 6355$ photospheric velocity in the phase range $-$5 to 5 d with respect to peak. However, the specific measurement of the photospheric velocity will change depending upon whether we employ a single- or a double-velocity component model. Historically this feature has been treated as a single component around peak, believed to be free of any HV components. However, our study, as well as the results from Silverman et al. (2015), suggests that HV components might be more common at these phases (later than $-$5 d) than previously thought, with one third of SN Ia spectra between $-$5 d and peak displaying a HVF. Therefore, we pose the question, if in these large spectroscopic studies we were to consider and account for HV Si ii $\lambda 6355$ absorption close to maximum light, what percentage of the HV${}_{\text{W}}$ classifications would be overturned in favour of a NV${}_{\text{W}}$ subtype?

In the top panel of Fig. 16, we compare the photospheric velocity from the single component fits against the ‘true’ photospheric velocities (i.e. single component fits for spectra without a HVF, and double component fits for spectra with a HVF). All spectra identified to have a HV component exhibit lower photospheric velocities in their two-component fits as would be expected. We indicate the Wang classification cut off velocity of 12000 km s^-1 by the vertical and horizontal dotted lines, highlighting the region in the bottom right in which the objects would be classified as HV${}_{\text{W}}$ with the single component model, but as NV${}_{\text{W}}$ by the double-component model (orange points). When considering the 85 objects from the low-bias sample which possess a spectrum later than $-$5 d, there are five SNe Ia (seven of the 150 SNe from the full sample with a spectrum later than $-$5 d), indicating that $26\pm^{14}_{11}$% of HV${}_{\text{W}}$ classifications ($24\pm^{11}_{8}$% for the full sample) before peak would be incorrect if we were to not consider HV components. These uncertainties represent 68% confidence intervals and are calculated as binomial with the Clopper-Pearson interval. Such large uncertainties are the result of low number statistics and while this leaves the percentages fairly unconstrained, this suggests that HVF components could cause a significant rate of false HV${}_{\text{W}}$ classifications in the 5 d before maximum light.

We similarly perform this analysis with respect to the Branch classification scheme (Branch et al., 2006, 2009), which divides up the SN Ia population based upon the equivalent widths of the Si ii $\lambda 6355$ and Si ii $\lambda 5972$ lines. Broad line SNe Ia (BL) exhibit pEW${}_{\text{Si 6355}}>105$ Å with a pEW${}_{\text{Si 5972}}<30$ Å. While we have no information on the Si ii $\lambda 5972$, we can examine the effects of HV components upon the measured Si ii $\lambda 6355$ pEW. The bottom panel of Fig. 16 displays the measured photospheric pEW coming from the single component fits against the ‘true’ photospheric pEWs (as above for the Wang analysis). The cutoff width at 105 Å is indicated by the dotted lines, with the orange points once again corresponding to those spectra that would receive an incorrect classification if we were to not consider the HV components. In the Branch scheme parameter space there are four SNe that are misclassified (six for the full sample), resulting in $20\pm^{13}_{9}$% of incorrect BL classifications ($18\pm^{9}_{7}$% for the full sample) with respect to the Si ii $\lambda 6355$ line. In their sample of BL SNe Ia, Yarbrough et al. (2023) also found the BL population to have higher Si ii $\lambda 6355$ velocities than their CN counterparts, further supporting this idea as the HV components would also cause an offset in measured single component velocities to higher values. These uncertainties represent 68% intervals as before for the Wang reclassifications, and once again the low number of datapoints results in loosely constrained percentages.

With generally smaller velocity separations in this phase range around peak, false positives are more probable and could have the inverse effect in both these cases. Incorrectly employing a double-component model to a Si ii $\lambda 6355$ feature with no HVF would artificially decrease the photospheric velocity and width, potentially resulting in HV${}_{\text{W}}$ SNe being classified as NV${}_{\text{W}}$ and Branch BL SNe being missed. Therefore, the misclassification rates calculated here represent upper limits.

While the specific origins and formation channels of these HV components of the Si ii $\lambda 6355$ feature remain unclear, they provide clear evidence of intermediate-mass element material in the upper ejecta of the majority of SNe Ia. Regardless of the explosion mechanism or progenitor channel, any model aiming to reproduce these features will require some fraction of silicon at the corresponding high velocities.

In Fig. 17 we present the measured velocity distributions for the PV and HV components for the HVF classified spectra in our sample, with the Si-abundance profiles from a number of explosion models from the Heidelberg Supernova Model Archive (HESMA; Kromer et al., 2017), arising from three different explosion mechanisms. In pink we plot the delayed-detonation (DDT) models (Röpke et al. 2012; Seitenzahl et al. 2013; Ohlmann et al. 2014) arising from an initial subsonic deflagration that later transitions to a supersonic detonation. In blue, we plot the double-detonation models (Sim et al. 2012; Gronow et al. 2020; Gronow et al. 2021), which involves a detonation in a helium shell atop the white dwarf, sending a shock-wave inwards which converges in - and subsequently detonates - the core. Finally, in yellow we show the silicon abundances from the pure subsonic deflagration models (Fink et al. 2014; Kromer et al. 2015).

As expected, the pure deflagration models are very homogeneous in nature, possessing very flat distributions of the elements in the model. As a result of the low energy they also do not produce material moving as velocities faster than around 15000 km s^-1, and therefore, are incapable of producing the HV - and many of the PV - components that we see in the DR2. The DDT and double-detonation models possess more layered ejecta structures, with peaks in the silicon abundance profiles in the 10000 to 15000 km s^-1 range. The majority of these models also exhibit silicon at velocities covering the full HVF velocity distribution. The formation of such HVFs is dependant upon the ionisation profile, which in turn is dependant upon not only the abundances, but also the density profile of the plasma. This will be the focus of a follow-up modelling study with a handful of well-sampled HVF objects (Harvey et al. in prep.).

We also note here that the models presented in Fig. 17 are designed to reproduce the overall SN Ia photometric and spectroscopic evolution coming from the primary white dwarf itself. These HVFs may instead be the result of environmental components such as CSM, which are not included in the models. Through the spectroscopic modelling of the HVFs exhibited by SN 1999ee, Mazzali et al. (2005b) explored the ejecta conditions needed to recreate the evolution. In terms of an abundance enhancement, they required a region of the upper ejecta to be dominated by silicon and calcium ($\sim$90 and $\sim$10 per cent respectively), which is hard to justify in terms of nuclear burning. Therefore, they concluded that the density enhancements in these regions to be a more plausible reason for the HVFs. In Mazzali et al. (2005a) they suggested the sweeping up of surrounding CSM to be responsible for the line-forming region higher up in the ejecta, pointing to a thick disk and/or high density CSM environment produced by a companion wind. They also proposed angular fluctuations as a potential cause of density enhancement, with HVFs being the result of three-dimensional effects. This was further supported by Tanaka et al. (2006) who investigated various density enhancement geometries and degrees of photospheric coverage to produce the observed HVF population. This multidimensional requirement is supported by the observed polarisation difference between the HV and PV components of the Ca NIR triplet in SN 2001el (Wang et al. 2003; Kasen et al. 2003) and in the Si ii $\lambda 6355$ of SN 2019ein/ZTF19aatlmbo (Patra et al., 2022)