© The Authors 2026

1. Introduction

Understood to be the thermonuclear explosions of white dwarfs enabled by some interaction with a binary companion, Type Ia supernovae (SNe Ia) are the most common transients we observe in the night sky. While famed for their use as cosmological distance indicators due to their standardisability via the empirical Phillips relation (Pskovskii 1977; Phillips 1993; Tripp 1998), extensive observations of these objects in the past few decades have uncovered a large spectrophotometric diversity leading to subdivision into a number of different subclasses (see Taubenberger 2017 for a review).

Despite our empirical understanding of thermonuclear SNe, there remain many open questions surrounding their progenitor systems and explosion mechanisms, as well as a number of unexplained observed properties. One such observation is the tendency to find high-velocity features (HVFs), which are secondary higher-velocity absorption lines that appear to form in a region of the ejecta detached from the photosphere (e.g. Gerardy et al. 2004; Mazzali et al. 2005a; Wang et al. 2009; Childress et al. 2013; Maguire et al. 2014). They are most commonly seen in the Ca II and Si II features as secondary absorption minima several thousands of km s⁻¹ to the blue of the photospheric features. The focus of this paper is the Si IIλ6355 line.

Several scenarios have been proposed in the literature to explain these features, including abundance and/or density enhancements in the outer ejecta, as well as changes in the ionisation state of the relevant elements. Appearing more frequently than in the Si IIλ6355 feature, previous studies have mostly focussed upon the HVF in the Ca II near-infrared (NIR) triplet. Spectropolarimetry of SN 2001el presented by Wang et al. (2003) demonstrated high levels of polarisation across the spectrum, with the HV Ca II component polarised differently, suggesting that this component is kinematically distinct. A geometrical study of the high-velocity (HV) ejecta of this object by Kasen et al. (2003) concluded that the HV calcium was produced by a spherically non-symmetric clumped shell of material. Thomas et al. (2004) studied the HV Ca II in SN 2000cx through both its NIR and H&K absorption features; they found that while 1D models are sufficient to explain the HVF in the NIR triplet, a 3D solution was required to produce the higher-velocity structure simultaneously in the Ca II H&K and NIR features. This solution was comprised of a clump of material along the line of sight, partially covering the photosphere. The HVF in SN 2003du was proposed by Gerardy et al. (2004) to be the result of interaction between the ejecta and the local circumstellar environment. Mulligan & Wheeler (2017) used SYN++ to investigate if HV Ca II NIR features could be explained by the presence of a shell of circumstellar material (CSM), and put limits on the mass of such a shell. Tanaka et al. (2006) performed three-dimensional modelling of different geometries of HV material, which was compared to six SNe Ia in Tanaka et al. (2008). They suggested that asymmetries or abundance enhancements in the outermost layers of the ejecta are likely required to explain the HVFs seen.

The search for HVFs in early SN Ia spectra by Mazzali et al. (2005a) found all objects to exhibit HV absorption in the Ca II NIR feature; some objects also displayed these features in the Si IIλ6355, leading to the conclusion that HVFs are ubiquitous across the SN Ia class. It was suggested that density enhancements instead of abundance enhancements are more likely the explanation for these features given the very high abundance enhancements required to match observations of SN 1999ee (Mazzali et al. 2005a,b). Silverman et al. (2015) found that Si IIλ6355 HVFs were significantly rarer than HVFs in the Ca II NIR triplet, and spectra must be obtained at the earliest epochs for the Si IIλ6355 HVFs to be detected. The connection of HVFs of Si IIλ6355 to other observables were also studied in Zhao et al. (2015). Quimby et al. (2006) argued that the ubiquity of these features disfavours line-of-sight effects such as filaments of material. High-cadence spectroscopy of the otherwise normal SN 2009ig by Marion et al. (2013) exhibited HV components to a number of features, Si II, Si III, S II, Ca II, and Fe II, allowing them to constrain well the transition from HVF to photospheric-velocity features (PVFs).

The post-maximum light spectra of SN 2012fr were compared to the Chandrasekhar mass deflagration model W7 (Nomoto et al. 1984) using the radiative-transfer code PHOENIX by Cain et al. (2018), but the origin of HVFs at early phases was not investigated. Modelling of individual very early spectra of 14 SNe Ia was performed in Ogawa et al. (2023) using the radiative-transfer code TARDIS (Kerzendorf & Sim 2014), but they did not specifically investigate the presence of HVFs using distinct density or abundance enhancements above the photosphere.

In Harvey et al. (2025) (hereafter H1) the Zwicky Transient Facility Cosmology Data Release II (ZTF Cosmo DR2) was searched for Si IIλ6355 HVFs; these features were identified in 75 SNe from the 307 that passed our sample cuts. From the phase distribution of the spectra in this sample, we calculated a rate of 78${\pm }_{17}^{12}$% per cent of spectra before −11 d (with respect to maximum light) to show some HV component to the Si IIλ6355. When comparing the distributions of various observables between the SNe found to exhibit HVFs and those that did not, we identified no difference between the two populations in terms of SALT light curve width parameter, x₁ (Guy et al. 2010), peak magnitude in the ZTFg band, the magnitude decline rate from peak to 15 d after in the g band (Δm_15, ZTFg), host galaxy stellar mass, or the host galaxy local g-r colour.

The aim of this paper is to explore the possibility of producing the observed HVF evolution found in a number of SNe Ia by the introduction of a density enhancement in the outer ejecta. We developed models for a number of well-sampled SNe Ia using the radiative-transfer code TARDIS, aiming simply to reproduce the photospheric features. We then ran a number of simulation grids adding in various density enhancements to the base models with the aim of reproducing the observed HVF evolution. The bulk of this work focuses on the Si IIλ6355; however, HVFs of other species are also discussed.

In Section 2 we describe the collation of our sample and in Section 3 we present the modelling method. The results of the modelling are analysed in Section 4, followed by a discussion of the implications of these results in Section 5.

2. Observations

For our SN Ia sample, for which we produced custom TARDIS models, we targeted events that exhibit clearly separated Si IIλ6355 HVFs and are photometrically and spectroscopically well sampled within the first two weeks after explosion. We conducted a literature search and found the following suitable events that have early spectra that display distinct Si IIλ6355 HVF and PVF: SN 1994D (Höflich 1995; Meikle et al. 1996; Patat et al. 1996), SN 2009ig (Foley et al. 2012; Marion et al. 2013; Chakradhari et al. 2019), SN 2012fr (Childress et al. 2013; Zhang et al. 2014; Contreras et al. 2018), SN 2018cnw (H1; Rigault et al. 2025), SN 2021fxy (DerKacy et al. 2023b), and SN 2021aefx (Hosseinzadeh et al. 2022; Ni et al. 2023). We considered this modest sample of six SNe Ia with four to five spectra per SN; each SN has a spectrum taken at least 11 d before maximum light. While this sample is biased, consisting of hand-picked objects from the literature, it provides a reasonable number of events with which to study the diversity of strong Si IIλ6355 HVFs and their related explosions.

We collated spectra for each object, covering the full evolution of the HV Si IIλ6355 component, aiming for a one day cadence where possible, resulting in 27 spectra for our six SNe Ia. The details of these spectra, including references, can be found in Table B.1. We supplemented the first spectrum of SN 2021fxy from DerKacy et al. (2023b), with two new spectra. The first was obtained as part of the ePESSTO+ collaboration (Smartt et al. 2015), with the ESO Faint Object Spectrograph and Camera v2 (EFOSC2; Buzzoni et al. 1984) on the European Southern Observatory’s (ESO) New Technology Telescope (NTT). It was reduced and calibrated using the pipeline described in Smartt et al. (2015). The second spectrum of SN 2021fxy was obtained with the Spectrograph for the Rapid Acquisition of Transients (SPRAT; Piascik et al. 2014) on the Liverpool Telescope (LT; Steele et al. 2004). The SPRAT spectrum was reduced using the pipeline of Barnsley et al. (2012) along with a custom Python pipeline (Prentice et al. 2018). While earlier spectra exist for SN 2012fr, as presented in Childress et al. (2013), these exhibit sole HV components without photospheric-velocity (PV) counterparts (or with very weak PV counterparts) which leads to complications with our modelling method as we are first required to create custom models to fit the PV components. Therefore, the first spectrum we model of SN 2012fr is at −11 d from peak. Earlier spectra of SN 2021aefx from Hosseinzadeh et al. (2022) were initially included for modelling; however, for these epochs we were unable to produce acceptable PV model fits. This is likely due to the domination from the HV components at very early times.

Each spectrum was flux calibrated to the available photometry, corrected for Milky Way extinction, de-redshifted, and scaled from flux to luminosity. The redshifts, distance moduli, Milky Way extinction values, time of first light and maximum light, as well as the bands used for flux calibration can be found in Table B.2. Due to the difficulties in estimating host galaxy extinction values, we did not correct for extinction arising in the host galaxies.

While most of these peak and first light dates were taken from the literature, the date of first light for both SN 2021fxy and SN 2018cnw were measured through fitting a power law to the early light-curve data. This fitting was performed on ZTF forced photometry data (Reusch 2020), with the prescription described in Miller et al. (2020) and allowing for a variable power law index. The listed phases for SN 2012fr were calculated using the maximum light date from Contreras et al. (2018), so as to use the same source for both first light and maximum light dates, and therefore differ slightly from those presented in Childress et al. (2013). For SN 1994D, no time of first light could be estimated because neither the date of first light or a measurement of the rise time could be found in the literature. As such, the initial explosion time estimate for the modelling of this object is drawn from the discovery date (MJD 49418) and then pushed backwards accordingly during the modelling. There also exists a discrepancy of one day between the modified Julian date (MJD) of the final SN 2021fxy spectrum used here, drawn from the public WISeREP archive (Yaron & Gal-Yam 2012) and the MJD presented in DerKacy et al. (2023b). With the time since explosion at ∼12 d this difference has a negligible effect upon the results presented here. The final corrected spectra are presented in Fig. 1.

