© The Authors 2025

1. Introduction

Broad-lined type Ic (type Ic-BL) supernovae (SNe) are an observationally rare class of SNe. Photometrically, these events are among the brightest SNe, reaching peak absolute magnitudes of around −18.5 mag after 10–15 days (Taddia et al. 2019). Spectroscopically, they are distinguished by the absence of hydrogen and helium absorption features and the presence of broad absorption features of Fe, Si, and Ca in their photospheric spectra (for reviews of SN classifications, see Gal-Yam et al. 2017; Filippenko 1997). These broad features are the result of the rapid expansion of the outer layers of the SN ejecta, with typical ejecta velocities being in the ∼15 000–20 000 km/s range at peak light (Srinivasaragavan et al. 2024; Taddia et al. 2019; Modjaz et al. 2016). Type Ic-BL SNe are believed to result from the death of a Wolf-Rayet star that has been stripped of its hydrogen and helium layers, either through episodic or eruptive mass loss, mass loss due to winds, or interaction in a binary system (e.g. Smith 2014; Crowther 2007; Mokiem et al. 2007).

Gamma-ray bursts (GRBs) are brief, intense pulses of gamma-rays and have been the subject of extensive theoretical and observational interest from the astronomy community for over five decades. The duration of the gamma-ray prompt emission has led to the separation of GRBs into long and short bursts (Kouveliotou et al. 1993). Long GRBs (those with prompt emission longer than 2 seconds) have been linked to the death of massive stars by the ‘collapsar’ model (Woosley 1993; MacFadyen & Woosley 1999). This model proposes that long-GRBs are the result of a jet launched by a compact object created when a massive star undergoes core-collapse. The deceleration of this jet in the ambient medium outside the former star is responsible for the GRB afterglow (Mészáros & Rees 1997), whilst dissipation of energy within the jet powers the prompt emission (Daigne & Mochkovitch 1998; Rees & Meszaros 1994).

The most notable prediction from the collapsar model is that a supernova should accompany nearly every long GRB, with the star’s outer layers being accelerated by deposition of energy from the central engine within the stellar material. This type of SN is known as an engine-driven SN. In 1998, the detection of GRB980425 and its associated type Ic-BL SN, SN1998bw, put the collapsar model at the forefront of long GRB science (Galama et al. 1998). In the years since the discovery of GRB980425-SN1998bw, more than 60 ‘GRB-SN associations’ have been detected (e.g. Finneran et al. (2025), Cano et al. (2017)). Nearly all SNe associated with GRBs are type Ic-BL SNe¹.

The association between GRBs and type Ic-BL SNe has been extensively investigated over the past 25 years (for a review of GRBs associated with SNe, hereafter GRB-SNe, see Cano et al. (2017); a review of GRB jets in SNe may be found in Corsi & Lazzati (2021)). Observations of long GRBs at low redshifts indicate that the majority of these events are associated with a supernova, with simulations showing that the same central engine can power both the GRB and the type Ic-BL SN (Barnes et al. 2018). In cases where an extensive follow-up campaign was conducted, SN-like emission has been ruled out for just five events at low redshift: GRB 060505 (Fynbo et al. 2006), GRB 060614 (Fynbo et al. 2006; Gal-Yam et al. 2006; Della Valle et al. 2006), GRB 111005A (MichałowskI et al. 2018; Tanga et al. 2018), GRB 211211A (Rastinejad et al. 2022; Troja et al. 2022), and GRB 230307A (Levan et al. 2024; Yang et al. 2024). Evidence of kilonova emission has been found for GRB230307A (Rastinejad et al. 2022; Troja et al. 2022) and GRB211211A (Levan et al. 2024; Yang et al. 2024), and this has accelerated the shift away from the typical long-short classification of GRBs.

In contrast with the abundance of nearby long-GRBs with an associated type Ic-BL SN, fewer than one in four type Ic-BL SNe have an associated GRB detection². We term these SNe ‘ordinary type Ic-BL SNe’ in this work. Despite this fact, type Ic-BL SNe appear to have very similar parameters regardless of whether or not they are associated with a GRB. Type Ic-BL SNe with GRBs synthesise ∼0.4 ± 0.2 M_⊙ of nickel, with ejecta of ∼6 ± 4 M_⊙ (e.g Cano et al. 2017), while ordinary type Ic-BL SNe synthesise ∼0.3 ± 0.2 M_⊙ of nickel, with ejecta of ∼4 ± 3 M_⊙ (e.g. Taddia et al. 2019; Prentice et al. 2016; Lyman et al. 2016; Taddia et al. 2015; Drout et al. 2011). These two populations also have similar SN kinetic energies (∼10⁵² ergs) (c.f Taddia et al. 2019; Cano et al. 2017).

These kinetic energies are extremely high relative to other SNe (Corsi & Lazzati 2021). This has been seen as evidence of an additional energy source within these events, as these energies cannot be achieved through gravitational collapse alone (Corsi & Lazzati 2021; Barnes et al. 2018). It has been inferred that the source of this additional energy may be a rapidly rotating accreting black hole or magnetar (an engine) that launches a GRB jet within the star (Corsi & Lazzati 2021). High-velocity material generated by the propagation of this jet through the stellar material may be responsible for the high kinetic energies of these SNe. For this reason they are sometimes referred to as engine-driven SNe. However, because GRB jet-like emission has not been observed for almost three quarters of type Ic-BL SNe, it is necessary to consider mechanisms by which a jet may be created by the progenitor of a type Ic-BL SN but go undetected.

The ultra-relativistic nature of GRB jets gives rise to the phenomenon of ‘orphan’ afterglows. These are GRB afterglows observed without a corresponding prompt gamma-ray component (Rhoads 1997; Mészáros et al. 1998). Observers whose line of sight lies outside the beaming cone of the GRB radiation will not see the prompt emission and will only observe the afterglow radiation once the jet begins to spread laterally as it decelerates. In contrast, the radiation emitted by a type Ic-BL SN is visible in all directions. It is therefore possible that many ordinary type Ic-BL SNe could have a contemporaneous GRB that is invisible because the observer is too far off-axis. Late-time radio observations suggest that there may be hidden GRBs in some type Ic-BL SNe (Corsi et al. 2023, 2016; Soderberg et al. 2006). At X-ray wavelengths, simulations performed using the Python package afterglowpy (Ryan et al. 2020) indicate that typical GRBs associated with type Ic-BL SNe may be invisible at observing angles above 25–30 degrees. To date, SN2020bvc is the only candidate off-axis afterglow associated with a type Ic-BL SN (Izzo et al. 2020).

Another justification for the lack of GRB emission in three quarters of type Ic-BL SN events may be the phenomenon of jet choking. As the GRB jet propagates through the outer layers of the progenitor, it may become stalled, a scenario known as a choked jet. If this happens, little to no GRB prompt emission will be produced (Nakar & Piran 2017). This process may transfer energy to the progenitor, producing a high-velocity hot cocoon around the jet (Nakar & Piran 2017; Ramirez-Ruiz et al. 2002). Such a cocoon has been inferred from spectroscopic observations of some GRB-SNe (Izzo et al. 2019). Such energy transfer may be responsible for the high-velocity features of type Ic-BL SNe.

To look for evidence of energy injection from a GRB jet in a type Ic-BL SN, Modjaz et al. (2016) measured the Fe II expansion velocities of a small sample of type Ic-BL SNe with and without an associated GRB. Their results showed that the mean expansion velocities of GRB-SNe are ∼6000 km/s higher than those of ordinary type Ic-BL SNe. However, we note that the average velocities of both populations found by Modjaz et al. (2016) are consistent within one sigma, due to large uncertainties on the average velocities of GRB-SNe. Subsequent theoretical work by Barnes et al. (2018) suggested that the deposition of jet energy within the SN ejecta does not significantly impact the SN spectra or feature velocities. For this reason, it is not clear why there should be velocity differences between ordinary type Ic-BL SNe and GRB-SNe. It is therefore possible that all type Ic-BL SNe harbour a GRB jet.

Previous studies of type Ic-BL SN velocities have relied on relatively small sample sizes (e.g. 11 GRB-SNe and 10 ordinary type Ic-BL SNe in Modjaz et al. 2016). With the rapid increase in the detection rate of type Ic-BL SNe enabled by all-sky survey telescopes such as the Zwicky Transient Facility (ZTF) (Dekany et al. 2020; Bellm et al. 2019), it has become possible to study a much larger dataset than ever before (e.g. 26 type Ic-BL SNe in Srinivasaragavan et al. (2024); 34 type Ic-BL SNe in Taddia et al. (2019)). The abundance of spectral data in online repositories allows for a quantitative investigation of velocity evolution in these SNe.

In this paper we investigate the spectral velocity evolution of the largest sample of ordinary type Ic-BL SNe and GRB-SNe collected to date, consisting of almost 900 spectra. A spline fitting method was adopted to measure the velocities of the Fe II and Si II absorption features within the type Ic-BL SN spectrum. The focus of the investigation is to determine whether the velocities of these events can help us distinguish between type Ic-BL SNe with and without GRBs. Section 2.1 describes the collection of sample data for this research, highlighting some of the challenges involved in the collation of a large heterogeneous dataset. Section 3 details the methodology and its application to the dataset. The results of this analysis are described in Sect. 4, with the discussion and conclusions presented in Sects. 5 and 6, respectively.

2. Sample selection and data collection

2.1. Sample selection

The SNe studied in this analysis were identified via an extensive search of the literature and online repositories. In particular, we made extensive use of the Weizmann Interactive Supernova Data Repository³ (WISeREP) (Yaron & Gal-Yam 2012) and the GRBSN webtool⁴ (Finneran et al. 2025). Our sample selection process is illustrated in Fig. 1.

Fig. 1. Flow of sample selection showing the data sources in orange, number of SNe at each stage of the process in blue, and filtering criteria in green. The number of each SN type at each stage of the process is shown. The search on WISeREP produced 8 GRB-SNe, 44 type Ic-BL SNe, and 130 type Ic SNe. Since the classifications reported on WISeREP may be unreliable, we performed a literature search for these SNe to confirm their types. For those SNe that could not be found in the literature, we classified them using SNID (Blondin & Tonry 2007). This step produced 52 confirmed type Ic-BL SNe and 8 GRB-SNe. Whilst searching the literature to confirm SN types, we identified an additional 1 GRB-SN and 61 type Ic-BL SNe. A total of 29 GRB-SNe were gathered from the GRBSN webtool, some of which had already been found in the literature and WISeREP searches. Combining data for events found in multiple data sources and imposing a minimum of three spectra from GRBSN or WISeREP resulted in a final sample of 61 type Ic-BL SNe and 13 GRB-SNe. We also obtained spectra for GRB180427A-SN2018fip from the ESO archive, bringing the number of spectra above the threshold of three and allowing for its inclusion in the final dataset.

We began by searching WISeREP for SNe of types Ic and Ic-BL with at least three spectral epochs and measured redshift. This allowed us to include SNe that may have been labelled as type Ic SNe prior to the common usage of the Ic-BL SN type in the early 2010s. This approach also ensures that SNe that have an ambiguous type are included in the initial sample. For example, Gangopadhyay et al. (2020) mention that SN2016P has spectral features whose widths are in between type Ic and type Ic-BL SNe, while Prentice et al. (2019) say that its velocities are more akin to type Ic SNe. We also cross-checked the list of type Ic-BL SNe provided by WISeREP with those on the Transient Name Server⁵ (TNS). This search is believed to be complete up to early April 2024.

The initial WISeREP search produced 8 GRB-SNe, 130 type Ic SNe and 44 type Ic-BL SNe. However, transient classifications reported on WISeREP may not be accurate, as they tend to be based on very early observations, often from only a single spectrum. Later observations may alter the classification of the SN, but this information may not be updated on WISeREP. This could lead to contamination of our sample with events that do not meet the criteria of a type Ic/Ic-BL SN. To confirm the classification provided by WISeREP, we followed two approaches. First, a literature search was performed for each SN. This method was used to confirm the SN type of 8 GRB-SNe and 43 type Ic-BL SNe.

For those SNe where no reliable classification could be found in the literature, we classified the SNe using the Supernova Identification (SNID) tool (Blondin & Tonry 2007). To improve the classification accuracy of SNID for type Ic-BL SNe, we included the stripped-envelope supernova (SESN) templates provided by the NYU Supernova Group⁶ (Modjaz et al. 2016; Liu et al. 2016).

During classification, we set the SN redshift to the value published on WISeREP. The age range of the SNID template spectra used for classification was restricted to t_peak − 10 and t_peak+45 days. Epochs larger than 45 days can show emission features that are relatively similar among different SN types and can lead to uncertain classification. The −10 day restriction limits the chance that a continuum dominated template spectrum is used during the fitting. This analysis identified 9 type Ic-BL SNe that (as of April 2024) have not yet been published.

The literature search and SNID classification of our WISeREP list resulted in a confirmed sample of 8 GRB-SNe and 52 type Ic-BL SNe. A separate literature search for type Ic-BL SNe identified an additional GRB-SN and 61 type Ic-BL SNe. These had not been identified in the initial WISeREP search because they were not listed as type Ic SNe or type Ic-BL SNe. A search of the GRBSN webtool turned up 29 GRB-SNe, 9 of which had already been found in the literature or on WISeREP. Combining events found in multiple data sources, at this point the sample consisted of 113 ordinary type Ic-BL SNe and 29 GRB-SNe.

This sample was further reduced by imposing a minimum of three spectral observations to allow for investigation of velocity evolution, and the need for a pre-determined redshift to take into account any possible cosmological effects. The final sample consists of 61 type Ic-BL SNe without observed GRBs and 13 GRB-SNe. A list of these SNe may be found in Table A.1.

2.2. Sample bias

