Carnegie Supernova Project. II. Near-infrared Spectral Diversity and Template of Type Ia Supernovae

With the goal of constructing spectral templates as a function of the color-stretch parameter (s_BV ; Burns et al. 2014) and phase, the following selection criteria were applied in this order:

• 1.  
  Spectra of SNe Ia that belong to one of the following peculiar subgroups: 02cx-like (Li et al. 2003), 02ic-like (Hamuy et al. 2003), 03fg-like (Howell et al. 2006), and Ca-strong (Galbany et al. 2019) were removed due to their peculiar spectral features. This step excludes 39 spectra of 14 SNe.
• 2.  
  Each SN must have a reliable time of the B-band maximum (), and each spectrum must have a phase that is less than or equal to 100 days. This step further excludes 83 spectra of 30 SNe Ia.
• 3.  
  Each spectrum was visually inspected for significant host-galaxy contamination by comparing with spectra in a similar parameter space. This final step further excludes eight spectra of four SNe.

Out of the initial 461 spectra of 142 SNe, 331 spectra of 94 SNe Ia remain after the three-step sample cut. Figure 1 illustrates the change of sample size after applying each selection criterion. All spectra that passed through the first two steps, 339 spectra of 98 SNe, are available in the CSP website. ³²  A log of these spectra is listed in Appendix A.

Zoom In Zoom Out Reset image size

and the s_BV of the sample SNe were obtained by fitting multiband CSP light curves using SNooPy with max_model. An arbitrary cutoff phase was set at 100 days relative to the B-band maximum, well past the phase range relevant for light-curve fitting.

The peculiar SNe Ia were identified based on their light curves and optical spectra. Note that spectra of the slow-declining 91T-like (Filippenko et al. 1992a; Phillips et al. 1992) and fast-declining 91bg-like (Filippenko et al. 1992b; Leibundgut et al. 1993) SNe Ia were kept in the sample intentionally for the following reasons: (1) Their NIR spectral features, while extreme, were shown to extend the properties of normal SNe Ia (e.g., Hoeflich et al. 2002; Hsiao et al. 2015; Phillips et al. 2022). (2) Studies have shown that 91bg-like and fast-declining SNe Ia in general are standardizable and do not impact Hubble constant measurements significantly (e.g., Burns et al. 2014; Hoogendam et al. 2022; Phillips et al. 2022; Yang et al. 2022). (3) The inclusion of these SNe will help broaden the s_BV baseline since they are on the edges of the distribution.

Only eight spectra show significant host contamination. These are usually SNe Ia situated close to the host-galaxy core and taken at late phases when the background flux is comparable to or higher than the SN itself. The effects are shallow spectral slope and more diluted spectral features especially toward the red end. These were easy to identify through visual inspection and comparison to other spectra at similar phases. They also showed up as outliers during the dimensionality reduction step (Section 3.2) if there were included.

The final sample of spectra that passed the selection criteria (331 spectra of 94 SNe) cover a wide range of heliocentric redshift, s_BV , and phase (see Figure 2). A full range in s_BV and phase is desirable as the goal here is to build a set of spectral templates as a function of s_BV and phase. A wide range in redshift is also useful as spectral features are shifted out of regions that are heavily affected by telluric absorption. The sample SNe cover a heliocentric redshift range of 0.003 ≤ z ≤ 0.087 with a median of 0.021 and a color-stretch range of 0.48 ≤ s_BV ≤ 1.30 with a median of 0.97. The r-band absolute magnitudes of the host galaxies are obtained using the Z-PEG (Le Borgne & Rocca-Volmerange 2002) SED software following the procedure described in Uddin et al.(2017).

Zoom In Zoom Out Reset image size

The majority of the SNe Ia (70%) in this sample were discovered by untargeted searches. This is the result of the survey strategy of CSP-II, designed to avoid bias toward bright and large host galaxies (Phillips et al. 2019). As shown in the bottom panels of Figure 2, SNe from untargeted searches extend out to higher redshifts and cover a wider range of the host-galaxy absolute magnitude compared to those from targeted searches. They also have a similar distribution in s_BV except for a possible excess at the bright and slowly declining end. The normalized distributions of the full CSP-II SNe Ia sample (combined “Physics” and “Cosmology” samples from Phillips et al. 2019) are also shown in Figure 2. The selected sample for the template is representative of the full CSP-II sample in terms of s_BV .

