The Twins Embedding of Type Ia Supernovae. I. The Diversity of Spectra at Maximum Light

Initial methods to standardize the luminosities of SNe Ia involved correcting the observed peak brightnesses of SNe Ia for correlations with the widths of their light curves (Phillips 1993) and their B − V colors at maximum light (Riess et al. 1996; Tripp 1998). Current cosmological analyses fit the light curves of each SN Ia using an empirical model of the time-evolving spectral energy distribution (SED). The most commonly used such model, SALT2 (Guy et al. 2007, 2010; Betoule et al. 2014), has one component c for the color and one component x₁ for the intrinsic diversity that effectively captures the width of the light curve. With the SALT2 model, the luminosities of SNe Ia can be estimated with an accuracy of ∼0.15 mag, and the distance to an individual SN Ia can be inferred with an accuracy of ∼8%.

Light-curve width and color have been shown to not be sufficient to capture all of the diversity of SNe Ia. Most notably, the inferred SALT2 luminosities of SNe Ia show differences of ∼0.1 mag when comparing SNe Ia in host galaxies with different masses, metallicities, colors, or star formation rates (Kelly et al. 2010; Sullivan et al. 2010; D’Andrea et al. 2011; Gupta et al. 2011; Childress et al. 2013; Hayden et al. 2013; Rigault et al. 2013, 2015, 2020; Roman et al. 2018). These observed differences with host-galaxy properties imply that there is additional unmodeled variability of SNe Ia.

One major open question is to determine the dimensionality of SNe Ia, or how many different modes of variability there are. This is a somewhat ill-posed question: SNe Ia are the result of highly complex explosions and undoubtedly require a parameter space with a very large number of dimensions to fully capture all of their intrinsic diversity. Certain physical components presumably have a very large effect on the observed spectral time series of SNe Ia, such as the composition of the white dwarf (Timmes et al. 2003), total ejecta mass (Scalzo et al. 2014), amount of ⁵⁶Ni produced (Arnett 1982), photosphere ejecta velocity (Foley et al. 2011), or asymmetry of the explosion (Maeda et al. 2011). Other physical components may only be observable with highly specialized observations with limited effect on the full spectral time series, perhaps including blueshifted Na i D lines (Phillips et al. 2013), the presence of unburned carbon (Thomas et al. 2011), or polarization (Wang & Wheeler 2008). For the purposes of cosmology, it is important to consider that we are only required to model components that have a non-negligible effect on the estimated intrinsic luminosity of SNe Ia at maximum light. Furthermore, all of these observational signatures of diversity are likely correlated and caused by a smaller number of underlying parameters of the explosion.

There have been many efforts to empirically identify modes of variability of SNe Ia other than light-curve width and color. Nugent et al. (1995) showed that ratios of the equivalent widths of the Si ii 5972 Å and Si ii 6355 Å lines map out a spectral sequence of SNe Ia. Branch et al. (2006) used these two lines to separate the spectra of SNe Ia into four subgroups and showed that at least two dimensions are required to parameterize the intrinsic diversity of SNe Ia. Wang et al. (2009) and Foley & Kasen (2011) showed that there is diversity in the velocity of the Si ii 6355 Å feature that is uncorrelated with the width of the light curve and that affects standardization. Further studies have collected thousands of spectra of SNe Ia and have confirmed that the intrinsic diversity of SNe Ia is multidimensional (Blondin et al. 2012; Silverman et al. 2012; Folatelli et al. 2013).

The Nearby Supernova Factory (SNfactory; Aldering et al. 2002) has collected spectrophotometric time series of hundreds of SNe Ia that have enabled a wide range of new analyses of the diversity of SNe Ia. Using this data set, Nordin et al. (2018) showed that there are at least two components of the intrinsic diversity of SNe Ia in the U band that affect standardization. Léget et al. (2020) showed that the diversity of 13 different spectral features can be explained with an underlying three-dimensional parameter space that they call SUGAR. Saunders et al. (2018) developed the SNEMO model using the full spectral time series from SNfactory directly and found that 15 linear components are required to parameterize the diversity of this data set.

The intrinsic diversity of SNe Ia is complicated to parameterize because a typical variation, such as a change in the expansion velocity of the photosphere, leads to highly nonlinear effects in the observed spectra and photometry. Models such as SALT2, SUGAR, or SNEMO are linear, meaning that they attempt to describe the time-evolving SED of an SN Ia at each wavelength and time as the sum of a set of linear components. Fitting a linear model to nonlinear phenomena will result in a model with many more (redundant) components than a similar nonlinear model. Sasdelli et al. (2016) used deep learning to model the spectra of SNe Ia and showed that a nonlinear four-dimensional parameter space can capture the intrinsic diversity of SNe Ia as well as a 15 component linear model. However, their analysis used the derivative of the spectrum rather than the spectrum directly, which removes information such as the brightness and color that are necessary for standardization. Note that the linearity of a model is distinct from the linearity of standardization using the parameters of that model. For example, Rubin et al. (2015) implement nonlinear standardization in terms of SALT2 parameters, but they are still restricted to using the linear SALT2 model that cannot capture nonlinear spectral variation.