Fig. 1. Spectral sequences of the six SNe Ia to be modelled. Fits to the Si IIλ6355 feature are displayed in black for each spectrum. In the case of a preferred two-component fit, the individual components are displayed in grey. The spectra are plotted in normalised luminosity and offset, as is the case with all spectra presented throughout this study.

We employed the fitting algorithm developed for the HVF search in the ZTF Cosmo DR2 (H1) to fit the Si IIλ6355 profiles for these 27 spectra. This algorithm comprised Markov chain Monte Carlo (MCMC) fitting with single and double-component models and selection between the models with the Bayesian Information Criterion. The best-fit model in each case can be seen plotted in black in Fig. 1. For each spectrum, we defined the continuum regions to the blue and to the red of the Si IIλ6355 feature manually, as automated methods may result in incorrect selections for such extreme HVFs as those seen in our sample.

In Fig. A.1 we show how these spectral series compare in terms of the component velocities, the velocity separation between the PVF and HVF (Δv), and the ratio of the HVF pseudo-equivalent widths (pEW) to the PVF pEW (R_HVF) to the HVF sample found for the ZTF Cosmo DR2. The evolution of the five newly added objects from this paper (SN 2018cnw, also known as ZTF18abauprj, was part of the ZTF Cosmo DR2 sample) generally support the findings from H1.

For the pEW ratio, R_HVF, a decrease is seen for individual objects, but no global decrease, matching the results of H1. SN 1994D is an exception; with the smallest velocity separation in our sample, the two components of SN 1994D are highly entangled, leading to a high level of degeneracy between the two components in the spectral model fits, which is in turn reflected in the measured uncertainties.

Harvey et al. (2025) saw a global decrease in velocity separation with time; however, they did not generally observe a decrease in velocity separation for individual objects. Constraints on the evolution of individual objects was limited in the H1 sample as six of the eight multi-epoch objects possessed only two spectra. In this sample we have fewer targets; however, each target has four to five spectra, allowing for tighter constraints upon the phase evolution of HVFs in individual SNe Ia. Here we see an increase in velocity separation with phase in five of the six SNe Ia, implying that the PV component velocity decreases faster than that of the HV component.

As the six SNe Ia explored in this work are well known literature targets with clear HVFs, they tend to have larger velocity separations than the majority of other HVF spectra from the ZTF Cosmo DR2 sample, and therefore are not representative of the full range of SN Ia HVFs found in nature.

3. Method

Our aim in this work was to investigate the origin of HV Si IIλ6355 features in very early SN Ia spectra. We initially found best matching models for the photospheric components of the spectra and their overall spectral shape, before investigating the impact of HV density enhancements to explain the presence of the HV Si IIλ6355 features. In Section 3.1, we firstly describe our TARDIS modelling, including our choice of input parameters, such as density and abundance profiles. We then present our methodology for introducing density enhancements in Section 3.2. In Section 3.4, we describe the use of neural network (NN) emulators combined with a MCMC to determine the TARDIS parameters for the density enhancements that best match the observed spectra.

3.1. TARDIS photospheric modelling

TARDIS¹ is an open-source radiative-transfer code that synthesises spectra from one-dimensional (1D) models of SN ejecta (Kerzendorf & Sim 2014). The five inputs to a TARDIS simulation are the time since explosion (t), the luminosity (L), the inner boundary velocity (v_ph), the density profile (ρ), and the abundance profile (A). The TARDIS modelling presented here employed the nebular treatment for ionisation, dilute-lte for excitation, dilute-blackbody for radiative rates, and the macroatom treatment (based upon the macro atom described in Lucy 2002) for line interaction (see Kerzendorf & Sim 2014 for further details). The profiles of the abundance and density are defined as functions of velocity, and divided into discrete shells. The time of explosion for each of the SNe Ia was initially taken as the time of first light (Table B.2). For each SN, the initial estimate of the explosion time was allowed to be pushed backwards where needed in the development of the photospheric models (PV models) as a control of the temperature to fine tune the match to the observed spectra. This potential difference between the explosion time and the time of first light arises from the time taken for the first photons to escape the ejecta, and is known as the dark phase (Piro & Nakar 2013). Initial estimates for the inner boundary velocities roughly followed the measurements of the photospheric velocities in the top panel of Fig. A.1. However, these were varied to improve matches where necessary.

To produce our best matching PV models, we initially removed the Si IIλ6355 HV components from our observed spectra by adding the best-fitting HV Gaussian (see Section 2) to each observed spectrum. While the mechanism forming these HV components likely has significant impacts on other regions of the spectrum, we focussed only upon the Si IIλ6355 and thus this simplistic treatment is justified.

3.1.1. Density profiles

To develop a custom model for each of the SNe Ia in our sample, we first chose a density profile. As many of the explosion models have density profiles that are very similar (e.g. see a comparison in Fig. 10 of Harvey et al. 2023), the choice of density profile here is somewhat arbitrary. Therefore, we chose a density profile of the form 10^{mv + d₀} with values of m and d₀ chosen so that the profile roughly resembles the density profile from the N100 delayed-detonation model (Seitenzahl et al. 2013) in the region of velocity, v > 10 000 km s⁻¹. Ejecta at velocities of < 10 000 km s⁻¹ is below the photosphere at the epochs considered in this work, and was therefore not constrained. The initial density profile as a function of velocity chosen for all the objects in this study is

$$\begin{array}{c}\hfill {\rho }_0(v)={10}^{mv+d_0},\end{array}$$(1)

with m = −1.93 × 10⁻⁴ and d₀ = 1.43. This density profile is defined in g cm⁻³ at 100 s post explosion. Throughout the development of the best matching PV model for each SN Ia, tweaks were made to the slope (m) and offset (d₀) of this density profile to more closely reproduce the evolution of the corresponding Si IIλ6355 profile. However, the final density profiles for any SN Ia in our sample did not differ significantly from the original form.

3.1.2. Abundance profiles

TARDIS operates with a photospheric approximation, emitting photon packets from a solid inner boundary to pass through the ejected material above. As time progresses with successive spectral observations, this photosphere recedes inwards to lower velocities, revealing previously unseen material. It is this principle that allows for the development of custom abundance profiles with the technique known as abundance tomography (Stehle et al. 2005). The first spectrum in a series is used to constrain the abundances of the outermost material, which is then kept constant in the modelling of following spectra, with each successive spectrum constraining a new shell of material below the previous photosphere.

The custom models were divided into shells separated by 250 km s⁻¹ with abundances composed of O, Mg, Si, S, Ca, Ti, Cr, Fe, and Ni. C was not included in these models for reasons that shall be discussed in 3.1.3. We populated the outermost material with elemental abundances that are uniform with velocity. The inner boundary of this layer of material varies from object to object, with a minimum value of 17 000 km s⁻¹ for SN 1994D and a maximum of 21 500 km s⁻¹ for SN 2009ig. While we did not place strict limits on the allowed abundances in this region, we aimed to follow the relative abundances seen in delayed-detonation models at these higher velocities; this means that it is dominated by unburnt material, with small amounts of intermediate-mass elements (IMEs) and little to no iron-group elements (IGEs). We initialised the model as 100% O and gradually introduce other species to fill out the various absorption features.

Each TARDIS simulation was run with 30 iterations to ensure convergence of the temperature profile, with this value being chosen after it proved to be sufficient in initial manual testing. For these initial grids, the final packet count was set to 10⁵ as this was found to produce spectra with high enough signal-to-noise (S/N) for general comparison to the observed spectra, while keeping the overall computation times reasonable; which was especially important for the ∼7000–9000 density enhancement simulations per object later on (see Section 3.3).

3.1.3. Output base PV models

Through abundance tomography we developed a PV model for each of the six SNe Ia. The synthetic spectra calculated by TARDIS can be seen in Fig. A.2. The goal of each PV model was to closely match the evolution of the Si II PVF, while also matching the overall shape of the spectra, governed principally by the temperature. As we were not focussed upon the relative abundances of the other elements, we were required to designate an element to act as a filler species to fill the normalised abundance in each shell to unity. This element should have minimal impact upon the temperature and opacity of the ejecta, and also must not produce any absorption lines in the region of the feature of interest as to avoid contamination. Therefore, we employed oxygen as this filler species, which results in overly strong O I 7777 Å features in all models. The Ca II NIR and H&K features were not a focus of our PV modelling. Once the best-fitting models were established based on the full Si IIλ6355 profile, we tested the plausibility of reproducing the HV calcium features simultaneously with the same density enhancement in the outer ejecta (see Section 5.2).

The sharp steps seen in the elemental density profiles in Figure A.2 are not strict constraints of the models, but rather artefacts of the modelling process. Regions of the ejecta were set to be uniform with velocity between the various photospheric velocities for the different spectra. Smoother, continuous curves such as those from the explosion models are more realistic, with our custom models acting as lower-resolution approximations.

As the earliest spectroscopic epochs trace the outermost ejecta, they are dominated by contributions from HVFs and it is for these earliest epochs that we find the poorest matches to our PV models. In the earliest epoch of SN 2009ig, Marion et al. (2013) identified HV components in the Fe IIλλ5018, 5169 lines, as well as the Si IIIλ4560 line, which causes a shoulder feature to the red of the Ca II H&K. These features can also be seen in the earliest epochs of SN 2021fxy and SN 2021aefx in Fig. 1.