Next, we examine the distribution of the sample spectra and SNe in phase−s_BV space (see Figure 3). Unsurprisingly, most spectra were taken near maximum light and belong to normal-bright SNe Ia. The sample spectra have a median s_BV of 0.98 and a median phase of 14 days. Figure 3 then identifies the data that would provide the most improvement for the next iteration of template building: NIR spectra of fast-declining SNe in general, particularly at late phases. These observations are naturally difficult as these SNe are rare, intrinsically faint, and fade more rapidly than normal-bright SNe.

Zoom In Zoom Out Reset image size

Our approach is to model the dependence of spectroscopic properties on phase and s_BV as a hypersurface in the phase−s_BV space (Section 3.3). The dearth of data in certain regions and the over-reliance on individual well-observed SNe (e.g., SN 2012fr has 25 spectra in the sample) may produce potential biases if the data represent the extremes in the respective bins. The best way to improve this is a larger sample that fills in the phase−s_BV space. To characterize this potential bias, we compare the templates built with and without the 25 spectra of SN 2012fr in Section 4.1.

Since there are no spectrophotometric standard stars in the NIR that can be observed by FIRE, flux calibration was done using telluric standards (Hsiao et al. 2019). Each standard is an A-type main-sequence star similar to Vega (A0 V star) observed at a similar airmass, usually immediately before or after the science observation. The flux calibration is performed simultaneously with the telluric correction given the B and V magnitudes of the A0 V star, following the procedures outlined by Vacca et al. (2003). Here, we examine the spectrophotometric accuracy of our sample NIR spectra calibrated using this method in terms of the relative colors. Note that the vast majority of the spectra were observed at parallactic angle to prevent the change of spectroscopic colors caused by the slit loss.

To assess the spectrophotometric accuracy after flux calibration, the synthetic broadband Y − J and J − H colors measured directly from the spectra are compared to those from the corresponding light curves. All spectra were corrected for Milky Way (MW) extinction using the Schlafly & Finkbeiner (2011) dust map and the Fitzpatrick (1999) law assuming R_V = 3.1. Spectroscopic colors were then measured from synthetic photometry using YJH filters. Due to the design of CSP-II, the photometric coverage has typically nightly cadence in the optical but has large gaps in the NIR (Phillips et al. 2019). Here, we took advantage of the simultaneous optical-to-NIR light-curve fitting of SNooPy, using the SNooPy models to interpolate the sparse data in the NIR. The MW extinction-corrected Y − J and J − H colors interpolated to the dates of the spectroscopic observations were then taken as the comparison photometric colors.

The comparison between the observed spectroscopic and photometric colors is shown in Figure 4. The rms of the difference between spectroscopic and photometric colors is around 0.2 mag in both Y − J and J − H, with the mean systematic difference within 0.1 mag. Furthermore, Hsiao et al. (2019) showed that the median color differences between 13 NIR spectra and well-sampled light curves are 0.08 mag in Y − J and 0.03 mag in J − H. This indicates that the spectrophotometric accuracy of this sample of NIR spectra (without the contribution of the interpolation uncertainty) is roughly at the 10%–20% level.

Zoom In Zoom Out Reset image size

Another way to examine the spectroscopic colors is to look at their evolution over time. Figure 5 shows the resulting time evolution of the spectroscopic Y − J and J − H colors. There is clear evidence that the color evolution depends heavily on the s_BV parameter, indicating that the spectroscopic colors are accurate enough to preserve such a signal. To avoid incurring further uncertainty from the interpolation of sparse data and the narrower phase coverage of SNooPy, we chose not to match the spectroscopic colors to the photometric colors in this work.

Zoom In Zoom Out Reset image size

As the entire data set is observed with a single instrument and configuration by design, it is possible to create a common rest-frame wavelength grid that is representative of the pixel grid and can minimize the interpolation uncertainties. First, an average wavelength grid is determined by taking the mean calibrated wavelength of the sample spectra at each pixel in the observer frame. This wavelength grid was then de-redshifted by the median redshift of the sample (z_median = 0.021). Finally, the sampling frequency of the grid was doubled by adding a mid-point between neighboring points, forming a rest-frame wavelength grid with 3673 wavelength points. The resulting wavelength interval ranges from ∼3 Å in the Y band to ∼9 Å in the K band.