The ordinary type Ic-BL SNe presented in this study were generally discovered by magnitude limited surveys such as the ZTF Bright Transient Survey (Fremling et al. 2020), which conducts a magnitude limited survey down to m_peak≤ 18.5 mag. This results in the sample containing a larger fraction of bright SNe than exists in reality within the volume sampled by the magnitude limit. This magnitude cut is made more severe by the requirement for these SNe to have spectroscopic confirmation (which is more likely for brighter SNe) and for there to be several spectra with sufficient signal-to-noise ratio (S/N) to measure a velocity (which requires a brighter SN, a longer observation, or a larger telescope).

There are some benefits to the abundance of bright SNe in the sample, in that it is often easier to fit the early light curve of bright SNe in order to determine their t₀, since there are likely more observations from surveys earlier in their evolution. This can help mitigate the issues around uncertainty on t₀, including its influence on parameter estimation.

Another aspect of including an excess of bright SNe is that the brightness of a type Ic-BL SN is directly related to the amount of nickel produced in the explosion (Arnett 1982). In GRB-associated SNe, simulations have shown that much of this Nickel production takes place on the jet-axis (Shankar et al. 2021), possibly making the SN brighter along this axis. Consequently, an over-abundance of bright, survey detected SNe might bias the sample towards type Ic-BL SNe with a polar viewing angle. This has implications for the interpretation of an absence of a detectable GRB jet in many of these explosions.

2.3. Redshift distribution

Figure 2 shows the redshift distributions of ordinary type Ic-BL SNe and GRB-SNe. Ordinary type Ic-BL SNe are detected up to z ∼ 0.2, while GRB-SNe range in redshift up to z ∼ 0.6. This is a result of differences in how these events are discovered. Ordinary type Ic-BL SNe tend to be discovered in magnitude limited surveys, which imposes a redshift limit on the SNe that can be detected. In contrast, GRB-SNe are normally discovered during optical follow-up of a GRB trigger. Knowing where and when to look for the SN allows us to detect type Ic-BL SNe associated with GRBs at greater distances than those SNe detected by surveys.

Fig. 2. Redshift distribution of the SNe used in this analysis. The most distant SNe are GRB-SNe as a result of targeted follow-up campaigns after the discovery of their associated GRB. In contrast, the bulk of the ordinary type Ic-BL SN population are discovered within redshift 0.1.

With the differences in redshift distributions of GRB-SNe and ordinary type Ic-BL SNe comes the potential for differences in the environmental conditions of these events, in particular metallicity. However, several studies indicate that the environments of long-GRBs and type Ic-BL SNe are very similar (Sanders et al. 2012b; Japelj et al. 2018; Modjaz et al. 2008). As a result, the difference in redshift distributions is not expected to significantly effect any conclusions drawn here. The impact of redshift on the measured spectral velocities has been accounted for by de-redshifting each spectrum before computing thevelocity.

2.4. Data download

The SN spectra and metadata were downloaded from WISeREP and the GRBSN webtool. We also downloaded spectra for GRB180427A-SN2018fip from the ESO Archive⁷. Spectra from WISeREP were downloaded using the WISeREP API⁸ (Müller-Bravo 2023). This is a third party tool that was modified to provide both the spectra and their metadata from WISeREP and was also used to partially automate the SNID classifications discussed previously.

Numerous instances of duplicate spectra were found in the data downloaded from WISeREP. This was common for telescopes that maintain their own database. WISeREP automatically uploads new spectra from these repositories. In some cases the original observers re-uploaded these spectra following publication of the full dataset, resulting in duplicate spectra. There was also significant overlap between GRBSN Webtool data and WISeREP data, because both tools contain data sourced from the SN literature.

To deal with duplicate and similar spectra, a visual examination was performed for all pairs of spectra with the same observation date. This step was also used to remove any near-identical spectral pairs. For example, when galaxy-subtracted spectra were available, these were retained instead of the original spectrum. If emission line removal or de-reddening of the spectrum was detected for one of the spectra in a pair, the spectrum with the highest S/N and least evidence of artefacts was retained. We found 232 duplicate spectra via this search. These spectra were excluded from velocity measurements. The final sample studied here includes 875 unique spectra, which is the largest to-date for type Ic-BL SNe.

In the case of GRB-SNe, the GRB trigger time was sourced from the GRBSN webtool and used as the explosion time (t₀) (the GRB trigger time is typically accurate to within a few seconds). For ordinary type Ic-BL SNe, the explosion time needs to be calculated from the rising of their light curve using a low-order polynomial fit. Thus, the precision of t₀ in these cases depends on how well sampled the light curve is. When available, explosion times were sourced from the literature. Otherwise, t₀ was estimated from the light curve using publicly available photometric data. If the resulting fit produced uncertainties larger than 10 days, or in the absence of accessible photometric data, the explosion time was assumed to be the halfway point between the last non-detection and the first detection, or 48 hours prior to the first detection if no constraining non-detection was available.

3. Methods

3.1. Visual examination

A visual examination was performed for all spectra. The purpose of this review was to remove corrupted data, as well as those that show nebular features or are dominated by continuum emission. Supernovae of type Ic-BL enter the nebular phase on a timescale of a few months. During this transition, the optical depth of the SN decreases due to its continued expansion, and the recombination of elements in the ejecta, rendering it transparent to optical radiation. As a result, the SN spectrum transitions from an absorption dominated regime to an emission dominated one (see e.g. Sahu et al. 2018). In a nebular spectrum, the minimum wavelength of absorption features may be shifted due to contamination by flux from adjacent emission lines, or the feature may disappear entirely. To avoid any contamination of the results for velocity with spectra from the nebular phase, each spectrum was examined for nebular features prior to fitting. This avoids any bias that may be introduced by filtering the spectra based on an arbitrary time when the nebular transition is expected to happen. Likewise, spectra that appear to be dominated by continuum are removed, since they have no clearly identifiable features to fit. This is often the case during the early evolution of a GRB-SN, when the contribution of the GRB afterglow is significant.

3.2. Redshift correction and de-reddening

The SNe analysed in this work span a wide range of redshifts. The expansion of the universe causes spectral features formed in the SN frame to appear at lower velocities in the observer frame. To account for this effect, we divide the observed wavelengths of our spectra by 1 + z, where z is the redshift of the source.

To perform this correction we use the reported redshift on WISeREP/GRBSN webtool. We initially assumed that the downloaded spectra had not been corrected for redshift; however, we then noticed that some spectra from these sources had already undergone redshift correction. If a spectrum that has already been redshift-corrected is corrected a second time, our measured velocities would appear larger than their true values. To try to avoid this issue we performed a visual examination of each spectrum. When we found evidence of prior redshift-correction, we flagged this spectrum to ensure that the correction was not applied a second time.

One marker we used to determine whether a spectrum had already been redshift-corrected was the 6563 Å Hα host-galaxy line. This line was chosen because it can be easily identified, and several of our SNe exhibit this line (∼40%, though we note that it may have been removed from some spectra prior to their publication). If this line was centred at its rest-wavelength in the downloaded spectrum, then the spectrum has already been corrected for redshift. Since this check was performed visually, it is only capable of catching larger discrepancies.

In cases where there was no Hα line from the SN host, we compared our redshift-corrected spectrum with the redshift-corrected spectrum in the published literature. When easily identified features were similarly positioned in both the published spectrum and our redshift-corrected spectrum, we considered our spectrum to have been properly redshift-corrected.

In some cases we were unable to verify whether a redshift correction had already been applied to a spectrum. This often happened due to a lack of obvious emission features or because there were no publications available as a cross-check. In these instances we applied the redshift correction ourselves before performing a visual examination of the spectral sequence. With this examination we ensured that similar spectral features remain aligned across the dataset throughout the spectral evolution. This allowed us to unearth some instances in which spectral features were shifted with respect to the same features in the rest of the spectral sequence. This generally indicates that this spectrum had already been redshift-corrected.

As the SN light propagates through its host galaxy and the Milky Way, it is scattered and absorbed by dust grains in a process known as reddening. A correction for this reddening effect may be applied to the spectra of an SN, for example by following the prescription of Cardelli et al. (1989). However, de-reddening typically scales the continuum of the spectrum, rather than affecting the minimum flux of absorption features, which is a function of the feature’s velocity. Liu et al. (2016) found that there was no difference between the de-reddened velocity and the velocity before reddening in their SN sample. In this analysis, no correction for reddening was performed, though it should be noted that some SNe may have been de-reddened prior to being uploaded to WISeREP.

3.3. Emission line removal

Emission lines from the host galaxy or the sky, as well as cosmic-ray artefacts, are present in some of our spectra. These may make it difficult to isolate the continuum, and can alter the location and shape of the absorption lines. Additionally, spectral smoothing is also impacted by these phenomena. Since the accurate determination of the wavelength at which the spectral features reach their minimum flux is critical for this analysis, any emission lines in the spectra needed to be removed prior to the fitting of the features. We note that we have not explicitly removed telluric features from our spectra. These absorption features are likely not an issue in the Fe II region of the spectrum but may contribute to the Si II feature. In some cases these lines may have been removed during the initial data reduction.

Typically, emission lines are removed during the initial processing of the observational data. However, for most of the SNe in this sample, the original files are not available, and it would be impractical to perform a full data reduction analysis for such a large dataset. For this reason, an interactive, interpolation-based emission line removal code, called emlineclipper⁹, was developed. This method can be applied directly to the reduced one-dimensional spectra.

The emlineclipper program first displays the input spectrum to the user, allowing them to choose boundaries that bracket the emission features that they wish to remove; this range is known as the ‘bounding region’ and is demarcated by the ‘bounding lines’. The code then fits a spline to the data in two 100 Å windows on each side of the bounding region. The residuals between the spline fit and the input spectrum are computed within these windows, followed by computation of the mean and standard deviation of the residual array. These values are then used to resample the spectrum by adding noise to the spline within the bounding region. Results from emission line removal using this code are shown in Fig. 3; the residuals are comparable in magnitude to the noise in the spectrum, and the resampled spectrum follows the evolution of the continuum near the emission lines very closely. The program can also handle cases of multiple nearby emission lines.

Fig. 3. Emission line removal using the emlineclipper Python code. Line selection was performed manually by selecting the bounding region around the emission lines. The data in a small range on each side of the emission line region were fit using a spline. The mean and standard deviation of the spline fit residuals in this region were used to resample the flux in the emission line region.

3.4. Smoothing and mitigation of error sources

Online repositories host spectra from a wide range of instruments spanning many decades of observations. Thus, the dataset contains spectra of varying resolutions and S/Ns. In order to increase the S/N, prior to fitting the spectra were smoothed using a Savitzky-Golay (SG) filter (Savitzky & Golay 1964). This increases the likelihood that the minimum wavelength is recovered for an absorption feature, especially if the input data is noisy. An SG filter has two tunable parameters: the degree of the polynomial and the window length. A cubic polynomial was used to reduce the risk of over-fitting the data, especially near the peaks. As there is no formal method to determine the ideal window length, a statistical analysis was conducted to determine the optimum window length (Finneran & Martin-Carrillo 2024). This analysis suggests that the optimum window length is between 2–5% of the spectrum length in resolution elements. This limit proved sufficient to avoid loss of information near the peaks when smoothing spectra in the majority of cases. The fine-tuning of smoothing parameters in this method typically introduces an additional uncertainty of ∼500–1000 km/s (Finneran & Martin-Carrillo 2024).

Despite efforts to reduce the distortion of the spectrum during smoothing, the results are somewhat sensitive to the exact choice of smoothing parameters. Additionally, the fitting of the absorption region may be sensitive to the choice of the endpoints of this region. These effects can be minimised by performing the velocity measurement for a large sample of Monte Carlo spectra and averaging the result.

In general, the spectral uncertainty arrays, which are required for the Monte Carlo analysis, are not provided by online repositories. These arrays are produced during the sky-subtraction phase of spectral analysis, and cannot be reconstructed without the original data. A method to produce pseudo-uncertainty arrays is described in Appendix B. One benefit of the Monte Carlo approach is the ability to simultaneously quantify the uncertainty introduced by the S/N of the data, the choice of the fitting region, and the uncertainty introduced by the fitted function.

3.5. Feature identification

Measuring the velocity of a spectral feature requires a high degree of confidence in the identification of the species responsible for the feature and hence the rest-frame wavelength used during the analysis. Rather than carrying out detailed modelling of the spectra for each SN, line identifications were made based on established patterns within spectra, in combination with comparison to the identifications in the literature where available.

Feature identification is important for each of the three features that dominate the type Ic-BL spectrum. The Fe II feature is well known to be a blend of at least three spectral lines, with rest-frame wavelengths of 4924 Å, 5018 Å and 5169 Å, (e.g. Modjaz et al. 2016). Depending on the chosen rest-frame wavelength, the Fe II velocity may differ by up to 15 000 km/s once the feature becomes de-blended later on in the evolution (c.f. Prentice et al. 2018). As it is not clear which line is responsible for the minimum in any given spectrum, the 5169 Å line was selected in order to maintain consistency with previous studies. Although Modjaz et al. (2016) developed a method intended to counteract the effects of Fe II blending, they found it to be inapplicable to the Si II feature. Similar methods to those used here have been applied to GRB-SNe and ordinary type Ic-BL SNe before, with the benefit that they can be applied to all of the spectral features that we wish to investigate. In cases where the Fe II feature becomes de-blended, the line that most smoothly continues the decaying velocity trend was followed; this choice was informed by visual examination of the spectral sequence. This avoids having a discontinuity in the velocity evolution as soon as the lines de-blend (e.g. Prentice.2018). Applying this identification consistently to all SNe in the sample shouldfacilitate comparison of the evolutionary trends between SNe without significant issue.

In the case of the Si II feature, the canonical rest-frame wavelength is 6355 Å. However, it has been proposed that blending with the Na I line may be an issue, as observed by Modjaz et al. (2016) for PTF10qts. There is also an ongoing debate as to whether the feature near 6100 Å, commonly identified with Si II, is in fact entirely due to silicon, or whether it is contaminated by other elements including Hα in absorption (Parrent et al. 2016). For type Ic-BL SNe the source of this hydrogen would likely be the material expelled by the Wolf-Rayet progenitor star. Once again consistency in identification is the key to studying the evolution.