While the PV model for SN 2021fxy matches well the overall shape of the −10.8 d, −8.0 d, and −7.0 d spectra, the simulated spectrum at the earliest epoch (−13.9 d) is overluminous in the blue compared to the observation. We see close agreement in the Si region at −13.9 d, −10.8 d, and −8.0 d, with the simulated PV feature slightly too shallow at the latest epoch.

For SN 2018cnw, while we were able to reproduce the photospheric components in the −12.7 d, −11.8 d, and −5.8 d epochs, we struggled to obtain a good match to the earliest spectrum at −15.1 d (the earliest spectrum in our sample) with our PV model. We suspect that at this early phase, the photosphere is at such high velocity that it is in close proximity to the likely density or abundance enhancement causing the HVF and this impacts the PVF strongly. As such, a good match with the PV model alone is not possible. For completeness we plot the best match we could achieve with the PV model for the −15.1 d spectrum along with the good fits to the other epochs in Fig. A.2.

As seen in Fig. A.1, the velocity evolution for SN 2012fr of the observed photospheric component is remarkably flat, and our PV model required a very gradual decrease in the inner velocity parameter to well match this evolution. While this makes for better velocity matches, our PV model produces photospheric Si IIλ6355 components that are overly strong in the −11.0 d and −9.6 d epochs as can be seen in Fig. A.2.

The density enhancements to be introduced to these models in Section 3.2 bring with them changes to the temperature (and therefore ionisation) structure of the material that may cause strong C IIλ6580 absorption and potentially degrade the match of the model to the observed Si IIλ6355 feature. As such we included no carbon in the PV models, instead filling the majority of this outermost layer with oxygen. Once we obtained a working model for the HVFs, a tuned amount of carbon could be introduced to the outer ejecta, as leftover unburnt material from the original white dwarf.

In Fig. 2 we present the species density profiles (Harvey et al. 2023) of the nine elements in our six PV models along with the profiles from the double-detonation (blue; Gronow et al. 2021) and delayed-detonation (red; Seitenzahl et al. 2013) models from the Heidelberg Supernova Model Archive (HESMA; Kromer et al. 2017). These species-density profiles were calculated as the product of the mass fraction profile for the species and the density profile of the model, and therefore can be used to compare models with different density profiles. From Fig. 2 it is clear to see that the PV models more or less populate the same region of the parameter space as the delayed-detonation models. Slight exceptions to this are the augmented levels of Ti in the models built to match SN 2009ig and SN 2012fr.

Fig. 2. Species densities for the PV models for the six SNe Ia in our sample. These species densities are calculated as the product of the abundance and density profiles and can be used to compare models with different density profiles. The faint red and faint blue profiles correspond to the HESMA delayed-detonation (Seitenzahl et al. 2013) and double-detonation (Gronow et al. 2021) models, respectively. The colours of the PV model profiles correspond to the colours in Fig. 1.

3.2. Density enhancements at high velocity

In this section, we describe how density enhancements were introduced into the best matching PV models obtained in Section 3.1 to attempt to match the HV components in the Si IIλ6355 feature. While enhanced abundances of the HVF species (e.g. Si, Ca) in the upper ejecta have previously been suggested to produce these HV components, there is little justification for abundance enhancements in these regions for individual species. The shell burning region in the double-detonation models can produce this double-peaked abundance profile structure for the IMEs, but it brings with it large amounts of IGEs in the upper ejecta, making for heavily reddened spectra, in disagreement with our sample of observed spectra. This can be seen in Fig. 2 for the double-detonation models shown in blue, where species densities of IGEs (Ti, Fe, Ni) stay much higher than the delayed-detonation models shown in red at velocities beyond ∼18 000 km s⁻¹. We performed initial testing to explore the possibility of reproducing HVFs through a Gaussian enhancement to the Si abundance relative to the other elements. Even in the regime of extreme abundance enhancements this proved insufficient to reproduce HVFs.

Previous studies have argued that enhancements in the density profile, rather than in the abundance profile, are necessary for HV formation (Mazzali et al. 2005b). During initial tests with density enhancements we observed significantly greater variation in the blue region of the synthesised Si profile than for the abundance enhancement testing, supporting the potential for density enhancements to be responsible for HVF formation as suggested in the literature. These density enhancements alter the temperature and ionisation profiles in such a way that they enable the necessary conditions for HVF formation, highlighting the critical role of density in driving the ionisation balance and spectral features of the outer layers.

The modelling work presented here therefore attempts to reproduce the HVF evolution through Gaussian enhancements to the density profile. These density enhancements are only required in the line of sight of the observer, and therefore are not required to be symmetric around the exploding star. Line-of-sight variation of such density enhancements produced by three-dimensional (3D) models may potentially explain the observed diversity found in H1. The choice of Gaussian shaped density enhancements was somewhat arbitrary. Investigation of the effect of the functional form of the density enhancement, such as using alternative profiles like Lorentzian, exponential, or more complex models, would be an interesting direction for future work and could provide deeper insight into the underlying physical processes.

In addition to our base PV models we introduced a simple Gaussian density enhancement to the base density profile ρ₀(v), giving a new density profile,

$$\begin{array}{c}\hfill \rho (v)={\rho }_0(v)+r{\rho }_0(b)exp\frac{-{(v-b)}^2}{2c^2},\end{array}$$(2)

where b is the velocity location of the centroid of the Gaussian, c is the enhancement width, and r is the amplitude of the enhancement measured as a multiple of the base density at the velocity b. Our aim was then to test the impact of various amplitudes, widths, and positions of this density enhancement on the output model spectra. The impact of the injected density enhancements upon the evolution of the Si IIλ6355 structure depends upon the fractional abundance of Si in this HV layer. As such, we also allowed this outermost Si abundance to vary within our models while remaining uniform with velocity.

3.3. Simulation grids

For each of the SNe Ia in our sample we defined a grid over which to run TARDIS simulations with the aforementioned density enhancements. These grids are four-dimensional, corresponding to the three parameters governing the density enhancement, as well as the silicon abundance in the outer ejecta. For each SN our grid spans six values for the amplitude ratio r (0.5, 1, 2, 4, 8, 16), ten widths c between 200 and 2000 km s⁻¹ in increments of 200 km s⁻¹, six outer silicon abundances X_Si from 1% to 6% of the total outer layer mass in steps of 1%, and five velocity values b separated by 1000 km s⁻¹, the range of which varies from object to object depending upon the average velocity of the object’s HV component. The choice of these values for the Gaussian parameters and Si abundance were somewhat arbitrary but chosen to roughly cover the extremes of the HV Si II components seen in our observed spectra.

Each spectral series was therefore investigated initially with a grid of 1800 models covering a range of plausible density enhancements. In Section 3.4, we introduce the use of NNs to artificially increase the resolution of the grids of the parameters of the density enhancement, using these 1800 models as the training sample.

3.4. Neural networks

Recent years have seen a number of studies implementing machine learning techniques to automate the process of spectroscopic modelling of SNe (Vogl et al. 2020; Chen et al. 2020, 2024; Kerzendorf et al. 2021; O’Brien et al. 2021; O’Brien et al. 2024; Magee et al. 2024). The general approach is to construct a NN–or a series of NNs–that can emulate the performance of a radiative-transfer code such as TARDIS in a fraction of a second; a speed increase factor of ∼10 000 from the ∼5 minutes required for per simulation. This emulator can then be used as a mathematical function in an MCMC framework to calculate the best TARDIS input parameters. An emulator technique was employed in this study as a tool to interpolate between the simulation grid points for our density enhancements at high velocity, artificially increasing the resolution of the grid.

As we were solely interested in the formation of the Si IIλ6355 PV and HV components, we restricted our investigation to the region enclosed between the two previously defined continuum regions for each spectrum (covering ∼5800–6200 Å depending upon the spectrum). Each simulation from the grid has some four-dimensional input vector corresponding to the four parameters of the grid (X_Si, r, b, c), which maps to an output spectral vector of length 75, where the region around the Si IIλ6355 in the synthetic spectrum has been linearly interpolated at 75 wavelength points. The spectral vector length of 75 was chosen so that the interpolated sections of the spectra would be sampled with wavelength separations of 4–7 Å, which was shown in H1 to be sufficient for the characterisation of Si IIλ6355 lines with HV components. For comparison, the NNs from Kerzendorf et al. (2021) mapped input vectors of 12 dimensions onto spectral vectors with 500 points covering the range 3400–7600 Å.

3.4.1. Training

We built and trained our emulator networks with TENSORFLOW (Abadi et al. 2015) and KERAS (Chollet et al. 2015). Our training sample for each network comprised the full 1800 models of the simulation grid for that SN. While the size of this training set was ∼50 times smaller than that of Kerzendorf et al. (2021), our input and output dimensions were much smaller (4 versus 12 and 75 versus 500, respectively), and the morphological diversity amongst our simulations was significantly reduced since we were only considering one spectral feature (Si IIλ6355) compared to a range of 3400–7600 Å in their work. We subsequently ran a further 180 simulations per epoch, randomly choosing interpolated input vectors from the simulation grid, which were then split in half to create the validation and test samples.