All sample spectra were de-redshifted to the rest frame and corrected for MW extinction using the Schlafly & Finkbeiner (2011) dust map and the Fitzpatrick (1999) law assuming R_V = 3.1. They were then interpolated onto the common wavelength grid described above. Note that, to prevent spurious flux values from dominating the PCA, sigma-clipping was performed before the interpolation, removing all flux values 5σ above or below the local mean with an average window size of 100 Å.

If PCA is performed directly on the observed spectra covering so wide of a wavelength range as the NIR, the principal components (PCs) would be dominated by regions of relatively high intrinsic flux (e.g., the z and Y bands) and spurious variations (e.g., the telluric regions and at the edges of the spectra). The model SEDs need to include the more subtle variations of the spectral features. To achieve this, dimensionality reduction was performed in seven wavelength regions. The wavelength regions were somewhat arbitrarily defined but designed to correspond to the zYJHK ³³  filters with some overlap. The overlapping region ranges from ∼90 Å between W1 and W2 to ∼370 Å between W6 and W7, with roughly 50 sampling points in each. The seven wavelength regions are presented graphically in the top panel of Figure 6 and tabulated in Table 1.

Zoom In Zoom Out Reset image size

Additional selection criteria were imposed on the sample spectra in each wavelength region:

• 1.  
  Each spectrum must have a median signal-to-noise ratio (S/N) ≥ 3.
• 2.  
  In the telluric regions, W4 and W6, spectra that are over- or under-corrected for telluric absorptions were excluded by visual inspection.
• 3.  
  In region W7 (K band), spectra with spurious flux values at the edge of the detector were excluded by visual inspection.

Separating a spectrum into wavelength regions has the added benefit that if the spectrum, for example, has poor telluric corrections, we do not have to discard other wavelength regions containing useful information. Approximately two-thirds of the sample have reliable telluric corrections, consistent with what we found in our SN II sample using the Paα feature (Davis et al. 2019).

The distributions of the S/N in each region are presented in the lower panels of Figure 6, and the final sample count and median S/N in each wavelength region are tabulated in Table 1. The S/N is highest in the Y band and decreases toward the blue and the red. As expected, the lowest sample count and S/N are found in the telluric regions, W4 and W6, with the median S/N staying above 10.

An SN Ia spectrum contains a large number of elements, but they do not vary independently. Spectral variations are largely captured by only a few parameters (e.g., Nugent et al. 1995). It is thus possible to drastically reduce the dimensionality of SN Ia spectroscopic data, and PCA was chosen as the statistical technique for this purpose. It is widely used for modeling the SN spectra (e.g., Hsiao et al. 2007; Davis et al. 2019; Shahbandeh et al. 2022). PCA produces an orthogonal set of eigenvectors or PCs ranked by their eigenvalues, which give a measure of the amount of data variation they capture. In this work, the matrix diagonalization was done using the Python package scikit-learn (Pedregosa et al. 2011).

As mentioned in Section 3.1, the data were separated into wavelength regions in order to capture subtle variations. In each wavelength region, flux elements were placed on a uniform wavelength grid and normalized to have the same integrated flux in the region before PCA. Our chosen light-curve parameter is the color-stretch parameter s_BV , since Δm₁₅(B) and SALT x1 do not reliably classify fast decliners (Burns et al. 2014). The PCA operation has no knowledge of the phase and s_BV labels for each spectrum nor the wavelength value for each flux element.

The resulting PCs are presented in Figure 7. Note that, in the figure, the PCs in the overlapped area between neighboring wavelength regions are not expected to be connected, since PCA is independently performed in each region. In regions W1–W6, it takes less than five PCs to capture >90% of the variations. Regions W2 and W3 (Y and J bands), in particular, only required two PCs each, indicating that the spectra are rather uniform. On the other hand, W7 (K band) requires six PCs to reach 90% of the variation, perhaps hinting at substantial diversity.

Zoom In Zoom Out Reset image size

Figure 7 also shows the reconstructed spectra of 2σ variation on top of the mean spectrum. Familiar spectral features in the NIR then emerge, such as the Ca ii IR triplet in W1 (z band), Mg ii in W2 (Y band), the H-band break in W5 (H band), and Co ii in W7 (K band). Another example of the captured information can be found in PC1 and PC3 of W1 and W2, which appear to describe the strength and velocity of the spectral features, respectively, in both regions.