The Ca II feature is also a triplet, consisting of lines at 8498 Å, 8542 Å, and 8662 Å (e.g. Rho et al. 2021). Following Rho et al. (2021), a mean wavelength of 8567 Å was used as the rest-wavelength for this feature. De-blending of this triplet is not commonly observed for type Ic-BL SNe. Indeed, as this triplet is in the near-infrared range of the spectrum it has rarely been studied for a large population of type Ic-BL SNe. Another complication with the Ca II feature is its proximity to the O I 7774 Å feature, with which it may be blended. Although we attempted to fit this feature, we were not able to obtain robust fits and thus omit these fits from the present analysis.

3.6. Velocity determination

The velocities of the three absorption features discussed above were determined using the method discussed in Appendix B. This method is very similar to that used by Silverman et al. (2012) for type Ia SNe. Similar methods have also been applied to SESNe, including those with associated GRBs (e.g. Liu et al. 2016; Patat et al. 2001). Selection of the fitting region was performed manually for all features as shown in Fig. 4. Once the observed wavelength has been determined, Eq. (B.1) is used to compute the expansion velocity. An example spline fit is shown in Fig. 5. The fit is performed on each of the Monte Carlo spectra and a histogram of velocities is generated. The median and 16th and 84th percentile uncertainties are quoted for all velocities presented in this analysis.

Fig. 4. Selection of the initial boundaries for the Fe II (5169 Å) and Si II (6355 Å) features in a spectrum of SN2020bvc. The ‘observed spectrum’ (grey) is the raw SN spectrum following redshift and emission line corrections. The ‘smoothed spectrum’ (red) is the SG-filtered observed spectrum. The hatched Fe and Si regions denote the user-selected regions in which to perform the spline fit for these features. The final boundaries of the fitting regions were allowed to vary by δ resolution elements with respect to the initial boundaries.

Fig. 5. Spline fit to the Si II feature of 2020bvc. The wavelength at which the spline fit reaches its minimum value is used in the Doppler formula to compute the expansion velocity of the feature. The ‘input spectrum’ (shown in grey) was generated by adding random noise to the SG-filtered observed spectrum (i.e. the reduced 1D spectrum from the source catalogue). The solid orange line represents the ‘smoothed spectrum’, created by applying an SG filter to the input spectrum. A solid red line illustrates the ‘spline fit’ applied to the smoothed spectrum between the wavelengths λblue and λred. The wavelength at which the spline fit reaches its minimum was taken as the observed wavelength of the feature and is denoted as λobs (purple line).

A potential source of error is introduced in the case of noisy spectra, where the smoothing by the SG filter may not be sufficient to produce a smooth feature with a clearly defined minimum. While more aggressive smoothing could be performed, this risks removing information from the spectrum. Instead, the smoothness of the cubic spline used in the fit may be adjusted by changing the number of knots used in the spline (hereafter ‘knot density’). An analysis of the stability and uncertainty of the measured velocity with different knot densities shows that the optimal knot density is ∼10% of the number of resolution elements within the limits of the feature (Finneran & Martin-Carrillo 2024). This value was adopted for the majority of spectra in this analysis.

Figure 6 shows a spectral sequence of SN2007d which was used to visually examine the spline fitting results. Due to the heterogeneity of the data, visual examination is an important step to confirm that the correct feature has been identified in all spectra for each SN. In cases where some misidentifications or poor fits were noticed, the fits were iteratively improved by: changing the delta parameter; selecting different wavelengths for the bounding maxima so that the fitting region was better constrained; refitting the feature using a spline with an increased knot density; or increasing the spline’s resolution.

Fig. 6. Spectral sequence of the type Ic-BL SN SN2007D. The observed spectra (grey) and smoothed spectra (multiple colours) are plotted as solid lines. Also shown are the determined minimum wavelengths of the Fe II (blue) and Si II (green) features, as well as the corresponding uncertainty on the minimum wavelength, which is influenced by the S/N of the spectrum. A trend of increasing minimum wavelength with time is visible for Fe II and Si II; this is a direct result of the deceleration of the SN ejecta. This plot was generated for each SN in the sample and was used to verify the feature selection and to classify nebular spectra. The last three spectra were classified as nebular, as they have a relatively flat continuum and a strong emission feature near 6300 Å. These spectra were obtained from WISeREP and come from Modjaz et al. (2014), Shivvers et al. (2019).

4. Results

The expansion velocities for each feature measured from the spectra can be found here. In all cases, the velocities are only reported for spectra younger than 60 days post explosion (rest-frame time), in order to avoid the nebular phase, where typical assumptions about the formation of the Fe II and Si II features, and homologous expansion, break down.

4.1. Trends in the velocity evolution of GRB-SNe and type Ic-BL SNe

4.1.1. Fe II

Figure 7 shows the Fe II velocity plotted against the rest-frame time since the SN explosion¹⁰. The velocity of this line was measured in at least one spectrum for 12 GRB-SNe and 59 typeIc-BL SNe. The absence of an Fe II velocity measurement for a particular spectrum can be caused by: cases of non-trivial line identification, either due to lack of reference line identification in the literature or multiple features in the region near Fe II; cases where the wavelength range of the SN spectrum did not include the Fe II region; cases where the Fe II feature did not have a well-defined minimum, either as a result of noise or there being no minimum value in this region of the spectrum; and cases in which the removal of an emission line altered the minimum of the feature so dramatically that clipping and smoothing were not feasible.

Fig. 7. Velocity evolution of the Fe II feature. Left: Plot of all SNe in the sample for which Fe II velocities were measured. All SNe follow a rapid decay at early times (< 20 days) before beginning to plateau. A continuum of velocities exists for ordinary type Ic-BL SNe and GRB-SNe. Right: Selection of SNe for which the velocity evolution has been fitted. The velocities have been scaled to highlight the overall trend. The scaling factor is computed from the fitted velocity at 15 days divided by the fitted velocity of GRB980425-SN1998bw at 15 days. The majority show a power law trend, with others showing broken power law trends. There are no evolutionary differences between the ordinary type Ic-BL SN and GRB-SN samples.

The left panel of Fig. 7 shows that velocities of both GRB-SNe and ordinary type Ic-BL SNe decline for 15–20 days prior to entering a plateau phase. There does not appear to be a single velocity that separates GRB-SNe from ordinary type Ic-BL SNe. There is significant spread in the velocities of both populations, which suggests that a continuum of events likely exists, with some GRB-SNe expanding more rapidly than ordinary type Ic-BL SNe. This is quantified in statistical terms in Sect. 4.4. Additionally, SNe that begin at high velocities tend to plateau at high velocities; this might indicate that the plateau is intrinsic to the expansion of a GRB-SN or ordinary type Ic-BL SN, or indeed to SNe in general, rather than being a distinguishing feature of either class.

Figure 7 also shows that some SNe (e.g. GRB980425-SN1998bw) exhibit a rise in velocity during their evolution around 30 days post explosion. It is not clear what has caused this fluctuation in velocity. For all SNe presented in our plots of velocity we verified from the spectral sequences that our fits give reasonable results for the locations of feature minima during this time period for all SNe. This type of behaviour has been observed in the analysis presented by Modjaz et al. (2016) (their Fig. 5) and Taddia et al. (2019) (their Fig. 15); however, no physical rationale for this was given by either paper. This trend is highlighted in this analysis because the uncertainties of the method are small with respect to these fluctuations.

The right panel of Fig. 7 reveals that a large fraction of the SN sample, both ordinary type Ic-BL SNe and GRB-SNe, follow power law decays. It also shows that the velocities, though scaled, still have some scatter, around 1000 km/s at t₀+10 days rest-frame. The scatter becomes larger after 15–20 days partially due to some SNe showing evolutions consistent with a broken power law decay. These decays tend to follow a steep-shallow evolution. Examples of likely broken power law decays are GRB980425-SN1998bw, SN2002ap and SN2009bb. The scatter at late times could also be due to the effects of line blending (see Prentice et al. 2018). The lack of early observations for many SNe makes it difficult to determine whether power laws emerge early in the evolution, or appear later.

4.1.2. Si II

Figure 8 shows the Si II velocity plotted against the rest-frame time since the SN explosion¹¹. The velocity of this line was measured in at least one spectrum for 11 GRB-SNe and 60 type Ic-BL SNe.

Fig. 8. Velocity evolution of the Si II feature. Same sample criteria and graphs as in Fig. 7. Left: Plot of all SNe in the sample for which Si II velocities were measured. Many ordinary type Ic-BL SNe and GRB-SNe evolve rapidly at early times (< 20 days). Following this, they then begin to plateau. A continuum of velocities exists for ordinary type Ic-BL SNe and GRB-SNe. GRB100316D-SN2010bh has a much higher velocity at all epochs studied than the rest of the sample. Right: Selection of SNe for which the velocity evolution has been fitted. The velocities have been scaled to highlight the overall trend. The scaling factor is computed from the fitted velocity at 15 days divided by the fitted velocity of GRB980425-SN1998bw at 15 days. Several ordinary type Ic-BL SNe and GRB-SNe can be fitted by power law decays. Supernovae that evolve more rapidly tend to have a higher initial velocity and vice-versa. Some type Ic-BL SNe and GRB-SNe show broken power law decays, with the breaks all occurring around 15 days.

The left panel of Fig. 8 shows a steady decline of the Si II velocity from 20 000–25 000 km/s to below 10 000 km/s in the first 30 days, though there are some SNe with velocities larger than this. Some SNe appear to show a plateau phase after ∼25 days, while others simply keep decaying to very low velocities. Similar to the velocity evolution of the Fe II feature, there does not seem to be a clear distinction between GRB-SNe and ordinary type Ic-BL SNe. This is investigated statistically in Sect. 4.4.

The right panel of Fig. 8 shows that, based on their Si II velocity, the SNe in the sample can be divided between those with broken power law evolutions and those with power law decays. Many of those with broken power law evolution follow a shallow-steep decay, unlike the broken power laws for Fe II. In the case of SNe with power law decays, their slopes appear to match the decay rates of the first segment of the broken power laws quite well. The breaks in the power laws appear to occur between 10–20 days post explosion in the rest-frame of the event.

There is some scatter in the velocity when viewed in log-log scale, this can be attributed to intrinsic differences between the SNe. The scatter is similar in magnitude to that found for Fe II over the same time periods, and appears to be similar throughout the evolution, though the effect of the log-log scale serves to visually enhance the scatter after 20 days when the velocity falls below 10 000 km/s.

Figure 8 includes four well-sampled high-velocity events: GRB100316D-SN2010bh, GRB980425-SN1998bw, GRB030329-SN2003dh and SN2014ad. The evolution of GRB980425-SN1998bw and SN2014ad seems quite similar, with both declining from a similar initial velocity and plateauing at around 20–22 days. However, the final plateau velocity of SN2014ad is larger by nearly a factor of two (∼18 000 km/s versus ∼9000 km/s). GRB100316D-SN2010bh appears to have a shallow velocity decrease before 10 days, followed by a steepening of its decay. GRB030329-SN2003dh achieves the highest Si II velocity (∼45 000 km/s), which rapidly declines to velocities that are more in line with the general population by day 10.

4.2. Parameterising the velocity evolution

Motivated by the evolutionary trends visible in the log-log velocities of Fe II and Si II, their velocity evolution was fit with either a power law or broken power law function. While this is the first quantitative analysis of velocity evolution for type Ic-BL SNe, similar analysis has been carried out for type Ia SNe by Zheng et al. (2017) and for type Ib SNe by Branch et al. (2002).

For the power law model, the function adopted was

$$\begin{array}{c}\hfill v=a{t}^b,\end{array}$$(1)

where t is the rest-frame time since explosion, a is the velocity of the SN one day after the explosion, and b is the slope of the velocity evolution; a negative b value implies a velocity that decreases over time.

The broken power law adopted is similar to those used in GRB afterglow fitting (e.g. Martin-Carrillo et al. 2014) and SN studies (e.g. Zheng et al. 2017). It can be expressed as

$$\begin{array}{c}\hfill v=A{[{\left(\frac{t}{t_b}\right)}^{-{\alpha }_{1}s}+{\left(\frac{t}{t_b}\right)}^{-{\alpha }_{2}s}]}^{-1/s},\end{array}$$(2)

where A is the SN velocity at the break time, α₁ and α₂ are the slopes of the power law segments, t_b is the rest-frame time of the break, and s is a parameter that controls the smoothness of the break.

While the power law fit was performed for all SNe with at least three velocity measurements, the broken power law fit was only performed if there were five or more velocity measurements, and only if the data showed clear deviation from the single power law decay. For each SN fitted with a broken power law model, an F-test was used to determine whether that model provided a better fit than the simple power law model, requiring a three sigma confidence level to accept the broken power law fit.

Marginalisation of the fit parameters was performed using the emcee Python package¹² (Foreman-Mackey et al. 2013). A basic check for convergence, based on the autocorrelation time of the walkers was performed for each fit, as recommended by the emcee documentation. Although the uncertainties on the computed feature velocities are asymmetric, the largest value was used during fitting as a conservative estimate of theuncertainty.

Similarly to Zheng et al. (2017), the value of the smoothness parameter, s, was fixed during the fitting process, since the fits seem to be only weakly dependent on it. In general, a value of s  =   ±  20 seemed to produce reasonable results for the majority of fits. However, in a few cases, the s parameter had to be changed to s  =   ±  5, indicating a smoother break. The sign of s dictates whether the power law evolution is steep-shallow or shallow-steep.

The quality of the fit was estimated based on the posterior distributions generated by emcee. Using this information, the data were cast into one of three samples: Gold, Silver, or Bronze. Figure 9 shows examples of SNe belonging to each of these samples, along with their best fits. The best-fit curves for each SN may be viewed in the figures presented here.