The validation sample was used throughout the training process to provide feedback on the performance of the network in the form of the mean squared error (MSE) between the predictions and the simulations. While the NNs were not explicitly trained upon the validation data, the resulting MSE evolution curves informed our choices of network hyperparameters and as such these spectra influenced the training process indirectly. The test data, however, were never seen by the NNs, enabling for an unbiased assessment of the ability of the NNs to accurately emulate the performance of TARDIS, with different density enhancements post-training.

We performed several preprocessing steps upon the training spectra and input labels before constructing our networks. Being the outputs of Monte Carlo simulations, the generated spectra possess noise which should not be captured by the NNs. As such, the training spectra were smoothed with a Savitzky–Golay filter from the scipy signal module (Virtanen et al. 2020), with a window length of 15 and a polynomial order of three. These smoothed spectra were subsequently interpolated to the 75 wavelength points spaced linearly across the Si IIλ6355 feature between the predefined continuum regions. We scaled the luminosity axis of the smoothed Si IIλ6355 feature to have the minimum point at y = 5 and the maximum at y = 10, before finally taking the logarithm with base 10. The input values were mapped to have values between 0 and 1, with these corresponding to the boundary values of the grid. For the parameter r we took the logarithm with base 2 before carrying out this mapping step.

Our initial choice of network architecture corresponded to the best architecture published in Kerzendorf et al. (2021). We subsequently tested the impact of changing layer widths and the number of layers, finding no loss of accuracy in reducing the depth of the network to a single hidden layer. Therefore, our chosen architecture for all NNs presented here consisted of two fully connected layers of width 200, with an input dimension of four (corresponding to the density enhancement parameters) and an output dimension of 75 (the emulated spectral vector). The activation function for each of these layers was chosen to be softplus, with the optimiser as nadam. A schematic view of the NN architecture is shown in Fig. 3. Similarly to Kerzendorf et al. (2021), the introduction of a dropout fraction did not improve the accuracy of the emulator predictions, and therefore was not implemented. Preliminary testing indicated that the MSE for the validation sample reached a minimum at around ∼1000 epochs, after which it gradually tended upwards. The MSE for the training sample, however, continued to decrease, while at a very slow rate. Such behaviour indicateed that after ∼1000 epochs we entered a regime of overfitting, and as such the final networks were trained for a total of 1000 training epochs, with a batch size of 4.

Fig. 3. Schematic view of the chosen NN architecture. The blue nodes on the left correspond to the four inputs governing the density enhancement and silicon abundance, which are fed into the input layer. The input layer and hidden layer are both comprised of 200 neurons with a softplus activation function and are represented by the grey nodes. Finally, the orange nodes correspond to the normalised luminosity outputs at the 75 wavelength points across the range of the Si IIλ6355 feature.

3.4.2. Performance

For each spectral epoch, we trained a NN to predict the impact of the density enhancement on the morphology of the Si IIλ6355 feature, therefore amounting to four or five networks per SN depending on the number of spectra for that object. For each of these NNs, we possessed a further 90 simulated spectra as our validation set and another 90 spectra as our test set; each corresponding to 5 per cent the size of the training set. The top row of panels in Fig. 4 displays the evolution of the MSE of the NN predictions against the validation data with training epoch. In all cases these curves follow the same evolution, with a sharp initial decrease followed by a slow decline to a plateau in the regions of ∼10⁻⁴.

The MSE curves for networks of SN 2012fr and SN 2009ig plateau at the highest values among the sample. We believe this to be a result of the grid parameter choices for these objects, as they possess the highest range of values for the density enhancement velocity location b (see Section 3.3). With the enhancement at higher velocities, the formed HVF is further detached from the photosphere and there exists more diversity among the simulation grid in terms of feature shape. This higher morphological diversity is characterised by the NNs with the same number of training spectra as for the less separated features, resulting in lower predicted accuracies. We see the converse effect for the NNs of SN 1994D, with the lowest MSEs as the corresponding grid has the lowest range of b values, reducing the diversity between the models and allowing for more accurate characterisation. Regardless of this difference, the predictions of all networks provided sufficient matches to the test data to proceed with the MCMC fitting.

All of the test set spectra from all the epochs for a given supernova were combined and ranked by the χ² difference between prediction and simulation. The median and worst-case predictions for each SN are presented in the second and third rows of Fig. 4, respectively. Once again the worst matches of the NN predictions to the simulated models of SN 2012fr and SN 2009ig were the poorest among the six SNe; however, they still described the simulated feature reasonably well.

Fig. 4. Performance summary of the 27 NNs (four or five per SN based on the number of observed spectral epochs). The evolution of the MSE is shown against the training epoch in the top panels. The middle and bottom panels show the median and worst-case predictions, respectively, from the simulated test data corresponding to each SN epoch. The grey lines are the unsmoothed TARDIS outputs and the colours lines are the NN predictions. The median and worst predictions are not necessarily at the same epoch for each SN.

3.5. MCMC

With the trained NNs reproducing the TARDIS simulation grids with a high degree of accuracy, we ran MCMC fits for each of the SNe Ia over all the corresponding observed spectral epochs using EMCEE (Foreman-Mackey et al. 2013). We imposed uniform priors upon the four input parameters between the boundaries of the grid (X_Si, r, b, and c). For a single spectral epoch, the log-likelihood function takes the form

$$\begin{array}{c}\hfill lnL(X_{\text{Si}},r,b,c)=-\frac{1}{2}{\displaystyle \sum {\left(\frac{y_{\text{obs}}-y_{\text{NN}}(X_{\text{Si}},r,b,c)}{y_{\text{obs,err}}}\right)}^2},\end{array}$$(3)

with y_obs and y_obs,err as the preprocessed luminosity values of the observed spectra and the corresponding uncertainties, respectively, and y_NN(X_Si, r, b, c) are the predictions from the relevant NNs. Having set the TARDIS simulation parameters in the development of the PV models, we reduced the dimensionality of the parameter space to the four parameters pertaining to density enhancement and the outer silicon abundance. While such assumptions enable probabilistic modelling, there was less flexibility to vary the simulated Si IIλ6355 profile and as such perfect matches of the models to the observed spectra were unlikely. To account for model mismatch we assumed conservative relative uncertainties on the observed spectra of 5%. The total log-likelihood for a single SN was calculated as the sum of these individual log-likelihoods for all spectral epochs investigated.

The initial value for all four parameters in all fits was set to be 0.5, i.e. the centre of the parameter space. We ran the MCMC fits with 32 walkers, for an initial burn-in period of 100 iterations, followed by the final chains with 1000 iterations. In four of the six cases these hyperparameter choices resulted in autocorrelation times of ∼1 indicating that the chains mixed well and explored the parameter space efficiently. The high autocorrelation times for the remaining two objects (SN 2018cnw and SN 2021fxy) appeared to be caused by the early sections of the chains. The burn-in period for these two objects was therefore increased to 200 iterations, resulting in autocorrelation times of ∼1.

As discussed in Section 3.1.3, the first epoch of SN 2018cnw at −15.1 d was not well fit by the PV model. In the MCMC fitting of the HV components of SN 2018cnw, we initially attempted to fit using all four epochs; however, due to the poor matching of the −15.1 d spectrum this produced a poor match to the Si IIλ6355 in all four epochs. Therefore, we repeated the fitting with the first epoch excluded, obtaining a close match to the remaining three spectral epochs. In the case of the first epoch of SN 2018cnw, the assumption that the density enhancement only affects the formation of the HV feature, and not the overall shape of the spectrum, breaks down and is what we suspect to be the cause of this discrepancy. Simultaneous fitting of the first epoch along with the other epochs would likely require leaving the TARDIS input parameters as free parameters, which is not computationally feasible. The results of this second fit (excluding the earliest spectrum) are used below.

4. Results

In this section, we discuss the outputs of our models for the observed Si IIλ6355 features in our sample of six SNe Ia. In particular, this focusses upon the contribution from the best-fitting HV Si IIλ6355 components. As discussed in Section 3.5, the best-fitting density enhancement parameters were calculated across all spectral epochs simultaneously since we required a single density enhancement to explain the HVFs and their evolution with time. We present the best NN fit to the Si IIλ6355 region of each SNe in Fig. 5. The original best-matching PV only model is also shown for comparison to highlight the contribution of the density enhancement to the feature.

Fig. 5. Best-fitting density enhancement model spectra (solid coloured lines) compared to observed spectra (grey) and the Si IIλ6355 region from the PV simulations (dotted). The faint regions of the model spectra correspond to the regions that are not constrained in the MCMC fitting. The asterisk by the −15.7 d spectrum of SN 2018cnw indicates that this spectrum was not used in the MCMC parameter inference (see Section 3.5).

The corner plots for the density enhancement parameters for each of the six MCMCs are shown in Fig. A.3. The values of these best-fitting Gaussian parameters for the density enhancement and the abundance of Si II at the velocity of the enhancement are given in Table 1. For SN 2018cnw, SN 2012fr and SN 2009ig, the MCMC contours are well constrained with Gaussian distributions for each of the parameters. However, for SN 2021fxy, SN 2021aefx, and SN 1994D, the MCMC contours are close to the boundaries of the simulation grid and are less well sampled, suggesting that they would prefer values that are higher or lower than the boundaries of the grid.

A comparison between the HVF and PVF velocities measured for the observed and simulated spectra is shown in Fig. 6. These measurements were performed using the code developed in H1. As in H1, we adopted a double-component fit for the simulated features only in cases where this was the preferred model, with Δv > 4000 km s⁻¹. In general, good agreement is seen suggesting that the models are broadly reproducing the data in the region of the Si IIλ6355 feature. In the sections below, we describe the main characteristics of the best-fitting models for the Si IIλ6355 feature region for each of the six SNe.