Each input spectrum can then be projected onto a PC in the multidimensional data space, providing an associated projection value. Mathematically, each input spectrum f _m of n flux elements is then represented by the sum of all PCs ξ _i weighted by the projections p_i,m :

The projections are shown in Figure 8 with respect to phase and s_BV . It is clear that the first several PCs in all regions have strong dependence on phase and s_BV . The dependence is the weakest in the telluric regions W4 and W6 because of the lower S/N of the input data and the general lack of strong features. W5 (H band) shows particularly strong dependence in the top four PCs, giving confidence that the most prominent spectral features in the NIR can be well predicted with light curves alone. This offers the motivation for the next step: modeling the dependence of spectral properties, represented by the PC projections, as a hypersurface in phase-s_BV space.

Zoom In Zoom Out Reset image size

In the previous step, the spectral data have been reduced in dimensionality such that each spectrum is represented by a few PC projections. Here, we will determine the dependence of the spectral properties on phase and s_BV . The method chosen to map the projection hypersurface in phase-s_BV space is GPR.

GPR has been shown to be successful at interpolating supernova light-curve data with intermittently missing data (e.g., Vincenzi et al. 2019; Pessi et al. 2022). The data characteristics are similar here (Figure 3). It also has the advantage of producing an uncertainty estimate to the hypersurface, which allows for the construction of not only the spectral templates but also the statistical flux uncertainties. In this work, GPR was performed using the Python package scikit-learn (Pedregosa et al. 2011). An example of the resulting hypersurface is shown in Figure 9, successfully capturing the spectral properties as a function of phase and s_BV .

Zoom In Zoom Out Reset image size

The kernel setup for the GPR in this work was a constant kernel multiplied by a radial-basis function (RBF) kernel plus a white noise kernel. ³⁴  The hyperparameters are optimized by maximizing the log-marginal likelihood. The details of the initial hyperparameter setting and optimal values can be found in Appendix C. Before performing GPR, the projection outliers were excluded by clipping spectra with projection values 5σ above or below the mean. The remaining projection values were then normalized to range from 0–1 to ensure consistency across all wavelength regions and PCs while maintaining a uniform kernel setup.

With the hypersurface in hand, one can obtain the modeled PC projection value and the associated uncertainties for each PC in each wavelength region, given the phase and s_BV . Once it was decided which PCs are to be included, template spectra were then constructed via reverse transformation to return the projections to flux space. Finally, the spectra from individual wavelength regions were merged by matching the integrated flux in the overlapping regions and using the weighted average flux. ³⁵  Now the question becomes: what is the minimum number of PCs necessary to describe the spectroscopic diversity and how to select them efficiently? This is addressed in Section 3.4.

An important feature of the new template is flux uncertainty estimation, indicating the confidence of the SED prediction. For example, the lack of data in certain parameter space translates to larger uncertainties in the template. The flux uncertainties were determined using procedures similar to the template spectra construction and a Monte Carlo approach. First, for each wavelength region and PC, the PC projection value was independently re-sampled 1000 times, assuming that the GPR hypersurface uncertainty has a Gaussian distribution. Crucially, we included all top 20 PCs for the uncertainty estimate, regardless of whether the PCs were included in the template construction. This is because those PCs not included in the template construction represent the observed SED variations that were not captured by the template. Each set of the re-sampled PC projections was then reverse transformed to the original flux space. The ratio of the template flux over the standard deviation of the re-sampled fluxes at each wavelength was then taken as the “S/N.” Finally, the S/Ns in the seven wavelength regions were merged by taking the weighted mean in the overlapping areas. An example of the template construction and the merging process can be found in Figure 10. Note that the choice of the number of iterations is based on the convergence of the median S/N.

Zoom In Zoom Out Reset image size

The PCA+GPR method we have presented here is robust, but has some limitations. First, it requires interpolations of the flux and their uncertainties onto a common wavelength grid prior to PCA. Second, PCA is a linear transformation and may require more dimensions to fully capture the variations in complex data distributions. Third, GPR results are sensitive to the kernel design (e.g., Stevance & Lee 2023). Furthermore, the dimension reduction and template generation parts do not communicate with each other, for example, one cannot improve the other iteratively. In Appendix B, we show that the neural network approach using cVAE can circumvent some of these drawbacks and achieve similar results. However, the spectral features in the telluric region are not well modeled with the cVAE approach compared with the current method.