Fig. 9. Sample of SNe representative of the Gold, Silver, and Bronze samples. This figure shows their velocities as circles and a random sample of Markov chain Monte Carlo power law and broken power law fits as solid lines. Left: Ordinary type Ic-BL SN 2016coi, a Gold sample SN. Centre: Ordinary type Ic-BL SN iPTF13ebw, a Silver sample SN due to a fit that did not converge for one of the parameters. Right: Ordinary type Ic-BL SN 2017ifh, a Bronze sample SN due to the scatter on the data points and the quality of the fit.

The Gold sample is composed of SNe whose velocity evolution fits show little scatter, have converged for all parameters¹³, and that have near-Gaussian posterior probability distributions for the fit parameters. Supernovae in the Gold sample must have a t₀ constrained via modelling or non-detections to within 3 days.

The Silver sample is composed of SNe whose velocity evolution fits have a less clear trend than SNe from the Goldsample, but that have a sufficient number of data points to produce a good fit. The requirements on scatter and error size are relaxed compared with the Gold sample. An SN may also belong to the Silver sample if the velocity evolution fit did not converge for all parameters but the corner plots show good marginalisation of the posterior parameter distributions. Supernovae in the Silver sample must have a t₀ constrained via modelling or non-detections to within 5 days.

The Bronze sample comprises fits that have a less clear velocity evolution than the Silver sample or very few data points. The requirements on scatter and error size are relaxed compared with the Silver sample. Supernovae may also be placed in the Bronze sample if their t₀ constrained via modelling or non-detections has an error of more than 5 days or if the t₀ time is estimated as 48 hours before the last non-detection.

A final Exclude sample was created for SNe that seemed to have no clear evolutionary trend in their data. Most often, this is due to a very low number of scattered velocities or unrealistic values for the scaling constants (a in Eq. (1)). This sample was not included in the analysis of the evolution parameters.

4.2.1. Fe II velocity evolution parameters

A total of 44 SNe satisfied our criterion for power law fits in the case of the Fe II velocities¹⁴. Table 1 shows the number of SNe assigned to each sample for the Fe II velocity evolution fits. More than half of the events are in the Gold or Silver samples.

As we demonstrate in Appendix C, shifts in t₀ alter the marginalised fit parameters, particularly the power law index. The potential magnitude of this effect grows as the uncertainty in t₀ increases. Among our samples, the Gold sample has the lowest uncertainty in t₀ (< 3 days). Statistical tests presented here are only performed for the Gold sample. This will reduce the likelihood that the conclusions drawn from these tests are biased due to uncertainty in t₀. All figures shown in this subsection show the data from all three samples (Gold, Silver and Bronze) to provide an overview of the whole population.

Figure 10 shows the Fe II decay index distributions of GRB-SNe and ordinary type Ic-BL SNe¹⁵. A sample of velocities for type Ic SNe taken from Modjaz et al. (2016) is shown for comparison purposes. Regardless of subtype, the Fe II velocities of all of these SNe exhibit similar decay indices, suggesting that their decay rates are intrinsic to (stripped-envelope) SNe, rather than being a distinguishing feature of type Ic-BL SNe or GRB-SNe. Since jet activity is not generally invoked in descriptions of type Ic SNe, the similarity of these distributions implies that jets may not play a significant role in determining the pre-break evolution of the Fe II velocity. The median decay indices of GRB-SNe and ordinary type Ic-BL SNe from the Gold sample are presented in Table 2. Their similarity supports the theory that GRB-SNe and ordinary type Ic-BL SNe undergo similar early decline.

Fig. 10. Distributions of the power law index for the Fe II velocity of GRB-SNe and ordinary type Ic-BL SNe. For broken power laws, the decay index of the first segment is plotted. Fit results for the sample of type Ic SNe from Modjaz et al. (2016) are shown for comparison purposes. Left: Histograms of Fe II decay indices including data from the Gold, Silver, and Bronze samples. Right: Histograms of Fe II decay indices including data from the Gold sample only. Allowing for the low population statistics, the decay index appears to be very similar for type Ic SNe, ordinary type Ic-BL SNe, and GRB-SNe. This indicates that the presence or absence of a central engine or a relativistic jet has no impact on the rate of velocity evolution of a supernova.

The observation of visually similar decay-index distributions motivated the use of a Kolmogorov–Smirnov (KS) test to interrogate the null-hypothesis that these decay indices are drawn from the same underlying distribution. This test was performed using the scipy Python library (Virtanen et al. 2020), using the ks_2samp function (method is based on Hodges 1958). A two-sample KS test compares the cumulative distribution functions (CDFs) of two discrete empirically observed distributions, computing the maximum difference between the two CDFs. This statistic is then compared to the KS distribution to determine the p-value. In this test, the null-hypothesis is typically that the two distributions are drawn from the same population, with the alternative hypothesis being that the distributions are drawn from distinct populations. In this work, we reject the null-hypothesis if p < 0.01.

A KS test was used to compare the Fe II decay index distributions of GRB-SNe and ordinary type Ic-BL SNe from the Gold sample. The results of this test are presented in Table 3. The test yielded a p-value of 0.85, and so we cannot reject the null-hypothesis that the Fe II decay indices of GRB-SNe and ordinary type Ic-BL SNe are drawn from the same underlying distribution. This conclusion is in line with the significant overlap of the median decay indices presented in Table 2. It should be noted that the Gold sample has a limited number of events (for Fe II there are 5 GRB-SNe and 13 ordinary type Ic-BL SNe), and that these results may be altered if additional observations are made in future.

There are two outliers among the SNe in Fig. 10, both of which are type Ic-BL SNe. SN2002ap belongs to the Bronze sample, and shows evidence of a very extreme power law decay at early times, as shown in Fig. 11, with a decay index of $-3_{-1}^{+1}$, based on a broken power law fit. Although there is clear evidence that SN2002ap is not a pure power law decay, the fit does not put stringent constrains on the break time (assuming a smooth break, s = −5). The velocity of the first data point is almost 40 000 km/s, which as shown in Fig. 7 is not particularly unusual. However, there is a decay of almost 10 000 km/s to the next data point less than 1 day later, which results in the high decay index. It should be noted that, although the high decay index of SN2002ap is abnormal compared with the rest of the sample, it is still consistent within its 2σ error. In the literature SN 2002ap was found to have a spectral shape similar to GRB980425-SN1998bw and SN1997ef (Gal-Yam et al. 2002; Mazzali et al. 2002). Mazzali et al. (2002) also found that the spectral evolution is more rapid than SN1997ef, and that its kinetic energy was similar to those of GRB-SNe, though no GRB was detected. Also, they found no evidence of asymmetry, suggesting a weak or absent jet. Evidence of helium in the spectrum was also found by Mazzali et al. (2002), making it interesting that the velocity evolved so rapidly, given the larger mass available within the ejecta to slow down the recession of thephotosphere.

Fig. 11. Velocities of Fe II for SN2002ap. The best fit is a broken power law, showing the steep decay at early times. The decay index is poorly constrained, and the break region has been poorly fit. As a result of these issues, SN2002ap was assigned to the Bronze sample.

The second outlier is SN2009bb, with a decay index of $-2_{-1}^{+1}$ (see Fig. 12). The large uncertainty on the decay index means that it is consistent with the rest of the sample at 2σ. This SN is a member of the Gold sample, and has well marginalised parameters, with the fit appearing to capture the break behaviour quite well. Similar to SN2002ap, this SN was fit with s = −5, as s = −20 produced an unrealistically sharp break. Notably SN2009bb is an example of an engine driven type Ic-BL SN (Pignata et al. 2011), but with no evidence of an aspherical explosion.

Fig. 12. Velocities of Fe II for SN2009bb. The best fit is a broken power law, showing the steep decay at early times. The decay indices and break region are well constrained. Consequently, SN2009bb was assigned to the Gold sample.

4.2.2. Si II velocity evolution parameters

The distributions of the Si II decay indices are shown in Fig. 13, with the median values for GRB-SNe and ordinary type Ic-BL SNe from the Gold sample shown in Table 2. The Si II decay indices of Gold sample GRB-SNe and ordinary type Ic-BL SNe show considerable overlap, and their median values are consistent within uncertainty. There is a larger difference between the decay index distributions when comparing data from the Gold, Silver and Bronze samples. To test the similarity of the distributions, a KS test was performed for the Gold sample. The results of this test are presented in Table 3. With a p-value of 0.21, there is insufficient evidence to reject the null-hypothesis, that the observed Si II decay indices come from the same underlying distribution. This strengthens the case that GRB-SNe and ordinary type Ic-BL SNe undergo similar velocity evolution at early times (pre-break). However, the population statistics are still very low at this time (for Si II there are 6 GRB-SNe and 11 ordinary type Ic-BL SNe) and so this result will need to be re-tested with a larger sample size in future.

Fig. 13. Same as Fig. 10 but for Si II. The decay indices of Gold sample GRB-SNe are similar to those of Gold sample ordinary type Ic-BL SNe. This could be an indication that the decay index of Si II is not affected by central engine or jet activity.

The lack of positive Si II decay indices, associated with increasing velocities over time, aligns with assumptions that Si II is linked to the photospheric evolution. In homologous expansion, the greatest velocities are found in material at the edge of the expansion. This means that as the ejecta becomes optically thin and the photosphere recedes, the velocity of the features decrease. If an increasing velocity was detected, it would mean that the photosphere is moving to higher velocities, which would be difficult to explain given the expansion and cooling of the ejecta (it could be true during the very early evolution however, see Liu et al. 2018).

4.3. Broken power law fits

As mentioned previously, a subset of well sampled SNe show velocity evolutions that can be better constrained using a broken power law model. While theoretical studies suggest that the evolution of the photospheric velocity should follow a power law (Branch et al. 2002), observation of broken power law evolution may indicate that the spectra are sensitive to the evolution of more than just the photospheric velocity. It may also indicate that for these SNe the commonly used lines are not appropriate tracers of the photospheric velocity. The results of the broken power law fits can be found here. The best-fit curves for each SN may be viewed in the figures presented here.

We also considered the possibility that a large uncertainty in t₀ measurements for our type Ic-BL SN sample may have been responsible for our observation of broken power law evolution. This is because shifting the velocity data in time may introduce conditions under which power law evolution may be best fit (erroneously) by a broken power law. We quantified this using simulated data from a broken power law, and then artificially shifting this data in time and re-fitting it using power law and broken power law fits. The best-fit was determined using an F-test. We did this for multiple shifts of t₀ (from −6 to +10 days) and concluded that errors in the fit parameters grow as uncertainty in t₀ increases. However, we also found that no amount of uncertainty on t₀ can cause a power law fit to become the preferred fit compared to a broken power law fit and vice-versa. It is based on these conclusions that we decided to perform statistical tests on our measured parameters using only the Gold sample, as its events have well-constrained t₀. We present our full analysis in Appendix C.

4.3.1. Fe II broken power law fits

There are five SNe consistent with broken power law evolution for the Fe II feature. Figure 14 shows the Fe II decay indices pre-break and post-break for these events. A solid black line indicates the set of points for which the pre-break and post-break Fe II decay indices are equal. Points above this line indicate that the event has a steeper pre-break velocity evolution than post-break. Four out of five SNe follow a steep-shallow evolution (only GRB100316D-SN2010bh follows a shallow-steep decay), and there is no distinction between GRB-SNe and ordinary type Ic-BL SNe in this regard. This late shallow evolution indicates that these SNe experience a near-plateau phase, which could potentially suggest that the photosphere is stationary relative to the velocity stratification of the expanding ejecta. In the case of GRB100316D-SN2010bh, the shallow-steep decay is mostly based on the final observation, which shows a greatly reduced Fe II velocity in comparison to the rest of the evolution.

Fig. 14. Decay indices of Fe II pre-break and post-break. A solid black line indicates the set of points for which the pre-break and post-break Fe II decay indices are equal. Points above this line indicate that the event has a steeper pre-break velocity evolution than post-break. The majority of Fe II features studied undergo a steep-shallow break.

A comparison between the Fe II velocity at break, A, and the break time, t_b, (see Fig. 15) shows that the break velocity is strongly correlated with the break time for 4 of the 5 broken power laws; with larger break times implying higher velocities at break. The outlier is SN2016coi, which shows a much lower break velocity than expected. A linear fit (excluding SN2016coi) shows a positive correlation (see Table 4 for details on the resulting parameters). This differs from the negative correlation between these parameters seen in the corner plots for the individual fits. It seems that the correlation between the break velocity and the time of the break shown in Fig. 15 is intrinsic to the population, rather than a result of the fitted models. However, a larger dataset is needed to confirm this result. Both ordinary type Ic-BL SNe and GRB-SNe lie along this correlation, so it cannot be used as a discriminating test without more data.

Fig. 15. Break velocity (A) versus rest-frame break time (tb) for Fe II. The black line shows the line of best fit. A later break time implies a higher break velocity. We note that we have excluded SN2016coi from this fit.

4.3.2. Si II broken power law fits

Figure 16 shows that all Si II features for which a broken power law was the best fit undergo a shallow-steep decay. There appears to be a larger spread in the post-break decay index compared with the pre-break index. Notably, GRB171205A-SN2017iuk has the steepest post-break power law and the most extreme change of index at the break.

Fig. 16. Same as Fig. 14 but for Si II feature. A solid black line indicates the set of points for which the pre-break and post-break Fe II decay indices are equal. Points below this line indicate that the event has a shallower pre-break velocity evolution than post-break. All SNe have a shallow-steep decay. There is no clustering of GRB-SNe or ordinary type Ic-BL SNe that can be used to distinguish the populations.