Using observations from SNfactory, Fakhouri et al. (2015, hereafter F15) introduced an alternative method of standardizing SNe Ia. When an underlying physical parameter of the explosion is varied, we would expect to see a relatively smooth sequence in the spectra of SNe Ia. The authors developed a method of estimating the “spectral distance” between any pair of SNe Ia and call pairs with low spectral distances “twins.” Standardization can then be done by inferring the luminosity of a new SN Ia directly from its set of twins. As long as one has a large enough reference sample to span the full range of diversity of SNe Ia, one can find SNe Ia with similar spectra to any new SN Ia. The exact functional form of how the spectrum is affected by changes in some underlying physical parameter of the explosion is irrelevant because the twins method only does local comparisons. The twins method does not provide a parameterization of SNe Ia. Rubin (2020) showed that the statistics of the twins pairings from F15 are consistent with the intrinsic variability of SNe Ia being described by an underlying three to five-dimensional parameter space, although they did not explicitly construct this parameter space.

In this work, we extend the twins methodology of Fakhouri et al. (2015) to develop a nonlinear parameterization of the spectral diversity of SNe Ia. Assuming that we have a sample of observed SNe Ia that continuously spans the full range of intrinsic diversity of SNe Ia, we can recover the underlying parameterization by identifying sequences of observed spectra where each spectrum in the sequence has a small spectral distance to its neighbors. We call this parameterization the “Twins Embedding.”

For this analysis, we use a large data set of high signal-to-noise spectra of SNe Ia from SNfactory that is described in Section 2. We perform a sequential analysis to decompose the variability of SNe Ia. First, we model the differential time evolution of the spectra of SNe Ia near maximum light in Section 3, and we estimate the spectra of all of our SNe Ia at maximum light. In Section 4, we introduce a second procedure that we call “Reading Between the Lines” (hereafter RBTL) to estimate the contributions to the spectra from distance uncertainties and dust extinction and produce dereddened spectra of SNe Ia that nominally have only intrinsic variability remaining. Finally, we perform a nonlinear decomposition of the remaining variability of these dereddened spectra in Section 5, and we produce a parameterization of the intrinsic diversity of SNe Ia that we call the “Twins Embedding.” In Section 6, we explore the properties of the Twins Embedding and show how it can be used to recover a wide range of previously studied indicators of intrinsic diversity of SNe Ia. In a companion article (Boone et al. 2021, hereafter Paper II), we show how the Twins Embedding can be used to improve the standardization of SNe Ia.

The RBTL model is implemented as follows. For each SN i, we begin with a spectrum at maximum light with associated uncertainties from Section 3. We represent the mean spectrum of an SN Ia at maximum light as f_mean(λ_k ). Each SN i then has a parameter Δm_i representing its difference in brightness compared to the mean spectrum in magnitudes and a parameter representing the coefficient of the extinction–color relation C(λ_k ) that best matches the SN’s spectrum to the mean function. We choose to use the extinction–color relation C(λ_k ) from Fitzpatrick (1999) with a fiducial R_V = 2.8. The modeled flux of the spectrum at maximum light of SN i can then be written as

We assume that the intrinsic dispersion of SNe Ia, η(λ_k ), is the same for all SNe Ia and uncorrelated in wavelength. For computational reasons, we implement this uncertainty as a fraction of the modeled flux. However, we interpret it in the following text as the corresponding difference in magnitudes. The total uncertainty of the spectrum at maximum light for an SN relative to the modeled spectrum is therefore modeled as

As in Section 3, we implement this model using the Stan modeling language (Carpenter et al. 2017), and we use Stan to obtain the MAP estimate of the posterior distribution. Finally, we apply the inverse of the magnitude and extinction corrections to obtain “dereddened spectra” f_dered.,i(λ_k ) for each of our spectra at maximum light:

The resulting values of the RBTL intrinsic dispersion η(λ_k ) can be found in Table 2, and the values of Δm and for each SN can be found in Table 3. The dereddened spectra at maximum light are available on the SNfactory website at https://snfactory.lbl.gov/snf/data/.

Table 3. Measurements of All of SNe Ia in Our Sample

	SALT2 Parameters		RBTL Parameters		Twins Embedding Coordinates