We considered three strategies for selecting the PCs for template construction:

• 1.  
  Use the first N PCs and keep N fixed for all wavelength regions.
• 2.  
  Use the first N PCs that capture a total variance that exceeds a fixed threshold in each wavelength region.
• 3.  
  Use only the PCs with GPR coefficient of determination (R²) scores that exceed a fixed threshold in each wavelength region.

The R² measures how well the GPR hypersurface predicts the PC projections of sample spectra and is defined as follows:

where is the PC projection of the i^th spectrum; is the GPR predicted PC projection given the same phase and s_BV values of the i^th spectrum; and is the mean of the PC projections of the sample spectra. A larger R² value means that the GPR hypersurface captures more of the data variance. A value of R² approaching 0 means that the predictions are close to the mean of projections and largely independent of the parameters.

To compare the three strategies quantitatively, we utilized the χ² test:

where is the observed spectrum flux of the i^th spectrum, is the predicted template spectrum flux given the same phase and s_BV values of the i^th spectrum, and the observed flux errors were used as the uncertainties σ_i . Note that the template uncertainties are not considered since they might be correlated with the observed flux uncertainties. Each used for the test was color-corrected to match the color of before the comparison. The degrees of freedom (dof) were taken as the number of flux points minus the number of PCs used. Excluding the spectra with poor telluric correction and spurious flux values at the K-band edge (see Section 3.1), 225 observed spectra among the sample were used for the χ² test.

Each wavelength region requires a different number of PCs to reach a total variance threshold. For example, the telluric regions W4 and W6, as well as W7 (K band) require more PCs to capture 95% of the data variance than the other regions; see panel (a) of Figure 11. As expected, the GPR R² scores are the highest for the top-ranking PCs, as shown in panel (b) of Figure 11. However, they are not ordered perfectly, revealing some lower-ranked PCs with strong dependence on phase and s_BV . Note that the GPR R² scores are generally lower in the telluric regions W4 and W6. While this result may be revealing real spectral behavior, it more likely points to the observational challenges in these regions.

Zoom In Zoom Out Reset image size

The performance of the three proposed PC selection strategies was then assessed via the χ² test. In panel (c) of Figure 11, we plot the average number of PCs required across wavelength regions versus the median reduced χ² for each strategy. Using a GPR R² score threshold (strategy 3) shows a clear advantage over the other two strategies, as it requires fewer PCs in order to reach the same accuracy in the SED prediction. For example, to reach a reduced χ² ∼ 1.43, the average number of PCs needed is 15 for strategy 1, 13 (total variance of ≥99.5%) for strategy 2, and eight (R² ≥ 0.2) for strategy 3. This trend is consistent if we take the template uncertainties into consideration for the calculation, but the reduced χ² values are consistently lower. Hence, we adopted strategy 3 with a threshold of R² ≥ 0.2. The number of PCs used in each region is listed in Table 1. The GPR hyperparameters of the selected PCs are all converged to the optimal values, except the RBF length scale for s_BV of PC1 in the W4 (telluric) region; see details in Appendix C. Note that even though some higher-ranked PCs are omitted when constructing the template flux, it is crucial to include these when estimating the uncertainties in the spectral template.

Having PC projections that are highly correlated with phase and s_BV would be ideal for constructing spectral templates. However, spectral features with low GRP R² scores could aid the search for a secondary parameter for improving SN Ia standardization or even point to interesting physics. For example, the early-phase NIR spectra of the transitional SN Ia iPTF13ebh have been shown to have strong C i features, most notably in the Y band (Hsiao et al. 2015). These spectra also have the most extreme PC projection values compared with other spectra at early phases, such as PC2, PC6, PC7, PC8, and PC10 in W2 (Y band; see Figure 8).

The template spectra constructed using the PCA + GPR method are able to replicate the spectral diversity in the sample spectra. In Figure 12, we highlight the spectral dependence on s_BV in three wavelength regions and epochs: pre-maximum in the Y band, and post-maximum in the H and K bands.

Zoom In Zoom Out Reset image size