As shown in Fig. 17, in contrast with the rest of the sample, SN2016coi and PTF12gzk show evidence of a third power law segment in their Si II velocity evolution. PTF12gzk was found to have a relativistic velocity component by Horesh et al. (2013), who suggested that it may be similar to GRB-SN events, with similar explosion parameters (Ben-Ami et al. 2012)¹⁶. In the case of SN2016coi, it has a similar radio and X-ray emission to SN2009bb and SN2012ap, two relativistic SNe; however the velocity of the shock-wave is sub-relativistic (Terreran et al. 2019). SN2016coi’s classification as a type Ic-BL SN (Prentice et al. 2018; Terreran et al. 2019) is robust, though evidence of helium in the early spectrum was also found (Prentice et al. 2018). As for PTF12gzk, this event appears to be a transitional object between type Ic-BL SNe and GRB-SNe.

Fig. 17. Broken power law fits to the Si II velocities of SN2016coi, PTF12gzk, and GRB060218-SN2006aj. All SNe evolve from a similar maximum velocity of 20 000 km/s, initially following a shallow decay before a break at around 18 days. A second break is visible at around 40 days for SN2016coi and PTF12gzk; there is no data beyond 20 days for GRB060218-SN2006aj, so it is not possible to prove or disprove the existence of a second break for this event.

Evidence for the third power law segment appears at t₀+40 days for both events. As such, the fits for these broken power laws were restricted to data prior to 40 days, and no attempt was made to fit the velocities beyond 40 days. Figure 16 shows that these SNe have a similar set of power law characteristics to GRB060218-SN2006aj. The break times of these SNe are all consistent within their uncertainties with a break at 18 days. We were unable to find spectra for GRB060218A-SN2006aj between ∼20 and 60 days on WISeREP or GRBSN, and so it was not possible to confirm the presence or absence of this third decay component in this GRB-SN. The third decay component could indicate the existence of a separate ejecta component that has not yet been studied in detail, and whose influence becomes apparent only after maximum light. The similarity with respect to GRB060218A-SN2006aj supports the conclusion that these events are similar to (some) GRB-SNe.

Figure 18 shows that there is a negative linear correlation between the break velocity and break time for the Si II feature. The slope and intercept values can be found in Table 4. The slope appears to be similar to the slope of the correlation visible in the corner plots of the posterior parameter distributions for the individual fits. Given this, it is difficult to know whether the correlation seen here is something intrinsic to the SNe considered in this analysis, or whether the result is strongly influenced by the behaviour of the fitted function. Nonetheless, all SNe with Si II velocity measurements exhibit this correlation, regardless of the type of event. Additionally, it can be seen from Fig. 16 that there are vastly different slopes on each side of the break for different events, so it is not clear why the break time should correlate so well with the velocity at the break. Finally, this correlation may also be influenced by the decision to fix the smoothness parameter during fitting. If future SNe exhibit broken power law evolution, this analysis may reveal whether the correlation is intrinsic, or a result of the modelling.

Fig. 18. Same as Fig. 15 but for Si II. The black line shows the line of best fit. A later break time implies a lower break velocity.

4.4. Comparing the t₀+15 day velocities of ordinary type Ic-BL SNe and GRB-SNe

If there is energy injection in GRB-SNe that is not present in (all) type Ic-BL SNe, then there should be observable differences in the velocities of the two populations. Based on the results described in Sect. 4.1, this section provides a more quantitative answer to this question.

For this analysis, the velocities computed from the (broken) power law fits at t₀+15 days were used, since many of the SNe in this sample have observations around this time, and thus their fits are well constrained¹⁷

In the case of those SNe with broken power law fits, the power law index was determined based on the break time of the broken power law. Thus, if an SN has a break after t₀+15 days, then the index of the first power law segment was used; conversely, if the break occurs before the chosen epoch, the index of the second power law segment was used.

In the following subsections, statistical tests on the velocity at t₀+15 days were computed using SNe from the Gold sample only. This decision is justified by the following arguments: The velocity evolution begins at t₀; if t₀ is not well known, comparisons based on time periods after t₀ are not valid. The SNe where the evolution is not clear may have poorly constrained fits at the chosen epoch, which would make it hard to compare their velocities. The SNe in the Gold sample have well sampled velocity evolution and strong t₀ constraints, increasing the level of confidence in the prediction of velocity from the fitted power laws for any given epoch.

4.4.1. Fe II

Figure 19 shows the distribution of Fe II velocities predicted from the power law fits to the GRB-SN, type Ic-BL SN and type Ic SN samples at t₀+15 days. The left panel of this figure shows the data from the Gold, Silver and Bronze samples. In this plot, the distribution of ordinary type Ic-BL SN velocities appears to peak at a lower velocity than that of GRB-SNe. Ordinary type Ic-BL SN velocities are intermediate between those of type Ic SNe and GRB-SNe. The GRB-SN with the maximum velocity is only around 5000 km/s faster than the fastest ordinary type Ic-BL SN.

Fig. 19. Histograms of Fe II velocities at t0+15 days. The histograms are shown with bins of 5000 km/s. Left: Histograms of Fe II velocities including data from the Gold, Silver, and Bronze samples. Right: Histograms of Fe II velocities including data from the Gold sample only.

The right panel of Figure 19 presents the Fe II velocity distributions at t₀+15 days for Gold sample SNe. In this case the GRB-SN and ordinary type Ic-BL SN peak velocities are quite similar, and an ordinary type Ic-BL SN is the highest velocity event. This similarity is also reflected by the median Fe II velocities of the Gold sample SNe shown in Table 5. The median velocities of GRB-SNe and ordinary type Ic-BL SNe are consistent within one sigma. It remains the case that type Ic SNe expand less rapidly than both ordinary type Ic-BL SNe and GRB-SNe.

A KS test was performed on the Gold sample to test the null-hypothesis that the GRB-SN and ordinary type Ic-BL SN velocities at t₀+15 days are drawn from the same underlying velocity distribution. The results of the KS test are shown in Table 6. The null-hypothesis (that the distributions of GRB-SN and ordinary type Ic-BL SN velocities come from the same underlying distributions) cannot be rejected based on the high p-value of 0.89. It is therefore not possible to identify a type Ic-BL SN with an associated GRB solely based on its Fe II velocity at t₀+15 days.

Figure 20 shows the relationship between the Fe II velocity and the Fe II decay index at t₀+15 days for SNe from the Gold and Silver samples¹⁸. These samples have t₀ constrained to within 5 days, and so predictions of their decay indices and velocities at t₀+15 days are likely reasonable. The velocities of GRB-SNe and ordinary type Ic-BL SNe show significant overlap, although some ordinary type Ic-BL SNe are located at lower velocities. There appears to be no correlation between decay index and velocity for ordinary type Ic-BL SNe or GRB-SNe at t₀+15 days. The three SNe presented here that are consistent with broken power law fits seem to have velocities consistent with the other, single power law SNe. These broken power law SNe have relatively shallow decay indices, indicating that they have all transitioned to the slow decay segment of their velocity evolution.

Fig. 20. Velocity versus decay index of Fe II at t0+15 days for GRB-SNe and ordinary type Ic-BL SNe from the Gold and Silver samples. Velocities have been computed from the power law (PL) and broken power law (BPL) fits. A sample of type Ic SNe from Modjaz et al. (2016) are shown for comparison purposes. Both ordinary type Ic-BL SNe and GRB-SNe have higher velocities than the type Ic SNe population, but there is no clear distinction between the GRB-SN and ordinary type Ic-BL SN velocities.

There is a population of low-velocity type Ic-BL SNe, with similar velocity to the type Ic SN population. Among these is SN2003bg, whose type is ambiguous. Some authors have proposed that it evolves from a BL-like SN to a type IIb SN, before becoming similar to a type Ibc SN in the nebular phase (Hamuy et al. 2009; Soderberg et al. 2006; Mazzali et al. 2009). Its spectral evolution is certainly very different to classical type Ic-BL SNe, with more narrow features than expected after the early epochs. This SN may be more similar to type Ic SNe, both in terms of its decay index and velocity at t₀+15 days. The other members of this group are: SN2018giu, which was classified using SNID as a type Ic-BL SN based on a single spectrum, so the classification may not be definitive; SN2017dcc, which both Taddia et al. (2019) and Prentice et al. (2019) found to be a type Ic-BL SN; and SN2009ca, which may be an SLSN according to Stritzinger et al. (2023). The fastest SNe are SN2014ad and GRB030329-SN2003dh; both of these SNe also have similar decay indices.

4.4.2. Si II

Figure 21 presents the distributions of velocity computed from the power law and broken power law fits to Si II velocities at t₀+15 days. When comparing data from all samples (left panel of Fig. 21), it appears that GRB-SNe exhibit higher Si II velocities than ordinary type Ic-BL SNe. However, given the uncertainties in t₀ among the Bronze and Silver samples, and the impact of such an uncertainty on the velocities computed from the fitted models, these velocities may not be very accurate. Thus, the right panel of Fig. 21 presents the Gold sample only, showing that the velocity distributions of GRB-SNe and type Ic-BL SNe are visually quite similar. Their similarity was investigated using a KS test, the results of which are shown in Table 6. The large p-value produced by this test (p = 0.46) means that we cannot reject the null-hypothesis that these two velocity samples are drawn from a common distribution. Additionally, there is significant overlap in the median velocities of both populations from the Gold sample, which are presented in Table 5. This result suggests that the Si II velocities of SNe with well determined t₀ cannot be used to identify a type Ic-BL SN associated with a GRB.

Fig. 21. Histograms of Si II velocities at t0+15 days. The histograms are shown with bins of 5000 km/s. Left: Histograms of Si II velocities including data from the Gold, Silver, and Bronze samples. Right: Histograms of Si II velocities including data from the Gold sample only.

Figure 22 shows the relationship between the Si II velocity and the Si II decay index at t₀+15 days for the Gold and Silver samples. It does not appear to be possible to identify GRB-associated type Ic-BL SNe based purely on velocity or slope within this plot. Once again the fastest SNe are SN2014ad and GRB030329-SN2003dh; both of these SNe also have similar decay indices. SN2003bg is once again the slowest SN in this sample, which could be due to the same reasons as listed in the case of the Fe II feature.

Fig. 22. Velocity versus decay index of Si II at t0+15 days for GRB-SNe and ordinary type Ic-BL SNe from the Gold and Silver samples. Velocities have been computed from the power law (PL) and broken power law (BPL) fits. GRB-SNe show similar velocities to some ordinary type Ic-BL SNe, but evolve more slowly.

Many of the SNe with broken power law fits show a shallow index Figure 22. Since the majority of the Si II feature fits show a shallow-steep decay, this indicates that many SNe have a break after 15 days. There is no evidence that SNe that follow broken power law evolution have higher velocities than SNe that are consistent with a single power law model.

4.5. Comparison of the Fe II and Si II decay indices

The range of the decay index distribution for ordinary typeIc-BL SNe is larger for both Fe II and Si II compared with the distributions for GRB-SNe. This could perhaps be intrinsic, but the more likely explanation is that there are a larger number of ordinary type Ic-BL SNe than GRB-SNe in this sample, and so it is more likely that events with a wider range of indices may be observed.

In Figure 23, SNe from the Gold, Silver and Bronze samples are presented, showing the decay indices of both features. The caveats around determining the decay index when t₀ is uncertain also influence this plot. However, it is still useful to include Silver and Bronze events here to understand the behaviour of the whole sample. Additionally, the impact of t₀ on decay index is likely to shift both the Fe II index and Si II index simultaneously, so this plot is still useful in determining which feature slows down more rapidly. The majority of the SNe show an Si II decay index above −1 and an Fe II index larger than −0.8. Within this cluster there is substantial overlap between GRB-SNe and ordinary type Ic-BL SNe. The velocity of the Fe II feature decays faster than the Si II feature velocity in just three cases: GRB171205A-SN2017iuk, SN2009bb and SN2002ap.

Fig. 23. Decay index of Si II versus the decay index of Fe II for SNe from the Gold, Silver, and Bronze samples. The decay indices have been taken from either the power law (PL) slope or the pre-break slope for a broken power law (BPL) fit. For the majority of SNe, the silicon velocity decreases more rapidly than the iron velocity. There appears to be no significant differences between GRB-SNe and ordinary type Ic-BL SNe in this regard. For broken power laws, the decay index of the first segment is plotted.

Izzo et al. (2019) propose that GRB171205A-SN2017iuk has a spectrum characterised by an excess of high velocity material (> 30 000 km/s in the first few days); with the high velocity material attributed to a cocoon around the GRB jet. They find that the early spectra are consistent with Si II velocities of up to 115 000 km/s one day after the burst, and declining to around 20 000 km/s at t₀+10 days in the rest frame of the GRB-SN. The spectra in our sample for this event begin at t₀+8 days, at which time the velocity is in a plateau at 10 000 km/s. The velocity then reduces in the following weeks. Wang et al. (2018) also fit the Si II velocity, and find a result similar to that presented here around the same epoch. As the spectrum is relatively featureless, the discrepancy between the velocities reported here and those of Izzo et al. (2019) is difficult to explain.

The Si II velocities of GRB171205A-SN2017iuk are best fit with a broken power law (Fig. 24). As this is a GRB-SN it has a well constrained t₀, so conclusions drawn from its decay indices are robust. In contrast with many SNe in the sample, the decay index is nearly flat for the first segment of the fit, with a value of $-0.1_{-0.3}^{+0.2}$. The presence of a cocoon was also suggested by D’Elia et al. (2018). If the cocoon is active during the first few days of the explosion, it is unlikely that the (relatively) late-time spectra presented here are influenced by its activity. Additionally, the velocity space of these observations is well below the 30 000 km/s that may be associated with cocoons.

Fig. 24. Broken power law fit to the Si II velocities of GRB171205A-SN2017iuk.

4.6. Co-evolution of Fe II and Si II

Figure 25 shows the difference between the Fe II and Si II velocity over time¹⁹. There is no clear difference between the GRB-SN and ordinary type Ic-BL SN populations at early times.However, from around 20 days, it appears that GRB-SNe show a larger difference in this metric.