SN Ia Name	x1	c		Δm	ξ1	ξ2	ξ3
SNF 20050728–006	0.46 ± 0.31	0.117 ± 0.035	0.276 ± 0.014	−0.094 ± 0.039	−0.937	−0.557	0.339
SNF 20050729–002	0.92 ± 0.34	−0.107 ± 0.038	−0.260 ± 0.018	0.128 ± 0.030	0.621	−0.536	0.152
SNF 20050821–007	0.42 ± 0.29	−0.043 ± 0.036	−0.205 ± 0.014	0.029 ± 0.040	−2.580	2.821	0.355
SNF 20050927–005	0.19 ± 0.35	0.018 ± 0.035	0.033 ± 0.014	0.301 ± 0.056	−0.176	−1.688	1.221
SNF 20051003–004	1.16 ± 0.16	−0.100 ± 0.030	−0.349 ± 0.013	0.071 ± 0.065	2.030	0.489	−1.463
SNF 20060511–014	−0.78 ± 0.18	−0.041 ± 0.035	−0.030 ± 0.013	0.037 ± 0.048	−0.396	−1.828	1.501
SNF 20060512–001	0.40 ± 0.14	−0.007 ± 0.030	0.046 ± 0.012	−0.096 ± 0.057	1.810	2.584	2.179
SNF 20060512–002	−0.26 ± 0.23	0.086 ± 0.033	0.079 ± 0.015	−0.197 ± 0.045	−1.993	1.080	−1.411
SNF 20060521–001	−1.63 ± 0.29	−0.111 ± 0.036	−0.057 ± 0.013	0.047 ± 0.034	1.238	−2.242	−0.497
SNF 20060521–008	−1.35 ± 0.32	0.122 ± 0.042	0.404 ± 0.016	−0.206 ± 0.041	−3.916	−0.466	−1.604
⋯

Note. For each SN Ia, we show its SALT2 fit parameters, its extracted RBTL extinction and magnitude residual Δm, and its coordinates in the Twins Embedding. The uncertainties on Δm contain contributions from both measurement uncertainties and peculiar velocities. The first 10 lines of this table are shown in this table. The full table is available in machine-readable version.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

The estimated spectra at maximum light of the 203 SNe in our sample both before and after dereddening are shown in Figure 4 along with the modeled intrinsic dispersion. The dereddened spectra show remarkable similarity, especially in wavelength regions away from the main absorption lines, with intrinsic dispersions of <0.10 mag at almost all wavelengths. Note that because we fit for Δm and , any mode of the intrinsic dispersion that affects the spectrum in a similar way to these components will not be captured in the recovered intrinsic dispersion. The very low recovered intrinsic dispersion of ∼0.02 mag between 6600 and 7200 Å implies that there is almost no uncorrelated dispersion at these wavelengths. As a result, the RBTL model effectively relies heavily on flux measurements at these wavelengths for standardization. Perhaps not coincidentally, in this spectral region, the opacity is dominated by electron scattering rather than line absorption. The existence of regions of low dispersion having some unique association with the physics of radiative transfer in SN Ia atmospheres further supports the rationale for the RBTL technique.

Zoom In Zoom Out Reset image size

The Ca ii H&K lines, the Ca ii IR triplet, and the Si ii 6355 Å feature are the locations in the spectra of SNe Ia with the largest intrinsic dispersions. Interestingly, there are several regions of the spectra that are blanketed by lines but that still show relatively low intrinsic dispersion. From 3900 Å to 5900 Å, the spectra of SNe Ia show a variety of absorption lines, but the intrinsic diversity is recovered to be ∼0.08 mag at most wavelengths, with the exception of a handful of stronger lines that introduce diversity at up to 0.13 mag.

This model is effectively using the intrinsic diversity of SNe Ia to weight each wavelength when determining how to fit for the brightness and extinction of each SN. Hence, the lines with strong diversity are deweighted, and the model effectively “reads between the lines” to estimate the brightness and extinction.

The process of recovering a low-dimensional nonlinear parameter space from high-dimensional observations is referred to as “manifold learning.” Manifold learning has been shown to be very effective for constructing nonlinear parameterizations of the spectra and light curves of astronomical objects (Richards et al. 2009, 2012; Daniel et al. 2011; Matijevič et al. 2012; Sasdelli et al. 2016). A major challenge with nonlinear dimensionality reduction is that any transformation of a given parameterization is an equally valid parameterization. Many different manifold learning techniques exist, each of which imposes different assumptions and constraints on the recovered parameterization. Our goal is not to find an “optimal” parameterization of SNe Ia (which is not well defined), but instead, one that is useful for cosmological applications.