The left panels of Figure 12 present the pre-maximum spectra in the Y band at around −10 days. At these early phases, the photosphere arises in the outer and intermediate layers of the ejecta. Faster declining and fainter SNe clearly show stronger Mg ii λ 1.0092 and λ1.0927 μm features. The additional absorption feature on the blue wing of Mg ii λ1.0927 μm is captured by the template and appears to be only present in subluminous SNe Ia. This feature has been identified to be unburned C i λ1.0693 μm in several fast decliners: SN 1999by (Hoeflich et al. 2002), iPTF13ebh (Hsiao et al. 2015), SN 2012ij (Li et al. 2022), and SN 2015bp (Wyatt et al. 2021). An alternative identification of He i λ1.0830 μm has also been proposed in the context of the double-detonation scenario (Boyle et al. 2017).

At around two weeks past maximum, the prominent H-band break is at its peak strength (e.g., Hsiao et al. 2013). The template shows drastic differences in the profile shapes between slow and fast-declining SNe Ia (middle panels of Figure 12). These variations have been investigated in previous studies. The strength of the H-band break at its peak varies strongly with the light-curve decline rate (Hsiao et al. 2013, 2015). The velocity of the blue edge of the H-band break measures the ⁵⁶Ni distribution and is also strongly correlated with s_BV (Ashall et al. 2019a, 2019b). In the region between 1.6 and 1.7 μm, the spectral local minimum in a slow decliner is in contrast to a local emission feature in a fast decliner, which may be the result of a photosphere well within the central Ni-rich region (Hoeflich et al. 2002). Any viable explosion scenario needs to be able to explain these trends.

At approximately one month past maximum, the photosphere has receded to reveal the iron-peak elements in the inner layers. The K band is dominated by multiple lines of Co ii at this epoch (e.g., Gall et al. 2012). Interestingly, the template only shows a mild dependence of the feature shapes on the light-curve shape (right panels of Figure 12). The observed spectra show a range of feature strength but a weaker correlation with s_BV compared to the H band. This perhaps points to a secondary parameter or contribution from a different spectral line that is not associated with the decay products of ⁵⁶ Ni.

In Figure 13, we present the time-series comparison between our new template, the template of Hsiao et al. (2007) updated to include the NIR (Hsiao 2009), and the observed spectra within our sample. Note that the Hsiao template gives the average SED that is designed to match a stretch s = 1 SN Ia. For the Hsiao template, the variation with light-curve shape is accounted for through matching the time on the stretched time axis (t_stretched = t/s_BV ), assuming that a fast decliner also has a more rapid spectral evolution. Within the training sample, three well-observed SNe were selected for comparison: iPTF13ebh (s_BV = 0.63) representing a fast decliner, SN 2013gy (s_BV = 0.91) representing an average object, and SN 2013hh (s_BV = 1.20) representing a slow decliner. As expected, the Hsiao template performs well when tasked to match an average SN Ia. On the other hand, the new template performs well for all s_BV and shows a clear advantage at the extremes. This result demonstrates that the spectral feature dependence on the light-curve shape cannot be solely described by the speed of the evolution.

Zoom In Zoom Out Reset image size

Ideally, the performance of the new template should be assessed using an independent sample. In the left panel of Figure 14, the new and Hsiao templates are compared to the sample of Marion et al. (2009). Note that the phases relative to V_max listed by Marion et al. (2009) were converted to phases relative to B_max assuming that V_max occurs 2 days after B_max. Also note that Δm₁₅(B) values were converted to s_BV using Equation (4) of Burns et al. (2014). The sources of Δm₁₅(B) values are tabulated in Table 2. Overall, the new template provides a better match than the Hsiao template especially for low and high s_BV values, such as SN 2005am (s_BV = 0.72). There are three notable disagreements for the NIR Ca ii triplet: SN 2003W at −6 days, SN 2004bv at −3 days, and SN 2001bf at 6 days. These could point to real diversity in the feature that is not captured by the template or observational artifacts near the short wavelength limit of the detector.

Zoom In Zoom Out Reset image size