Fig. 25. Difference between the Fe II and Si II velocities over time for ordinary type Ic-BL SNe and GRB-SNe from the Gold, Silver, and Bronze samples. Left: Absolute velocity differences for each SN. Right: Velocity differences scaled by the velocity difference for GRB980425-SN1998bw at t0+15 days.

The Fe II lines are thought to be a good tracer of the photospheric velocity (e.g. Branch et al. 2002). Figure 25 shows that the Fe II (photospheric) velocity is larger than the velocity in the Si II forming region for the majority of SNe, regardless of their type; with the difference normally in the range of 0–10 000 km/s. As discussed earlier, some SNe have uncertain t₀ measurements. This will not influence the y-axis distribution of this plot since we are comparing the directly measured velocities of Fe II and Si II, whose difference is independent of t₀. The only noticeable effect would be a movement of the velocity difference in time by at most a few days. Of course, this scatter may occur in either direction, and so overall this plot is still useful for gaining an overview of the evolution of the two populations.

Figure 25 suggests that the Fe II feature often maintains a high velocity relative to the Si II feature at late times, and also that this velocity difference is relatively constant throughout the evolution. This is in tension with the idea that the Fe II line traces the photospheric velocity. Since the Si II feature must be formed in a region above the photosphere, its velocity should be higher than that of Fe II, assuming that Fe II is a good tracer of photospheric velocity. The relative difference between the features could be due to misidentification of the underlying feature responsible for the absorption in the Fe II region, likely due to the effects of blending. It has been shown by Prentice et al. (2018) that large velocity differences can emerge between the velocity of a blended feature and the velocities of its component lines at late times when de-blending occurs. The choice to use 5169 Å as the rest-wavelength of the Fe II feature could explain the large velocity discrepancies seen at late time. Alternatively it may be the case that Fe II is not a suitable proxy for photospheric velocity, as suggested by Prentice et al. (2018).

As shown in Fig. 25, for many SNe there is an early period prior to 10 days when the Fe II velocity is declining more rapidly than the Si II velocity; this manifests as a decrease from a delta of around 20 000 km/s at early times to below 10 000 km/s between 10–20 days. Following this, the difference between the two velocities rises again. For some SNe this velocity difference exceeds 20 000 km/s after 30 days, while for others a range of around 5000–10 000 km/s is maintained out to day 60. The behaviour of both the GRB-SNe and the ordinary type Ic-BL SNe appears to be very similar in this regard.

Some SNe in Fig. 25 show unusual evolution. GRB130702A shows an increasing delta between Fe II and Si II out to around 35 days. This appears to be driven by a sharp decline in the Si II velocity around this time. On the opposite side of the distribution around this time is GRB100316D-SN2010bh, whose sharp decline in Fe II velocity results in a negative delta at this time.

5. Discussion

5.1. Comparisons to previous analyses

In order to verify that our method gives reliable results, we compared the Fe II and Si II velocities of a subset of our SNe to those in the literature. The SNe in this subset are: SN1997ef (data from Patat et al. 2001), GRB980425-SN1998bw (data from Modjaz et al. 2016; Patat et al. 2001) and GRB130702A (data from D’Elia et al. 2015). These SNe were chosen because they have well-studied velocity evolution. The spectral sequences of these SNe may be found here. Figure 26 shows the results of this investigation. The Si II velocity measurements of SN1997ef and GRB980425-SN1998bw appear to be reliable, as they show strong similarities to the literature values in magnitude and evolutionary trend. This is to be expected, since both the method applied in this analysis and the methods used in the literature rely on the determination of the wavelength of minimum flux within a spectral feature. We conclude that the estimates of Si II velocities for type Ic-BL SNe presented here are reliable.

In contrast with the results for Si II velocities, the picture for Fe II velocities is mixed. In recent years, a method of velocity determination proposed by Modjaz et al. (2016) has become commonplace for measuring the Fe II velocities of type Ic-BL SNe. This method blue-shifts and convolves a type Ic SN template spectrum in order to try and match a type Ic-BL SN spectrum. To determine the Fe II velocity of the type Ic-BL SN, the fitted blue-shift velocity is added to the velocity of the Fe II 5169 Å line in the template spectrum. This is intended to deal with the issue of line-blending within the Fe II feature, which can influence the results of velocity measurement methods that assume that the minimum flux of the Fe II feature is due to the Fe II 5169 Å line (as our method does).

We observe in Fig. 26 that the Fe II velocities of GRB980425-SN1998bw are approximately 12 000 km/s higher when measured from the minimum of this feature using our method than they are when measured using the method of Modjaz et al. (2016). The magnitude of this discrepancy suggests that Fe II 4924 Å or 5018 Å may be the dominant lines in this feature, and it is their velocities that are being (erroneously) inferred from the minimum flux wavelength using our method. We investigate the strengths and weaknesses of our method and the method of Modjaz et al. (2016) in much greater detail in Finneran & Martin-Carrillo (2024).

Figure 26 also shows the Fe II velocities of GRB130702A-SN2013dx determined by D’Elia et al. (2015) using the Modjaz et al. (2016) method. In this case, there is excellent agreement between our result and those of D’Elia et al. (2015). This implies that the minimum of this Fe II feature is strongly influenced by Fe II 5169 Å.

The comparisons presented in Fig. 26 suggest that the Fe II velocities presented here may not always reflect the Fe II 5169 Å velocity, but rather the average velocity of the Fe II feature. Notably, however, Fig. 26 shows that the velocity evolution obtained by both methods is similar, and so we conclude that the evolutionary trends presented here are robust.

Modjaz et al. (2016) studied the Fe II velocity evolution of 11 GRB-SNe and 10 type Ic-BL SNe. For ordinary type Ic-BL SNe, this sample represents a six-fold increase in the number of type Ic-BL SNe with measured Fe II velocity, compared to Modjaz et al. (2016). This is in part due to the explosion in type Ic-BL SN classification rates since 2018 with the advent of ZTF (Dekany et al. 2020). Differences between the GRB-SNe in the sample and those of the Modjaz et al. (2016) sample can be explained by the selection criteria: only those SNe that had at least three epochs of spectral data in online repositories are presented here.

Qualitatively, the evolutionary trend for Fe II feature velocity seen by Modjaz et al. (2016) is confirmed. However, the quantitative analysis reveals that there is no velocity that divides GRB-SNe from ordinary type Ic-BL SNe, nor are the medianvelocities of the two groups at t₀+15 days statistically different (the opposite conclusion to that of Modjaz et al. 2016).

The statistical results presented here were based on SNe in the Gold sample, in order to mitigate the effects of uncertainty in t₀ on the modelled parameters and the velocities. It should be noted that Modjaz et al. (2016) do not account for or quantify the uncertainty on t₀ in their sample, even though it will influence their comparison of velocities at t₀+15 days, by shifting the observed velocities in time. This is one possible explanation for our conclusion regarding velocity differences differing from that of Modjaz et al. (2016). This difference may also be due to the different methods used to determine velocity (template fitting vs. spline fitting) or it could be due to issues with line blending.

Overall, the range of velocities found here and in Modjaz et al. (2016) seems to be quite similar, with all SNe lying between 5000–50 000 km/s. In contrast with the results of Modjaz et al. (2016), the velocities of GRB100316D-SN2010bh are similar to GRB980425 -SN1998bw, though it is still the highest velocity GRB-SN. The remaining GRB-SNe are surrounded by ordinary type Ic-BL SNe, in a grouping that evolves from 25 000–35 000 km/s to 5000–20 000 km/s. This could be evidence of jet activity in at least some type Ic-BL SNe, or it could mean that jet activity has no impact on the magnitude or evolution of the velocities.

SN2014ad seems to be an outlier among ordinary type Ic-BL SNe; its velocities and their evolution are similar to those of GRB980425-SN1998bw and GRB100316D-SN2010bh. Sahu et al. (2018) corroborate this finding,observing that SN2014ad is characterised by higher expansion velocities than typical GRB-SNe and ordinary type Ic-BL SNe (measured from the Si II line); they also concluded that parameters derived from its light curve make it similar to GRB-SNe.

The sample of Si II features presented here is the largest of its kind for GRB-SNe and ordinary type Ic-BL SNe. As such, no direct population-level comparisons can be drawn with previous research. Statistical results for the Si II feature are in line with those found for Fe II and suggest that the velocities and decay indices of GRB-SNe are similar to those of ordinary type Ic-BL SNe. Overall, the behaviour of the Si II features is consistent with the idea that GRB-SNe and ordinary type Ic-BL SNe arise from a population of similar progenitors.

5.2. Jet activity in type Ic-BL SNe

It appears that neither the velocity nor its evolution are distinguishing characteristics between GRB-SNe and ordinary type Ic-BL SNe, at least for SNe with well-constrained t₀ values. Previous studies have ruled out GRBs (or put very stringent constraints on the phase space of the explosion) for several SNe in this sample (e.g. Soderberg et al. 2006; Corsi et al. 2016, 2023). Despite this, there is significant overlap between the ordinary type Ic-BL SNe and GRB-SNe studied here, including for some that are known not to harbour a GRB. This overlap exists for all lines studied here and therefore cannot purely be attributed to issues with line-identification or blending, since these issues do not affect every event or every feature in the spectrum.

Fig. 26. Comparison of SN velocities derived from the minima of absorption features in this work (orange circles) with SN velocities from the literature (green squares). Top left: Fe II velocities of GRB980425-SN1998bw. The literature data is from Modjaz et al. (2016), who use a custom velocity determination method in which a blue-shifted and convolved type Ic template spectrum is fit to the type Ic-BL spectrum in order to determine the Fe II velocity. This is intended to handle the issue of rest-wavelength assignment for a heavily blended feature. In this instance, there is a clear disparity between the velocities measured by both methods. This type of disparity was observed by Modjaz et al. (2016). A possible explanation for this is that Fe II 5169 Å may not be the dominant Fe line in this feature for this SN. Consequently, using the wavelength of the minimum flux of the feature with Eq. (B.1) leads to an over-estimate of the velocity. The discrepancy in this case is approximately 12 000 km/s, which suggests that Fe II 4924 Å or 5018 Å may be the dominant lines in this feature, and it is their velocities that are being (erroneously) inferred from the minimum flux wavelength. Top right: Fe II velocities of GRB130702A-SN2013dx. The literature data is from D’Elia et al. (2015), who use the method proposed by Modjaz et al. (2016). A much better agreement between methods is observed in this case, implying that the minimum of this Fe feature is strongly influenced by Fe II 5169 Å. Bottom left: Si II velocities of SN1997ef. Literature data are from Patat et al. (2001), who obtained the Si II velocity from the minimum flux wavelength of the Si II 6355 Å feature. As expected, the results of both methods show excellent agreement. Bottom right: Si II velocities of GRB980425-SN1998bw. Literature data are from Patat et al. (2001), who obtained the Si II velocity from the minimum flux wavelength of the Si II 6355 Å feature. As expected, the results of both methods show excellent agreement.

Jetted emission is inferred to power the afterglow light curves of GRB-SNe; thus, the SNe observed in association with GRBs must contain a jet. To form a GRB the jet must have sufficient energy to break out of the stellar envelope without choking (e.g. Corsi & Lazzati 2021). It has been suggested that energy injection from a central engine may be the source of the high velocities observed in GRB-SNe. However, models of the spectra of GRB-SNe suggest that the velocity increases due to a jet are small, and are present at early times (Barnes et al. 2018). The results of this study suggest that jet energy injection does not have a significant influence on the velocity of the ejecta. Thus, it is not possible to conclude that jets are present in all type Ic-BL SNe, or infer the frequency with which these jets may become choked and be seen as ordinary type Ic-BL SNe.

5.3. Factors influencing the velocity evolution

The velocities of both absorption features studied show evidence of a plateau phase beginning around t₀+20 days in the rest frame of the event. Li et al. (2023) showed that a similar plateau observed for SN2020bvc may be modelled assuming an r-process component within the ejecta, rather than the classical model of a SN powered by the decay of nickel produced by the SN shock or a GRB jet. However, a study by Anand et al. (2024) ruled out r-process enrichment in a sample of 25 type Ic-BL SNe discovered by ZTF. Given the weight of this large sample it is difficult to conclude that the plateaus are all due to r-process elements. If r-process SNe are indeed rare then plateaus are not anticipated for many SNe. A possible explanation for some of these plateaus may then be that the decay indices of these SNe are impacted by the uncertainty on t₀. However, some GRB-SNe exhibit this plateau also, so this cannot be the only explanation.

It is also possible that plateaus may be a natural result of the underlying ejecta’s density profile. Liu et al. (2018) presents the theoretical framework for the evolution of the photospheric velocity. Out of all the density profiles studied in this work, only the broken power law model shows evidence of a rapid decline in velocity followed by a slower decline phase. In this model there is an inner core where the density falls off more slowly with increasing radius, and an outer shell where the density declines more rapidly. Liu et al. (2018) points out that the actual photospheric radius evolves more slowly than their theoretical prediction; it is then possible that plateaus could be a consequence of this.

Mixing of the radioactive nickel in the ejecta can also impact the photospheric velocity evolution and the magnitude of those velocities. Moriya et al. (2020) performed simulations of SESN progenitors, finding that the greater the degree of mixing in the ejecta, the smoother and more rapid the decline in the photospheric velocity is. They found evidence of plateaus in thephotospheric velocity around 15 days prior to the bolometric luminosity peak. There is no evidence of early plateaus such as these in this sample, and so qualitatively it appears that many of the SNe studied here are highly mixed, in line with the results of studies performed on SN light curves.

5.4. Physical interpretation of power law and broken power law fits

A power law decline is anticipated for the photospheric velocity, and likely also for absorption features formed near thephotosphere (e.g. Branch et al. 2002). The majority of the SNe in this sample can be fit with a power law decline during their photospheric evolution. Differences between the decay indices for the Fe II and Si II lines may be put down to differences in the optical depths for these lines, though in most cases they appear to follow similar declining trends, with the Si II feature appearing to decline more rapidly than the Fe II feature in most cases.