For this analysis, we choose to construct a parameterization of SNe Ia that preserves the spectral distances between twin SNe as described in F15. We define the spectral distance γ_ij between two SNe Ia labeled i and j as

Our goal is to construct a parameterization of SNe Ia where the Euclidean distance between the coordinates of any two twin SNe Ia is equal to the spectral distance between those two SNe Ia. We do not require that the spectral distances of non-twins be preserved.

To accomplish these goals, we make use of the Isomap algorithm from Tenenbaum & Silva (2000). This algorithm is designed to embed points from a high-dimensional space into a low-dimensional one while preserving the distances between nearby points in the high-dimensional space. By using the spectral distance of Equation (11) as the distance measure for the Isomap algorithm, we will therefore generate an embedding that preserves the distances between twin SNe Ia.

The Isomap algorithm proceeds as follows. First, for each SN Ia, we find its K nearest neighbors using spectral distances. For each pair of SNe Ia in the sample, we then find the shortest path between the two SNe Ia passing through pairs of neighbors. We calculate the “geodesic distance” for each pair of SNe Ia as the sum of spectral distances between neighbors along each step in this path. After computing the geodesic distances between each pair of SNe Ia, we obtain a distance matrix. To generate an embedding, we center the distance matrix and compute its eigendecomposition (see Tenenbaum & Silva 2000 for details). For a D-dimensional embedding, we keep the D eigenvectors of the centered distance matrix with the largest eigenvalues. Each eigenvector then corresponds to a “component” of the embedding and captures a distinct mode of variability of SNe Ia. The values of the eigenvectors are the “coordinates” of each SN Ia within the embedding.

The Isomap algorithm has two parameters that must be set to produce an embedding: the number of neighbors K, and the dimensionality D of the embedding. As will be discussed in Section 5.5, the number of neighbors K does not have a major impact on the resulting embedding, and for a range of different values tested between ∼6 and 50, we obtain nearly identical embeddings. We choose to use 10 neighbors for the rest of our analysis. We will discuss the dimensionality of SNe Ia in Section 5.3.

The Isomap algorithm generates a low-dimensional embedding of SNe Ia, but it does not produce a model of the spectrum of an SN Ia given its coordinates in the embedding. Instead, to do this we use GP regression to model the spectra as a function of the Isomap coordinates. The details of this model are described in the Appendix. We use a separate GP model for each wavelength, and we optimize the hyperparameters of each GP independently.

For some of the SNe in our sample, the uncertainty on the spectrum at maximum light is comparable to the recovered intrinsic dispersion η(λ_k ). The Isomap algorithm is unable to take this uncertainty into account, so if spectra with large uncertainties are included in our analysis, we could confuse intrinsic diversity with the uncertainty in our estimate of the spectrum at maximum light. To mitigate this, we remove any SNe at this stage of the analysis whose uncertainties on the spectrum at maximum light are large compared to the intrinsic dispersion. We choose to require that the total measurement variance of the estimate of the spectrum at maximum light be less than 10% of the total intrinsic variance of SNe Ia. Choosing this fractional variance threshold requires somewhat of a trade-off: with better-measured spectra, we could potentially recover more components of the intrinsic dispersion, but a stricter threshold reduces the number of SNe Ia in the sample and thus our ability to reconstruct the parameter space. Out of the original sample of 203 SNe Ia, 173 SNe Ia have uncertainty on their spectrum at maximum light that passes this stringent fractional variance threshold. As will be shown in Section 5.5, whether this cut is applied has almost no effect on the embedding coordinates for the SNe Ia that pass it.

Note that in the RBTL analysis we included the estimates of the uncertainties for the spectra at maximum light, so the spectra that were cut due to the fractional variance threshold will not have a major impact on the RBTL analysis. We retrained the RBTL model on only the spectra that pass the uncertainty of the spectrum at the maximum light requirement and found that the changes were negligible (the estimated brightnesses change by <0.005 mag). As a result, we choose to use the RBTL model trained on all of the spectra for further analysis.

We use the Isomap + GP model to investigate the dimensionality of SNe Ia. For a given number of Isomap components D, we examined what fraction of the variance is explained by the model. The results of this procedure are shown in Figure 5. We find that the first three components of the Isomap + GP model each explains a significant fraction (49.3, 26.2, and 11.1% respectively) of the total variance, and together, they explain 86.6% of the total variance. Measurement uncertainty accounts for 2.9% of the remaining variance for this sample, so our model explains 89.2% of the intrinsic variance when this is taken into account. Additional components beyond the third do not explain a significant amount of the remaining variance.

Zoom In Zoom Out Reset image size