Fig. 6. Comparison of the velocities from the best-fitting models and the observed spectra. The solid symbols represent the PV components; the open symbols are their HV counterparts. The circles and triangles correspond to the observations and simulations, respectively; the crosses display the measurements from the PV model spectra in Fig. A.2.

4.1. SN 2021fxy

The early spectra and light curves of SN 2021fxy were studied in DerKacy et al. (2023b). They estimated using the blue edge of the Si IIλ6355 absorption trough of the first spectrum (−13.9 d) that the line-forming region of Si II extends to at least −28 000 km s⁻¹ but that the Si IIλ6355 feature has mostly faded by the −7.0 d spectrum.

As seen in Fig. 5, the best-fitting density enhancement model provides a good match to the Si IIλ6355 in the four spectra of SN 2021fxy. However, as in the PV model alone (see Section 3.1.3), the PV component strength is still underproduced in the last spectrum at −7.0 d. The HV component from the best-fitting model matches closely to the observed spectrum in all four epochs.

In Fig. A.3, the MCMC contours are seen to overlap the upper boundaries of the grid in the X_Si and c parameters, indicating that a more extended density enhancement combined with a higher Si abundance would be marginally preferred for SN 2021fxy. The width of the density enhancement for SN 2021fxy is already the highest for our sample and the Si abundance is the joint largest with SN 1994D. The velocities measured using Gaussian fits (as described in Section 2) to the best-fitting models and the observed spectra are in excellent agreement for SN 2021fxy, again suggesting that the model is a good match to the data (Fig. 6).

4.2. SN 2018cnw

Identified as having a HV component in H1 in three of the four available epochs, SN 2018cnw exhibits the largest velocity separation measured in the ZTF Cosmology DR2. With four epochs from −15.1 to −5.8 d covering ∼9 d (the largest phase range in this sample) the HV component is seen to fade away completely to leave the solitary photospheric component.

As discussed in Section 3.5, the earliest spectrum of SN 2018cnw was excluded from the MCMC fitting. As expected there is therefore a significant difference seen in the velocities obtained from the Gaussian fits to the observed and best-fitting models for this first epoch in Fig. 6. However, the velocities measured from the following epochs are in good agreement with the fits to the observed data. The reproduction of the Si IIλ6355 structure in these later three epochs shown in Fig. 5 are seen to be good matches and the MCMC contours are also well constrained (Fig. A.3).

4.3. SN 2012fr

With the largest velocity separation in the sample, and higher than any object found in the ZTF DR2 (H1), SN 2012fr represents the extreme end of the HVF population. An extensive observational study of the HV features of Si II and Ca II was first presented in Childress et al. (2013). The distinct nature of the HVF for SN 2012fr was also confirmed in other studies (Zhang et al. 2014; Contreras et al. 2018). For the spectra investigated in this work, we measured the HV component to have a mean separation of 8300 km s⁻¹ from its PV counterpart, making for a double-component profile that can be easily identified by eye. Our spectral series for SN 2012fr spans ∼5 days from −11 to −6 days with respect to peak. These spectra demonstrate the full evolution of a HVF, from starting off the dominant component, rapidly weakening to be of equal strength as the PV component, and then to fading away completely. We find consistent Si IIλ6355 HV and PV component velocities compared to the values of Childress et al. (2013) of ∼22 000 to 20 000 km s⁻¹ for the HVF and ∼13 000 to 12 000 km s⁻¹ for the PVF, over the same phase range as studied here.

In Fig. 7, we show the fractional abundance of singly ionised silicon (Si⁺) for both the PV model and the density enhancement model across multiple epochs. The introduction of the density enhancement leads to elevated temperatures throughout the ejecta, up to and including, the region of the enhancement. These higher temperatures suppress recombination, resulting in lower Si⁺ abundances around ∼20 000 km s⁻¹, just below the core of the density enhancement. At higher velocities, within the enhancement itself, the increased electron density drives a sharp rise in the recombination rate, leading to significantly higher Si⁺ abundances in the density enhancement model. This region of enhanced Si⁺ is responsible for producing the HVFs observed in the synthetic spectra.

Fig. 7. Top: Density profiles from the PV model (grey) and the density enhancement model (green) for SN 2012fr. Bottom: Relative fraction of silicon material in the singly ionised state for the PV model (grey) and the density enhancement model (green) across the first four epochs of SN 2012fr. The black vertical lines correspond to the measured HVF velocities from the synthesised spectra. The shaded region represents the velocity range centred at b with width 2c (i.e. from b − c to b + c), where b and c are the centroid and width of the Gaussian density enhancement, respectively, as given in Table 1.

In our best-fitting density enhancement model, we see the velocity of the PV component pushed slightly higher compared to that in the PV model, resulting in the small offset we see in Fig. 6, where the PV velocities measured from the best-fitting models are slightly higher than the observed spectra. We similarly see a velocity difference between the observed and simulated HV components.

Our best-fitting model with a density enhancement for 2012fr closely matches the observed HVF evolution of the object (Fig. 5), and consists of a low outer abundance of Si (the lowest in our sample) with a strong (highest in our sample) and wide density enhancement at a large separation from the photosphere. This result implies that the highly separated HVFs find their origins in fairly extreme density enhancements, which is reflected in the fact that they are so rare.

Being such an extreme object, we also find the strongest parameter correlations in the fitting parameters posteriors. As seen in the corresponding panels of Fig A.3, we find the outer silicon abundance (X_Si) to correlate positively with the enhancement centroid velocity (b) and negatively with the enhancement strength (r). As the density enhancement moved out to higher velocities where there is less material, we required a higher concentration of silicon to form the HVF. Conversely, as the concentration of silicon increased, we required the density enhancements to be less extreme in amplitude.

4.4. SN 2021aefx

SN 2021aefx has been the focus of a number of studies at both early (Ashall et al. 2022; Hosseinzadeh et al. 2022; Ni et al. 2023) and late times, including extensively with James Webb Space Telescope (JWST), in both imaging (Chen et al. 2023) and spectroscopy (Kwok et al. 2023; DerKacy et al. 2023a; Blondin et al. 2023; Ashall et al. 2024). A prominent HVF of the Si IIλ6355 was identified in the earliest spectra and the feature was seen to evolve rapidly with time.

The posterior contours overlap the lower grid boundary for the width of the Gaussian density enhancement parameter c, resulting in a density enhancement confined to a very narrow range of velocities; significantly narrower than the other events. While the best density enhancement model fills out the two-component Si IIλ6355 feature well (Fig. 5), the specific shape of the absorption structure is not well matched, with the MCMC fit tending towards a more entangled PV-HV pair that resembles more of a single broad component, compared to the fairly distinct components seen in Fig. 1. This results in a two-component fit being preferred in only the first epoch for the density enhancement model, as opposed to the first three epochs for the observations. This can be seen in the SN 2021aefx panel of Fig. 6 for which we only possess a HVF velocity measurement at −13.8 d.

4.5. SN 2009ig

We measured the highest HVF velocities in our sample for SN 2009ig, in agreement with high velocities measured in previous studies of this SN (Foley et al. 2012; Marion et al. 2013; Chakradhari et al. 2019). However, as the PV component velocities are also at the higher end of the range, we do not see as extreme a velocity separation as we saw for SN 2012fr (Fig A.1). Due to these high velocities, the best density enhancement model identified by the MCMC fitting has the highest density enhancement velocity (b) among all the objects investigated.

While the first epoch of SN 2009ig is at a similar phase to the first epoch of SN 2018cnw and has a similar inner boundary velocity, we did not encounter the same issues with regards to matching the Si IIλ6355 velocity in the PV model, nor did we require the exclusion of the first spectral epoch in the MCMC fitting. This is likely due to the density enhancement being at such high velocities, causing less disruption to the rest of the spectrum. Although the density enhancement fills out the velocity extent of the feature, the fitting of the simulated features with the HVF classification code prefers a double-component classification only in the earliest epoch. This can be seen in Fig. 6. The single-component classification for the simulated features at the remaining four epochs indicates that the simulated line profiles closely resemble single broad Gaussians.

4.6. SN 1994D

Patat et al. (1996) considered the Si IIλ6355 of SN 1994D as a single component in their spectra but noted that there was a very steep decline in the velocities across the earliest spectra (−11 to −8 d from peak), which we now interpret as the phases where an additional component from a HVF is present. SN 1994D exhibits the smallest velocity separation in our sample, appearing to bridge the gap between broad-line SNe Ia and those with clearly resolved HV components. With a mean HVF velocity of only ∼18 000 km s⁻¹ over the epochs investigated, we extended the uniform abundance HV shell in our models inwards to 17 000 km s⁻¹. The velocity range for the density enhancement chosen for the simulation grid ran from 18 000–22 000 km s⁻¹.

While the best-fitting density enhancement model matches well the evolution of the SN 1994D Si IIλ6355 feature in Fig. 5, the MCMC results cluster up against the boundaries of the grid. The parameter contours in Fig. A.3 overlap the upper limit on the outer silicon abundance (X_Si), as well as the lower limit on the density enhancement strength (r). In contrast to the extreme density enhancements required for the highly separated HVFs, higher silicon abundances with very slight density enhancements appear to give rise to the less separated HVFs.

These small density enhancements are of the same scale as the general variation we see in the angle-averaged density profiles from the simulations from the Heidelberg Supernova Model Archive (HESMA), and are therefore likely to be very common. The regularity of such density enhancements agrees with the fact that these smaller velocity separations were found to be the most common in H1.