Interpretation of the broken power law fits is more difficult. The analysis presented in Appendix C strengthens the case for these broken power law fits, showing that uncertainty in t₀ is not responsible for the observation of broken power law fits. Broken power law evolution has not been studied for a large sample of type Ic-BL SNe before, and there has been limited evidence for this type of evolution in other SN populations. The Si II velocities of type Ia SNe were investigated by Zheng et al. (2017), who found some limited evidence of broken power law evolution. However, it is unlikely that this was due to the presence of jets, since type Ia SNe are not typically associated with jet-generated emission. It is therefore possible that the broken power laws seen here are not due to the jet. As mentioned in the preceding section, they may be also be linked to the density profile of the ejecta; however, many GRB-SNe have been found to reside in relatively constant density, ISM-like media (e.g. Cano et al. 2017, and references therin). In addition, these fits cannot purely be due to identification or blending issues, as they are observed for both the Fe II feature (often heavily blended) and the Si II feature (usually not heavily blended) in a handful of SNe.

It is possible that these fits reveal the presence of additional ejecta components, whose velocity evolution may differ due to differences in temperature, composition, or orientation relative to the observer. Broken power laws may also be linked to aspherical explosion geometry. However, in this case it is difficult to explain why they are not observed in all GRB-SNe. One possibility is that the cocoon formed by the GRB jet as it tunnels through the star may be responsible for this break. However, the times of the breaks, and their velocities, are not consistent with the proposed cocoon properties (e.g Izzo et al. 2019). Additionally, there is only one observation of cocoon emission for any of the events in this sample: GRB171205A-SN2017iuk (Izzo et al. 2019), which does not stand out significantly from the overall population. For all Si II features, the broken power law fits follow a shallow-steep decay, whereas for Fe II the decays tended to be steep-shallow. It is difficult to explain why different velocity components would be present for different lines in the same SN. Only SN2016coi showed a similar broken power law evolution for both Fe II and Si II, though the break times and slopes are not compatible.

6. Conclusions

This paper has presented measurements of the Fe II and Si II expansion velocities of a large sample of type Ic-BL SNe. These events were drawn from online archives and are divided between a population of type Ic-BL SNe with an associated GRB detection and ordinary type Ic-BL SNe where no GRB component was found. This sample of 61 ordinary type Ic-BL SNe and 13 GRB-SNe represents the largest sample of velocity measurements for these types of events to date. The velocities have been measured using a spline fitting method.

The presence of a GRB component has no discernible impact on the expansion velocity of a type Ic-BL SN. For both the Fe II and Si II features, there is no separation between the two populations in velocity-time space. This has been confirmed statistically and by visual examination of the data for Fe II and Si II. This observation is in line with expectations from theory that show that the interaction of the jet with the stellar mantle should have little impact on the SN observables. For this reason the analysis presented here is not capable of ruling out a GRB component for type Ic-BL SNe where no such component has been observed.

The velocity evolution of a type Ic-BL SN may be fit with either a power law or broken power law model. A single power law fit is in agreement with the expectations for a single component photosphere, while those with broken power law fits may indicate the presence of a second (or even third) velocity component in the ejecta of some type Ic-BL SNe. This could be evidence of an aspherical explosion. There is no particular preference for broken power laws among type Ic-BL SNe with an associated GRB, meaning that this observation cannot be used to rule out jets in ordinary type Ic-BL SNe. The decay slopes of the Fe II and Si II features are also identical for both populations. This is likely further evidence that the jet has no direct impact on the expansion velocity of the SN.

There is some evidence that broken power law fits may be related to engine driven explosions. Both PTF12gzk and SN2016coi show strong similarities to GRB060218-SN2006aj in terms of their velocity and velocity evolution. Emission from a GRB cocoon is expected to peak within ∼1 day of the explosion (Izzo et al. 2019). Since breaks are observed at around 15 days post explosion, this raises the question of what the origin of this additional velocity component could be.

These results demonstrate that there is structure within the velocity evolution that had not been identified prior to this work. It is likely that a greater number of broken power law fits may be found in data with a higher sampling rate. If all type Ic-BL SNe harbour a GRB jet, the impacts of this may only be found at very early times (< 10 days), which are not probed in detail by this study.

Both of these issues can be addressed by changes to observing strategies, with greater emphasis placed on rapid classification and multi-wavelength follow-up of new type Ic-BL SNe. Reliable classifications are ideally needed within hours to days rather than weeks, as is currently the case for type Ic-BL SNe without a GRB detection. With upcoming facilities such as the Vera Rubin Observatory (Ivezić et al. 2019) expected to observe many more type Ic-BL SNe, rapid classification will become even more important. To this end, the authors are investigating the possibility of classifying new type Ic-BL SNe using a machine learning algorithm applied to the pre-peak light curve (Cotter et al., in prep.).

In tandem with improving classification speed and accuracy, efforts should be made to perform multi-wavelength follow-up of all type Ic-BL SNe, particularly at X-ray and radio wavelengths. This will allow the frequency of off-axis GRB jets in these explosions to be constrained and ordinary type Ic-BL SNe and GRB-SNe to be more clearly differentiated.

Data availability

Additional tables and figures are available online at https://zenodo.org/records/15825201.

Acknowledgments