For comparison purposes, we also perform a linear PCA decomposition of the same data. We find that the linear model requires seven components to explain as much of the variance as our three-component nonlinear model, which supports our claim that there is a significant amount of nonlinear variability in the spectra of SNe Ia.

We show the unexplained dispersion as a function of wavelength in Figure 6. After the RBTL procedure, the residual dispersion is ∼0.3 mag for the Ca ii and Si ii features. In contrast, with the three-component Isomap + GP model, there is less than 0.1 mag of residual dispersion at the wavelengths associated with these features, and the dispersion is less than 0.05 mag at almost all wavelengths. Again, adding components beyond the third has a negligible effect on the intrinsic dispersion.

Zoom In Zoom Out Reset image size

Alternatively, we look at how well the spectral distances from Equation (11) are preserved in the embedding. We find that the Euclidean distances between two SNe Ia in the Isomap embedding have Pearson correlations with the F15 spectral distances of 0.77 for one component, 0.90 for two components, and 0.95 for either three or four components with very little improvement if additional components are added. With at least three components, we are therefore able to accurately preserve the spectral distances (and thus twin pairings) of F15.

Sasdelli et al. (2016) previously showed that deep learning can produce a nonlinear four-dimensional representation of SNe Ia. However, their analysis modeled the derivatives of the spectra in wavelength rather than modeling the spectra directly and did not include an explicit model of how spectra vary near maximum light or a means of measuring the brightness and color of a spectrum (which is necessary for cosmological applications). Because we explicitly model the extrinsic contributions to the spectrum, our model is able to explain significantly more of the variance in the spectra: with three components, we can explain 86.6% of the intrinsic variance compared to 82% with four components for the analysis of Sasdelli et al. (2016). Furthermore, the distances between two SNe Ia in the Twins Embedding capture the spectral distances of F15 while they have no meaning in the analysis of Sasdelli et al. (2016).

We can use GPs to predict the spectrum of an SN Ia given its coordinates from the Isomap decomposition. To test how well our model would perform on new observations, we generate leave-one-out predictions where we condition the GP on spectra from all of the SNe Ia in our sample except for one SN Ia and predict the spectrum at the Isomap coordinates of the remaining SN Ia. Results of this procedure for three SNe Ia are shown in Figure 7. In this plot, we show the results for a “normal” SN Ia, a 91T-like SN Ia (Filippenko et al. 1992b), and a 91bg-like SN Ia (Filippenko et al. 1992a). We find that the three-component model is able to predict accurate spectra for the full range of SNe Ia including 91T-like and 91bg-like SNe Ia.

Zoom In Zoom Out Reset image size

There is still a small amount of residual dispersion after applying any of these models. To study this, we performed leave-one-out predictions for all of the SNe Ia in our sample and evaluated the correlation matrix of the residuals. The results of this procedure are shown in Figure 8. For the base RBTL model (left panel of Figure 8), we see a strong off-diagonal structure in the correlation matrix implying that the residuals at different wavelengths are highly correlated. In contrast, for the three-component Isomap + GP model (right panel of Figure 8), there is very little correlation between different wavelengths. This implies that the remaining variance is mostly uncorrelated across wavelengths and explains why adding additional components does little to improve the model.

Zoom In Zoom Out Reset image size

We conclude that a three-dimensional embedding is sufficient to explain the vast majority of the intrinsic diversity of SNe Ia at maximum light (along with extrinsic contributions from dust and brightness removed in Section 4). As this embedding was effectively constructed using pairs of twin SNe Ia, we refer to it hereafter as the “Twins Embedding.” We label the three components of the Twins Embedding ξ₁, ξ₂, and ξ₃. The RBTL parameters and coordinates of each SN Ia in the Twins Embedding can be found in Table 3.

As described in Section 5.1, nonlinear dimensionality reduction is challenging because there is in general no unique solution. We investigated how stable the Twins Embedding is by making various modifications to the input data set and applying our algorithms to generate variants on the embedding. In general, there is no guarantee that the recovered embedding will be axis-aligned with our original one, especially if we modify the binning of the spectrum since that will change the weights of different regions. To perform a robust quantitative comparison, we fit a GP model to predict the coordinates in the original embedding from the coordinates in some alternative embedding. We use the same GP model that was used to predict the flux of each spectrum across the Twins Embedding, the details of which can be found in the Appendix. For each variant, we calculate the fraction of variance that is explained by the GP model. The results of this procedure are shown in Figure 9.

Zoom In Zoom Out Reset image size

Changing the number of neighbors for the Isomap algorithm has very little effect on the recovered embedding for values between ∼5 and 50 neighbors. For fewer than five neighbors, we find that we are unable to recover the third component. This is likely due to the fact that using a small number of neighbors adds noise to the distance matrix. Removing the uncertainty cuts described in Section 5.2 has almost no effect on the recovered embedding. Similarly, randomly discarding 20% or 50% of the SNe Ia in the sample does not noticeably affect the embedding.