As discussed in Section 2.2, the small number of NIR spectra of subluminous SNe Ia, especially at the late phase, is a limiting factor for the template at this particular phase space. Here, the template is compared to the spectra of two subluminous SNe from the literature: SN 1999by (s_BV = 0.42; Hoeflich et al. 2002) and SN 2015bp (s_BV = 0.67; Wyatt et al. 2021), shown in the right panel of Figure 14. These spectra were not included in our training sample. The spectra of SN 2015bp at various phases, including at more than 1 month past maximum, are well represented by the new template. This is despite the time gap of 50 days in the four FIRE spectra of SN 2015bp included in the training sample. Because of the limited range in the light-curve shape of our sample, reproducing the spectral features of SN 1999by with such a low s_BV requires extrapolation on the hypersurface with GPR. Nevertheless, the spectral features of the template match exceptionally well to those of SN 1999by. Note especially that the C i features at early times were reproduced and are stronger than those of SN 2015bp. The Hsiao template simply fails at these low s_BV values.

We also checked the broadband colors of the new template using the SNooPy model. The broadband colors of the new spectral template and the light-curve model of SNooPy are in general agreement within 0.2 mag for all colors and s_BV values (Figure 15). The large discrepancies at the early phases may be due to the lack of observed data. The general agreement gives assurance that the merging of the seven wavelength regions at the last step of the template construction preserves the broadband color information. Note, however, that the spectral template should always be color-matched to the observed photometric colors of the SN Ia in question during the light-curve fitting process.

Zoom In Zoom Out Reset image size

Another test we did is to address the potential training bias toward SN 2012fr, the best-observed SN in the sample, with 25 spectra out of the total 331 sample spectra. Following the methodology of leave-one-out cross validation, we constructed another template without SN 2012fr following the same methodology to compare with the template built with the whole sample. The two templates were color-corrected to have matching broadband colors. The maximum flux difference is around 6% in the region near 0.85 μm, where SN 2012fr presents the detached high-velocity Ca ii feature at early times. In general, the differences between the two templates are negligible, with a median flux difference less than 1%.

We utilized the definitions of K-correction uncertainties from Hsiao et al. (2007) to evaluate the performance of the new template. K-corrections obtained using the template spectra (K_temp) were compared to those from observed spectra of our sample (K_obs) for different combinations of filters and redshifts.

The rms and the mean of the K-correction differences were then taken as the statistical and systematic uncertainties, respectively. The K-correction uncertainties were also compared to those computed using the Nugent (Nugent et al. 2002) and Hsiao (Hsiao et al. 2007) templates. The computations were done using the SNooPy package, while the template spectra were color-corrected to match the broadband YJH colors of the observed spectra before the computation.

Note that these uncertainty measurements were designed solely to evaluate how well the template can predict the spectral features of the observed spectra. The broadband colors of the template were corrected to match those of the observed spectra, as it is common practice in light-curve fitting to correct the template to the observed colors. These estimates do not include other effects, such as the uncertainty in the relative flux calibration discussed in Section 2.2.

In the most ideal case, one should assess the performance of the new template using a validation set of observed spectra that is completely independent of the template training set. In the NIR, the sample size is severely restricted by the telluric region, as discussed in Section 3.1. There is currently no spectral sample with reliable telluric regions besides the FIRE sample presented here, and after inspecting the telluric corrections, our sample was reduced by roughly one-third. Given this difficulty, we chose to evaluate its performance using a cross-validation technique described below.

First, a subsample with reliable telluric regions was chosen to form a “validation” pool. Next, we randomly selected 66 spectra from this pool to form the validation set, which constitutes 20% of the entire sample. The remaining 80% of the entire sample was then used as the training set. The training set was used to construct a template with the procedures outlined in Section 3, while the validation set was used to estimate the K-correction uncertainties. The process of splitting the sample, constructing the template, and estimating the K-correction uncertainties was repeated for 1000 iterations. The resulting uncertainty estimates converge rather quickly, generally after 50 iterations, while the 1000 iterations yielded Gaussian-like distributions.

Figure 16 presents the K-correction uncertainties as a function of redshift for several cases of single-filter K-corrections at low redshifts and cross-filter K-corrections at high redshifts. The mean and standard deviation of the 1000 iterations are shown as the solid curves and shaded regions, respectively. The observer-frame filters for the high-redshift cases are the Wide Field Instrument (WFI) F129, F159, and F184 filters of the RST (Hounsell et al. 2018). ³⁷  The remaining filters are of RetroCam on the du Pont Telescope, measured by Rheault et al. (2014).