In fitting the simulated Si IIλ6355 features, a single-component model was preferred at the earliest epoch, with the simulated line profile more Gaussian in shape than in the observed spectrum. In the two following epochs the classification code prefers a double-component model, as seen for the observed spectra.

5. Discussion

In this section, we discuss the combined results of our TARDIS modelling of six SNe Ia showing HV components in their Si IIλ6355 features in very early time spectra (−14 to −6 d with respect to maximum light). In Section 5.1, we discuss the properties of the density enhancements required in our study to explain the diversity of the observed Si IIλ6355 HV features. In Section 5.2, we altered the outer Ca abundance of our best-fitting models to test whether the Ca HVFs can be simultaneously formed by the same density enhancement and in Section 5.3, we discuss our results in the context of the hydrodynamical delayed-detonation (Seitenzahl et al. 2013) and double-detonation models of Gronow et al. (2021).

5.1. Density enhancement morphology and origin

In Fig. 8, we show the density profiles for each of our six SNe with the density enhancement that best matches the observed HV Si II features. There is seen to be significant variation in the positions of the density enhancements, their widths, and their strengths. As shown in Table 1, the centroid of the enhancement varies from ∼18 700 km s⁻¹ for SN 1994D to 26 100 km s⁻¹ for SN 2009ig. The centroid of the enhancement is seen to broadly correlate with the measured HV component velocities shown in Fig. A.1, where SN 1994D had the lowest HV component velocities and SN 2009ig had the highest.

Fig. 8. Best-fitting enhanced density profiles for each of the six supernovae. The PV model density profiles for the targets are shown as the coloured dotted lines; the grey solid and dashed lines correspond respectively to the angle-averaged density profiles from the delayed-detonation and double-detonation simulations from HESMA.

The strength of the density enhancement also varies between the SNe Ia in the sample, with the strongest density enhancements typically required to produce the more extreme velocity separations. As the highly detached HV components are formed far above the photosphere, they have densities far lower than the photospheric density and as such require these fairly extreme enhancements to produce line-forming regions. The converse affect can be seen for SN 1994D where the HV line-forming region has densities much closer to that of the photospheric region and as a result only a slight density enhancement is required. The stronger density enhancements are likely far less prevalent in nature, and therefore explain the scarcity of events such as SN 2012fr.

These density enhancements constitute additional material introduced atop the base density profile, and a mechanism is therefore required to provide this additional mass. Firstly, we consider the extreme-velocity material (above the density enhancement region) within the supernovae and whether these enhancements may be the bunching up of this faster moving material, relocating it to slower velocities. While potentially plausible for the weakest undulations – as seen in SN 1994D – the majority of the enhancements found here invoke masses significantly larger than the masses found at velocities above the density enhancement in literature models. This implies that these enhancements are not the result of the compacting of faster moving material.

Another potential candidate for the source of this additional mass is CSM surrounding the system. This material, however, moves at non-explosive velocities and we would therefore require some mechanism to accelerate the CSM to the velocities of the density enhancements (∼20 000–30 000 km s⁻¹). CSM is therefore also unlikely to provide this additional mass.

The most promising candidate for this additional mass is likely to be the supernova itself. Whether linked to a particular explosion mechanism, or produced through some asymmetry of the progenitor system, differences in the distribution of kinetic energy from the explosion could impact the distribution of material at high velocities. We therefore propose these density enhancements to be generated by the supernova explosions themselves.

5.2. Calcium HVFs

With best-fitting density enhancement models for the Si IIλ6355 HVFs determined, we turned our focus to the HV components of the Ca II H&K and NIR features. To test if the current enhanced density profiles can simultaneously reproduce the Si IIλ6355 and Ca II HVFs, we ran a small grid of 20 models for each SN Ia varying the outer calcium abundance – in the same region as for the outer silicon abundance – over the range 0.001–5% with a logarithmic spacing. Increasing the calcium abundance yielded stronger Ca II H&K and NIR features; however, the Ca II features formed at velocities far lower than the observed Ca II HV components. The synthetic spectra from the models with highest calcium abundances (5%) can be seen in Fig. 9 over the H&K and NIR regions for all epochs with wavelength coverage in at least one of these regions in the corresponding observed spectrum.

Fig. 9. Ca II H&K and NIR features synthesised by the best-fitting density enhancement profiles with Ca abundances as defined in the PV models (dotted) and with 5% outer calcium (solid). The coloured vertical lines correspond to the NIR wavelengths at the inner velocity boundary of the outer ejecta shell.

As is clear from Fig. 9, we get strong Ca II NIR formation at the inner edge of this HV region in which we increased the Ca abundance (indicated by the vertical coloured lines). In our observed spectra we do not find Ca II NIR components with velocities this low, and therefore require significantly less calcium around the lower boundary of this region, while retaining the silicon abundance constrained by the MCMC fitting. In some of the earliest epochs in the simulated spectra, we also see some higher-velocity Ca II NIR absorption coming from the density enhancement region. However, these sit thousands of km s⁻¹ lower than the observed Ca II NIR HV components. Therefore, we require low Ca abundances in the vicinity of the Si IIλ6355 density enhancement and suggest that the observed Ca II NIR HVFs form from a secondary density enhancement higher up in the ejecta.

From this we speculate on the potential for three line-forming regions: the photosphere, the density enhancement for the Si IIλ6355 HVFs (Si-dominated density enhancement), and another higher-velocity density enhancement for the Ca II HVFs (Ca-dominated density enhancement). As we have never observed any triple Si IIλ6355 components (with the third component aligned with the Ca II HVF) we postulate that there is a drop off in the silicon abundance in the outermost ejecta where we find this density enhancement. We therefore would expect the ratio of Ca/Si to decrease as we move from high velocities into the lower-velocity material where the Si becomes more dominant.

If this is the case, we expect photospheric components of Ca and Si to follow very similar velocity evolutions, with the HVFs of the two species separated by several thousand km s⁻¹. Depending upon the Ca abundance in the vicinity of the Si-dominated density enhancement, we might see a third component of the Ca II NIR sitting between the typically seen PV and HV components, tracing the evolution of the Si IIλ6355 HVF. This idea of these three line-forming regions is visualised in Fig. 10 where we plot the velocity measurements for the PV and HV components of the Si IIλ6355 and Ca II NIR features of SN 2021aefx and SN 2012fr from Ni et al. (2023) and Silverman et al. (2015) respectively. We use the velocities from these studies because they cover a broader phase range than our measurements, as well as having measurements of the Ca II NIR triplet. These literature measurements are in close agreement with the velocity measurements made in this study.

Fig. 10. Velocity measurements of the PV and HV components of the Si IIλ6355 (orange circles) and Ca II NIR (blue triangles) features in SN 2021aefx and SN 2012fr from Ni et al. (2023) and Silverman et al. (2015), respectively. The open points are components classified as HVFs in the literature; the solid points are those classified as PVFs. The lines join together points that we propose to be forming in the same line-forming regions, labelled as the photosphere, Si-dominated density enhancement (DE), and Ca-dominated DE. The phase ranges covered by the modelling in this paper are shaded in grey. The phases for the SN 2012fr measurements are given relative to the date of maximum used throughout this work, and therefore differ slightly from the values provided in Silverman et al. (2015).

The Ca II NIR photospheric feature in SN 2021aefx identified by Ni et al. (2023) (solid blue triangles) appears to closely trace the Si IIλ6355 HVF up until −12 d at which point it drops by ∼8000 km s⁻¹ in two days to then follow the evolution of the Si IIλ6355 PVF. We instead propose the presence of three distinct Ca II NIR components: one at the photosphere which appears at ∼ − 10 d, a second formed by the Si-dominated density enhancement derived in the MCMC fitting above, and a third coming from a Ca-dominated density enhancement (corresponding to the HV component identified in Ni et al. (2023)). The lack of evidence of a third silicon component implies a drop off in the silicon abundance at the highest velocities.

SN 2012fr similarly shows a late-forming photospheric Ca II NIR component that aligns in velocity space with the photospheric Si IIλ6355, and a HV Ca II NIR component several thousand km s⁻¹ faster than the Si IIλ6355 HVF. However, Silverman et al. (2015) do not identify a Ca II NIR component that traces the evolution of the Si IIλ6355 HVF, implying a lower calcium abundance for this object in the region of the Si-dominated density enhancement.

5.3. Comparison to hydrodynamical models

Based on our argument in Section 5.1 that the density enhancements arise from the explosions themselves, and our hypothesis of three line-forming regions described in Section 5.2, we investigated if some common explosion models could provide plausible matches in the appropriate regions of the outer ejecta.

Delayed-detonations exhibit a single explosion site in the core, producing a high concentration of IGEs in the core transitioning through the IMEs to unburnt material moving out through the ejecta. The Ca abundance profile peaks at lower velocities than Si, however with significant overlap, resulting in similar PVF velocities for the two species. As seen, however, the observed Ca HVFs are significantly faster than their Si counterparts, which is difficult to explain with the Ca peaking at lower velocities as in this model. The double-detonations instead have two explosion sites, in the core and in the shell. As such these models produce this IGE to IME gradient from the core outwards, and from the shell inwards, resulting in double-peaked abundance profiles, with the Ca outer peak sitting at higher velocities than that of Si. This behaviour can be seen in Fig. 11 where we plot the angle-averaged ratio of the abundances of Ca to Si for delayed-detonation and double-detonation explosion models from Seitenzahl et al. (2013) and Gronow et al. (2021), respectively. The two explosion mechanisms largely overlap at lower velocities; however, they exhibit large separation as we move into the outer ejecta. The prediction drawn here of a low Ca/Si ratio increasing towards higher velocities generally matches the abundance ratios produced by the double-detonation models. The Si IIλ6355 and Ca II NIR lines have different conditions for line formation, as such we emphasise that while these modelling results imply this increasing Ca/Si ratio, specific modelling of the HVFs of both features simultaneously would be required to constrain the amount of silicon that could be present in the region of the Ca II HVF and thus the range of the Ca/Si ratio.