Our sample was selected to have very high signal to noise, so resampling the uncertainties of each spectrum has little effect on the recovered embedding. We find that we are able to degrade the spectrum down to an S/N of ∼5 per 1000 km s⁻¹ before we start to see a degradation in the recovered embedding. Interestingly, we are able to bin the spectrum down to ∼10000 km s⁻¹ before we see significant changes in the embedding, which corresponds to a very low-resolution spectrum (R ∼ 30).

Finally, we investigated how changing the wavelength range of our spectra from the baseline of 3300–8600 Å affects the recovered spectrum. Decreasing the cutoff wavelength at the red end has almost no effect on the embedding, but increasing the cutoff at the blue end does, with a significant degradation in the third component for cutoffs above ∼4000 Å. As discussed in detail in Nordin et al. (2018), there is a large amount of variability in the U band, particularly near the Ca ii H&K lines, and not having observations of this region of the spectrum significantly degrades the performance of our model.

We now examine what effects each of the three components of the Twins Embedding has on the spectrum of an SN Ia. For each component, we calculate the median spectrum in 20 evenly populated bins of the coordinates for that component. We plot these median spectra in Figure 10. We find that each component captures a different smooth, nonlinear spectral sequence.

Zoom In Zoom Out Reset image size

The first component of the Twins Embedding (ξ₁) primarily affects the pseudo-equivalent widths (pEWs) of the Ca ii features. This component maps out a full spectral sequence of the changing pEW of the Ca ii H&K feature, and a similar spectral sequence is seen for the Ca ii IR triplet. Additionally, ξ₁ shows a spectral sequence for the Si ii 6355 Å feature. This component is the only one that has a major impact on the emission profile of this feature and appears to be capturing diversity in both the optical depth and line velocity of the Si ii 6355 Å feature. ξ₁ has very little effect on the depth of the Si ii 5972 Å feature that is typically associated with the width of the light curve.

The second component (ξ₂) primarily affects the pEWs of the Si ii lines. A spectral sequence is identified in the Si ii 4130 Å, Si ii 5972 Å, and 6355 Å features. Interestingly, this component directly changes the depths of these lines without having a major effect on the line velocities themselves. ξ₂ also has a large effect around 3700 Å that appears to be due to a set of Si ii lines at rest-frame wavelengths between 3853 and 3863 Å. This component does not appear to affect the Ca ii H&K absorption typically associated with these wavelengths, and by comparing with ξ₁, we can see that the Twins Embedding cleanly separates two modes of variability in this wavelength range. ξ₂ has a very complex spectral sequence for redder wavelengths, spanning both the O i and Ca ii features.

The third component ξ₃ primarily affects the ejecta velocity profiles of the SN. This component identifies a spectral sequence in the velocity of the S ii doublet feature near 5400 Å and in the velocity profile of the Si ii 6355 Å feature. Interestingly, ξ₃ seems to capture an overall shift in velocity in the emission profile: as this component is varied, we see the full absorption profile shift in wavelength for both the S ii doublet and the Si ii 6355 Å feature. This contrasts with ξ₁, which also alters the observed line velocity of these emission features by increasing the optical depth of these elements at larger velocities. Similar effects are seen for the velocities of many other lines in the spectra. ξ₃ does not have a major effect on any of the line depths in the spectrum. Interestingly, we see a sequence in the pseudo-continuum level around the Si ii 5972 Å feature that does not affect the pEW of this feature.

If the Twins Embedding truly captures all of the observable intrinsic diversity of SNe Ia, then we would expect to be able to recover any indicator of intrinsic diversity of SNe Ia from some transformation of the Twins Embedding. We tested this hypothesis by training GP models to predict a wide variety of indicators from the Twins Embedding, and then by measuring the fraction of the variance in those indicators that is explained by the GP models.

First, we extracted the pEWs and line velocities of the Ca ii H&K feature and several of the Si ii features using the methodology described in Chotard et al. (2011) on each of the dereddened spectra. For SNe Ia that were included in the U-band analysis of Nordin et al. (2018), we include their spectral indicators uNi, uTi, uSi, and uCa at maximum light, and we also include the uCa measurement prepeak and uTi measurement postpeak that the authors show affect standardization. For SNe Ia that were included in the SUGAR analysis of Léget et al. (2020), we include their SUGAR coordinates. Similarly, for SNe Ia that were included in the SNEMO analysis of Saunders et al. (2018), we include their SNEMO7 coordinates. Finally, we compare our embedding to the SALT2 x₁ values measured as described in Section 2. As SALT2 is not able to accurately model all of the peculiar SNe Ia that we have included in our sample, for this comparison, we reject SNe Ia with bad SALT2 fits where the normalized median absolute deviation (NMAD) of the SALT2 model residuals is more than 0.12 mag or where more than 20% of the SALT2 model residuals have an amplitude of more than 0.2 mag.