Zoom In Zoom Out Reset image size

As shown in the upper panels of Figure 16, the statistical K-correction uncertainties of the new template are significantly reduced compared to the previous templates. For example, at z = 0.05, improvements of 0.09, 0.03, and 0.02 mag in the Y, J, and H bands, respectively, can be gained by switching from the Hsiao template to the new one. Note that the result of Boldt et al. (2014) that the smallest uncertainties can be found in the Y band is not reproduced here. From a larger sample, we found that all YJH bands yield similar uncertainties. For cross-filter K-corrections, a minimum would be present at a redshift where the observed and rest-frame filters are aligned (Kim et al. 1996). On average, there is a 0.02 mag improvement in statistical uncertainty compared to the Hsiao template in cross-filter K-corrections.

The lower panels of Figure 16 demonstrate that the systematic K-correction uncertainties are essentially eliminated by adopting the new template. This is true for all of the redshifts and filter combinations considered here, and represents as 0.1 mag improvements in some cases. The improvements in the Y band are the most evident. Since the broadband colors of the template spectra were matched to the observed spectra before the comparison, these uncertainty estimates indicate the ability of the template to predict the spectral feature profile shapes.

Here, we explore the impact of adopting the new SN Ia NIR spectral template on cosmological analyses, both at low (e.g., Hubble constant measurements) and high (e.g., dark energy measurements) redshifts.

We first examine the impact on the NIR peak magnitude of adopting the new template on the low-redshift end in single-filter K-corrections. The results are presented in comparison to the Hsiao template in Figure 17. Note that we chose to make the calculations at B-band maximum for convenience, and the analysis yielded a rough estimate of potential magnitude change around peak magnitude. Both templates are scaled and color-corrected to match the YJH magnitudes of the SNooPy light-curve template, which is a function of s_BV . Recall that the Hsiao template at peak is a single spectrum and does not account for the variation of spectral features. The spread of roughly 0.1 mag in K-corrections at z = 0.1 results entirely from color corrections.

Zoom In Zoom Out Reset image size

On the other hand, the new template at peak is a series of spectra that are a function of s_BV . Figure 17 shows the difference between the two templates in terms of the K-correction as a function of redshift. It essentially reflects the fact that the diversity is not only in the broadband colors but also in the spectral features. Not accounting for the variation of spectral features as a function of s_BV can cause errors in K-corrections as large as 0.3 mag in Y and 0.1 mag in J at z = 0.1. Also note that the two templates show the largest difference in slow decliners (high s_BV values). The prominent H-band feature is at maximum strength around 12 days after the maximum. Both slow and fast decliners show K-correction differences as large as 0.1 mag at z = 0.1.

For cosmological analyses using distant SNe Ia, the redshift evolution of the SN Ia population is an important factor to consider. Based on previous high-redshift surveys, the SN Ia population at high redshifts tends to comprise slower decliners and have higher optical light-curve stretch values on average (e.g., Howell et al. 2007; Nicolas et al. 2021). This would introduce an additional systematic shift in the K-corrections if the template SEDs do not reflect the spectral feature variations with light-curve shape.

We simulate such a scenario, considering a shifting s_BV distribution at three different redshifts. In the left panels of Figure 18, three populations of SNe Ia are presented at three different redshifts, simulating the distribution of SN Ia light-curve shapes observed locally and at high redshifts (Howell et al. 2007; Riess et al. 2007). The mean of the s_BV distributions shifts by as much as 5% from the local sample to z = 1.12. We then calculated K-corrections to the rest-frame Y band and examined the difference between the new and Hsiao templates at B-band maximum.

Zoom In Zoom Out Reset image size

The resulting K-correction differences from the corresponding s_BV distribution are presented in the right panels of Figure 18. The observer-frame filters for the two high-redshift bins, H and K, align well with the rest-frame Y band at each redshift. The effective wavelength of the de-redshifted H (z = 0.56) and K (z = 1.12) bands match that of the rest-frame Y band within 0.02 μm. Even so, the K-correction differences are systematic and significant, with means of 0.10, 0.01, and 0.07 mag at z = 0.03, 0.56, and 1.12, respectively. The results again highlight the importance of accounting for the spectral feature variations as a function of light-curve shape (e.g., Jones et al. 2022).