Fig. 11. Ca/Si ratios as a function of velocity of the outer ejecta for the delayed-detonation explosion models of Seitenzahl et al. (2013) shown in pink and the double-detonation explosion models of Gronow et al. (2021) shown in orange.

From the intrinsic asymmetry of the double-detonation mechanism, we would expect variation in the density profiles from these models from various lines of sight. In Fig. 12 we present in grey the line-of-sight density profiles for 100 different lines of sight for five 3D double-detonation models (M0803, M0905, M1005, M1010 and M1105) from Gronow et al. (2021) with varying core and shell masses. These models show density enhancements of varying strengths at velocities ∼10 000–15 000 km s⁻¹, with some variation seen between the lines of sign and between the models. However, these density enhancements lie at velocities far below those inferred for our six SNe from the modelling performed here of ∼19 000–26 000 km s⁻¹.

Fig. 12. Density profiles from 100 lines of sight for five double-detonation models (grey). In pink are the same profiles scaled up to have 300% of the kinetic energy. The vertical coloured lines correspond to the peaks of the density enhancements (b values) derived from the modelling for our six SNe; the colours match those for these objects used throughout this work.

To investigate what would be required for a closer match with the data, we scaled these density profiles to new kinetic energies using the following equations of Hachinger et al. (2009) and Ashall et al. (2016):

$$\begin{array}{c}\hfill \rho \prime ={\rho }_0{\left(\frac{E\prime }{E_0}\right)}^{-\frac{3}{2}}{\left(\frac{M\prime }{M_0}\right)}^{\frac{5}{2}},\end{array}$$(4)

$$\begin{array}{c}\hfill v\prime =v_0{\left(\frac{E\prime }{E_0}\right)}^{\frac{1}{2}}{\left(\frac{M\prime }{M_0}\right)}^{-\frac{1}{2}}.\end{array}$$(5)

Here ρ is the density profile, E is the kinetic energy, v is the velocity profile, and M is the mass (taken here to be constant). The subscript 0 denotes the initial values with the scaled values labelled with dashes. The pink lines in Fig. 12 display these line-of-sight density profiles with 300% the kinetic energy of the base profiles. While the density enhancements from these scaled profiles align roughly with the locations of the inferred density enhancements from our modelling of each SNe Ia (vertical coloured lines in Fig. 12), the requirement of a 300% increase in kinetic energy is unrealistic and leads us to conclude that the double-detonation mechanism alone cannot explain the origin of HVFs.

Due to the more symmetric nature of the delayed-detonation mechanism, we find far less variation in the density profiles from different lines of sight, with the angle-averaged profiles representing well the general density profile. We do not see such clear density bumps in these delayed-detonation angle-averaged density profiles and as such the delayed-detonation mechanism alone is also insufficient to explain the origin of HVFs.

5.4. The impacts of physics treatments

As outlined in Section 3.1, our simulations employ the nebular ionisation and dilute-lte excitation approximations, as implemented in TARDIS. These treatments approximate the complex physical processes governing ionisation and excitation in the SN ejecta, reducing computation times and enabling parameter space exploration. Since they introduce some limitations, it is important to assess whether they significantly affect our key conclusion: that a density enhancement in the outer ejecta is required to reproduce the observed HVFs.

5.4.1. Excitation

The dilute-lte excitation approximation assumes that level populations follow Boltzmann distributions at a local radiation temperature, scaled by a dilution factor, which is updated as part of the iterative MCMC fitting procedure. Kerzendorf & Sim (2014) compared the dilute-lte treatment in TARDIS to full non-local thermodynamic equilibrium (NLTE) calculations by examining departure coefficients, which quantify deviations of level populations from local thermodynamic equilibrium (LTE). They found that for the levels responsible for the Si IIλ6355 feature, the difference between the departure coefficients from dilute-lte and the NLTE calculations was very small in the region around the photosphere, resulting in little change to the strength or shape of the formed line.

As we move higher into the ejecta and the densities drop, the approximation of LTE begins to break down as the radiation field decouples from the plasma. This is reflected by the marginal deviation of the departure coefficients from the two treatments several thousand km s⁻¹ above the photosphere. While the level populations in the region of the photosphere are relatively unaffected by the choice of excitation treatment, this becomes slightly more important in the HVF region. The differences would primarily affect the precise shape and depth of the Si IIλ6355 feature rather than eliminating it altogether. In particular, a full-NLTE excitation treatment may slightly alter the absorption strength or velocity width of the HVF, which would in turn shift the optimal values of the density enhancement parameters. These changes, however, are expected to be modest, and therefore the overall conclusion – that a density enhancement in the outer ejecta is required to reproduce the observed HVF – remains robust even under a more complete excitation treatment.

5.4.2. Ionisation

The nebular ionisation treatment is based on the modified nebular approximation of Mazzali & Lucy (1993). It provides a computationally inexpensive approximation to NLTE ionisation without solving the full system of statistical equilibrium equations.

To assess its accuracy, we refer to the detailed radiative-transfer code comparison study (Blondin et al. 2022), which included TARDIS (using the nebular ionisation and dilute-lte excitation modes) and the 1D full-NLTE code CMFGEN (Hillier & Miller 1998; Hillier & Dessart 2012). Among other metrics, they compared the ionisation fractions of various elements at multiple epochs. In Fig. 13, we present the ratio of singly ionised silicon predicted by TARDIS to that of CMFGEN from a low-luminosity (∼0.1 M_⊙ of ⁵⁶Ni) Chandrasekhar mass delayed-detonation model, which serves as a proxy for the deviation from full-NLTE ionisation. The times at which the CMFGEN calculations were performed do not align perfectly with those of the TARDIS simulations, and we therefore take the closest epoch in each case. We see larger deviations from NLTE ionisation balance closer to the photosphere, with this discrepancy becoming more pronounced at later epochs. In the simulations up to 15 d post explosion, we find close agreement (within 1%) at velocities above 20 000 km s⁻¹ (the region of the density enhancements). This similarly holds true for the other low-luminosity model investigated.

Fig. 13. Ratio of the relative fraction of singly ionised silicon between the TARDIS and NLTE CMFGEN simulations of the low-luminosity delayed-detonation model DDC25 (Blondin et al. 2013) as a function of velocity. The phases in brackets refer to the epochs of the CMFGEN models, while those outside the brackets correspond to the TARDIS times since explosion.

The same ionisation balance comparison for higher-luminosity models (⁵⁶Ni masses of ∼0.6 M_⊙) yields a different result, however. In these models we find that the TARDIS ionisation balance differs greatly from the NLTE calculations. This deviation is likely in part due to the nebular treatment, as at early phases the densities are many orders of magnitude above the regime of nebular theory and the ejecta is optically thick to the ionising UV photons. The photospheric approximation within TARDIS as the approximation of an opaque photosphere may also contribute to this deviation as it makes TARDIS agnostic to the presence of radioactive material deeper in the ejecta, which may affect the radiation field in ways not captured by the blackbody boundary condition.

We therefore caution that the nebular ionisation treatment may not be sufficiently accurate to determine whether density enhancements are required for the formation of HVFs. We strongly encourage specific modelling of Si IIλ6355 and Ca II HVFs in full-NLTE to assess whether such enhancements are indeed necessary, or if NLTE effects alone can produce sufficient quantities of detached singly ionised material to reproduce the features observed here.

6. Conclusion

Modelling of HVFs in the Si IIλ6355 feature for six well-studied literature SNe has been performed with the spectral-synthesis code TARDIS by way of introducing Gaussian density enhancements in the upper ejecta. A simplistic custom abundance profile was derived for each object through abundance tomography to match the PV features, with a HV layer with uniform abundances. A grid of simulations with various density enhancements in the HV layer were then performed, and NNs trained upon the Si IIλ6355 structure arising from the resulting grids. These networks were then used in an MCMC framework to infer the best-fitting density enhancement for each target.

The extreme separations seen in the features of SN 2012fr were matched closely, with the quality of the fitting tending to decrease as the separation decreased and the features became more entangled. This is likely due to the fact that in cases of large velocity separation, the density enhancement sits much farther from the photosphere, and therefore the impact of the density enhancement upon the PV components and overall spectrum is smaller, making the base PV model more accurate. The SNe with larger velocity separations required stronger density enhancements, whereas the smaller separations were the result of weaker enhancements. We suggest that these extreme enhancements are harder to create as an explanation for the scarcity of HVFs with large separations from their PV counterparts.

By varying the Ca abundance in this outer layer of material we investigated the possibility of forming the Si and Ca HVFs from a single density enhancement. As expected the velocities of the resulting Ca HVFs are too low, leading us to speculate that there may be a secondary density enhancement higher up in the ejecta. This would make for three distinct line-forming regions: the photosphere and then two density enhancements. We never see a third Si component formed by this second density enhancement, and only occasionally do we identify Ca components forming at the first density enhancement, suggesting that we require small Ca abundances compared to silicon at the photosphere which increases into the outer ejecta as the Si abundance drops. This corresponds to an increasing Ca/Si ratio with increasing velocity. This pattern in Ca/Si can be seen in the models of the double-detonation mechanism, but not in the delayed-detonation mechanism. Further investigation of this ratio would be interesting, with simultaneous modelling of Si IIλ6355 and Ca II HVFs being required to place meaningful constraints.