We use the same GP model that was used to predict the flux of each spectrum over the Twins Embedding. The details of this GP model can be found in the Appendix. When available, we include the measurement uncertainties of the indicators in the GP model. In Figure 11, we show the results of this procedure. We find that the Twins Embedding captures the majority of the diversity in all of these indicators, with 68%–94% of the variance explained depending on the indicator. We verified with leave-one-out GP predictions that we are able to perform out-of-sample predictions of each of these indicators within the quoted accuracies.

Zoom In Zoom Out Reset image size

Adding an additional component to the Twins Embedding does not significantly improve the ability to recover any of these indicators while removing the third component ξ₃ significantly decreases performance, further confirming that three dimensions are required to capture the diversity of SNe Ia. Interestingly, with transformations of the Twins Embedding, we are able to recover the majority of the variance in each of the components of the SNEMO7 model of Saunders et al. (2018). Given the location of an SN Ia in the three-dimensional Twins Embedding, we can therefore predict its location in the six-dimensional SNEMO7 parameter space with high accuracy. Our results imply that the six components of SNEMO7 are not all independent and can be captured by a three-dimensional submanifold. SNEMO7 is a linear model, and this is consistent with the results that we found for a simple linear model in Section 5.3. Note that the SNEMO7 parameter space was constructed using full spectral time series of SNe Ia while the Twins Embedding was constructed using only information at maximum light. These results then suggest that the vast majority of the information in the spectral time series of SNe Ia is captured at maximum light.

In our comparison to the SUGAR model (Léget et al. 2020), we find that we are able to reproduce the first two components of the SUGAR model with high accuracy and the majority of the variance of the third SUGAR component. Hence, the SUGAR parameter space is very similar to the Twins Embedding despite having been constructed very differently (a linear decomposition of spectral features). We reran our algorithm to calculate the fraction of explained variance for each indicator using a transformation of the SUGAR parameter space rather than a transformation of the Twins Embedding. We find that both parameterizations of SNe Ia are able to explain similar fractions of the variance for most of the indicators that we examined. However, the SUGAR model performs significantly worse on the U-band features of Nordin et al. (2018). For example, it is only able to recover 42% of the variance of the uNi indicator at maximum light compared to 71% for the Twins Embedding using the same set of SNe Ia. This can be explained by the fact that the SUGAR model was trained using only a specific set of spectral features, and is potentially missing information that is not captured in those spectral features. On the other hand, the Twins Embedding was trained using the entire spectrum at maximum light. Furthermore, the SUGAR analysis removed “outlier SNe Ia” that are not well described by a linear model while the Twins Embedding includes all of these SNe Ia because it is not restricted to linear variability.

The spectra of SNe Ia near maximum light are often labeled in the literature using the “Branch classification scheme” (Branch et al. 2006, hereafter B06). In this classification scheme, the spectra of SNe Ia are subdivided into four classes based on the pEWs of the Si ii 5972 and 6355 Å absorption features. As has been previously shown for large samples of spectra of SNe Ia (e.g., Blondin et al. 2012), the Branch classifications do not identify distinct subtypes of SNe Ia, and we find a continuous distribution across all of the label boundaries.

We compare the Branch classifications to our Twins Embedding in Figure 12. We find that the first two components of the Twins Embedding (ξ₁ and ξ₂) cleanly separate the different Branch classes from each other: with a set of cuts in the Twins Embedding, we would be able to recreate the Branch classification scheme nearly perfectly. There is some overlap at the border of each class due to uncertainties in the pEW measurements used for the Branch classification. Interestingly, while core normal SNe Ia lie in a tight cluster in the parameter space used for Branch classifications, they are spread over a fairly large region of the Twins Embedding.

Zoom In Zoom Out Reset image size