The authors would like to thank the anonymous reviewer for their thorough feedback that has greatly improved this manuscript. GF and AMC acknowledge support from the UCD Ad Astra programme. LC and AMC acknowledge support from the Irish Research Council Postgraduate Scholarship No GOIPG/2022/1008. This work has made extensive use of WISeREP (available at: https://www.wiserep.org) as one source for our type Ic-BL spectra. In order to access data from WISeREP we used the WISeREP API, created by Tomás E. Müller Bravo; we thank them for providing a modified code from which we built our data-collection pipeline. GRB-SN data was gathered from the GRBSN webtool, which can be found at https://grbsn.watchertelescope.ie/. We also cross-referenced the list of type Ic and type Ic-BL SNe using the Transient Name Server (available at: https://www.wis-tns.org/). We have made use of the Rochester Astronomy List of Bright Supernovae to gather data about the detection dates of OGLE-2013-SN-134 and LSQ14bef. We also make use of the Central Bureau forAstronomical telegrams (available at: http://tamkin2.eps.harvard.edu/cbet/RecentCBETs.html) to locate some legacy data for SNe.

References

Appendix A: Supernova sample

Appendix B: Measuring the velocities of spectral features using splines

An absorption feature is characterised by its width, depth (strength) and the wavelength at which it reaches its minimum flux. Absorption features in SNe are generated above the photosphere, when photons from the photosphere are absorbed by electrons bound within the atoms of the ejecta, forming a spectral line. The lines observed in SNe do not appear at only one wavelength value, instead they show a broadened profile across a range of wavelengths. This occurs due to the bulk motion of the particles of the gas relative to the photons. This leads to several narrow lines side by side in the spectrum, producing a broadened line. If the rest-wavelength of a spectral line can be identified, then the expansion velocity of the gas can be measured using the Doppler shift formula for wavelength:

$$\begin{array}{c}\hfill v=\frac{1-\beta }{1+\beta }\times c,\end{array}$$(B.1)

where c is the speed of light and β is given by

$$\begin{array}{c}\hfill \beta =\frac{{\lambda }_{\mathit{obs}}^2}{{\lambda }_{\mathit{true}}^2}.\end{array}$$(B.2)

Here, λ_obs is the wavelength of the absorption feature measured from the spectrum, and λ_true is the wavelength of the same line emitted from a gas at rest with respect to the observer.

Silverman et al. (2012) proposed a method of determining the Doppler-shifted wavelength of an absorption feature in the spectrum of a type Ia SN. The basic tenet of this method is to fit a smooth function to the spectrum in the region of an absorption feature. The wavelength is determined from the point at which the fit reaches its minimum. This method has also been used for SESN in the past (Liu et al. 2016), and a similar idea was used to measure the Si II 6355 Å feature for GRB130702A-SN2013dx (D’Elia et al. 2015), a type Ic-BL SN with an associated GRB.

For this analysis, a method similar to those in the literature was implemented in Python. The spectrum is first smoothed using a Savitzky-Golay filter as described in Sect. 3.4. The user can adjust the smoothing window length, and a comparison of the result is displayed graphically.

In order to estimate the error on the velocity, a Monte Carlo analysis was conducted. Once the smoothed spectrum is generated, it is used to generate an uncertainty array for the analysis, in a similar manner to Liu et al. (2016). First, the residuals between the input and smoothed spectra are calculated. Then, a noise array is generated, where the noise at each is drawn from a Gaussian distribution with zero mean and standard deviation computed from the residual array. Samples from the noise array are added to the smoothed spectrum to produce 1000 mock spectra.

Prior to fitting, the user selects the region of interest in the spectrum, in the case of type Ic-BL SN spectra, the features are: Fe II at 5169 Åand Si II at 6355 ÅṪhe red and blue boundaries of the features are selected manually on the spectrum by clicking the plot. The ultimate location of the boundaries is determined from the peaks in the smoothed mock spectra that are nearest to the user-selected points. A delta value can also be specified, which defines the region (in resolution elements) in which the code searches for a local maximum relative to the user’s input. The value of δ is kept low, typically 30-50 resolution elements and is chosen to ensure the selection of the most appropriate local maximum for the feature in question (see Fig. 5). This approach facilitates setting an approximate starting boundary, allowing the wavelength to vary based on the exact maxima in the Monte Carlo spectra. In any case, the method is relatively insensitive to the choice of the bounding region, a conclusion that was also made by Silverman et al. (2012).

Once the boundaries are determined for a given spectral feature, a cubic spline is fit to each smoothed mock spectrum (see Fig. 5), producing a distribution of velocities. The median velocity and the 16th and 84th percentiles are taken as the estimate for velocity and its error respectively. The cubic spline is agnostic to the shape of the complex absorption and emission processes that may shape a spectral feature. For example, the absorption in a P-Cygni profile would not produce the same feature as pure absorption would. The wavelength of the spline fit at its minimum is converted to a Doppler velocity using Eq. B.1. Figure 5 shows a detailed example of this method applied to one of the SNe analysed in this paper.

Appendix C: Robustness of broken power law fits to uncertainties in explosion time

The broken power law fits presented in this paper were motivated by the observation of apparent breaks in the slope of the velocity evolution for both Si II and Fe II, and their adoption as the best-fit to the observations is supported by the results of an F-test when compared to a power law fit. However, it is possible that the model parameters, and potentially the best-fitting model, may change if the SN explosion epoch, t₀, is uncertain. This is not the case for GRB-SNe, as their t₀ is very well constrained by the trigger time of gamma-ray observing satellites derived from the prompt gamma-ray emission. However, for the majority of the ordinary type Ic-BL SNe in this sample, there is some uncertainty in t₀, and the level of uncertainty is not uniform across the sample. It is therefore important to understand the impact that uncertainty in t₀ may have on the marginalised parameters of the models.

Fig. C.1. Simulated broken power law and synthetic data used to test robustness of broken-power law fits. The circular points have been sampled from the broken power law, represented by the solid line. The parameters of the broken power law can be found in Table C.1. A 500 km/s Gaussian scatter has been added to the sample data points to make the data more realistic.

Fig. C.2. Violin plots showing distributions of broken power law parameters for a given shift in t0. The dotted lines show the simulated values of these parameters for a given t0 shift. All values are subject to variation due to these t0 shifts. Top: Distributions for the Si II velocity at the break, A. Second from top: Distributions for the pre-break power law index, α1. Second from bottom: Distributions for the post-break power law index, α2. Bottom: Distributions for the break time, tb.

To estimate the impact of this source of uncertainty, a simulated dataset was created using the broken power law given in equation 2. Table C.1 shows the model parameters used to generate the simulated data. These values are the average best-fit values for the velocity evolution of the Si II feature in ordinary type Ic-BL SNe (all of these SNe are in the Gold sample). Figure C.1 shows the broken power law model and a sample of 7 data points taken from this curve between t₀+8 and t₀+32 days with regular spacing. A 500 km/s Gaussian scatter was added to make the data more representative of real observations.

From this simulated dataset, multiple datasets were generated in which the data were shifted along the time axis by between -6 and +10 days (in increments of 1 day). This simulates the potential uncertainty on the t₀ value, with a shift of -1 day representing a true t₀ that occurs 1 day after the t₀ of the non-shifted data-set. These datasets were then fit using the same Markov chain Monte Carlo fitting procedure used to compute velocity evolution parameters in previous analyses. Both power law and broken power law models were fit to these datasets, and an F-test was then used to determine which model was the best fit. Based on a 2-σ confidence level, the p-values of the F-tests indicate that this range of uncertainty in t₀ does not alter the best-fit from a broken power law to a power law. This result suggests that the observation of broken power laws is robust. The parameters arising from this analysis are shown in Table C.2.

The parameters presented in Table C.2 show that although the observation of broken power laws is robust, some parameters of the best-fit broken power law are subject to significant variance under changes in t₀. Figure C.2 shows how the best-fit parameters and their distributions are altered as t₀ is shifted. This figure shows that A is consistently over-estimated and t_b consistently under-estimated compared to their values used when generating the model dataset, no matter the t₀ shift applied. This may have been caused by the introduction of Gaussian noise to the simulated dataset. Table C.2 confirms that the values for a time-shift of 0 days are consistent within their uncertainties with those of the parameters used to generate this data.

The value of A is robust to changes in t₀, and the value of t_b changes as expected when the data are shifted in time. However, the values of the power law indices are not robust to shifting the t₀. Both α₁ and α₂ are reduced (relative to the 0-day shift value) when data are shifted forward in time, i.e. the power laws become steeper. The opposite is true when the data are shifted backwards in time. This effect is more pronounced for α₂. However, there is no evidence that α₁ and α₂ converge towards each other, which would lead to a power law becoming the best-fit. This is in line with the results of the F-tests.

It can be inferred that if data generated from a power law was used in such analysis, a similar change in b (the power law index) would result. This has implications for all conclusions drawn from the power law indices used in this paper. It is for this reason that statistical tests used to draw conclusions from quantities dependent on the power law index are performed using only the Gold sample SNe, which have the smallest uncertainty on t₀. The GRB-SNe are not impacted by this effect, because they have highly constrained t₀ values.

All Tables

All Figures

Fig. 1. Flow of sample selection showing the data sources in orange, number of SNe at each stage of the process in blue, and filtering criteria in green. The number of each SN type at each stage of the process is shown. The search on WISeREP produced 8 GRB-SNe, 44 type Ic-BL SNe, and 130 type Ic SNe. Since the classifications reported on WISeREP may be unreliable, we performed a literature search for these SNe to confirm their types. For those SNe that could not be found in the literature, we classified them using SNID (Blondin & Tonry 2007). This step produced 52 confirmed type Ic-BL SNe and 8 GRB-SNe. Whilst searching the literature to confirm SN types, we identified an additional 1 GRB-SN and 61 type Ic-BL SNe. A total of 29 GRB-SNe were gathered from the GRBSN webtool, some of which had already been found in the literature and WISeREP searches. Combining data for events found in multiple data sources and imposing a minimum of three spectra from GRBSN or WISeREP resulted in a final sample of 61 type Ic-BL SNe and 13 GRB-SNe. We also obtained spectra for GRB180427A-SN2018fip from the ESO archive, bringing the number of spectra above the threshold of three and allowing for its inclusion in the final dataset.

In the text

	Fig. 2. Redshift distribution of the SNe used in this analysis. The most distant SNe are GRB-SNe as a result of targeted follow-up campaigns after the discovery of their associated GRB. In contrast, the bulk of the ordinary type Ic-BL SN population are discovered within redshift 0.1.
In the text

	Fig. 3. Emission line removal using the emlineclipper Python code. Line selection was performed manually by selecting the bounding region around the emission lines. The data in a small range on each side of the emission line region were fit using a spline. The mean and standard deviation of the spline fit residuals in this region were used to resample the flux in the emission line region.
In the text

Fig. 4. Selection of the initial boundaries for the Fe II (5169 Å) and Si II (6355 Å) features in a spectrum of SN2020bvc. The ‘observed spectrum’ (grey) is the raw SN spectrum following redshift and emission line corrections. The ‘smoothed spectrum’ (red) is the SG-filtered observed spectrum. The hatched Fe and Si regions denote the user-selected regions in which to perform the spline fit for these features. The final boundaries of the fitting regions were allowed to vary by δ resolution elements with respect to the initial boundaries.

In the text

Fig. 5. Spline fit to the Si II feature of 2020bvc. The wavelength at which the spline fit reaches its minimum value is used in the Doppler formula to compute the expansion velocity of the feature. The ‘input spectrum’ (shown in grey) was generated by adding random noise to the SG-filtered observed spectrum (i.e. the reduced 1D spectrum from the source catalogue). The solid orange line represents the ‘smoothed spectrum’, created by applying an SG filter to the input spectrum. A solid red line illustrates the ‘spline fit’ applied to the smoothed spectrum between the wavelengths λblue and λred. The wavelength at which the spline fit reaches its minimum was taken as the observed wavelength of the feature and is denoted as λobs (purple line).

In the text

Fig. 6. Spectral sequence of the type Ic-BL SN SN2007D. The observed spectra (grey) and smoothed spectra (multiple colours) are plotted as solid lines. Also shown are the determined minimum wavelengths of the Fe II (blue) and Si II (green) features, as well as the corresponding uncertainty on the minimum wavelength, which is influenced by the S/N of the spectrum. A trend of increasing minimum wavelength with time is visible for Fe II and Si II; this is a direct result of the deceleration of the SN ejecta. This plot was generated for each SN in the sample and was used to verify the feature selection and to classify nebular spectra. The last three spectra were classified as nebular, as they have a relatively flat continuum and a strong emission feature near 6300 Å. These spectra were obtained from WISeREP and come from Modjaz et al. (2014), Shivvers et al. (2019).

In the text

Fig. 7. Velocity evolution of the Fe II feature. Left: Plot of all SNe in the sample for which Fe II velocities were measured. All SNe follow a rapid decay at early times (< 20 days) before beginning to plateau. A continuum of velocities exists for ordinary type Ic-BL SNe and GRB-SNe. Right: Selection of SNe for which the velocity evolution has been fitted. The velocities have been scaled to highlight the overall trend. The scaling factor is computed from the fitted velocity at 15 days divided by the fitted velocity of GRB980425-SN1998bw at 15 days. The majority show a power law trend, with others showing broken power law trends. There are no evolutionary differences between the ordinary type Ic-BL SN and GRB-SN samples.

In the text

Fig. 8. Velocity evolution of the Si II feature. Same sample criteria and graphs as in Fig. 7. Left: Plot of all SNe in the sample for which Si II velocities were measured. Many ordinary type Ic-BL SNe and GRB-SNe evolve rapidly at early times (< 20 days). Following this, they then begin to plateau. A continuum of velocities exists for ordinary type Ic-BL SNe and GRB-SNe. GRB100316D-SN2010bh has a much higher velocity at all epochs studied than the rest of the sample. Right: Selection of SNe for which the velocity evolution has been fitted. The velocities have been scaled to highlight the overall trend. The scaling factor is computed from the fitted velocity at 15 days divided by the fitted velocity of GRB980425-SN1998bw at 15 days. Several ordinary type Ic-BL SNe and GRB-SNe can be fitted by power law decays. Supernovae that evolve more rapidly tend to have a higher initial velocity and vice-versa. Some type Ic-BL SNe and GRB-SNe show broken power law decays, with the breaks all occurring around 15 days.

In the text

Fig. 9. Sample of SNe representative of the Gold, Silver, and Bronze samples. This figure shows their velocities as circles and a random sample of Markov chain Monte Carlo power law and broken power law fits as solid lines. Left: Ordinary type Ic-BL SN 2016coi, a Gold sample SN. Centre: Ordinary type Ic-BL SN iPTF13ebw, a Silver sample SN due to a fit that did not converge for one of the parameters. Right: Ordinary type Ic-BL SN 2017ifh, a Bronze sample SN due to the scatter on the data points and the quality of the fit.

In the text

Fig. 10. Distributions of the power law index for the Fe II velocity of GRB-SNe and ordinary type Ic-BL SNe. For broken power laws, the decay index of the first segment is plotted. Fit results for the sample of type Ic SNe from Modjaz et al. (2016) are shown for comparison purposes. Left: Histograms of Fe II decay indices including data from the Gold, Silver, and Bronze samples. Right: Histograms of Fe II decay indices including data from the Gold sample only. Allowing for the low population statistics, the decay index appears to be very similar for type Ic SNe, ordinary type Ic-BL SNe, and GRB-SNe. This indicates that the presence or absence of a central engine or a relativistic jet has no impact on the rate of velocity evolution of a supernova.

In the text

	Fig. 11. Velocities of Fe II for SN2002ap. The best fit is a broken power law, showing the steep decay at early times. The decay index is poorly constrained, and the break region has been poorly fit. As a result of these issues, SN2002ap was assigned to the Bronze sample.
In the text

	Fig. 12. Velocities of Fe II for SN2009bb. The best fit is a broken power law, showing the steep decay at early times. The decay indices and break region are well constrained. Consequently, SN2009bb was assigned to the Gold sample.
In the text

	Fig. 13. Same as Fig. 10 but for Si II. The decay indices of Gold sample GRB-SNe are similar to those of Gold sample ordinary type Ic-BL SNe. This could be an indication that the decay index of Si II is not affected by central engine or jet activity.
In the text

	Fig. 14. Decay indices of Fe II pre-break and post-break. A solid black line indicates the set of points for which the pre-break and post-break Fe II decay indices are equal. Points above this line indicate that the event has a steeper pre-break velocity evolution than post-break. The majority of Fe II features studied undergo a steep-shallow break.
In the text

	Fig. 15. Break velocity (A) versus rest-frame break time (tb) for Fe II. The black line shows the line of best fit. A later break time implies a higher break velocity. We note that we have excluded SN2016coi from this fit.
In the text

Fig. 16. Same as Fig. 14 but for Si II feature. A solid black line indicates the set of points for which the pre-break and post-break Fe II decay indices are equal. Points below this line indicate that the event has a shallower pre-break velocity evolution than post-break. All SNe have a shallow-steep decay. There is no clustering of GRB-SNe or ordinary type Ic-BL SNe that can be used to distinguish the populations.

In the text

Fig. 17. Broken power law fits to the Si II velocities of SN2016coi, PTF12gzk, and GRB060218-SN2006aj. All SNe evolve from a similar maximum velocity of 20 000 km/s, initially following a shallow decay before a break at around 18 days. A second break is visible at around 40 days for SN2016coi and PTF12gzk; there is no data beyond 20 days for GRB060218-SN2006aj, so it is not possible to prove or disprove the existence of a second break for this event.

In the text

	Fig. 18. Same as Fig. 15 but for Si II. The black line shows the line of best fit. A later break time implies a lower break velocity.
In the text

	Fig. 19. Histograms of Fe II velocities at t0+15 days. The histograms are shown with bins of 5000 km/s. Left: Histograms of Fe II velocities including data from the Gold, Silver, and Bronze samples. Right: Histograms of Fe II velocities including data from the Gold sample only.
In the text

Fig. 20. Velocity versus decay index of Fe II at t0+15 days for GRB-SNe and ordinary type Ic-BL SNe from the Gold and Silver samples. Velocities have been computed from the power law (PL) and broken power law (BPL) fits. A sample of type Ic SNe from Modjaz et al. (2016) are shown for comparison purposes. Both ordinary type Ic-BL SNe and GRB-SNe have higher velocities than the type Ic SNe population, but there is no clear distinction between the GRB-SN and ordinary type Ic-BL SN velocities.

In the text

	Fig. 21. Histograms of Si II velocities at t0+15 days. The histograms are shown with bins of 5000 km/s. Left: Histograms of Si II velocities including data from the Gold, Silver, and Bronze samples. Right: Histograms of Si II velocities including data from the Gold sample only.
In the text

	Fig. 22. Velocity versus decay index of Si II at t0+15 days for GRB-SNe and ordinary type Ic-BL SNe from the Gold and Silver samples. Velocities have been computed from the power law (PL) and broken power law (BPL) fits. GRB-SNe show similar velocities to some ordinary type Ic-BL SNe, but evolve more slowly.
In the text

Fig. 23. Decay index of Si II versus the decay index of Fe II for SNe from the Gold, Silver, and Bronze samples. The decay indices have been taken from either the power law (PL) slope or the pre-break slope for a broken power law (BPL) fit. For the majority of SNe, the silicon velocity decreases more rapidly than the iron velocity. There appears to be no significant differences between GRB-SNe and ordinary type Ic-BL SNe in this regard. For broken power laws, the decay index of the first segment is plotted.

In the text

	Fig. 24. Broken power law fit to the Si II velocities of GRB171205A-SN2017iuk.
In the text

	Fig. 25. Difference between the Fe II and Si II velocities over time for ordinary type Ic-BL SNe and GRB-SNe from the Gold, Silver, and Bronze samples. Left: Absolute velocity differences for each SN. Right: Velocity differences scaled by the velocity difference for GRB980425-SN1998bw at t0+15 days.
In the text

Fig. 26. Comparison of SN velocities derived from the minima of absorption features in this work (orange circles) with SN velocities from the literature (green squares). Top left: Fe II velocities of GRB980425-SN1998bw. The literature data is from Modjaz et al. (2016), who use a custom velocity determination method in which a blue-shifted and convolved type Ic template spectrum is fit to the type Ic-BL spectrum in order to determine the Fe II velocity. This is intended to handle the issue of rest-wavelength assignment for a heavily blended feature. In this instance, there is a clear disparity between the velocities measured by both methods. This type of disparity was observed by Modjaz et al. (2016). A possible explanation for this is that Fe II 5169 Å may not be the dominant Fe line in this feature for this SN. Consequently, using the wavelength of the minimum flux of the feature with Eq. (B.1) leads to an over-estimate of the velocity. The discrepancy in this case is approximately 12 000 km/s, which suggests that Fe II 4924 Å or 5018 Å may be the dominant lines in this feature, and it is their velocities that are being (erroneously) inferred from the minimum flux wavelength. Top right: Fe II velocities of GRB130702A-SN2013dx. The literature data is from D’Elia et al. (2015), who use the method proposed by Modjaz et al. (2016). A much better agreement between methods is observed in this case, implying that the minimum of this Fe feature is strongly influenced by Fe II 5169 Å. Bottom left: Si II velocities of SN1997ef. Literature data are from Patat et al. (2001), who obtained the Si II velocity from the minimum flux wavelength of the Si II 6355 Å feature. As expected, the results of both methods show excellent agreement. Bottom right: Si II velocities of GRB980425-SN1998bw. Literature data are from Patat et al. (2001), who obtained the Si II velocity from the minimum flux wavelength of the Si II 6355 Å feature. As expected, the results of both methods show excellent agreement.

In the text

	Fig. C.1. Simulated broken power law and synthetic data used to test robustness of broken-power law fits. The circular points have been sampled from the broken power law, represented by the solid line. The parameters of the broken power law can be found in Table C.1. A 500 km/s Gaussian scatter has been added to the sample data points to make the data more realistic.
In the text

Fig. C.2. Violin plots showing distributions of broken power law parameters for a given shift in t0. The dotted lines show the simulated values of these parameters for a given t0 shift. All values are subject to variation due to these t0 shifts. Top: Distributions for the Si II velocity at the break, A. Second from top: Distributions for the pre-break power law index, α1. Second from bottom: Distributions for the post-break power law index, α2. Bottom: Distributions for the break time, tb.

In the text