While we find the density profiles of various double-detonation models from various lines of sight to produce a range of density enhancements, these occur at velocities far below those derived in this work. To match the velocities, the double-detonation models would require three times the kinetic energy, and therefore the double-detonation mechanism alone cannot explain the HVFs. Delayed detonations of Chandrasekhar mass models also fail to reproduce this HVFs. Further work is required both to obtain and model a larger sample of well-studied SNe Ia and to explore the outer regions of the ejecta in further multi-dimensional hydrodynamical explosion models.

Currently there is no obviously favoured explosion mechanism that can describe the formation of HVFs. Based on the findings of this work, future spectroscopic modelling studies should aim to simultaneously reproduce the evolution of Si II and Ca II HVFs. Further investigation will allow us to constrain the absolute and relative abundances of these two key species, providing tighter constraints upon the composition of the outer ejecta. From an explosion modelling perspective, greater emphasis should be placed on understanding the HV ejecta and the variation we see between different models and mechanisms in the HV regime. A comprehensive mechanism capable of explaining HVF formation is essential to account for the wide range of targets in which they are observed. Updated modelling including NLTE effects as detailed in Section 5.4 are also important to confirm our results. Furthermore, from a data collection standpoint, early-time spectra are critical for the rapid identification of HVFs, providing the necessary foundation for subsequent spectroscopic follow-up and detailed modelling.

Acknowledgments

The research conducted in this publication was funded by the Irish Research Council under grant number GOIPG/2020/1387. K.M. acknowledges support from EU H2020 ERC grant no. 758638. K.M. and T.E.M.B. are supported by Horizon Europe grant no. 101125877. A.H. acknowledges support by the Klaus Tschira Foundation. A.H. is a Fellow of the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). J.P.A. acknowledges funding by ANID, Millennium Science Initiative, ICN12_009. T.-W.C. acknowledges the Yushan Fellow Program by the Ministry of Education, Taiwan for the financial support (MOE-111-YSFMS-0008-001-P1). Based on observations collected at the European Southern Observatory under ESO programme 106.216C.006/012. This work made use of the Heidelberg Supernova Model Archive (HESMA), https://hesma.h-its.org.

References

Appendix A: Figures

Fig. A.1. Comparison of the HVF SNe Ia sample collected here and the ZTF Cosmo DR2 HVF sample (grey) from H1. Here we plot the velocity evolution of the PV (solid) and HV (open) components (top), the separation between them in velocity space (middle), and the ratio of their pEWs (bottom). The same code has been used to fit both samples. The colours for each object match those in Fig. 1 and used throughout. The solid and dot-dashed lines in the velocity panel correspond respectively to the fits to the PV and HV components from the DR2 sample (H1).

Fig. A.2. Photospheric feature models for our six SNe Ia. The synthetic model spectra are shown in colour for each object with the observed spectra from Fig. 1 in grey. To aid comparison between the PV models and the observed Si IIλ6355 PV components, the HV components have been removed from the observed spectra by adding the best-fitting HV component Gaussian fit to the spectrum.

Fig. A.3. Corner plots of the parameter posteriors inferred by the density enhancement MCMC runs for the six SNe. The grid boundaries are indicated as grey shaded regions in panels where they are visible. The solid and dotted lines in the histogram panels correspond to the means and standard deviations, respectively.

Appendix B: Tables

All Tables

All Figures

	Fig. 1. Spectral sequences of the six SNe Ia to be modelled. Fits to the Si IIλ6355 feature are displayed in black for each spectrum. In the case of a preferred two-component fit, the individual components are displayed in grey. The spectra are plotted in normalised luminosity and offset, as is the case with all spectra presented throughout this study.
In the text

Fig. 2. Species densities for the PV models for the six SNe Ia in our sample. These species densities are calculated as the product of the abundance and density profiles and can be used to compare models with different density profiles. The faint red and faint blue profiles correspond to the HESMA delayed-detonation (Seitenzahl et al. 2013) and double-detonation (Gronow et al. 2021) models, respectively. The colours of the PV model profiles correspond to the colours in Fig. 1.

In the text

Fig. 3. Schematic view of the chosen NN architecture. The blue nodes on the left correspond to the four inputs governing the density enhancement and silicon abundance, which are fed into the input layer. The input layer and hidden layer are both comprised of 200 neurons with a softplus activation function and are represented by the grey nodes. Finally, the orange nodes correspond to the normalised luminosity outputs at the 75 wavelength points across the range of the Si IIλ6355 feature.

In the text

Fig. 4. Performance summary of the 27 NNs (four or five per SN based on the number of observed spectral epochs). The evolution of the MSE is shown against the training epoch in the top panels. The middle and bottom panels show the median and worst-case predictions, respectively, from the simulated test data corresponding to each SN epoch. The grey lines are the unsmoothed TARDIS outputs and the colours lines are the NN predictions. The median and worst predictions are not necessarily at the same epoch for each SN.

In the text

Fig. 5. Best-fitting density enhancement model spectra (solid coloured lines) compared to observed spectra (grey) and the Si IIλ6355 region from the PV simulations (dotted). The faint regions of the model spectra correspond to the regions that are not constrained in the MCMC fitting. The asterisk by the −15.7 d spectrum of SN 2018cnw indicates that this spectrum was not used in the MCMC parameter inference (see Section 3.5).

In the text

	Fig. 6. Comparison of the velocities from the best-fitting models and the observed spectra. The solid symbols represent the PV components; the open symbols are their HV counterparts. The circles and triangles correspond to the observations and simulations, respectively; the crosses display the measurements from the PV model spectra in Fig. A.2.
In the text

Fig. 7. Top: Density profiles from the PV model (grey) and the density enhancement model (green) for SN 2012fr. Bottom: Relative fraction of silicon material in the singly ionised state for the PV model (grey) and the density enhancement model (green) across the first four epochs of SN 2012fr. The black vertical lines correspond to the measured HVF velocities from the synthesised spectra. The shaded region represents the velocity range centred at b with width 2c (i.e. from b − c to b + c), where b and c are the centroid and width of the Gaussian density enhancement, respectively, as given in Table 1.

In the text

	Fig. 8. Best-fitting enhanced density profiles for each of the six supernovae. The PV model density profiles for the targets are shown as the coloured dotted lines; the grey solid and dashed lines correspond respectively to the angle-averaged density profiles from the delayed-detonation and double-detonation simulations from HESMA.
In the text

	Fig. 9. Ca II H&K and NIR features synthesised by the best-fitting density enhancement profiles with Ca abundances as defined in the PV models (dotted) and with 5% outer calcium (solid). The coloured vertical lines correspond to the NIR wavelengths at the inner velocity boundary of the outer ejecta shell.
In the text

Fig. 10. Velocity measurements of the PV and HV components of the Si IIλ6355 (orange circles) and Ca II NIR (blue triangles) features in SN 2021aefx and SN 2012fr from Ni et al. (2023) and Silverman et al. (2015), respectively. The open points are components classified as HVFs in the literature; the solid points are those classified as PVFs. The lines join together points that we propose to be forming in the same line-forming regions, labelled as the photosphere, Si-dominated density enhancement (DE), and Ca-dominated DE. The phase ranges covered by the modelling in this paper are shaded in grey. The phases for the SN 2012fr measurements are given relative to the date of maximum used throughout this work, and therefore differ slightly from the values provided in Silverman et al. (2015).

In the text

	Fig. 11. Ca/Si ratios as a function of velocity of the outer ejecta for the delayed-detonation explosion models of Seitenzahl et al. (2013) shown in pink and the double-detonation explosion models of Gronow et al. (2021) shown in orange.
In the text

	Fig. 12. Density profiles from 100 lines of sight for five double-detonation models (grey). In pink are the same profiles scaled up to have 300% of the kinetic energy. The vertical coloured lines correspond to the peaks of the density enhancements (b values) derived from the modelling for our six SNe; the colours match those for these objects used throughout this work.
In the text

	Fig. 13. Ratio of the relative fraction of singly ionised silicon between the TARDIS and NLTE CMFGEN simulations of the low-luminosity delayed-detonation model DDC25 (Blondin et al. 2013) as a function of velocity. The phases in brackets refer to the epochs of the CMFGEN models, while those outside the brackets correspond to the TARDIS times since explosion.
In the text

Fig. A.1. Comparison of the HVF SNe Ia sample collected here and the ZTF Cosmo DR2 HVF sample (grey) from H1. Here we plot the velocity evolution of the PV (solid) and HV (open) components (top), the separation between them in velocity space (middle), and the ratio of their pEWs (bottom). The same code has been used to fit both samples. The colours for each object match those in Fig. 1 and used throughout. The solid and dot-dashed lines in the velocity panel correspond respectively to the fits to the PV and HV components from the DR2 sample (H1).

In the text

	Fig. A.2. Photospheric feature models for our six SNe Ia. The synthetic model spectra are shown in colour for each object with the observed spectra from Fig. 1 in grey. To aid comparison between the PV models and the observed Si IIλ6355 PV components, the HV components have been removed from the observed spectra by adding the best-fitting HV component Gaussian fit to the spectrum.
In the text

	Fig. A.3. Corner plots of the parameter posteriors inferred by the density enhancement MCMC runs for the six SNe. The grid boundaries are indicated as grey shaded regions in panels where they are visible. The solid and dotted lines in the histogram panels correspond to the means and standard deviations, respectively.
In the text