The Twins Embedding thus implies that there is significant spectral diversity among core normal SNe Ia. To probe this claim, we calculated the median spectrum for core normal SNe Ia in ten equally populated bins of ξ₁. The results of this procedure are shown in Figure 13. All of these median spectra have similar pEWs for the Si ii 6355 Å absorption feature and have a small pEW for the Si ii 5972 Å feature, per the definition of a core normal SN. However, the pseudo-continuum level near the Si ii 5972 Å feature varies quite dramatically between the different median spectra. As analyses such as B06 focused only on the pEW and ignore the pseudo-continuum level, they were not able to identify this difference. We also see large differences for these spectra near 3800 Å, which we associate with a set of Si ii lines between 3853 and 3863 Å. These differences were not seen in B06 due to a lack of spectral coverage. Finally, we see a spectral sequence for an absorption feature near 8000 Å. This difference was identified in B06 for SN 2001el and was suggested to be due to high-velocity Ca ii. Hence, with only the first two components, the Twins Embedding is able to reproduce the Branch classification scheme and is able to identify additional diversity that is not solely limited to the pEWs and velocities of different lines.

Zoom In Zoom Out Reset image size

One major open question about SNe Ia is whether there exist subtypes of SNe Ia, perhaps from different progenitor channels. To look for signs of such subtypes, we examine the distribution of the SNe Ia in the Twins Embedding as seen in, e.g., Figure 12. For this sample of SNe Ia, we do not see any evidence of large-scale bimodality in the recovered parameter space: the core of the parameter space is well populated and appears to be continuously filled. Note that we constructed the Twins Embedding using nonlinear dimensionality reduction, with the goal of preserving the spectral distances between “twin” SNe Ia from F15. The distances between points in the embedding should be meaningful because we are capturing the vast majority of the variance, as shown in Section 5.3, but the exact shape of the distribution of SNe Ia in the embedding will depend on the algorithm used and should not be overinterpreted.

The Twins Embedding includes all SNe Ia, including ones that are typically referred to as “peculiar.” We identified where the peculiar subclasses of 91bg-like (Filippenko et al. 1992a), 91T-like (Filippenko et al. 1992b), and 02cx-like SNe Ia (Li et al. 2003) are located in the Twins Embedding. We use the labeling of peculiar SNe Ia for the SNfactory sample from Lin et al. (2021). The results of this procedure can be seen in Figure 14.

Zoom In Zoom Out Reset image size

We find that all three of these kinds of peculiar SNe Ia lie at the edges of the Twins Embedding. They are easily separable from the rest of the sample using their coordinates in the Twins Embedding. There are seven 91T-like SNe Ia in our sample that are all tightly clustered at large values of ξ₁. The dereddened spectra of the 91T-like SNe Ia are all nearly identical, which explains why they lie in such a small region of the Twins Embedding. The three 91bg-like SNe Ia lie at very low values of ξ₂ with a significant spread in ξ₁. The two 02cx-like SNe Ia in our sample are separated from the rest of the sample in all of the Twins Embedding components.

One way to probe whether different SNe Ia are from the same physical processes is to test whether we can find a set of other SNe Ia whose spectra form a continuous sequence from the first spectrum to the second. If such a sequence can be found, then it suggests that the observed differences between those two SNe Ia are due to some continuous underlying physical parameter rather than being due to discrete processes. By selecting SNe Ia that lie along a path in the Twins Embedding, we can identify such a sequence between any two SNe Ia. In panel (A) of Figure 15, we show an example of this procedure by identifying a set of eight SNe Ia that are evenly spaced between the Twins Embedding coordinates of PTF 11kjn (a 91bg-like SN Ia) and SNF 20080522–000 (a 91T-like SN Ia). At each step in this sequence, the adjacent SNe Ia are “twins” according to the definition of F15, except for the initial pairing of PTF 11kjn to SN 2013bs whose “twinness percentile” of 29% is slightly above the threshold of 20% defined in F15 due to the low density of SNe Ia in this region of the parameter space. This sequence of SNe Ia is consistent at all wavelengths: the spectral features evolve continuously when traversing the sequence of spectra, although the changes in the spectral features are highly nonlinear. We are able to produce similar sequences between any two 91bg-like, 91T-like, or “normal” SNe Ia. This suggests that all of these SNe Ia may come from the same underlying process and that their diversity could be explained by variations in some continuous underlying parameters.

Zoom In Zoom Out Reset image size

On the other hand, we find that we are unable to identify continuous sequences of SNe Ia between 02cx-like SNe Ia and the rest of the sample. An example of an attempt at such a sequence is shown in panel (B) of Figure 15. SN 2005cc has large spectral differences relative to the next point in the sequence, SNF 20060512–001, and this pair falls into the 64th percentile of F15 twinness, which is worse than a random matching. The spectral differences that are seen do not follow the established sequence from the other SNe Ia, especially between 3800 and 5000 Å. Similar results are found for the other 02cx-like SN Ia in our sample, SN 2011ay. This could be suggestive of a different underlying process, but it is difficult to make definitive conclusions with our small sample of only two 02cx-like SNe Ia because we could simply be lacking intermediate SNe Ia in our sample.