# SUGAR: An improved empirical model of Type Ia Supernovae based on spectral features

P.F. Léget, E. Gangler, F. Mondon, G. Aldering, P. Antilogus, C. Aragon, S.

Bailey, C. Baltay, K. Barbary, S. Bongard, et al.

► To cite this version:

P.F. Léget, E. Gangler, F. Mondon, G. Aldering, P. Antilogus, et al. SUGAR: An improved empirical model of Type Ia Supernovae based on spectral features. Astronomy and Astrophysics - A&A, 2020, 636, pp.A46. 10.1051/0004-6361/201834954. hal-02371622

HAL Id: hal-02371622

https://hal.science/hal-02371622

Submitted on 5 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# SUGAR: An improved empirical model of Type Ia Supernovae based on spectral features

P.F. Léget, E. Gangler, F. Mondon, G. Aldering, P. Antilogus, C. Aragon, S. Bailey, C. Baltay, K. Barbary, S. Bongard, et al.

► To cite this version:

P.F. Léget, E. Gangler, F. Mondon, G. Aldering, P. Antilogus, et al. SUGAR: An improved empirical model of Type Ia Supernovae based on spectral features. Astronomy & Astrophysics - A&A, 2020, 636, pp.A46. 10.1051/0004-6361/201834954. hal-02371622

HAL Id: hal-02371622

https://hal.science/hal-02371622v1

Submitted on 5 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# SUGAR: An improved empirical model of Type Ia supernovae based on spectral features*

P.-F. Léget¹,²,³, E. Gangler¹, F. Mondon¹, G. Aldering⁴, P. Antilogus³, C. Aragon⁴, S. Bailey⁴, C. Baltay⁵,
K. Barbary⁴, S. Bongard³, K. Boone⁴,⁶, C. Buton⁷, N. Chotard⁷, Y. Copin⁷, S. Dixon⁴, P. Fagrelius⁴,⁶, U. Feindt⁸,
D. Fouchez⁹, B. Hayden⁴, W. Hillebrandt¹⁰, A. Kim⁴, M. Kowalski¹¹,¹², D. Kuesters¹¹, S. Lombardo¹¹, Q. Lin¹³,
J. Nordin¹¹, R. Pain³, E. Pecontal¹⁴, R. Pereira⁷, S. Perlmutter⁴,⁶, K. A. Ponder¹⁵, M. V. Pruzhinskaya¹,¹⁶,
D. Rabinowitz⁵, M. Rigault¹, K. Runge⁴, D. Rubin⁴,¹⁷, C. Saunders³, L.-P. Says¹, G. Smadja⁷, C. Sofiatti⁴,⁶,
N. Suzuki⁴,²⁰, S. Taubenberger¹⁰,¹⁸, C. Tao⁹,¹³, and R. C. Thomas¹⁹ (The Nearby Supernova Factory)

¹ Université Clermont Auvergne, CNRS/IN2P3, Laboratoire de Physique de Clermont, 63000 Clermont-Ferrand, France
e-mail: pierrefrancois.leget@gmail.com

² Kavli Institute for Particle Astrophysics and Cosmology, Department of Physics, Stanford University, Stanford, CA 94305, USA

³ LPNHE, CNRS/IN2P3, Sorbonne Université, Paris Diderot, Laboratoire de Physique Nucléaire et de Hautes Énergies, 75005 Paris, France

⁴ Physics Division, Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley CA 94720, USA

⁵ Department of Physics, Yale University, New Haven CT 06250-8121, USA

⁶ Department of Physics, University of California Berkeley, 366 LeConte Hall MC 7300, Berkeley, CA 94720-7300, USA

⁷ Université de Lyon, 69622, Lyon, France; Université de Lyon 1, Villeurbanne; CNRS/IN2P3, Institut de Physique des Deux Infinis, Lyon, France

⁸ The Oskar Klein Centre, Department of Physics, AlbaNova, Stockholm University, 106 91 Stockholm, Sweden

⁹ Aix Marseille Université, CNRS/IN2P3, CPPM, UMR 7346, 13288 Marseille, France

¹⁰ Max-Planck Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany

¹¹ Institut fur Physik, Humboldt-Universitat zu Berlin, Newtonstr. 15 12489, Berlin

¹² Deutsches Elektronen-Synchrotron, 15735 Zeuthen, Germany

¹³ Tsinghua Center for Astrophysics, Tsinghua University, Beijing 100084, PR China

¹⁴ Centre de Recherche Astronomique de Lyon, Université Lyon 1, 9 avenue Charles André, 69561 Saint Genis Laval, France

¹⁵ Berkeley Center for Cosmological Physics, University of California Berkeley, 341 Campbell Hall, Berkeley, CA 94720, USA

¹⁶ Lomonosov Moscow State University, Sternberg Astronomical Institute, Universitetsky pr. 13, Moscow 119234, Russia

¹⁷ Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

¹⁸ European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching, Germany

¹⁹ Computational Cosmology Center, Computational Research Division, Lawrence Berkeley National Laboratory, 1 Cyclotron Road MS 50B-4206, Berkeley, CA 94720, USA

²⁰ Kavli Institute for the Physics and Mathematics of the Universe, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8583, Japan

Received 21 December 2018 / Accepted 16 September 2019

## ABSTRACT

**Context.** Type Ia supernovae (SNe Ia) are widely used to measure the expansion of the Universe. Improving distance measurements of SNe Ia is one technique to better constrain the acceleration of expansion and determine its physical nature.

**Aims.** This document develops a new SNe Ia spectral energy distribution (SED) model, called the SUPernova Generator And Reconstructer (SUGAR), which improves the spectral description of SNe Ia, and consequently could improve the distance measurements.

**Methods.** This model was constructed from SNe Ia spectral properties and spectrophotometric data from the Nearby Supernova Factory collaboration. In a first step, a principal component analysis-like method was used on spectral features measured at maximum light, which allowed us to extract the intrinsic properties of SNe Ia. Next, the intrinsic properties were used to extract the average extinction curve. Third, an interpolation using Gaussian processes facilitated using data taken at different epochs during the lifetime of an SN Ia and then projecting the data on a fixed time grid. Finally, the three steps were combined to build the SED model as a function of time and wavelength. This is the SUGAR model.

**Results.** The main advancement in SUGAR is the addition of two additional parameters to characterize SNe Ia variability. The first is tied to the properties of SNe Ia ejecta velocity and the second correlates with their calcium lines. The addition of these parameters, as well as the high quality of the Nearby Supernova Factory data, makes SUGAR an accurate and efficient model for describing the spectra of normal SNe Ia as they brighten and fade.

**Conclusions.** The performance of this model makes it an excellent SED model for experiments like the Zwicky Transient Facility, the Large Synoptic Survey Telescope, or the Wide Field Infrared Survey Telescope.

**Key words.** supernovae: general – cosmology: observations

<footnote>* The data are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/636/A46</footnote>
# SUGAR: An improved empirical model of Type Ia supernovae based on spectral features*

P.-F. Léget¹,²,³, E. Gangler¹, F. Mondon¹, G. Aldering⁴, P. Antilogus³, C. Aragon⁴, S. Bailey⁴, C. Baltay⁵,
K. Barbary⁴, S. Bongard³, K. Boone⁴,⁶, C. Buton⁷, N. Chotard⁷, Y. Copin⁷, S. Dixon⁴, P. Fagrelius⁴,⁶, U. Feindt⁸,
D. Fouchez⁹, B. Hayden⁴, W. Hillebrandt¹⁰, A. Kim⁴, M. Kowalski¹¹,¹², D. Kuesters¹¹, S. Lombardo¹¹, Q. Lin¹³,
J. Nordin¹¹, R. Pain³, E. Pecontal¹⁴, R. Pereira⁷, S. Perlmutter⁴,⁶, K. A. Ponder¹⁵, M. V. Pruzhinskaya¹,¹⁶,
D. Rabinowitz⁵, M. Rigault¹, K. Runge⁴, D. Rubin⁴,¹⁷, C. Saunders³, L.-P. Says¹, G. Smadja⁷, C. Sofiatti⁴,⁶,
N. Suzuki⁴,²⁰, S. Taubenberger¹⁰,¹⁸, C. Tao⁹,¹³, and R. C. Thomas¹⁹ (The Nearby Supernova Factory)

¹ Université Clermont Auvergne, CNRS/IN2P3, Laboratoire de Physique de Clermont, 63000 Clermont-Ferrand, France
e-mail: pierrefrancois.leget@gmail.com

² Kavli Institute for Particle Astrophysics and Cosmology, Department of Physics, Stanford University, Stanford, CA 94305, USA

³ LPNHE, CNRS/IN2P3, Sorbonne Université, Paris Diderot, Laboratoire de Physique Nucléaire et de Hautes Énergies, 75005 Paris, France

⁴ Physics Division, Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley CA 94720, USA

⁵ Department of Physics, Yale University, New Haven CT 06250-8121, USA

⁶ Department of Physics, University of California Berkeley, 366 LeConte Hall MC 7300, Berkeley, CA 94720-7300, USA

⁷ Université de Lyon, 69622, Lyon, France; Université de Lyon 1, Villeurbanne; CNRS/IN2P3, Institut de Physique des Deux Infinis, Lyon, France

⁸ The Oskar Klein Centre, Department of Physics, AlbaNova, Stockholm University, 106 91 Stockholm, Sweden

⁹ Aix Marseille Université, CNRS/IN2P3, CPPM, UMR 7346, 13288 Marseille, France

¹⁰ Max-Planck Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85748 Garching, Germany

¹¹ Institut fur Physik, Humboldt-Universitat zu Berlin, Newtonstr. 15 12489, Berlin

¹² Deutsches Elektronen-Synchrotron, 15735 Zeuthen, Germany

¹³ Tsinghua Center for Astrophysics, Tsinghua University, Beijing 100084, PR China

¹⁴ Centre de Recherche Astronomique de Lyon, Université Lyon 1, 9 avenue Charles André, 69561 Saint Genis Laval, France

¹⁵ Berkeley Center for Cosmological Physics, University of California Berkeley, 341 Campbell Hall, Berkeley, CA 94720, USA

¹⁶ Lomonosov Moscow State University, Sternberg Astronomical Institute, Universitetsky pr. 13, Moscow 119234, Russia

¹⁷ Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

¹⁸ European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching, Germany

¹⁹ Computational Cosmology Center, Computational Research Division, Lawrence Berkeley National Laboratory, 1 Cyclotron Road MS 50B-4206, Berkeley, CA 94720, USA

²⁰ Kavli Institute for the Physics and Mathematics of the Universe, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8583, Japan

Received 21 December 2018 / Accepted 16 September 2019

## ABSTRACT

**Context.** Type Ia supernovae (SNe Ia) are widely used to measure the expansion of the Universe. Improving distance measurements of SNe Ia is one technique to better constrain the acceleration of expansion and determine its physical nature.

**Aims.** This document develops a new SNe Ia spectral energy distribution (SED) model, called the SUPernova Generator And Reconstructor (SUGAR), which improves the spectral description of SNe Ia, and consequently could improve the distance measurements.

**Methods.** This model was constructed from SNe Ia spectral properties and spectrophotometric data from the Nearby Supernova Factory collaboration. In a first step, a principal component analysis-like method was used on spectral features measured at maximum light, which allowed us to extract the intrinsic properties of SNe Ia. Next, the intrinsic properties were used to extract the average extinction curve. Third, an interpolation using Gaussian processes facilitated using data taken at different epochs during the lifetime of an SN Ia and then projecting the data on a fixed time grid. Finally, the three steps were combined to build the SED model as a function of time and wavelength. This is the SUGAR model.

**Results.** The main advancement in SUGAR is the addition of two additional parameters to characterize SNe Ia variability. The first is tied to the properties of SNe Ia ejecta velocity and the second correlates with their calcium lines. The addition of these parameters, as well as the high quality of the Nearby Supernova Factory data, makes SUGAR an accurate and efficient model for describing the spectra of normal SNe Ia as they brighten and fade.

**Conclusions.** The performance of this model makes it an excellent SED model for experiments like the Zwicky Transient Facility, the Large Synoptic Survey Telescope, or the Wide Field Infrared Survey Telescope.

**Key words.** supernovae: general – cosmology: observations

<footnote>* The data are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/636/A46</footnote>
# 1. Introduction

Type Ia supernovae (SNe Ia) are excellent cosmological probes: they are very luminous objects that are visible up to a redshift of $z \sim 2$ (Guillochon et al. 2017), and their luminosity dispersion is naturally low and can be further reduced by an appropriate standardization process. As precise distance indicators, comparing their luminosity and their redshift allowed Perlmutter & Aldering (1998), Perlmutter et al. (1999), Riess et al. (1998), and Schmidt et al. (1998) to demonstrate that the expansion of the Universe is accelerating; a feature that gave rise to the dark energy paradigm. This result, obtained with a small sample of SNe Ia, has since been repeatedly confirmed with larger samples (Astier et al. 2006; Guy et al. 2010; Suzuki et al. 2012; Rest et al. 2014; Betoule et al. 2014; Scolnic et al. 2018), and in combination with other cosmological probes like the cosmological microwave background (CMB; Planck Collaboration XIII 2016), baryon acoustic oscillations (BAO; Delubac et al. 2015), and cosmic shear (Troxel et al. 2018) led to the so-called concordance model, the flat-$\Lambda$ cold dark matter model.

Improving our description of SNe Ia as a cosmological probe is needed, not only to discriminate among alternate dark energy models (Copeland et al. 2006), but also in order to address existing tensions, such as the $3.4\sigma$ difference between the values of the Hubble constant $H_0$ from the local measurement by Riess et al. (2016) and the cosmological fit from the Planck Collaboration XIII (2016). The current uncertainty budget due to limited SN Ia statistics will be greatly improved by current surveys like the Zwicky Transient Facility (ZTF; Bellm 2014) or next generation surveys like the Large Synoptic Survey Telescope (LSST) or the Wide Field Infrared Survey Telescope (WFIRST) (LSST Dark Energy Science Collaboration 2012; Spergel et al. 2015). Since, systematic and statistical uncertainties are already of the same order of magnitude (Betoule et al. 2014; Scolnic et al. 2018), a better understanding of systematics will be needed to improve the accuracy of SNe Ia as a cosmology probe.

Part of the systematic uncertainty comes from the standardization process; the observed flux of the SN Ia has to be corrected for variations observed from one object to another. Two main contributions to this variation have been observed. The first one, seen by Rust (1974) and Pskovskii (1977, 1984), is the correlation between the peak luminosity of a SNe Ia and the light curve decrease time, the so-called brighter-slower effect. Various parametrizations have been proposed for this effect, the most commonly used being $\Delta m_{15}$ (Phillips 1993), stretch (Perlmutter et al. 1997), or $X_1$ (Guy et al. 2007). The second contribution, observed first by Hamuy et al. (1995) and Riess et al. (1996), is that the peak luminosity depends on color, the so-called brighter-bluer effect. This effect can be explained by the presence of dust in varying quantities along the line of sight that is possibly combined with an intrinsic color-brightness correlation after the stretch effect has been taken into account. Most standardization techniques in photometry thus rely on a stretch and color relation known as the Tripp (1998) relation. The SALT2 model (Guy et al. 2007) is one such standardization method, and has over the years become the reference in cosmological analysis. Despite many attempts, no consensus has yet emerged on how to go beyond a two-component model to describe SNe Ia light curves. However, the observation of a step in the standardized luminosity with respect to the host mass by Kelly et al. (2010) and Sullivan et al. (2010) has led to the inclusion of a corrective term in the subsequent cosmological analysis.

This corrective term linked to the environment hints that more fundamental properties of SNe Ia physics are not captured by the stretch-color standardization scheme, and that there is room for improvement.

However, the most obvious indication that the parameterization used is insufficient comes from the observation of a residual dispersion of SNe Ia luminosity around the Hubble diagram after standardization. In order to obtain statistically coherent results, Perlmutter et al. (1999) and Astier et al. (2006) introduced an intrinsic dispersion in luminosity as an additional uncertainty in the fit. Betoule et al. (2014) estimate the dispersion value at 0.11 mag using a SALT2 standardization, for an observed total dispersion of 0.16 mag. Since the intrinsic dispersion is due to unmodeled SN Ia variations that may depend on the redshift, this may result in a significant error on the extraction of cosmological parameters (Rigault et al. 2018). The precision of SNe Ia as cosmological probes therefore depends on this intrinsic dispersion in luminosity. However, the intrinsic dispersion strongly depends on the assumptions about measurement uncertainties. A way to characterize the overall accuracy obtained by a given standardization method is to use the weighted root mean square (wRMS) metric on the standardized magnitude (e.g., Blondin et al. 2011). Thus, in order to reduce the effects of unmodeled variations of SNe Ia, an improved standardization procedure must also involve reducing the wRMS.

Most efforts to improve the standardization of SNe Ia involve the search for a new parameter correlated with the intrinsic luminosity. After the discovery of the mass-step, many efforts focused on describing the environmental effects of the SNe Ia (Kelly et al. 2010; Sullivan et al. 2010; Rigault et al. 2013, 2015, 2018; Roman et al. 2018). A complementary approach consists of directly looking for this new parameter from the analysis of the light curves or observed spectra. For example, Mandel et al. (2017) proposed that the color of SNe Ia is a mixture of intrinsic color and extinction by dust, and proposed a Bayesian model to separate these two components. For their part, Chotard et al. (2011) showed that once the Ca II H&K variability is taken into account, the color law is compatible with a Cardelli et al. (1989) extinction curve. Moreover, Mandel et al. (2014) highlighted the dependence of the extinction curve on the minima of P Cygni profiles. Thus the use of additional variables related to spectral indicators seems to be a promising avenue to describe the variability of SNe Ia.

The use of spectral indicators to standardize SNe Ia has a long history (Nugent et al. 1995; Arsenijevic et al. 2008; Bailey et al. 2009; Wang et al. 2009; Foley & Kasen 2011; Chotard et al. 2011). As example, using only the flux ratio $R_{642/443}$, Bailey et al. (2009) is able to get a significatively better standardization in B-band in comparison to the classical Tripp (1998) relation. The method of using spectral information is also suggested by the analysis of Fakhouri et al. (2015), which shows that the best method to measure distance with SNe Ia is with spectroscopic twins SNe Ia. Another recent study by Nordin et al. (2018) has shown that the use of spectral information from the UV part of the spectra improves distance measurement compared to the Tripp (1998) relation. However, those methods do not currently lead to a spectral energy distribution (SED) model, which is necessary for cosmological analyses on purely photometric data, such as LSST (LSST Dark Energy Science Collaboration 2012). Thus we propose here a full SED model which will be based on spectral indicators generalizing the procedure originally developed in
# 1. Introduction

Type Ia supernovae (SNe Ia) are excellent cosmological probes: they are very luminous objects that are visible up to a redshift of $z \sim 2$ (Guillochon et al. 2017), and their luminosity dispersion is naturally low and can be further reduced by an appropriate standardization process. As precise distance indicators, comparing their luminosity and their redshift allowed Perlmutter & Aldering (1998), Perlmutter et al. (1999), Riess et al. (1998), and Schmidt et al. (1998) to demonstrate that the expansion of the Universe is accelerating; a feature that gave rise to the dark energy paradigm. This result, obtained with a small sample of SNe Ia, has since been repeatedly confirmed with larger samples (Astier et al. 2006; Guy et al. 2010; Suzuki et al. 2012; Rest et al. 2014; Betoule et al. 2014; Scolnic et al. 2018), and in combination with other cosmological probes like the cosmological microwave background (CMB; Planck Collaboration XIII 2016), baryon acoustic oscillations (BAO; Delubac et al. 2015), and cosmic shear (Troxel et al. 2018) led to the so-called concordance model, the flat-$\Lambda$ cold dark matter model.

Improving our description of SNe Ia as a cosmological probe is needed, not only to discriminate among alternate dark energy models (Copeland et al. 2006), but also in order to address existing tensions, such as the $3.4\sigma$ difference between the values of the Hubble constant $H_0$ from the local measurement by Riess et al. (2016) and the cosmological fit from the Planck Collaboration XIII (2016). The current uncertainty budget due to limited SN Ia statistics will be greatly improved by current surveys like the Zwicky Transient Facility (ZTF; Bellm 2014) or next generation surveys like the Large Synoptic Survey Telescope (LSST) or the Wide Field Infrared Survey Telescope (WFIRST) (LSST Dark Energy Science Collaboration 2012; Spergel et al. 2015). Since, systematic and statistical uncertainties are already of the same order of magnitude (Betoule et al. 2014; Scolnic et al. 2018), a better understanding of systematics will be needed to improve the accuracy of SNe Ia as a cosmology probe.

Part of the systematic uncertainty comes from the standardization process; the observed flux of the SN Ia has to be corrected for variations observed from one object to another. Two main contributions to this variation have been observed. The first one, seen by Rust (1974) and Pskovskii (1977, 1984), is the correlation between the peak luminosity of a SNe Ia and the light curve decrease time, the so-called brighter-slower effect. Various parametrizations have been proposed for this effect, the most commonly used being $\Delta m_{15}$ (Phillips 1993), stretch (Perlmutter et al. 1997), or $X_1$ (Guy et al. 2007). The second contribution, observed first by Hamuy et al. (1995) and Riess et al. (1996), is that the peak luminosity depends on color, the so-called brighter-bluer effect. This effect can be explained by the presence of dust in varying quantities along the line of sight that is possibly combined with an intrinsic color-brightness correlation after the stretch effect has been taken into account. Most standardization techniques in photometry thus rely on a stretch and color relation known as the Tripp (1998) relation. The SALT2 model (Guy et al. 2007) is one such standardization method, and has over the years become the reference in cosmological analysis. Despite many attempts, no consensus has yet emerged on how to go beyond a two-component model to describe SNe Ia light curves. However, the observation of a step in the standardized luminosity with respect to the host mass by Kelly et al. (2010) and Sullivan et al. (2010) has led to the inclusion of a corrective term in the subsequent cosmological analysis.

This corrective term linked to the environment hints that more fundamental properties of SNe Ia physics are not captured by the stretch-color standardization scheme, and that there is room for improvement.

However, the most obvious indication that the parameterization used is insufficient comes from the observation of a residual dispersion of SNe Ia luminosity around the Hubble diagram after standardization. In order to obtain statistically coherent results, Perlmutter et al. (1999) and Astier et al. (2006) introduced an intrinsic dispersion in luminosity as an additional uncertainty in the fit. Betoule et al. (2014) estimate the dispersion value at 0.11 mag using a SALT2 standardization, for an observed total dispersion of 0.16 mag. Since the intrinsic dispersion is due to unmodeled SN Ia variations that may depend on the redshift, this may result in a significant error on the extraction of cosmological parameters (Rigault et al. 2018). The precision of SNe Ia as cosmological probes therefore depends on this intrinsic dispersion in luminosity. However, the intrinsic dispersion strongly depends on the assumptions about measurement uncertainties. A way to characterize the overall accuracy obtained by a given standardization method is to use the weighted root mean square (wRMS) metric on the standardized magnitude (e.g., Blondin et al. 2011). Thus, in order to reduce the effects of unmodeled variations of SNe Ia, an improved standardization procedure must also involve reducing the wRMS.

Most efforts to improve the standardization of SNe Ia involve the search for a new parameter correlated with the intrinsic luminosity. After the discovery of the mass-step, many efforts focused on describing the environmental effects of the SNe Ia (Kelly et al. 2010; Sullivan et al. 2010; Rigault et al. 2013, 2015, 2018; Roman et al. 2018). A complementary approach consists of directly looking for this new parameter from the analysis of the light curves or observed spectra. For example, Mandel et al. (2017) proposed that the color of SNe Ia is a mixture of intrinsic color and extinction by dust, and proposed a Bayesian model to separate these two components. For their part, Chotard et al. (2011) showed that once the Ca II H&K variability is taken into account, the color law is compatible with a Cardelli et al. (1989) extinction curve. Moreover, Mandel et al. (2014) highlighted the dependence of the extinction curve on the minima of P Cygni profiles. Thus the use of additional variables related to spectral indicators seems to be a promising avenue to describe the variability of SNe Ia.

The use of spectral indicators to standardize SNe Ia has a long history (Nugent et al. 1995; Arsenijevic et al. 2008; Bailey et al. 2009; Wang et al. 2009; Foley & Kasen 2011; Chotard et al. 2011). As example, using only the flux ratio $R_{642/443}$, Bailey et al. (2009) is able to get a significatively better standardization in B-band in comparison to the classical Tripp (1998) relation. The method of using spectral information is also suggested by the analysis of Fakhouri et al. (2015), which shows that the best method to measure distance with SNe Ia is with spectroscopic twins SNe Ia. Another recent study by Nordin et al. (2018) has shown that the use of spectral information from the UV part of the spectra improves distance measurement compared to the Tripp (1998) relation. However, those methods do not currently lead to a spectral energy distribution (SED) model, which is necessary for cosmological analyses on purely photometric data, such as LSST (LSST Dark Energy Science Collaboration 2012). Thus we propose here a full SED model which will be based on spectral indicators generalizing the procedure originally developed in
Chotard et al. (2011). This method uses spectral features, which allows for the addition of more than one intrinsic parameter and they offer a possible way to separate intrinsic properties from extrinsic properties. From it, we build a full SED model which may be used for purely photometric surveys.

The aim of this study is to revisit the parametrization of SNe Ia SED in light of the SNfactory spectrophotometric dataset (Aldering et al. 2002), and to seek new sources of variability by statistical analysis. Our model is trained using spectral indicators, derived around maximum light in *B*-band, as features to describe the model. They provide both a reduced dimensionality description of spectra, and a description which is linked to the physics of the explosion. In addition, we select indicators insensitive to reddening in order to decouple the characterization of the reddening from effects purely linked to the intrinsic part of the explosion. This new SNe Ia SED model is named the Supernova Generator And Reconstructor (SUGAR) model.

Another approach to develop a new SED model was undertaken in parallel by Saunders et al. (2018), using the same dataset and based on a “SALT2 like” framework. Saunders et al. (2018) generalized a strategy that was originally proposed by Guy et al. (2007) and which was first performed on broad band photometry from the SNfactory dataset in Kim et al. (2013). In brief, Saunders et al. (2018) did a principal component analysis (PCA)-like study on interpolated spectral time series in order to go beyond the classical Tripp (1998) relation. Saunders et al. (2018) is able to significatively improve the SED description with respect to the SALT2 model. The SNEMO model developed in Saunders et al. (2018) differs from the SUGAR model in that SNEMO attempts to find principal components purely based on spectral time series variability in the relative luminosities at each wavelength and as a function of time, whereas SUGAR uses spectral features to try to find a compact description of such spectral time series. However, both models share some technical details in the model training and were developed in common.

This paper is organized as follows: Sect. 2 presents how the spectrophotometric time series of the SNfactory were obtained, Sect. 3 focuses on the intermediary data used in building the SUGAR model (spectral indicators at maximum light, dimensionality reduction through factor analysis, derivation of extinction parameters and time interpolation). Section 4 describes the SUGAR model: its formalism, training, and components. Section 5 presents the performance of the model for fitting spectral time series using SUGAR model. These results are compared to the performance achieved by SALT2 for the same data. Finally, Sect. 6 discusses adding additional components to the SUGAR model and some technical choices that were made in the training of the SUGAR model. The appendices describe details of the mathematical implementation of the SUGAR model. The SUGAR template is available online<fnref n="1" />, while the data used to train SUGAR are available online<fnref n="2" />.

## 2. Spectrophotometric time series from SNfactory

This analysis is based on 171 SNe Ia obtained by the SNfactory collaboration beginning in 2004 (Aldering et al. 2002) with the SuperNova Integral Field Spectrograph (SNIFS, Lantz et al. 2004) installed on the University of Hawaii 2.2-m telescope (Mauna Kea). SNIFS is a fully integrated instrument optimized for semi-automated observations of point sources on a structured

background over an extended optical window at moderate spectral resolution. SNIFS has a fully-filled 6.4" × 6.4" spectroscopic field-of-view subdivided into a grid of 15 × 15 contiguous square spatial elements (spaxels). The dual-channel spectrograph simultaneously covers 3200–5200 Å (*B*-channel) and 5100–10 000 Å (*R*-channel) with 2.8 and 3.2 Å resolution, respectively. The data reduction of the *x*, *y*, *λ* data cubes is summarized by Aldering et al. (2006) and updated in Sect. 2.1 of Scalzo et al. (2010). The flux calibration is developed in Sect. 2.2 of Pereira et al. (2013) based on the atmospheric extinction derived in Buton et al. (2013). In addition, observations are obtained at the SN Ia location at least one year after the explosion to serve as a final reference to enable subtraction of the underlying host and the host subtraction, as described in Bongard et al. (2011). For every SN Ia followed, the SNfactory creates a spectrophotometric time series, typically composed of ~14 epochs, with the first spectrum taken on average three days before maximum light in *B*-band (Bailey et al. 2009; Chotard et al. 2011). The sample of 171 SNe Ia contains the objects with good final references, objects that passed quality cuts suggested by Guy et al. (2010), and is restricted to objects with at least one observation in a time window of ±2.5 days around maximum light in *B*-band. The size of this window is kept identical with respect to the study of Chotard et al. (2011) and is discussed in Léget (2016). After flux calibration, host-galaxy subtraction and correction for Milky Way extinction, the flux of the 171 spectra is integrated in synthetic top-hat filters defined in Pereira et al. (2013) and for reference a SALT2 fit is applied with the model from Betoule et al. (2014) in order to obtain the $X_1, C$, and $m_B$ parameters. The spectra of the 171 SNe Ia are transformed to the rest frame with a fiducial cosmology and an arbitrary Hubble constant. It is allowing to respect blinding of future cosmological analysis using this dataset. For the spectral analysis here, the spectra are rebinned at 1500 km s⁻<fnref n="1" /> between 3254 and 8649 Å (197 bins per spectra) for computational efficiency while still resolving spectral features and converted to the absolute AB magnitude system. The training is done on 113 SNe Ia observed before 2010 and the validation is done on 58 SNe Ia observed after 2010.

## 3. Derived data

Any empirical SN Ia modeling must solve three problems<fnref n="3" />: Choose what features to model, account for color, and deal with the data sampling. For the features modeling of SUGAR, we used spectral features at maximum light to describe the intrinsic part of the SED. This is described successively in Sects. 3.1 and 3.2. Section 3.1 describes how the spectral features are selected and measured. Section 3.2 explains how the spectral features are projected onto a new basis that allows us to work in an orthogonal basis for the SUGAR training. To estimate the average color curve of SNe Ia, we generalize the method of Chotard et al. (2011). This is described in Sect. 3.3. To deal with the data sampling, we project the observed spectra onto a common time grid so that the SED model can be calculated on the same time grid. This is done in Sect. 3.4 using the Gaussian process method. The following sections describe those intermediate steps that determine the full SUGAR SED model in Sect. 4. We use the spectral features derived and color curve parameters

<footnote n="1">¹ http://supernovae.in2p3.fr/sugar</footnote>

<footnote n="2">² https://snfactory.lbl.gov/sugar</footnote>

<footnote n="3">³ In general, empirical SNIa modeling must also handle spectroscopic and photometric data that are not observed on the same day or by the same instruments. However, this is not a problem with the SNfactory dataset.</footnote>
Chotard et al. (2011). This method uses spectral features, which allows for the addition of more than one intrinsic parameter and they offer a possible way to separate intrinsic properties from extrinsic properties. From it, we build a full SED model which may be used for purely photometric surveys.

The aim of this study is to revisit the parametrization of SNe Ia SED in light of the SNfactory spectrophotometric dataset (Aldering et al. 2002), and to seek new sources of variability by statistical analysis. Our model is trained using spectral indicators, derived around maximum light in *B*-band, as features to describe the model. They provide both a reduced dimensionality description of spectra, and a description which is linked to the physics of the explosion. In addition, we select indicators insensitive to reddening in order to decouple the characterization of the reddening from effects purely linked to the intrinsic part of the explosion. This new SNe Ia SED model is named the Supernova Generator And Reconstructor (SUGAR) model.

Another approach to develop a new SED model was undertaken in parallel by Saunders et al. (2018), using the same dataset and based on a “SALT2 like” framework. Saunders et al. (2018) generalized a strategy that was originally proposed by Guy et al. (2007) and which was first performed on broad band photometry from the SNfactory dataset in Kim et al. (2013). In brief, Saunders et al. (2018) did a principal component analysis (PCA)-like study on interpolated spectral time series in order to go beyond the classical Tripp (1998) relation. Saunders et al. (2018) is able to significatively improve the SED description with respect to the SALT2 model. The SNEMO model developed in Saunders et al. (2018) differs from the SUGAR model in that SNEMO attempts to find principal components purely based on spectral time series variability in the relative luminosities at each wavelength and as a function of time, whereas SUGAR uses spectral features to try to find a compact description of such spectral time series. However, both models share some technical details in the model training and were developed in common.

This paper is organized as follows: Sect. 2 presents how the spectrophotometric time series of the SNfactory were obtained, Sect. 3 focuses on the intermediary data used in building the SUGAR model (spectral indicators at maximum light, dimensionality reduction through factor analysis, derivation of extinction parameters and time interpolation). Section 4 describes the SUGAR model: its formalism, training, and components. Section 5 presents the performance of the model for fitting spectral time series using SUGAR model. These results are compared to the performance achieved by SALT2 for the same data. Finally, Sect. 6 discusses adding additional components to the SUGAR model and some technical choices that were made in the training of the SUGAR model. The appendices describe details of the mathematical implementation of the SUGAR model. The SUGAR template is available online<fnref n="1" />, while the data used to train SUGAR are available online<fnref n="2" />.

## 2. Spectrophotometric time series from SNfactory

This analysis is based on 171 SNe Ia obtained by the SNfactory collaboration beginning in 2004 (Aldering et al. 2002) with the SuperNova Integral Field Spectrograph (SNIFS, Lantz et al. 2004) installed on the University of Hawaii 2.2-m telescope (Mauna Kea). SNIFS is a fully integrated instrument optimized for semi-automated observations of point sources on a structured

background over an extended optical window at moderate spectral resolution. SNIFS has a fully-filled 6.4" × 6.4" spectroscopic field-of-view subdivided into a grid of 15 × 15 contiguous square spatial elements (spaxels). The dual-channel spectrograph simultaneously covers 3200–5200 Å (*B*-channel) and 5100–10 000 Å (*R*-channel) with 2.8 and 3.2 Å resolution, respectively. The data reduction of the *x*, *y*, *λ* data cubes is summarized by Aldering et al. (2006) and updated in Sect. 2.1 of Scalzo et al. (2010). The flux calibration is developed in Sect. 2.2 of Pereira et al. (2013) based on the atmospheric extinction derived in Buton et al. (2013). In addition, observations are obtained at the SN Ia location at least one year after the explosion to serve as a final reference to enable subtraction of the underlying host and the host subtraction, as described in Bongard et al. (2011). For every SN Ia followed, the SNfactory creates a spectrophotometric time series, typically composed of ~14 epochs, with the first spectrum taken on average three days before maximum light in *B*-band (Bailey et al. 2009; Chotard et al. 2011). The sample of 171 SNe Ia contains the objects with good final references, objects that passed quality cuts suggested by Guy et al. (2010), and is restricted to objects with at least one observation in a time window of ±2.5 days around maximum light in *B*-band. The size of this window is kept identical with respect to the study of Chotard et al. (2011) and is discussed in Léget (2016). After flux calibration, host-galaxy subtraction and correction for Milky Way extinction, the flux of the 171 spectra is integrated in synthetic top-hat filters defined in Pereira et al. (2013) and for reference a SALT2 fit is applied with the model from Betoule et al. (2014) in order to obtain the $X_1, C$, and $m_B$ parameters. The spectra of the 171 SNe Ia are transformed to the rest frame with a fiducial cosmology and an arbitrary Hubble constant. It is allowing to respect blinding of future cosmological analysis using this dataset. For the spectral analysis here, the spectra are rebinned at 1500 km s⁻<fnref n="1" /> between 3254 and 8649 Å (197 bins per spectra) for computational efficiency while still resolving spectral features and converted to the absolute AB magnitude system. The training is done on 113 SNe Ia observed before 2010 and the validation is done on 58 SNe Ia observed after 2010.

## 3. Derived data

Any empirical SN Ia modeling must solve three problems<fnref n="3" />: Choose what features to model, account for color, and deal with the data sampling. For the features modeling of SUGAR, we used spectral features at maximum light to describe the intrinsic part of the SED. This is described successively in Sects. 3.1 and 3.2. Section 3.1 describes how the spectral features are selected and measured. Section 3.2 explains how the spectral features are projected onto a new basis that allows us to work in an orthogonal basis for the SUGAR training. To estimate the average color curve of SNe Ia, we generalize the method of Chotard et al. (2011). This is described in Sect. 3.3. To deal with the data sampling, we project the observed spectra onto a common time grid so that the SED model can be calculated on the same time grid. This is done in Sect. 3.4 using the Gaussian process method. The following sections describe those intermediate steps that determine the full SUGAR SED model in Sect. 4. We use the spectral features derived and color curve parameters

<footnote n="1">¹ http://supernovae.in2p3.fr/sugar</footnote>

<footnote n="2">² https://snfactory.lbl.gov/sugar</footnote>

<footnote n="3">³ In general, empirical SNIa modeling must also handle spectroscopic and photometric data that are not observed on the same day or by the same instruments. However, this is not a problem with the SNfactory dataset.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Pseudo-equivalent widths of absorption lines at maximum light as well as minima of P-Cygni profiles at maximum of light which are used in our analysis. They are represented on the spectrum of SN-training-77 at a phase of 1 day after maximum brightness.</figcaption></figure>

derived at maximum light combined with the interpolated spectra to infer the full SED at all epochs.

### 3.1. Spectral indicators at maximum light

Spectral indicators are metrics of empirical features of the input spectrum such as equivalent widths or line velocities. They offer an efficient characterization of spectral variability by representing the underlying spectral complexity with a few key numbers, and play a role in nonlinear dimensionality reduction of the original data. These two characteristics, interpretability and simplicity, make them ideal for describing the intrinsic part of SNe Ia SED, and consequently derive the extrinsic part in the same spirit as Chotard et al. (2011). This decoupling restricts the set of spectral indicators to pseudo-equivalent widths and the wavelengths of P-Cygni profile minima. As the spectral indicators evolve with phase, we select the spectrum closest to maximum light in B-band, if it is within a time window of ±2.5 days. This window size was chosen to optimize the trade-off between the total number of SNe Ia in the sample and the potential loss in precision due to time evolution.

Near maximum light and within the spectral range covered by the SNfactory, it is possible to systematically obtain the 13 spectral indicators represented in Fig. 1. They are the nine pseudo-equivalent widths of Ca II H&K, Si II λ4131, Mg II, Fe λ4800, S II W, Si II λ5972, Si II λ6355, O I λ7773⁴, and Ca II IR features, as well as four minima of P Cygni profiles (Si II λ4131, S II λ5454, S II λ5640, and Si II λ6355). For the pseudo-equivalent widths, we rely only on well-defined troughs that are present in all SNe Ia in the sample, and consider line blends as a whole. Some of the possible minima were also discarded either because the corresponding lines form a complex mixture or because the corresponding trough is too shallow to accurately define a minimum for some SNe Ia. This is why the Si II λ5972 feature velocity is rejected. The spectral indicators and their uncertainties were automatically derived from the spectra at maximum light following the procedure of Chotard et al. (2011), which is described in detail in Appendix A.

## 3.2. Factor analysis on spectral features

### 3.2.1. Factor analysis model

The space defined by the 13 selected spectral indicators has too high dimensionality to efficiently train a model. Additionally, some of the spectral indicators are correlated, and therefore contain redundant information. Most of the model variation can be captured in a reduced number of dimensions. Principal component analysis (Pearson 1901) is one of the methods to implement such a dimensionality reduction. It consists of diagonalizing the covariance matrix of the sample and projecting the data into the resulting eigenvectors basis. In this new basis, the variables are uncorrelated, and an approximation of the input data is found by neglecting the dimensions corresponding to the smallest eigenvalues. This method has been employed in the case of SNe Ia by Guy et al. (2007), Kim et al. (2013), and Sasdelli et al. (2015). However, in the case considered here, some directions are dominated by noise, so their eigenvectors would align along the direction of measurement errors rather than the intrinsic sample variance. To solve this problem, we employ a variant of PCA, Factor Analysis (Spearman 1904, 1927), which has the advantage of taking into account the variance caused by measurement uncertainties. This technique decomposes the observables into two terms, one representing the explanatory factors and one representing the noise affecting each variable. This is expressed by the following relation (Ghahramani & Hinton 1997):

$$x_i = \Lambda q_i + u_i, \quad (1)$$

where $x_i$ is the spectral feature vector with $s$ components for the $i$th SN Ia. $q_i$ is the explanatory factor vector of dimension $m \le s$, and is linearly related to $x_i$ by the matrix $\Lambda$ of dimensions $m \times s$. $u_i$ is the noise with variance $\Psi$ of dimensions $s \times s$. In order to fix the normalization of $\Lambda$ and $q_i$, we assume that the $q_i$ are drawn from a centered normal distribution:

$$P(q_i) \sim N(0, I). \quad (2)$$

In the framework of the PCA, and up to a normalization, the matrix $\Lambda$ and the $q_i$ are respectively equivalent to the eigenvector matrix and the projections into the new basis. Factor analysis thus consists of determining the matrices $\Lambda$ and $\Psi$ that maximize the likelihood, under the assumption that the $x_i$ are distributed according to a normal distribution with both intrinsic scatter and measurement noise:

$$P(x_i) \sim N(0, \Lambda\Lambda^T + \Psi). \quad (3)$$

To estimate $\Lambda$, $\Psi$, and the explanatory factors $q_i$, Ghahramani & Hinton (1997) propose a solution based on an expectation-maximization algorithm where $q_i$ and $\Lambda$ are estimated iteratively. Here, unlike in conventional factor analysis, a reliable estimate of the spectral indicator measurement error is provided. For each SN Ia $i$, $\Psi_i$ is the known diagonal matrix that contains the squared errors of the spectral indicators, hence, it is not necessary to fit for a global $\Psi$. We adapted the aforementioned expectation-maximization algorithm to take this change into account. Additional details can be found in Appendix B.

Finally, a prescription is needed to normalize each individual variable. The amplitude of the spectral indicator variation, expressed in Å, is not a good indicator of their impact on the spectral shape. As an example, Branch et al. (2006) show that the weak Si II λ5972 line can play a significant role in subclassing SNe Ia. As a consequence, the input data are normalized to

<footnote n="4">⁴ We refer to this as O I λ7773 but note it overlaps with Mg II λ7812.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Pseudo-equivalent widths of absorption lines at maximum light as well as minima of P-Cygni profiles at maximum of light which are used in our analysis. They are represented on the spectrum of SN-training-77 at a phase of 1 day after maximum brightness.</figcaption></figure>

derived at maximum light combined with the interpolated spectra to infer the full SED at all epochs.

### 3.1. Spectral indicators at maximum light

Spectral indicators are metrics of empirical features of the input spectrum such as equivalent widths or line velocities. They offer an efficient characterization of spectral variability by representing the underlying spectral complexity with a few key numbers, and play a role in nonlinear dimensionality reduction of the original data. These two characteristics, interpretability and simplicity, make them ideal for describing the intrinsic part of SNe Ia SED, and consequently derive the extrinsic part in the same spirit as Chotard et al. (2011). This decoupling restricts the set of spectral indicators to pseudo-equivalent widths and the wavelengths of P-Cygni profile minima. As the spectral indicators evolve with phase, we select the spectrum closest to maximum light in B-band, if it is within a time window of ±2.5 days. This window size was chosen to optimize the trade-off between the total number of SNe Ia in the sample and the potential loss in precision due to time evolution.

Near maximum light and within the spectral range covered by the SNfactory, it is possible to systematically obtain the 13 spectral indicators represented in Fig. 1. They are the nine pseudo-equivalent widths of Ca II H&K, Si II λ4131, Mg II, Fe λ4800, S II W, Si II λ5972, Si II λ6355, O I λ7773$^4$, and Ca II IR features, as well as four minima of P Cygni profiles (Si II λ4131, S II λ5454, S II λ5640, and Si II λ6355). For the pseudo-equivalent widths, we rely only on well-defined troughs that are present in all SNe Ia in the sample, and consider line blends as a whole. Some of the possible minima were also discarded either because the corresponding lines form a complex mixture or because the corresponding trough is too shallow to accurately define a minimum for some SNe Ia. This is why the Si II λ5972 feature velocity is rejected. The spectral indicators and their uncertainties were automatically derived from the spectra at maximum light following the procedure of Chotard et al. (2011), which is described in detail in Appendix A.

## 3.2. Factor analysis on spectral features

### 3.2.1. Factor analysis model

The space defined by the 13 selected spectral indicators has too high dimensionality to efficiently train a model. Additionally, some of the spectral indicators are correlated, and therefore contain redundant information. Most of the model variation can be captured in a reduced number of dimensions. Principal component analysis (Pearson 1901) is one of the methods to implement such a dimensionality reduction. It consists of diagonalizing the covariance matrix of the sample and projecting the data into the resulting eigenvectors basis. In this new basis, the variables are uncorrelated, and an approximation of the input data is found by neglecting the dimensions corresponding to the smallest eigenvalues. This method has been employed in the case of SNe Ia by Guy et al. (2007), Kim et al. (2013), and Sasdelli et al. (2015). However, in the case considered here, some directions are dominated by noise, so their eigenvectors would align along the direction of measurement errors rather than the intrinsic sample variance. To solve this problem, we employ a variant of PCA, Factor Analysis (Spearman 1904, 1927), which has the advantage of taking into account the variance caused by measurement uncertainties. This technique decomposes the observables into two terms, one representing the explanatory factors and one representing the noise affecting each variable. This is expressed by the following relation (Ghahramani & Hinton 1997):

$$x_i = \Lambda q_i + u_i, \quad (1)$$

where $x_i$ is the spectral feature vector with $s$ components for the $i$th SN Ia. $q_i$ is the explanatory factor vector of dimension $m \le s$, and is linearly related to $x_i$ by the matrix $\Lambda$ of dimensions $m \times s$. $u_i$ is the noise with variance $\Psi$ of dimensions $s \times s$. In order to fix the normalization of $\Lambda$ and $q_i$, we assume that the $q_i$ are drawn from a centered normal distribution:

$$P(q_i) \sim N(0, I). \quad (2)$$

In the framework of the PCA, and up to a normalization, the matrix $\Lambda$ and the $q_i$ are respectively equivalent to the eigenvector matrix and the projections into the new basis. Factor analysis thus consists of determining the matrices $\Lambda$ and $\Psi$ that maximize the likelihood, under the assumption that the $x_i$ are distributed according to a normal distribution with both intrinsic scatter and measurement noise:

$$P(x_i) \sim N(0, \Lambda\Lambda^T + \Psi). \quad (3)$$

To estimate $\Lambda$, $\Psi$, and the explanatory factors $q_i$, Ghahramani & Hinton (1997) propose a solution based on an expectation-maximization algorithm where $q_i$ and $\Lambda$ are estimated iteratively. Here, unlike in conventional factor analysis, a reliable estimate of the spectral indicator measurement error is provided. For each SN Ia $i$, $\Psi_i$ is the known diagonal matrix that contains the squared errors of the spectral indicators, hence, it is not necessary to fit for a global $\Psi$. We adapted the aforementioned expectation-maximization algorithm to take this change into account. Additional details can be found in Appendix B.

Finally, a prescription is needed to normalize each individual variable. The amplitude of the spectral indicator variation, expressed in Å, is not a good indicator of their impact on the spectral shape. As an example, Branch et al. (2006) show that the weak Si II λ5972 line can play a significant role in subclassing SNe Ia. As a consequence, the input data are normalized to

<footnote n="4">⁴ We refer to this as O I λ7773 but note it overlaps with Mg II λ7812.</footnote>
unit variance prior to the factor analysis, so that no spectral indicator is favored a priori. Within this framework, the eigenvalues of $\Lambda\Lambda^T$ represent the variance explained by each factor, and $s - \text{Tr}(\Lambda\Lambda^T)$ is the variance coming from the noise, where here $s = 13$ is the number of spectral indicators and Tr is the trace operator.

### 3.2.2. Outlier rejection

Outliers affect the sample variance of any population. As we care most about a correct description of the bulk of SNe Ia, it is desirable to identify and remove these outliers. Amongst them, SNe Ia of type SN1991T (Filippenko et al. 1992a) or SN1991bg (Filippenko et al. 1992b) are known to have different spectral and photometric behavior from other SNe Ia. However, basing an outlier rejection on this empirical identification has two issues. The first is that within the spectral indicator space, some of those subtypes may not appear as distinct subclasses, and do not offer objective grounds for rejection. The second issue is that the attribution of a given SN Ia to one of these subclasses by SNID (Blondin & Tonry 2007) may provide inconsistent results depending on the epoch considered for the identification. In order to apply a self-contained criterion for defining an outlier, we thus resorted to a $\chi^2$-based definition for identifying outliers: this quantity can be interpreted as the squared distance to the center of the distribution normalized by the natural dispersion. For SN Ia $i$, it is given by

$$ \chi_i^2 = \mathbf{x}_i^T (\Lambda\Lambda^T + \mathbf{\Psi}_i) \mathbf{x}_i. \quad (4) $$

In the case of a Gaussian distribution, $\chi_i^2$ should follow a $\chi^2$-distribution with 13 degrees of freedom. An iterative cut at $3\sigma$ on the value of $\chi_i^2$ rejects 8 SNe Ia from the sample. A visual inspection carried out on each of these 8 SNe Ia, as well as a sub-classification made using SNID, show that four of those objects exhibit shallow silicon features as defined by Branch et al. (2006), two of them have high-velocity silicon features as defined by Wang et al. (2009), one is within the broad line subcategory, and the remaining one has a failed estimation of one of the spectral indicators.

We studied the influence of the $3\sigma$ cut by carrying out the analysis successively with and without the cut. Even though this had no major effect on the direction of the vectors, the cut is applied to the rest of the analysis to avoid training the model on outliers.

### 3.2.3. Factor analysis results on spectral features

The expectation-maximization factor analysis (EM-FA) algorithm described in the previous sections is applied to the 13 spectral features measured at maximum light. The relative weight of the eigenvalues on the total variance is shown in Fig. 2, and the correlations between the eigenvectors, and the spectral indicators are presented in Fig. 3. The significance of the correlation is expressed in units of $\sigma$, meaning that the null hypothesis, that is 0 correlation is rejected at $\sigma$ confidence level. Figure 2 shows that two vectors dominate the variability in the space of spectral indicators. The first vector is linked with coherent variations of the pseudo-equivalent widths, except the widths of the Ca II H&K and S II W lines as seen in Fig. 3. The second vector is anticorrelated with these same two features, and describes a coherent variation of the velocities. Since the first vector is very similar to the pseudo-equivalent width of Si II $\lambda4131$, which is shown by Arsenijevic et al. (2008) to be highly correlated with

<figure><img src="image_3.png" /><figcaption>Fig. 2. Relative importance of the eigenvalues associated with the eigenvectors. They are ordered in decreasing order of variance, and the total amounts to 100% once the contribution of the noise (indicated by a red line) is taken into account. We can see that the first three vectors dominate the variability of this space.</figcaption></figure>

the stretch parameter, it is therefore natural to assume that this vector represents the main source of intrinsic variability already known for SNe Ia. The second vector is mostly driven by the line velocities, pEW(Ca II H&K), and pEW(S II W) lines, with smaller dependencies from the other pseudo-equivalent widths. Interestingly, while investigating the spectral diversity of SNe Ia beyond stretch and color, Chotard et al. (2011) demonstrated the role of pEW(Ca II H&K), while Wang et al. (2009) focused on the role of the speed of Si II $\lambda6355$. Our second vector unifies these two approaches. The interpretation of the next vectors is less straightforward. As their rank increases, they describe less of the sample variance, and thus exhibit lower correlations with the original data. The third vector is the last of the eigenvectors to show strong correlation with one of the spectral indicators, pEW Si II $\lambda5972$. The correlations with the SALT2 parameters are presented in Fig. 4. As expected from the correlations with the pseudo-equivalent widths, $q_1$ is strongly correlated with $X_1$. Also, $q_3$ exhibits a significant correlation with $X_1$, which may be an indication that the stretch is driven by two different parameters. Remaining factors are only weakly correlated with $X_1$. The correlations between SALT2 color and any of the factors are close to zero. The SALT2 absolute magnitude cosmology residual, $\Delta\mu_B$, are mainly correlated with $q_1$ and $q_3$, presumably due to the correlation of both with $X_1$.

Within this framework, the eigenvalues of $\Lambda\Lambda^T$ represent the variance explained by each factor, and $s - \text{Tr}(\Lambda\Lambda^T)$ is the variance coming from the noise, where here $s = 13$ is the number of spectral indicators and Tr is the trace operator. Noise, in this situation, represents 18% of the variability observed, so only the first two vectors clearly outweigh the noise and it is legitimate to ask whether or not the other vectors are sensitive to statistical fluctuations. For the rest of the analysis, we keep the first three factors for the final training and discuss this choice in the Sect. 6.2, by looking at the impact on the final model.

### 3.3. Extinction curve estimation

In this section, we estimate the average extinction curve by generalizing the procedure described in Chotard et al. (2011). For this,
unit variance prior to the factor analysis, so that no spectral indicator is favored a priori. Within this framework, the eigenvalues of $\Lambda\Lambda^T$ represent the variance explained by each factor, and $s - \text{Tr}(\Lambda\Lambda^T)$ is the variance coming from the noise, where here $s = 13$ is the number of spectral indicators and Tr is the trace operator.

### 3.2.2. Outlier rejection

Outliers affect the sample variance of any population. As we care most about a correct description of the bulk of SNe Ia, it is desirable to identify and remove these outliers. Amongst them, SNe Ia of type SN1991T (Filippenko et al. 1992a) or SN1991bg (Filippenko et al. 1992b) are known to have different spectral and photometric behavior from other SNe Ia. However, basing an outlier rejection on this empirical identification has two issues. The first is that within the spectral indicator space, some of those subtypes may not appear as distinct subclasses, and do not offer objective grounds for rejection. The second issue is that the attribution of a given SN Ia to one of these subclasses by SNID (Blondin & Tonry 2007) may provide inconsistent results depending on the epoch considered for the identification. In order to apply a self-contained criterion for defining an outlier, we thus resorted to a $\chi^2$-based definition for identifying outliers: this quantity can be interpreted as the squared distance to the center of the distribution normalized by the natural dispersion. For SN Ia $i$, it is given by

$$ \chi_i^2 = \mathbf{x}_i^T (\Lambda\Lambda^T + \mathbf{\Psi}_i) \mathbf{x}_i. \quad (4) $$

In the case of a Gaussian distribution, $\chi_i^2$ should follow a $\chi^2$-distribution with 13 degrees of freedom. An iterative cut at $3\sigma$ on the value of $\chi_i^2$ rejects 8 SNe Ia from the sample. A visual inspection carried out on each of these 8 SNe Ia, as well as a sub-classification made using SNID, show that four of those objects exhibit shallow silicon features as defined by Branch et al. (2006), two of them have high-velocity silicon features as defined by Wang et al. (2009), one is within the broad line subcategory, and the remaining one has a failed estimation of one of the spectral indicators.

We studied the influence of the $3\sigma$ cut by carrying out the analysis successively with and without the cut. Even though this had no major effect on the direction of the vectors, the cut is applied to the rest of the analysis to avoid training the model on outliers.

### 3.2.3. Factor analysis results on spectral features

The expectation-maximization factor analysis (EM-FA) algorithm described in the previous sections is applied to the 13 spectral features measured at maximum light. The relative weight of the eigenvalues on the total variance is shown in Fig. 2, and the correlations between the eigenvectors, and the spectral indicators are presented in Fig. 3. The significance of the correlation is expressed in units of $\sigma$, meaning that the null hypothesis, that is 0 correlation is rejected at $\sigma$ confidence level. Figure 2 shows that two vectors dominate the variability in the space of spectral indicators. The first vector is linked with coherent variations of the pseudo-equivalent widths, except the widths of the Ca II H&K and S II W lines as seen in Fig. 3. The second vector is anticorrelated with these same two features, and describes a coherent variation of the velocities. Since the first vector is very similar to the pseudo-equivalent width of Si II $\lambda4131$, which is shown by Arsenijevic et al. (2008) to be highly correlated with

<figure><img src="image_3.png" /><figcaption>Fig. 2. Relative importance of the eigenvalues associated with the eigenvectors. They are ordered in decreasing order of variance, and the total amounts to 100% once the contribution of the noise (indicated by a red line) is taken into account. We can see that the first three vectors dominate the variability of this space.</figcaption></figure>

the stretch parameter, it is therefore natural to assume that this vector represents the main source of intrinsic variability already known for SNe Ia. The second vector is mostly driven by the line velocities, pEW(Ca II H&K), and pEW(S II W) lines, with smaller dependencies from the other pseudo-equivalent widths. Interestingly, while investigating the spectral diversity of SNe Ia beyond stretch and color, Chotard et al. (2011) demonstrated the role of pEW(Ca II H&K), while Wang et al. (2009) focused on the role of the speed of Si II $\lambda6355$. Our second vector unifies these two approaches. The interpretation of the next vectors is less straightforward. As their rank increases, they describe less of the sample variance, and thus exhibit lower correlations with the original data. The third vector is the last of the eigenvectors to show strong correlation with one of the spectral indicators, pEW Si II $\lambda5972$. The correlations with the SALT2 parameters are presented in Fig. 4. As expected from the correlations with the pseudo-equivalent widths, $q_1$ is strongly correlated with $X_1$. Also, $q_3$ exhibits a significant correlation with $X_1$, which may be an indication that the stretch is driven by two different parameters. Remaining factors are only weakly correlated with $X_1$. The correlations between SALT2 color and any of the factors are close to zero. The SALT2 absolute magnitude cosmology residual, $\Delta\mu_B$, are mainly correlated with $q_1$ and $q_3$, presumably due to the correlation of both with $X_1$.

Within this framework, the eigenvalues of $\Lambda\Lambda^T$ represent the variance explained by each factor, and $s - \text{Tr}(\Lambda\Lambda^T)$ is the variance coming from the noise, where here $s = 13$ is the number of spectral indicators and Tr is the trace operator. Noise, in this situation, represents 18% of the variability observed, so only the first two vectors clearly outweigh the noise and it is legitimate to ask whether or not the other vectors are sensitive to statistical fluctuations. For the rest of the analysis, we keep the first three factors for the final training and discuss this choice in the Sect. 6.2, by looking at the impact on the final model.

### 3.3. Extinction curve estimation

In this section, we estimate the average extinction curve by generalizing the procedure described in Chotard et al. (2011). For this,
<figure><img src="image_4.png" /><figcaption>Fig. 3. Pearson correlation coefficients between the spectral indicators and the first five factors. The eccentricities of the ellipses indicate the magnitude of the correlation, and the color the importance of the correlation in units of σ. The vectors one and two correspond respectively to a global correlation of the pseudo-equivalent widths and of the line velocities. Only the first three vectors display correlations with a significance higher than 5σ.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Pearson correlation coefficients between the SALT2 parameters and the first five vectors. The eccentricities of the ellipses indicate the magnitude of the correlation, and the color the importance of the correlation in units of σ.</figcaption></figure>

a model of the SED at maximum light is derived. This model allows us to deduce the average extinction curve and the procedure to estimate it: the results are described in the following.

### 3.3.1. Empirical description

The SED model at maximum light presented here is a generalization of the model presented in Chotard et al. (2011), in which the

authors removed intrinsic variability tied to spectral features and measured the remaining spectral variation, which was found to be consistent with an extinction curve like that for the Galaxy. The major difference is that the intrinsic description is here motivated by the factors derived in Sect. 3.2 rather than being described only by pEW(Si II λ4131) and pEW(Ca II H&K) as was done in Chotard et al. (2011). To model the SED we define for the SN Ia i the vector $\mathbf{x}_i \equiv \{h_{i,1}, h_{i,2}, h_{i,3}\}$, which is the true value of the measured factor $\mathbf{q}_i = \{q_{i,1}, q_{i,2}, q_{i,3}\}$. We propose that these are related to the intrinsic absolute magnitude by:

$$M_i(\lambda) = M_0(\lambda) + \sum_{j=1}^{j=3} h_{i,j} \alpha_j(\lambda) + A_{\lambda_0, i} \gamma(\lambda), \quad (5)$$

where $M_0(\lambda)$ is the average spectrum in absolute magnitude, $\alpha_j(\lambda)$ is the spectrum related to the factor $h_{i,j}$, where j is the factor index running from one to three (the choice of the number of components is discussed in Sect. 6.2). $A_{\lambda_0, i}$ is a free parameter that can be interpreted as the absorption due to extinction at a reference wavelength $\lambda_0$ (set here at the median wavelength, 5305 Å, without loss of generality). $\gamma(\lambda)$ is an arbitrary function that represents the effect of the extinction, and we have not set any prior on its shape. Once the function $\gamma(\lambda)$ is fixed, it is possible to fit the total-to-selective extinction ratio of the Cardelli et al. (1989) law, $R_V$, and the absorption in V-band, $A_V$ to it; this is described in the next sub-section. The model parameters here are $M_0(\lambda)$, the $\alpha_j(\lambda)$, the $\gamma(\lambda)$, the $h_{i,j}$, and the $A_{\lambda_0, i}$. It is understood that we are modeling the SED at specific wavelengths where we have measurements. For notational efficiency we rewrite Eq. (5) as:

$$\mathbf{M}_i = A \mathbf{h}_i, \quad \text{where} \quad (6)$$

$$\mathbf{A} = (\mathbf{M}_0, \gamma, \alpha_1, \alpha_2, \alpha_3), \quad (7)$$

$$\mathbf{h}_i^T = (1, A_{\lambda_0, i}, h_{i,1}, h_{i,2}, h_{i,3}). \quad (8)$$
<figure><img src="image_4.png" /><figcaption>Fig. 3. Pearson correlation coefficients between the spectral indicators and the first five factors. The eccentricities of the ellipses indicate the magnitude of the correlation, and the color the importance of the correlation in units of σ. The vectors one and two correspond respectively to a global correlation of the pseudo-equivalent widths and of the line velocities. Only the first three vectors display correlations with a significance higher than 5σ.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Pearson correlation coefficients between the SALT2 parameters and the first five vectors. The eccentricities of the ellipses indicate the magnitude of the correlation, and the color the importance of the correlation in units of σ.</figcaption></figure>

a model of the SED at maximum light is derived. This model allows us to deduce the average extinction curve and the procedure to estimate it: the results are described in the following.

### 3.3.1. Empirical description

The SED model at maximum light presented here is a generalization of the model presented in Chotard et al. (2011), in which the

authors removed intrinsic variability tied to spectral features and measured the remaining spectral variation, which was found to be consistent with an extinction curve like that for the Galaxy. The major difference is that the intrinsic description is here motivated by the factors derived in Sect. 3.2 rather than being described only by pEW(Si II λ4131) and pEW(Ca II H&K) as was done in Chotard et al. (2011). To model the SED we define for the SN Ia the vector $\mathbf{x}_i \equiv \{h_{i,1}, h_{i,2}, h_{i,3}\}$, which is the true value of the measured factor $\mathbf{q}_i = \{q_{i,1}, q_{i,2}, q_{i,3}\}$. We propose that these are related to the intrinsic absolute magnitude by:

$$M_i(\lambda) = M_0(\lambda) + \sum_{j=1}^{j=3} h_{i,j} \alpha_j(\lambda) + A_{\lambda_0, i} \gamma(\lambda), \quad (5)$$

where $M_0(\lambda)$ is the average spectrum in absolute magnitude, $\alpha_j(\lambda)$ is the spectrum related to the factor $h_{i,j}$, where j is the factor index running from one to three (the choice of the number of components is discussed in Sect. 6.2). $A_{\lambda_0, i}$ is a free parameter that can be interpreted as the absorption due to extinction at a reference wavelength $\lambda_0$ (set here at the median wavelength, 5305 Å, without loss of generality). $\gamma(\lambda)$ is an arbitrary function that represents the effect of the extinction, and we have not set any prior on its shape. Once the function $\gamma(\lambda)$ is fixed, it is possible to fit the total-to-selective extinction ratio of the Cardelli et al. (1989) law, $R_V$, and the absorption in V-band, $A_V$ to it; this is described in the next sub-section. The model parameters here are $M_0(\lambda)$, the $\alpha_j(\lambda)$, the $\gamma(\lambda)$, the $h_{i,j}$, and the $A_{\lambda_0, i}$. It is understood that we are modeling the SED at specific wavelengths where we have measurements. For notational efficiency we rewrite Eq. (5) as:

$$\mathbf{M}_i = A \mathbf{h}_i, \quad \text{where} \quad (6)$$

$$\mathbf{A} = (\mathbf{M}_0, \gamma, \alpha_1, \alpha_2, \alpha_3), \quad (7)$$

$$\mathbf{h}_i^T = (1, A_{\lambda_0, i}, h_{i,1}, h_{i,2}, h_{i,3}). \quad (8)$$
In Chotard et al. (2011), the parameters $\alpha_j(\lambda)$ were computed in sequential order, and the reddening curve $\gamma(\lambda)$ as well as the extinction $A_{\lambda, i}$ were determined in a second step. We improve on this procedure by applying a global fit for all parameters at once. This is described below, and allows us to derive a global average $R_V$ for all SNe Ia, and an $A_V$ for each SN Ia.

### 3.3.2. Fitting the model at maximum light and $R_V$

In this section we explain the framework for fitting the free parameters of Eq. (5) and how we then determine the $A_V$ and $R_V$ parameters. Fitting the parameters of the absolute SED of Eq. (5) is done using an orthogonal distance regression. Estimating these parameters within the framework of an orthogonal distance regression amounts to minimizing the following $\chi^2$:

$$ \chi^2 = \sum_i (\mathbf{M}_i^{\text{obs}} - \mathbf{A}\mathbf{h}_i)^T \mathbf{W}_{\mathbf{M}_i} (\mathbf{M}_i^{\text{obs}} - \mathbf{A}\mathbf{h}_i) + (\mathbf{q}_i - \mathbf{x}_i)^T \mathbf{W}_{\mathbf{q}_i} (\mathbf{q}_i - \mathbf{x}_i), \quad (9) $$

where $\mathbf{M}_i^{\text{obs}}$ is the observed spectrum of the SN Ia $i$ in absolute AB magnitude and $\mathbf{W}_{\mathbf{M}_i}$ is the weight matrix of $\mathbf{M}_i^{\text{obs}}$, defined as:

$$ \mathbf{W}_{\mathbf{M}_i} = \begin{bmatrix} \cdots & 0 \\ 0 & \sigma_{\lambda i}^2 & \cdots \end{bmatrix} + (\sigma_{\text{cal}}^2 + \sigma_z^2) \begin{bmatrix} 1 & \cdots & 1 \\ \vdots & \ddots & \vdots \\ 1 & \cdots & 1 \end{bmatrix} + \mathbf{D} \quad (10) $$

where $\sigma_{\lambda i}$ is the uncertainty of the SN Ia $i$ at the wavelength $\lambda$ derived from spectral error, $\sigma_{\text{cal}}$ the per-spectrum calibration uncertainty, taken to be 0.03 mag, $\sigma_z$ the combination of the redshift error and the uncertainty of 300 km s$^{-1}$ due to peculiar velocities, and $\mathbf{D}$ is the dispersion matrix. The dispersion matrix is fit to deal with the remaining variability which is not described by the three factors and the extinction. It includes an estimate of the remaining chromatic and achromatic (gray) dispersions. The estimation of $\mathbf{D}$ is described in Appendix D. $\mathbf{W}_{\mathbf{q}_i}$ is the weight matrix of $\mathbf{q}_i$ that comes from projecting the spectral feature uncertainties into the factor sub-space, and is defined as:

$$ \mathbf{W}_{\mathbf{q}_i} = \Lambda^T \mathbf{\Psi}_i^{-1} \Lambda, \quad (11) $$

where $\Lambda$ and $\Psi_i$ are the same as in Sect. 3.2. The unknown parameters of the model are the matrix $\mathbf{A}$ and the vector $\mathbf{h}_i$. These are found by minimizing Eq. (9) using an expectation-minimization algorithm described in Appendix C.

Once the parameters that minimize Eq. (9) have been found, we can compare $\gamma_\lambda$ with an extinction curve to see if it is compatible. It is straightforward to directly estimate the global mean $R_V$ from $\gamma_\lambda$ by minimizing

$$ \chi^2_{R_V} = \left[ \gamma - \left( a + \frac{1}{R_V} b \right) \right]^T \left[ \gamma - \left( a + \frac{1}{R_V} b \right) \right], \quad (12) $$

where the $a$ and $b$ vectors are the extinction curve coefficients defined in Cardelli et al. (1989). Minimizing $\chi^2_{R_V}$ as a function of $R_V$ then gives:

$$ R_V = \frac{b^T b}{\gamma^T b - a^T b}. \quad (13) $$

We note that the value we quote for $R_V$ within the context of our model depends on the assumptions adopted in the evaluation of the dispersion matrix (discussed in detail in Sect. 6.1) and therefore does not necessarily correspond to the true mean properties of dust. This uncertainty associated with separating SN Ia color behavior from dust behavior is common to all SN Ia fitting methods. Once the value of the global mean $R_V$ is fixed, we determine the absorption in the V-band, $A_V$ for each SN Ia. This amounts to minimizing

$$ \chi^2_{A_V} = \sum_i \left( \delta \tilde{\mathbf{M}}_i - A_{V,i} \left( \mathbf{a} + \frac{1}{R_V} \mathbf{b} \right) \right)^T \mathbf{W}_{\mathbf{M}_i} \times \left( \delta \tilde{\mathbf{M}}_i - A_{V,i} \left( \mathbf{a} + \frac{1}{R_V} \mathbf{b} \right) \right), \quad (14) $$

where $A_{V,i}$ is the absorption in the V-band for the SN Ia $i$ and $\delta \tilde{\mathbf{M}}_i$ are the residuals corrected for the intrinsic variabilities

$$ \delta \tilde{\mathbf{M}}_i = \mathbf{M}_i^{\text{obs}} - \mathbf{M}_0 - \sum_j q_{i,j} \alpha_j. \quad (15) $$

Taking the derivative of $\chi^2_{A_V}$ by $A_V$ we find:

$$ A_{V,i} = \frac{\left( \mathbf{a} + \frac{1}{R_V} \mathbf{b} \right)^T \mathbf{W}_{\mathbf{M}_i} \delta \tilde{\mathbf{M}}_i}{\left( \mathbf{a} + \frac{1}{R_V} \mathbf{b} \right)^T \mathbf{W}_{\mathbf{M}_i} \left( \mathbf{a} + \frac{1}{R_V} \mathbf{b} \right)}. \quad (16) $$

### 3.3.3. Results

The model above is trained on the 105 SNe Ia that remain in the training sample. The vectors $\alpha_j$ for $j = 1, j = 2$, and $j = 3$, as well as their effect on the average spectrum $\mathbf{M}_0$, are presented in Fig. 5. To quantify the overall impact of each vector in units of magnitude, we compute the average rms of $h_j \alpha_j$, the deviation from the average spectrum $\mathbf{M}_0$.

The vector $\alpha_1$ presented in Fig. 5, has an average impact of 0.17 mag when multiplied by the standard deviation of $q_1$. This is the expected amplitude for a stretch effect. The structure of this vector is associated with variations of the line depths, consistent with the correlation of the factor $q_1$ with both the pseudo-equivalent widths and the stretch (cf. Figs. 3 and 4).

The vector $\alpha_2$ shown in Fig. 5, which is correlated with velocities, has a much weaker impact of 0.05 mag on the average spectrum $\mathbf{M}_0$. This is consistent with the fact that this effect has not yet been detected on purely photometric data. Moreover, the associated variability is centered in localized structures such as the regions of the Ca II H&K, Si II $\lambda$6355, and Ca II IR lines. A closer scrutiny shows that the variability especially affects the bluer edge of the features, leading to an overall effect on the line velocities. This is as expected from the correlations of $q_2$ and the minima of P-Cygni profile shown in Fig. 3.

The vector $\alpha_3$ shown in Fig. 5, has an impact of 0.05 mag on the average spectrum $\mathbf{M}_0$. Like the first component $q_1$, $q_3$ is correlated with the stretch (cf. Fig. 4), but the corresponding vector $\alpha_3$ is less structured than $\alpha_1$ in the optical band. However, the expected correlations with the spectral features are consistent with the behavior of $q_3$ (cf. Fig. 3). Moreover, the most prominent features affect the extreme UV part of the spectrum as well as the Ca II H&K and Ca II IR regions. For the latter, brighter SNe Ia exhibit a stronger trough in the higher velocity part of the blends, which could link this vector to the presence of high-velocity calcium structure in ejectas.

In summary, the analysis of the SED components at maximum light confirms that the stretch has the dominant effect on
In Chotard et al. (2011), the parameters $\alpha_j(\lambda)$ were computed in sequential order, and the reddening curve $\gamma(\lambda)$ as well as the extinction $A_{\lambda, i}$ were determined in a second step. We improve on this procedure by applying a global fit for all parameters at once. This is described below, and allows us to derive a global average $R_V$ for all SNe Ia, and an $A_V$ for each SN Ia.

### 3.3.2. Fitting the model at maximum light and $R_V$

In this section we explain the framework for fitting the free parameters of Eq. (5) and how we then determine the $A_V$ and $R_V$ parameters. Fitting the parameters of the absolute SED of Eq. (5) is done using an orthogonal distance regression. Estimating these parameters within the framework of an orthogonal distance regression amounts to minimizing the following $\chi^2$:

$$ \chi^2 = \sum_i (\mathbf{M}_i^{\text{obs}} - \mathbf{A}\mathbf{h}_i)^T \mathbf{W}_{\mathbf{M}_i} (\mathbf{M}_i^{\text{obs}} - \mathbf{A}\mathbf{h}_i) + (\mathbf{q}_i - \mathbf{x}_i)^T \mathbf{W}_{\mathbf{q}_i} (\mathbf{q}_i - \mathbf{x}_i), \quad (9) $$

where $\mathbf{M}_i^{\text{obs}}$ is the observed spectrum of the SN Ia $i$ in absolute AB magnitude and $\mathbf{W}_{\mathbf{M}_i}$ is the weight matrix of $\mathbf{M}_i^{\text{obs}}$, defined as:

$$ \mathbf{W}_{\mathbf{M}_i} = \begin{bmatrix} \cdots & 0 \\ 0 & \sigma_{\lambda i}^2 & \cdots \\ \vdots & \vdots & \vdots \\ 0 & \cdots & 1 \end{bmatrix} + (\sigma_{\text{cal}}^2 + \sigma_z^2) \begin{bmatrix} 1 & \cdots & 1 \\ \vdots & \ddots & \vdots \\ 1 & \cdots & 1 \end{bmatrix} + \mathbf{D} \quad (10) $$

where $\sigma_{\lambda i}$ is the uncertainty of the SN Ia $i$ at the wavelength $\lambda$ derived from spectral error, $\sigma_{\text{cal}}$ the per-spectrum calibration uncertainty, taken to be 0.03 mag, $\sigma_z$ the combination of the redshift error and the uncertainty of 300 km s$^{-1}$ due to peculiar velocities, and $\mathbf{D}$ is the dispersion matrix. The dispersion matrix is fit to deal with the remaining variability which is not described by the three factors and the extinction. It includes an estimate of the remaining chromatic and achromatic (gray) dispersions. The estimation of $\mathbf{D}$ is described in Appendix D. $\mathbf{W}_{\mathbf{q}_i}$ is the weight matrix of $\mathbf{q}_i$ that comes from projecting the spectral feature uncertainties into the factor sub-space, and is defined as:

$$ \mathbf{W}_{\mathbf{q}_i} = \Lambda^T \mathbf{\Psi}_i^{-1} \Lambda, \quad (11) $$

where $\Lambda$ and $\Psi_i$ are the same as in Sect. 3.2. The unknown parameters of the model are the matrix $\mathbf{A}$ and the vector $\mathbf{h}_i$. These are found by minimizing Eq. (9) using an expectation-minimization algorithm described in Appendix C.

Once the parameters that minimize Eq. (9) have been found, we can compare $\gamma_\lambda$ with an extinction curve to see if it is compatible. It is straightforward to directly estimate the global mean $R_V$ from $\gamma_\lambda$ by minimizing

$$ \chi^2_{R_V} = \left[ \gamma - \left( a + \frac{1}{R_V} b \right) \right]^T \left[ \gamma - \left( a + \frac{1}{R_V} b \right) \right], \quad (12) $$

where the $a$ and $b$ vectors are the extinction curve coefficients defined in Cardelli et al. (1989). Minimizing $\chi^2_{R_V}$ as a function of $R_V$ then gives:

$$ R_V = \frac{b^T b}{\gamma^T b - a^T b}. \quad (13) $$

We note that the value we quote for $R_V$ within the context of our model depends on the assumptions adopted in the evaluation of the dispersion matrix (discussed in detail in Sect. 6.1) and therefore does not necessarily correspond to the true mean properties of dust. This uncertainty associated with separating SN Ia color behavior from dust behavior is common to all SN Ia fitting methods. Once the value of the global mean $R_V$ is fixed, we determine the absorption in the V-band, $A_V$ for each SN Ia. This amounts to minimizing

$$ \chi^2_{A_V} = \sum_i \left( \delta \tilde{\mathbf{M}}_i - A_{V,i} \left( \mathbf{a} + \frac{1}{R_V} \mathbf{b} \right) \right)^T \mathbf{W}_{\mathbf{M}_i} \times \left( \delta \tilde{\mathbf{M}}_i - A_{V,i} \left( \mathbf{a} + \frac{1}{R_V} \mathbf{b} \right) \right), \quad (14) $$

where $A_{V,i}$ is the absorption in the V-band for the SN Ia $i$ and $\delta \tilde{\mathbf{M}}_i$ are the residuals corrected for the intrinsic variabilities

$$ \delta \tilde{\mathbf{M}}_i = \mathbf{M}_i^{\text{obs}} - \mathbf{M}_0 - \sum_j q_{i,j} \alpha_j. \quad (15) $$

Taking the derivative of $\chi^2_{A_V}$ by $A_V$ we find:

$$ A_{V,i} = \frac{\left( \mathbf{a} + \frac{1}{R_V} \mathbf{b} \right)^T \mathbf{W}_{\mathbf{M}_i} \delta \tilde{\mathbf{M}}_i}{\left( \mathbf{a} + \frac{1}{R_V} \mathbf{b} \right)^T \mathbf{W}_{\mathbf{M}_i} \left( \mathbf{a} + \frac{1}{R_V} \mathbf{b} \right)}. \quad (16) $$

### 3.3.3. Results

The model above is trained on the 105 SNe Ia that remain in the training sample. The vectors $\alpha_j$ for $j = 1, j = 2$, and $j = 3$, as well as their effect on the average spectrum $\mathbf{M}_0$, are presented in Fig. 5. To quantify the overall impact of each vector in units of magnitude, we compute the average rms of $h_j \alpha_j$, the deviation from the average spectrum $\mathbf{M}_0$.

The vector $\alpha_1$ presented in Fig. 5, has an average impact of 0.17 mag when multiplied by the standard deviation of $q_1$. This is the expected amplitude for a stretch effect. The structure of this vector is associated with variations of the line depths, consistent with the correlation of the factor $q_1$ with both the pseudo-equivalent widths and the stretch (cf. Figs. 3 and 4).

The vector $\alpha_2$ shown in Fig. 5, which is correlated with velocities, has a much weaker impact of 0.05 mag on the average spectrum $\mathbf{M}_0$. This is consistent with the fact that this effect has not yet been detected on purely photometric data. Moreover, the associated variability is centered in localized structures such as the regions of the Ca II H&K, Si II $\lambda$6355, and Ca II IR lines. A closer scrutiny shows that the variability especially affects the bluer edge of the features, leading to an overall effect on the line velocities. This is as expected from the correlations of $q_2$ and the minima of P-Cygni profile shown in Fig. 3.

The vector $\alpha_3$ shown in Fig. 5, has an impact of 0.05 mag on the average spectrum $\mathbf{M}_0$. Like the first component $q_1$, $q_3$ is correlated with the stretch (cf. Fig. 4), but the corresponding vector $\alpha_3$ is less structured than $\alpha_1$ in the optical band. However, the expected correlations with the spectral features are consistent with the behavior of $q_3$ (cf. Fig. 3). Moreover, the most prominent features affect the extreme UV part of the spectrum as well as the Ca II H&K and Ca II IR regions. For the latter, brighter SNe Ia exhibit a stronger trough in the higher velocity part of the blends, which could link this vector to the presence of high-velocity calcium structure in ejectas.

In summary, the analysis of the SED components at maximum light confirms that the stretch has the dominant effect on
<figure><img src="image_6.png" /></figure>

**Fig. 5.** Impact of each component on the average spectrum. Top panel, top: average spectrum $M_0$ and the predicted effect of a variation of $q_1$ by $\pm 1\sigma$. Bottom: corresponding $\alpha_1$ vector. Middle panel, top: average spectrum $M_0$ and the predicted effect of a variation of $q_2$ by $\pm 1\sigma$. Bottom: corresponding $\alpha_2$ vector. Bottom panel, top: average spectrum $M_0$ and the predicted effect of a variation of $q_2$ by $\pm 1\sigma$. Bottom: corresponding $\alpha_2$ vector.

magnitudes as expected, but one also has to take into account other variabilities which are difficult to detect in photometric bands since they are linked to localized features such as velocities.

The $\gamma(\lambda)$ curve is shown in Fig. 6 and is consistent with a dust extinction curve. This shows that three eigenvectors provide a sufficient description of the intrinsic variability from the purpose of deriving a color law. The fit of $\gamma(\lambda)$ to a Cardelli et al. (1989)

<figure><img src="image_7.png" /></figure>

**Fig. 6.** Extinction curve fit result. Top panel: empirical average extinction curve, $\gamma(\lambda)$, compared to the best fit by a Cardelli et al. (1989) law, which gives an $R_V = 2.6$. The shaded blue area represent a $\pm 0.5$ range on $R_V$, in order to illustrate the typical range of value that are measured for SN Ia in the literature. The black dot on this curve indicates the wavelength used on the lower panel graph. Lower panel: residuals after correction for the intrinsic behavior at 4012 Å as a function of $A_{\lambda 0}$. Each dot corresponds to a single SN Ia. The $\gamma_{4012\AA}$ slope between $M$ and $A_{\lambda 0}$ is indicated by the red continuous line. Dashed represent the Cardelli et al. (1989) law for $R_V = 2.6$. The shaded blue area represent a $\pm 0.5$ range on $R_V$, in order to illustrate the typical range of value that are measured for SN Ia in the literature.

law gives a value of $R_V = 2.6$. The Cardelli et al. (1989) law fit coincides remarkably well with $\gamma(\lambda)$, except in the UV, where it differs slightly. The $R_V$ value is mainly driven by the reddest SNe Ia. Indeed, as can be seen in the lower panel of Fig. 6, three SNe Ia have the dominant contribution on the final value of $R_V$ and removing any of these from the sample would significantly alter the result, hence the rather large value of the uncertainty.

We made the choice of fitting an average extinction curve for the whole sample, but there are indications from observations that extinction curves exhibit some diversity (Amanullah et al. 2015). Moreover, any such extinction variation, or intrinsic color variation (e.g. Foley & Kasen 2011; Polin et al. 2019), that the model at maximum is unable to capture goes into the dispersion matrix. However, in order to keep the model simple, we keep an average value of $R_V$ and we leave $R_V$ variation analysis for the future.

Figure 7 presents the set of spectra before any correction and the residuals resulting from the difference between the observed spectrum and Eq. (5). In the top panel, the effect of extinction is clearly visible for spectra before any correction. After correction, the residuals between the models and data are mainly shifted in magnitude with respect to one another and no effect on color appears to remain. This is the signature expected for a mostly gray residual offset.

The amplitude of the correction made by the three factors $q_1$, $q_2$, $q_3$, and the extinction curve, given by the wRMS as a function of wavelength, is presented in the lower panel of Fig. 7. Before any correction, the dispersion in the B-band is around 0.4 mag, as expected. Spectra before any correction exhibit both
<figure><img src="image_6.png" /></figure>

**Fig. 5.** Impact of each component on the average spectrum. Top panel, top: average spectrum $M_0$ and the predicted effect of a variation of $q_1$ by $\pm 1\sigma$. Bottom: corresponding $\alpha_1$ vector. Middle panel, top: average spectrum $M_0$ and the predicted effect of a variation of $q_2$ by $\pm 1\sigma$. Bottom: corresponding $\alpha_2$ vector. Bottom panel, top: average spectrum $M_0$ and the predicted effect of a variation of $q_2$ by $\pm 1\sigma$. Bottom: corresponding $\alpha_2$ vector.

magnitudes as expected, but one also has to take into account other variabilities which are difficult to detect in photometric bands since they are linked to localized features such as velocities.

The $\gamma(\lambda)$ curve is shown in Fig. 6 and is consistent with a dust extinction curve. This shows that three eigenvectors provide a sufficient description of the intrinsic variability from the purpose of deriving a color law. The fit of $\gamma(\lambda)$ to a Cardelli et al. (1989)

<figure><img src="image_7.png" /></figure>

**Fig. 6.** Extinction curve fit result. Top panel: empirical average extinction curve, $\gamma(\lambda)$, compared to the best fit by a Cardelli et al. (1989) law, which gives an $R_V = 2.6$. The shaded blue area represent a $\pm 0.5$ range on $R_V$, in order to illustrate the typical range of value that are measured for SN Ia in the literature. The black dot on this curve indicates the wavelength used on the lower panel graph. Lower panel: residuals after correction for the intrinsic behavior at 4012 Å as a function of $A_{\lambda 0}$. Each dot corresponds to a single SN Ia. The $\gamma_{4012\AA}$ slope between $M$ and $A_{\lambda 0}$ is indicated by the red continuous line. Dashed represent the Cardelli et al. (1989) law for $R_V = 2.6$. The shaded blue area represent a $\pm 0.5$ range on $R_V$, in order to illustrate the typical range of value that are measured for SN Ia in the literature.

law gives a value of $R_V = 2.6$. The Cardelli et al. (1989) law fit coincides remarkably well with $\gamma(\lambda)$, except in the UV, where it differs slightly. The $R_V$ value is mainly driven by the reddest SNe Ia. Indeed, as can be seen in the lower panel of Fig. 6, three SNe Ia have the dominant contribution on the final value of $R_V$ and removing any of these from the sample would significantly alter the result, hence the rather large value of the uncertainty.

We made the choice of fitting an average extinction curve for the whole sample, but there are indications from observations that extinction curves exhibit some diversity (Amanullah et al. 2015). Moreover, any such extinction variation, or intrinsic color variation (e.g. Foley & Kasen 2011; Polin et al. 2019), that the model at maximum is unable to capture goes into the dispersion matrix. However, in order to keep the model simple, we keep an average value of $R_V$ and we leave $R_V$ variation analysis for the future.

Figure 7 presents the set of spectra before any correction and the residuals resulting from the difference between the observed spectrum and Eq. (5). In the top panel, the effect of extinction is clearly visible for spectra before any correction. After correction, the residuals between the models and data are mainly shifted in magnitude with respect to one another and no effect on color appears to remain. This is the signature expected for a mostly gray residual offset.

The amplitude of the correction made by the three factors $q_1$, $q_2$, $q_3$, and the extinction curve, given by the wRMS as a function of wavelength, is presented in the lower panel of Fig. 7. Before any correction, the dispersion in the B-band is around 0.4 mag, as expected. Spectra before any correction exhibit both
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Spectra before and after applying the model at maximum light. Top panel: 105 spectra in the absolute AB magnitude before correction (upper set of spectra) and after correction for intrinsic properties and color (lower set of spectra). An arbitrary offset has been applied for each set. The color code indicates the value of $A_{\lambda_0}$. Lower panel: wRMS of the residuals for the uncorrected (in red) and corrected (in blue) spectra.</figcaption></figure>

localized structure coming from intrinsic variabilities and a slowly increasing dispersion in the bluer part of the spectrum, which is the signature of extinction. Once spectra are corrected by the three factors $q_1, q_2, q_3$, and the extinction curve, the wRMS of the residuals drops to a floor of 0.1 mag. However, localized structures remain, concentrated in UV, Ca II H&K, Si II λ6355, O I λ7773, and Ca II IR regions. The amplitude of the baseline is compatible with the value of the gray intrinsic dispersion added in cosmological fits for nearby supernovae (Betoule et al. 2014). The fact that a part of the residual is gray is also visible in the dispersion matrix as can be seen in Fig. 8. Indeed, the high correlation across all wavelengths is consistent with gray fluctuation. However, some effects observed in the dispersion matrix are not due to the gray effect, but to the unmodeled variabilities not captured by SUGAR in certain spectral zones, which explains some features observed in the dispersion matrix and in the Fig. 7.

### 3.4. Time interpolation using Gaussian process

Now, we want to extend our model from maximum light to the full spectral time series. Spectra from the SNfactory are taken at different epochs in the development of each SN Ia. The treatment of these observations taken at different times is an important step to handle before constructing the full SED model. In order to work with a fixed time grid, it is necessary to use a method of interpolation. One solution is to use Gaussian processes, as was done for SNe Ia in works like Kim et al. (2013), Fakhouri et al. (2015), and Saunders et al. (2018). A complete review of Gaussian processes can be found in Rasmussen & Williams (2006). The Gaussian process implementation presented here is based on this review. However, we have developed a specific implementation in order to take into account wavelength dependencies

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Square root of the diagonal of the dispersion matrix D (top) and the correlation matrix corresponding to the dispersion matrix D (bottom). The shading represents the degree of correlation, as given by the Pearson correlation coefficient, $\rho$.</figcaption></figure>

and to accelerate computation of the interpolation. These are described below.

A Gaussian process is a generalization of the Gaussian probability distribution, and is a non-parametric way to interpolate data. The main assumption of Gaussian processes is that observed data are the realization of Gaussian random fields that are characterized by an average function and a correlation function. Therefore, the distribution of SNe Ia magnitude $m$ at a given phase $t$ and a given wavelength $\lambda$ follows this relation:

$$m(t, \lambda) \sim N(m_0(t, \lambda), K_{\lambda}), \quad (17)$$

where $m_0(t, \lambda)$ is the average function of the SNe Ia SED, and $K_{\lambda}$ corresponds to the correlation matrix (commonly called the kernel) that describes time correlations between different epochs. In the absence of an explosion model that would yield the analytical form of $K_{\lambda}$, we choose a squared exponential kernel, to which the measurements and calibration uncertainties are added:

$$K_{\lambda}(t_i, t_j) = \phi_{\lambda}^2 \exp \left[ -\frac{1}{2} \left( \frac{t_i - t_j}{l_{\lambda}} \right)^2 \right] + \sigma_{\text{gray}}^2 \delta'_{ij}, \quad (18)$$

where $\phi_{\lambda}^2$ represent the spectral variance around the average function for the wavelength $\lambda$, $l_{\lambda}$ corresponds to the temporal correlation length for the wavelength $\lambda$, $t_i$ and $t_j$ correspond to the phase of observation $i$ $j$, and $\sigma_{\text{gray}}$ is the gray error taken to be 0.03 mag. All wavelengths are treated independently of each other, that is the interpolation is performed on $y'_{\lambda n}$, the nth SN Ia light curve for each wavelength $\lambda$ from which the average function has been subtracted. For a given set of global hyperparameters ($\phi_{\lambda}, l_{\lambda}$), the interpolation $y'_{\lambda n}$ and the covariance matrix of uncertainties on the interpolation, $\text{cov}(y'_{\lambda n})$, are given by:

$$y'_{\lambda n} = K_{\lambda n}(t', t)^T (K_{\lambda n}(t, t) + \sigma_{\lambda}^2 I)^{-1} y_{\lambda n} \quad (19)$$

$$\text{cov}(y'_{\lambda n}) = K_{\lambda n}(t', t') - K_{\lambda n}(t, t)^T (K_{\lambda n}(t, t) + \sigma_{\lambda}^2 I)^{-1} K_{\lambda n}(t', t), \quad (20)$$

where the vector $\sigma_{\lambda}^2$ is the error on the magnitude for the wavelength $\lambda$ and for all observed phases. $t'$ is the new time grid and
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Spectra before and after applying the model at maximum light. Top panel: 105 spectra in the absolute AB magnitude before correction (upper set of spectra) and after correction for intrinsic properties and color (lower set of spectra). An arbitrary offset has been applied for each set. The color code indicates the value of $A_{\lambda_0}$. Lower panel: wRMS of the residuals for the uncorrected (in red) and corrected (in blue) spectra.</figcaption></figure>

localized structure coming from intrinsic variabilities and a slowly increasing dispersion in the bluer part of the spectrum, which is the signature of extinction. Once spectra are corrected by the three factors $q_1, q_2, q_3$, and the extinction curve, the wRMS of the residuals drops to a floor of 0.1 mag. However, localized structures remain, concentrated in UV, Ca II H&K, Si II λ6355, O I λ7773, and Ca II IR regions. The amplitude of the baseline is compatible with the value of the gray intrinsic dispersion added in cosmological fits for nearby supernovae (Betoule et al. 2014). The fact that a part of the residual is gray is also visible in the dispersion matrix as can be seen in Fig. 8. Indeed, the high correlation across all wavelengths is consistent with gray fluctuation. However, some effects observed in the dispersion matrix are not due to the gray effect, but to the unmodeled variabilities not captured by SUGAR in certain spectral zones, which explains some features observed in the dispersion matrix and in the Fig. 7.

### 3.4. Time interpolation using Gaussian process

Now, we want to extend our model from maximum light to the full spectral time series. Spectra from the SNfactory are taken at different epochs in the development of each SN Ia. The treatment of these observations taken at different times is an important step to handle before constructing the full SED model. In order to work with a fixed time grid, it is necessary to use a method of interpolation. One solution is to use Gaussian processes, as was done for SNe Ia in works like Kim et al. (2013), Fakhouri et al. (2015), and Saunders et al. (2018). A complete review of Gaussian processes can be found in Rasmussen & Williams (2006). The Gaussian process implementation presented here is based on this review. However, we have developed a specific implementation in order to take into account wavelength dependencies

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Square root of the diagonal of the dispersion matrix **D** (top) and the correlation matrix corresponding to the dispersion matrix **D** (bottom). The shading represents the degree of correlation, as given by the Pearson correlation coefficient, $\rho$.</figcaption></figure>

and to accelerate computation of the interpolation. These are described below.

A Gaussian process is a generalization of the Gaussian probability distribution, and is a non-parametric way to interpolate data. The main assumption of Gaussian processes is that observed data are the realization of Gaussian random fields that are characterized by an average function and a correlation function. Therefore, the distribution of SNe Ia magnitude $m$ at a given phase $t$ and a given wavelength $\lambda$ follows this relation:

$$m(t, \lambda) \sim N(m_0(t, \lambda), \mathbf{K}_{\lambda}), \quad (17)$$

where $m_0(t, \lambda)$ is the average function of the SNe Ia SED, and $K_{\lambda}$ corresponds to the correlation matrix (commonly called the kernel) that describes time correlations between different epochs. In the absence of an explosion model that would yield the analytical form of $K_{\lambda}$, we choose a squared exponential kernel, to which the measurements and calibration uncertainties are added:

$$K_{\lambda}(t_i, t_j) = \phi_{\lambda}^2 \exp \left[ -\frac{1}{2} \left( \frac{t_i - t_j}{l_{\lambda}} \right)^2 \right] + \sigma_{\text{gray}}^2 \delta'_{ij}, \quad (18)$$

where $\phi_{\lambda}^2$ represent the spectral variance around the average function for the wavelength $\lambda$, $l_{\lambda}$ corresponds to the temporal correlation length for the wavelength $\lambda$, $t_i$ and $t_j$ correspond to the phase of observation $i$ $j$, and $\sigma_{\text{gray}}$ is the gray error taken to be 0.03 mag. All wavelengths are treated independently of each other, that is the interpolation is performed on $y'_{\lambda n}$, the nth SN Ia light curve for each wavelength $\lambda$ from which the average function has been subtracted. For a given set of global hyperparameters ($\phi_{\lambda}, l_{\lambda}$), the interpolation $y'_{\lambda n}$ and the covariance matrix of uncertainties on the interpolation, $\text{cov}(y'_{\lambda n})$, are given by:

$$y'_{\lambda n} = \mathbf{K}_{\lambda n}(t', t)^T (\mathbf{K}_{\lambda n}(t, t) + \sigma_{\lambda}^2 \mathbf{I})^{-1} y_{\lambda n} \quad (19)$$

$$\text{cov}(y'_{\lambda n}) = \mathbf{K}_{\lambda n}(t', t') - \mathbf{K}_{\lambda n}(t, t)^T (\mathbf{K}_{\lambda n}(t, t) + \sigma_{\lambda}^2 \mathbf{I})^{-1} \mathbf{K}_{\lambda n}(t', t), \quad (20)$$

where the vector $\sigma_{\lambda}^2$ is the error on the magnitude for the wavelength $\lambda$ and for all observed phases. $t'$ is the new time grid and
t the observed phases of the SN Ia. Gaussian processes assume that the observed data are distributed around an average function. This function is a necessary input, because in a region with no data points, the Gaussian process interpolation converges to this average function (after ~3 temporal correlation lengths). The calculation of this average function is performed as follows. The training spectra are binned according to phase from maximum light in B-band, ranging from -12 to +48 days. Each phase bin spans two days and contains at least three spectra. The weighted average spectrum is computed for each phase bin and given as a first approximation of the average function. However, due to some lower signal-to-noise spectra, this version of the average function is not smooth. A Savitsky-Gollay filter is used to smooth the average function, which is then used for the Gaussian process interpolation.

What remains is to determine the pair of hyperparameters $\theta_{\lambda} = (\phi_{\lambda}, l_{\lambda})$ that represent, respectively, the amplitude of the fluctuation around the average function and the temporal correlation length. To estimate the hyperparameters for each wavelength, we maximize the product of all individual likelihoods:

$$ \mathcal{L}_{\lambda} = \prod_{n} \frac{1}{2\pi^{\frac{N}{2}}} \frac{1}{|\mathbf{K}_{\lambda n}(\mathbf{t}, \mathbf{t}) + \sigma_{\lambda}^{2} \mathbf{I}|^{\frac{1}{2}}} \times \exp\left(-\frac{1}{2} \mathbf{y}_{\lambda n}^{T}(\mathbf{K}_{\lambda n}(\mathbf{t}, \mathbf{t}) + \sigma_{\lambda}^{2} \mathbf{I})^{-1} \mathbf{y}_{\lambda n}\right). \quad (21) $$

This differs from previous work for two main reasons: the wavelength dependence of the hyperparameters, and their estimation using all the SNe Ia. These two new features are justified mainly by SNe Ia physics and computation time. Since SNe Ia are standardizable candles, we treat them as realizations of the same Gaussian random field, that is we assume that they share the same set of hyperparameters. In this case, we can use all available SNe Ia to estimate these hyperparameters. The wavelength dependence of the hyperparameters is due to as example dust extinction, intrinsic variability, the second maximum in infrared, etc. It follows that the standard deviation around the average function and the temporal correlation length should vary across wavelength. Moreover, the wavelength dependence allows us to be very efficient in terms of computation time. Indeed, it sped up matrix inversion from the classical $O(N_{\lambda}^3 \times N_t^3)$ to $O(N_{\lambda} \times N_t^3)$, where $N_{\lambda}$ is the number of bin in wavelength (197 here), and $N_t$ the number of observed phases (~14 in average). On a Mac Book Pro with a 2.9 GHz Intel Core i5 processor and 8 GB of RAM, finding all hyperparameters ($2 \times 197$), and computing the interpolation, and the pull distribution of residuals took under 6 min.

Results of hyperparameter adjustment in terms of wavelength are shown in Fig. 9. First, it can be seen that the structure of hyperparameters as a function of wavelength is not random. Indeed, for the parameter $\phi_{\lambda}$, its value varies from 0.3 mag to more than 0.6 mag: this number is difficult to interpret directly because it is the amplitude around an average value at all phases of the data. Nevertheless, the increase in the ultraviolet and blue wavelengths can be explained by the dispersion caused by dust extinction of the host-galaxies and a larger intrinsic dispersion in this wavelength range (Maguire et al. 2012; Nordin et al. 2018). There are also structures in the peaks of this intrinsic dispersion that correspond to the regions of Ca II H&K, Si II $\lambda$6355, and Ca II IR lines. The correlation length, $l_{\lambda}$, has structure in the silicon and calcium regions and remains globally stable with no features between 6000 Å and 8000 Å (average value of seven days). Once the hyperparameters have been estimated, the interpolations for each SN Ia are performed onto a new phase grid. The

<figure><img src="image_10.png" /><figcaption>Fig. 9. Results of the Gaussian process interpolation. From top to bottom: amplitude of the kernel, correlation length and standard deviation of the pull distribution in terms of wavelength. The gray areas represent the presence of the absorption line that were used in our original factor analysis.</figcaption></figure>

phase grid chosen covers a range between -12 and +48 days around maximum brightness in B-band, divided into three day bins. Finally, to check if both the predictions (interpolation) and associated errors are in agreement with what is expected from a Gaussian process, we compute the pull distribution of the residuals. Because we assume that our data are realizations of a Gaussian random fields, the pull distribution should follow a centered unit normal distribution. For each wavelength we fit a Gaussian to the pull distribution and calculate its standard deviation. The results are shown in the bottom of Fig. 9. The pull standard deviation is on average lower than 1 with an average value of 0.8. These interpolation results are used to establish the SED model that is described in the following sections.

# 4. The SUGAR model

## 4.1. The model

The SUGAR model assumes that SED variation at any epoch may be described by spectral features measured at maximum. The SUGAR SED model is based on a linear combination of factors derived in Sect. 3.2.3 and an extinction curve taken as a Cardelli et al. (1989) law derived in Sect. 3.3. To model the SUGAR SED, we propose a model similar to the one we constructed for maximum light. This SUGAR model is based on the three factors for each SN Ia $i$, combined into the vector $\mathbf{x}_i = \{h_{i,1}, h_{i,2}, h_{i,3}\}$, which in the model is the true value of the measured factor $\mathbf{q}_i = \{q_{i,1}, q_{i,2}, q_{i,3}\}$, and are related to the intrinsic magnitude by:

$$ M_i(t, \lambda) = M_0(t, \lambda) + \sum_{j=1}^{j=3} h_{i,j} \alpha_j(t, \lambda) + A_{V,i} f(\lambda, R_V) + \Delta M_{\text{gray}, i}, \quad (22) $$

where $M_0(t, \lambda)$ is the magnitude of the average spectral time series, $\alpha_j(t, \lambda)$ is the intrinsic variation related to the factor $h_{i,j}$; as discussed in Sect. 3.2.3, the numbers of factors is set to three. As at maximum, $h_{i,j}$ represents the true value of the measured factor $q_{i,j}$ derived in Sect. 3.2. $A_{V,i}$ is the extinction in the V-band and $f(\lambda, R_V)$ is the Cardelli et al. (1989) law where $R_V$ is the extinction ratio that is derived in Sect. 3.3. Finally, the term
t the observed phases of the SN Ia. Gaussian processes assume that the observed data are distributed around an average function. This function is a necessary input, because in a region with no data points, the Gaussian process interpolation converges to this average function (after ~3 temporal correlation lengths). The calculation of this average function is performed as follows. The training spectra are binned according to phase from maximum light in B-band, ranging from -12 to +48 days. Each phase bin spans two days and contains at least three spectra. The weighted average spectrum is computed for each phase bin and given as a first approximation of the average function. However, due to some lower signal-to-noise spectra, this version of the average function is not smooth. A Savitsky-Gollay filter is used to smooth the average function, which is then used for the Gaussian process interpolation.

What remains is to determine the pair of hyperparameters $\theta_{\lambda} = (\phi_{\lambda}, l_{\lambda})$ that represent, respectively, the amplitude of the fluctuation around the average function and the temporal correlation length. To estimate the hyperparameters for each wavelength, we maximize the product of all individual likelihoods:

$$ \mathcal{L}_{\lambda} = \prod_{n} \frac{1}{2\pi^{\frac{N}{2}}} \frac{1}{|\mathbf{K}_{\lambda n}(\mathbf{t}, \mathbf{t}) + \sigma_{\lambda}^{2} \mathbf{I}|^{\frac{1}{2}}} \times \exp\left(-\frac{1}{2} \mathbf{y}_{\lambda n}^{T}(\mathbf{K}_{\lambda n}(\mathbf{t}, \mathbf{t}) + \sigma_{\lambda}^{2} \mathbf{I})^{-1} \mathbf{y}_{\lambda n}\right). \quad (21) $$

This differs from previous work for two main reasons: the wavelength dependence of the hyperparameters, and their estimation using all the SNe Ia. These two new features are justified mainly by SNe Ia physics and computation time. Since SNe Ia are standardizable candles, we treat them as realizations of the same Gaussian random field, that is we assume that they share the same set of hyperparameters. In this case, we can use all available SNe Ia to estimate these hyperparameters. The wavelength dependence of the hyperparameters is due to as example dust extinction, intrinsic variability, the second maximum in infrared, etc. It follows that the standard deviation around the average function and the temporal correlation length should vary across wavelength. Moreover, the wavelength dependence allows us to be very efficient in terms of computation time. Indeed, it sped up matrix inversion from the classical $O(N_{\lambda}^3 \times N_t^3)$ to $O(N_{\lambda} \times N_t^3)$, where $N_{\lambda}$ is the number of bin in wavelength (197 here), and $N_t$ the number of observed phases (~14 in average). On a Mac Book Pro with a 2.9 GHz Intel Core i5 processor and 8 GB of RAM, finding all hyperparameters ($2 \times 197$), and computing the interpolation, and the pull distribution of residuals took under 6 min.

Results of hyperparameter adjustment in terms of wavelength are shown in Fig. 9. First, it can be seen that the structure of hyperparameters as a function of wavelength is not random. Indeed, for the parameter $\phi_{\lambda}$, its value varies from 0.3 mag to more than 0.6 mag: this number is difficult to interpret directly because it is the amplitude around an average value at all phases of the data. Nevertheless, the increase in the ultraviolet and blue wavelengths can be explained by the dispersion caused by dust extinction of the host-galaxies and a larger intrinsic dispersion in this wavelength range (Maguire et al. 2012; Nordin et al. 2018). There are also structures in the peaks of this intrinsic dispersion that correspond to the regions of Ca II H&K, Si II $\lambda$6355, and Ca II IR lines. The correlation length, $l_{\lambda}$, has structure in the silicon and calcium regions and remains globally stable with no features between 6000 Å and 8000 Å (average value of seven days). Once the hyperparameters have been estimated, the interpolations for each SN Ia are performed onto a new phase grid. The

<figure><img src="image_10.png" /><figcaption>Fig. 9. Results of the Gaussian process interpolation. From top to bottom: amplitude of the kernel, correlation length and standard deviation of the pull distribution in terms of wavelength. The gray areas represent the presence of the absorption line that were used in our original factor analysis.</figcaption></figure>

phase grid chosen covers a range between -12 and +48 days around maximum brightness in B-band, divided into three day bins. Finally, to check if both the predictions (interpolation) and associated errors are in agreement with what is expected from a Gaussian process, we compute the pull distribution of the residuals. Because we assume that our data are realizations of a Gaussian random fields, the pull distribution should follow a centered unit normal distribution. For each wavelength we fit a Gaussian to the pull distribution and calculate its standard deviation. The results are shown in the bottom of Fig. 9. The pull standard deviation is on average lower than 1 with an average value of 0.8. These interpolation results are used to establish the SED model that is described in the following sections.

# 4. The SUGAR model

## 4.1. The model

The SUGAR model assumes that SED variation at any epoch may be described by spectral features measured at maximum. The SUGAR SED model is based on a linear combination of factors derived in Sect. 3.2.3 and an extinction curve taken as a Cardelli et al. (1989) law derived in Sect. 3.3. To model the SUGAR SED, we propose a model similar to the one we constructed for maximum light. This SUGAR model is based on the three factors for each SN Ia $i$, combined into the vector $\mathbf{x}_i = \{h_{i,1}, h_{i,2}, h_{i,3}\}$, which in the model is the true value of the measured factor $\mathbf{q}_i = \{q_{i,1}, q_{i,2}, q_{i,3}\}$, and are related to the intrinsic magnitude by:

$$ M_i(t, \lambda) = M_0(t, \lambda) + \sum_{j=1}^{j=3} h_{i,j} \alpha_j(t, \lambda) + A_{V,i} f(\lambda, R_V) + \Delta M_{\text{gray}, i}, \quad (22) $$

where $M_0(t, \lambda)$ is the magnitude of the average spectral time series, $\alpha_j(t, \lambda)$ is the intrinsic variation related to the factor $h_{i,j}$; as discussed in Sect. 3.2.3, the numbers of factors is set to three. As at maximum, $h_{i,j}$ represents the true value of the measured factor $q_{i,j}$ derived in Sect. 3.2. $A_{V,i}$ is the extinction in the V-band and $f(\lambda, R_V)$ is the Cardelli et al. (1989) law where $R_V$ is the extinction ratio that is derived in Sect. 3.3. Finally, the term
$\Delta M_{\text{gray } i}$ is a gray and time independent term comparable to the parameter $X_0$ of SALT2, which allows us to work independently of the distance. Indeed, the effect of the distance on magnitude is equivalent to an additive constant independent of wavelength. Moreover, the parameter $\Delta M_{\text{gray } i}$ also contains the average of the spectral time series residuals with the SUGAR model. The addition of $\Delta M_{\text{gray } i}$ differs from what was done when determining the extinction curve in Sect. 3.3 and this choice is discussed in Sect. 6.1. We model the SED at the same wavelengths and epochs where the Gaussian process interpolation is done. For notational efficiency we rewrite Eq. (22) as:

$$ M_i = Ah_i, \quad \text{where} \tag{23} $$

$$ A = (M_0, f(R_V), 1, \alpha_1, \alpha_2, \alpha_3), \tag{24} $$

$$ h_i^T = (1, A_{V,i}, \Delta M_{\text{gray } i}, h_{i,1}, h_{i,2}, h_{i,3}). \tag{25} $$

The unknown parameters of the SUGAR model are therefore the average spectrum $M_0(t, \lambda)$, the intrinsic coefficients $\alpha_j(t, \lambda)$, and the gray offset $\Delta M_{\text{gray } i}$. Each SN Ia is parametrized during the training by the $q_{i,j}$ determined in Sect. 3.2.3 and $A_{V,i}$ derived in Sect. 3.3. The average spectral time series $M_0(t, \lambda)$, the intrinsic coefficients $\alpha_j(t, \lambda)$, and the gray offset $\Delta M_{\text{gray } i}$ are determined in the same way as the modeled described at maximum light in the Sect. 3.3. The main differences with Sect. 3.3 are that the model now depends explicitly on the phases of observation, the extinction curve is fixed at the maximum light value (same $R_V$), and the model is trained on the Gaussian process interpolation of spectra derived in Sect. 3.4 instead of measured spectra. The procedure is described below.

## 4.2. Fitting the SUGAR model

We use an orthogonal distance regression method to estimate the parameters of the SUGAR model. This is a version of the method described in Sect. 3.3, modified to include the Gaussian process interpolation, the known extinction parameters, and the gray offset. To estimate the SUGAR parameters, we minimize the following $\chi^2$:

$$ \chi^2 = \sum_i (\mathbf{M}_i^{\text{gp}} - A\mathbf{h}_i)^T \mathbf{W}_i (\mathbf{M}_i^{\text{gp}} - A\mathbf{h}_i) + (\mathbf{q}_i - \mathbf{x}_i)^T \mathbf{W}_i (\mathbf{q}_i - \mathbf{x}_i), \tag{26} $$

where the vector $\mathbf{M}_i^{\text{gp}}$ are the Gaussian process interpolations of the observed spectra (Sect. 3.4) which is in the following vector format:

$$ \mathbf{M}_i^{\text{gp}} = (\mathbf{y}'_{3340 \AA} \quad \dots \quad \mathbf{y}'_{\lambda, i} \quad \dots \quad \mathbf{y}'_{8580 \AA, i}), \tag{27} $$

where $\mathbf{y}'_{\lambda, i}$ is the light curve at wavelength $\lambda$ interpolated by Gaussian processes onto new time phase as described in Sect. 3.4. As in Sect. 3.3, the uncertainties affecting the $\mathbf{x}_i$ are propagated through the weight matrix $\mathbf{W}_{\mathbf{x}_i}$. The $\mathbf{h}_i$ plays a role analogous to the $\mathbf{q}_i$ of Sect. 3.2 and are estimated like the other free parameters of the model. Finally, the matrix $\mathbf{W}_M$ is the weight matrix of the SN Ia $i$ and follows this matrix format:

$$ \mathbf{W}_M = \begin{pmatrix} \operatorname{cov}(\mathbf{y}'_{3340 \AA, i}) & & & & & 0 \\ & \ddots & & & & \\ & & \operatorname{cov}(\mathbf{y}'_{\lambda, i}) & & & \\ & & & \ddots & & \\ 0 & & & & \operatorname{cov}(\mathbf{y}'_{8580 \AA, i}) & \\ & & & & & \end{pmatrix}^{-1}, \tag{28} $$

where $\operatorname{cov}(\mathbf{y}'_{\lambda, i})$ is the covariance of interpolation uncertainties from the Gaussian process interpolations determined in Sect. 3.4. The interpolation is done only in terms of the phase relative to B-band maximum, this makes the weight matrix block diagonal for each SN Ia. In this analysis, we do not add any covariance in terms of wavelength coming from the calibration uncertainties, for reasons of computation time and because they are expected to be small. These block diagonal matrices allowed us to accelerate the minimization algorithm. To estimate the 16 653 free parameters of the SUGAR model, it took just less than 1 h on a Mac Book Pro with a 2.9 GHz Intel Core i5 processor and 8 GB of RAM, thanks to sparse linear algebra. We note that while the covariance terms were neglected, the uncertainty at each wavelength was taken into account during the training of the Gaussian processes (these are the blocks $\operatorname{cov}(\mathbf{y}'_{\lambda, i})$). The unknown parameters of the model (the matrix $A$ and vectors $\mathbf{h}_i$) are computed by minimizing Eq. (26) using an expectation and minimization algorithm. This algorithm, described in Appendix C, avoids degeneracies between parameters.

## 4.3. The SUGAR components

The SUGAR model is described by a set of spectral time series components ($M_0, \alpha_1, \alpha_2, \alpha_3$) at three day intervals. The SED of a given SN Ia is then obtained by a linear combination of those components employing the ($q_1, q_2, q_3$) factors describing the supernova. For illustrative purposes, we present a spectral time series in six day intervals as well as their integration in synthetic top-hat filter system comprised of five bands with the following wavelength ranges: $\hat{U}$ [3300.0–3978.0] Å; $\hat{B}$ [3978.0–4795.3] Å; $\hat{V}$ [4795.3–5780.6] Å; $\hat{R}$ [5780.6–6968.3] Å; $\hat{I}$ [6968.3–8400.0] Å. The diacritic hat serves as a reminder that these are not standard Johnson-Cousins filters. The effect of the vectors $\alpha_j$ on the average spectrum $M_0$, on the average light curves (obtained by integration of $M_0$ in the $\hat{U}$, $\hat{B}$, $\hat{V}$, $\hat{R}$, $\hat{I}$ bands) and on the average colors $\hat{U} - \hat{B}$, $\hat{B} - \hat{V}$, $\hat{V} - \hat{R}$, $\hat{R} - \hat{I}$ is presented in Figs. 10 and 11 for $\alpha_1$, Figs. 12 and 13 for $\alpha_2$, and Figs. 14 and 15 for $\alpha_3$. The upper and lower contours correspond to a variation of $\pm 1\sigma$ of the associated $q_j$ parameter. The $\alpha_j(t = 0)$ components obtained when training the color law with only data at maximum or with the full time series are almost identical: all conclusions from Sect. 3.3.3 remain valid and we focus here on the temporal behavior.

Of the three factors, $q_1$ has the strongest impact: this can be seen from the time series presented in Fig. 10 and on the broad band light curves presented in Fig. 11. This is associated with a stretch effect, visible in $\hat{U}$, $\hat{B}$, and $\hat{V}$ bands as an enlargement of the variation band as we move away from maximum. Between -12 days and +30 days, the $q_1$ factor has a strong influence on spectroscopic details. The brightest SNe Ia exhibit shallower troughs in their spectral features: this is especially visible for Si II $\lambda$4131, Si II $\lambda$5972, and Ca II IR, but also valid for most of the lines. This effect is more pronounced for SNe Ia at earlier phases, and fades over time. After +30 days, the factor $q_1$ is less sensitive to localized features, and shows a relative enhancement of the optical $\hat{B}$ to $\hat{R}$ bands with respect to $\hat{U}$ and $\hat{I}$. Interestingly, between +12 days and +24 days, the brightest SNe Ia are bluer than at maximum or at later phases. This is true in absolute value, but also relative to the average SN Ia. This shows that the SNe Ia color space is driven by a variable intrinsic component in addition to reddening by dust. This intrinsic color variation is linked with the position of the second peak in $\hat{I}$, which appears later for brighter SNe Ia.
$\Delta M_{\text{gray } i}$ is a gray and time independent term comparable to the parameter $X_0$ of SALT2, which allows us to work independently of the distance. Indeed, the effect of the distance on magnitude is equivalent to an additive constant independent of wavelength. Moreover, the parameter $\Delta M_{\text{gray } i}$ also contains the average of the spectral time series residuals with the SUGAR model. The addition of $\Delta M_{\text{gray } i}$ differs from what was done when determining the extinction curve in Sect. 3.3 and this choice is discussed in Sect. 6.1. We model the SED at the same wavelengths and epochs where the Gaussian process interpolation is done. For notational efficiency we rewrite Eq. (22) as:

$$ M_i = Ah_i, \quad \text{where} \tag{23} $$

$$ A = (M_0, f(R_V), 1, \alpha_1, \alpha_2, \alpha_3), \tag{24} $$

$$ h_i^T = (1, A_{V,i}, \Delta M_{\text{gray } i}, h_{i,1}, h_{i,2}, h_{i,3}). \tag{25} $$

The unknown parameters of the SUGAR model are therefore the average spectrum $M_0(t, \lambda)$, the intrinsic coefficients $\alpha_j(t, \lambda)$, and the gray offset $\Delta M_{\text{gray } i}$. Each SN Ia is parametrized during the training by the $q_{i,j}$ determined in Sect. 3.2.3 and $A_{V,i}$ derived in Sect. 3.3. The average spectral time series $M_0(t, \lambda)$, the intrinsic coefficients $\alpha_j(t, \lambda)$, and the gray offset $\Delta M_{\text{gray } i}$ are determined in the same way as the modeled described at maximum light in the Sect. 3.3. The main differences with Sect. 3.3 are that the model now depends explicitly on the phases of observation, the extinction curve is fixed at the maximum light value (same $R_V$), and the model is trained on the Gaussian process interpolation of spectra derived in Sect. 3.4 instead of measured spectra. The procedure is described below.

## 4.2. Fitting the SUGAR model

We use an orthogonal distance regression method to estimate the parameters of the SUGAR model. This is a version of the method described in Sect. 3.3, modified to include the Gaussian process interpolation, the known extinction parameters, and the gray offset. To estimate the SUGAR parameters, we minimize the following $\chi^2$:

$$ \chi^2 = \sum_i (\mathbf{M}_i^{\text{gp}} - A\mathbf{h}_i)^T \mathbf{W}_i (\mathbf{M}_i^{\text{gp}} - A\mathbf{h}_i) + (\mathbf{q}_i - \mathbf{x}_i)^T \mathbf{W}_i (\mathbf{q}_i - \mathbf{x}_i), \tag{26} $$

where the vector $\mathbf{M}_i^{\text{gp}}$ are the Gaussian process interpolations of the observed spectra (Sect. 3.4) which is in the following vector format:

$$ \mathbf{M}_i^{\text{gp}} = (\mathbf{y}'_{3340 \AA} \quad \dots \quad \mathbf{y}'_{\lambda, i} \quad \dots \quad \mathbf{y}'_{8580 \AA} \quad \dots), \tag{27} $$

where $\mathbf{y}'_{\lambda, i}$ is the light curve at wavelength $\lambda$ interpolated by Gaussian processes onto new time phase as described in Sect. 3.4. As in Sect. 3.3, the uncertainties affecting the $\mathbf{x}_i$ are propagated through the weight matrix $\mathbf{W}_{\mathbf{x}_i}$. The $\mathbf{h}_i$ plays a role analogous to the $\mathbf{q}_i$ of Sect. 3.2 and are estimated like the other free parameters of the model. Finally, the matrix $\mathbf{W}_M$ is the weight matrix of the SN Ia $i$ and follows this matrix format:

$$ \mathbf{W}_M = \begin{pmatrix} \text{cov}(\mathbf{y}'_{3340 \AA} \mathbf{x}_i) & & & & & 0 \\ & \ddots & & & & \\ & & \text{cov}(\mathbf{y}'_{\lambda \mathbf{x}_i}) & & & \\ & & & \ddots & & \\ 0 & & & & \text{cov}(\mathbf{y}'_{8580 \AA} \mathbf{x}_i) & \\ & & & & & \end{pmatrix}^{-1}, \tag{28} $$

where $\text{cov}(\mathbf{y}'_{\lambda \mathbf{x}_i})$ is the covariance of interpolation uncertainties from the Gaussian process interpolations determined in Sect. 3.4. The interpolation is done only in terms of the phase relative to B-band maximum, this makes the weight matrix block diagonal for each SN Ia. In this analysis, we do not add any covariance in terms of wavelength coming from the calibration uncertainties, for reasons of computation time and because they are expected to be small. These block diagonal matrices allowed us to accelerate the minimization algorithm. To estimate the 16 653 free parameters of the SUGAR model, it took just less than 1 h on a Mac Book Pro with a 2.9 GHz Intel Core i5 processor and 8 GB of RAM, thanks to sparse linear algebra. We note that while the covariance terms were neglected, the uncertainty at each wavelength was taken into account during the training of the Gaussian processes (these are the blocks $\text{cov}(\mathbf{y}'_{\lambda \mathbf{x}_i})$). The unknown parameters of the model (the matrix $A$ and vectors $\mathbf{h}_i$) are computed by minimizing Eq. (26) using an expectation and minimization algorithm. This algorithm, described in Appendix C, avoids degeneracies between parameters.

## 4.3. The SUGAR components

The SUGAR model is described by a set of spectral time series components ($M_0, \alpha_1, \alpha_2, \alpha_3$) at three day intervals. The SED of a given SN Ia is then obtained by a linear combination of those components employing the ($q_1, q_2, q_3$) factors describing the supernova. For illustrative purposes, we present a spectral time series in six day intervals as well as their integration in synthetic top-hat filter system comprised of five bands with the following wavelength ranges: $\hat{U}$ [3300.0–3978.0] Å; $\hat{B}$ [3978.0–4795.3] Å; $\hat{V}$ [4795.3–5780.6] Å; $\hat{R}$ [5780.6–6968.3] Å; $\hat{I}$ [6968.3–8400.0] Å. The diacritic hat serves as a reminder that these are not standard Johnson-Cousins filters. The effect of the vectors $\alpha_j$ on the average spectrum $M_0$, on the average light curves (obtained by integration of $M_0$ in the $\hat{U}$, $\hat{B}$, $\hat{V}$, $\hat{R}$, $\hat{I}$ bands) and on the average colors $\hat{U} - \hat{B}$, $\hat{B} - \hat{V}$, $\hat{V} - \hat{R}$, $\hat{R} - \hat{I}$ is presented in Figs. 10 and 11 for $\alpha_1$, Figs. 12 and 13 for $\alpha_2$, and Figs. 14 and 15 for $\alpha_3$. The upper and lower contours correspond to a variation of $\pm 1\sigma$ of the associated $q_j$ parameter. The $\alpha_j(t = 0)$ components obtained when training the color law with only data at maximum or with the full time series are almost identical: all conclusions from Sect. 3.3.3 remain valid and we focus here on the temporal behavior.

Of the three factors, $q_1$ has the strongest impact: this can be seen from the time series presented in Fig. 10 and on the broad band light curves presented in Fig. 11. This is associated with a stretch effect, visible in $\hat{U}$, $\hat{B}$, and $\hat{V}$ bands as an enlargement of the variation band as we move away from maximum. Between −12 days and +30 days, the $q_1$ factor has a strong influence on spectroscopic details. The brightest SNe Ia exhibit shallower troughs in their spectral features: this is especially visible for Si II $\lambda$4131, Si II $\lambda$5972, and Ca II IR, but also valid for most of the lines. This effect is more pronounced for SNe Ia at earlier phases, and fades over time. After +30 days, the factor $q_1$ is less sensitive to localized features, and shows a relative enhancement of the optical $\hat{B}$ to $\hat{R}$ bands with respect to $\hat{U}$ and $\hat{I}$. Interestingly, between +12 days and +24 days, the brightest SNe Ia are bluer than at maximum or at later phases. This is true in absolute value, but also relative to the average SN Ia. This shows that the SNe Ia color space is driven by a variable intrinsic component in addition to reddening by dust. This intrinsic color variation is linked with the position of the second peak in $\hat{I}$, which appears later for brighter SNe Ia.
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Average spectral time series and the effect of a variation of $q_1$ by $\pm 1\sigma$. Each phase is separated by a constant magnitude offset.</figcaption></figure>

As at maximum, the effect of the $q_2$ factor for the whole time series (Fig. 12) is mainly localized around specific spectral features, with little impact on broad band light curves (Fig. 13), with the exception of $\hat{I}$-band. This factor associates higher line velocities around maximum with deeper absorption troughs, strongly visible before maximum and up to +18 days in the Ca II H&K region and at all phases in the Ca II IR region. The net effect described by $q_2$ is that higher velocities are associated with slightly dimmer SNe Ia. At later phases velocity effects are still observable and at all phases after -6 days, high-velocity SNe Ia have localized differences in the 4700-5100 Å spectral region. In addition, they are bluer in $\hat{V}$, $\hat{R}$, and $\hat{I}$ after maximum, with a maximal effect on the $\hat{I}$-band at around +12 days. A slight stretch effect is visible in $\hat{U}$, and $\hat{B}$, and also results in a later phase for the second maximum in $\hat{I}$ for the brightest SNe Ia in $\hat{B}$, as was mentioned for $q_1$. The shape of variations of the light curve in $\hat{I}$ are very different for each factor, indicating that this band can help reconstruct the $q_2$ factor when using photometric data only. Furthermore, while the $q_2$ factor has a small impact on the variations of individual light curves, it has a sizable influence on colors, especially $\hat{R}-\hat{I}$, $\hat{U}-\hat{B}$, and after +25 days, on $\hat{B}-\hat{V}$. In addition to the $\hat{I}$-band, the different color variation pattern in $\hat{B}-\hat{V}$ at late phases offers a way of disentangling the $q_1$ and $q_2$ factors when dealing with photometric data.

As seen in Fig. 14, the influence of $q_3$ is minimal around maximum light, similar to the results from Sect. 3.3. It grows at other phases and appears as a stretch effect on the light curves presented in Fig. 15. Unlike $q_1$, the stretch described by $q_3$ shows little correlation with the magnitude at maximum. While $q_3$ has a larger impact on the light curves than $q_2$, the reverse is true for individual colors. At late phases, $q_3$ exhibits a brighter-redder correlation that might be employed to distinguish this vector from $q_1$ when using photometric data only.

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Average light curves and color evolution with phase in synthetic bands, and the predicted effect of a variation of $q_1$ by $\pm 1\sigma$.</figcaption></figure>

The influence of $q_3$ on spectral structures is mainly visible in calcium regions, although variations can also be observed in other regions. Around Ca II IR, between -6 days and maximum, brighter SNe Ia exhibit a deeper absorption trough at high calcium velocity. This high-velocity calcium feature is visible for all $q_3$ values at -12 days, but fades faster for dimmer SNe Ia. Meanwhile the lower velocity counterpart is only visible in $q_3$ for the dimmest supernovae at -12 days and appears later for brighter SN Ia.

# 5. SUGAR performance

## 5.1. Fitting spectra with SUGAR

The measurement of spectral indicators is very sensitive to the signal-to-noise ratio of the spectrum, and at high redshift, measuring the indicators in the observer's reference frame is not possible. However, once the model has been trained, the parameters $q_1$, $q_2$, $q_3$, $A_V$, and $\Delta M_{\text{gray}}$ for a given SN Ia can be directly estimated from its spectral time series. Estimating the parameters $q_1$, $q_2$, $q_3$, $A_V$, and $\Delta M_{\text{gray}}$ is done directly on spectra by minimizing the following $\chi_i^2$:

$$ \chi_i^2 = (\mathbf{M}'_i - \mathbf{A}'\mathbf{x}_i)^T \mathbf{W}\mathbf{M}'_i (\mathbf{M}'_i - \mathbf{A}'\mathbf{x}_i), \quad (29) $$

where $\mathbf{M}'_i$ is the vector containing the entire spectral time series of SN Ia $i$, ordered in the same way as for the SUGAR training:

$$ \mathbf{M}'_i = (\mathbf{y}_{3340\AA} \cdots \mathbf{y}_{\lambda_i} \cdots \mathbf{y}_{8580\AA})_i. \quad (30) $$

where $\mathbf{y}_{\lambda_i}$ corresponds to the light curve at wavelength $\lambda$ for the observed phases. The model is projected onto observation phases using cubic spline interpolation. For notation efficiency, it is thus ordered as:

$$ \mathbf{A}' = (\mathbf{M}'_0, 1, f(\lambda, R_V), \alpha_{1'}, \alpha_{2'}, \alpha_{3'}). \quad (31) $$

The vector $\mathbf{x}_i$ contains the parameters that describe the SN Ia $i$ according to the SUGAR model:

$$ \mathbf{x}_i^T = (1, \Delta M_{\text{gray}} i, A_V i, q_{i,1}, q_{i,2}, q_{i,3}). \quad (32) $$
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Average spectral time series and the effect of a variation of $q_1$ by $\pm 1\sigma$. Each phase is separated by a constant magnitude offset.</figcaption></figure>

As at maximum, the effect of the $q_2$ factor for the whole time series (Fig. 12) is mainly localized around specific spectral features, with little impact on broad band light curves (Fig. 13), with the exception of $\hat{I}$-band. This factor associates higher line velocities around maximum with deeper absorption troughs, strongly visible before maximum and up to +18 days in the Ca II H&K region and at all phases in the Ca II IR region. The net effect described by $q_2$ is that higher velocities are associated with slightly dimmer SNe Ia. At later phases velocity effects are still observable and at all phases after -6 days, high-velocity SNe Ia have localized differences in the 4700 - 5100 Å spectral region. In addition, they are bluer in $\hat{V}$, $\hat{R}$, and $\hat{I}$ after maximum, with a maximal effect on the $\hat{I}$-band at around +12 days. A slight stretch effect is visible in $\hat{U}$, and $\hat{B}$, and also results in a later phase for the second maximum in $\hat{I}$ for the brightest SNe Ia in $\hat{B}$, as was mentioned for $q_1$. The shape of variations of the light curve in $\hat{I}$ are very different for each factor, indicating that this band can help reconstruct the $q_2$ factor when using photometric data only. Furthermore, while the $q_2$ factor has a small impact on the variations of individual light curves, it has a sizable influence on colors, especially $\hat{R} - \hat{I}$, $\hat{U} - \hat{B}$, and after +25 days, on $\hat{B} - \hat{V}$. In addition to the $\hat{I}$-band, the different color variation pattern in $\hat{B} - \hat{V}$ at late phases offers a way of disentangling the $q_1$ and $q_2$ factors when dealing with photometric data.

As seen in Fig. 14, the influence of $q_3$ is minimal around maximum light, similar to the results from Sect. 3.3. It grows at other phases and appears as a stretch effect on the light curves presented in Fig. 15. Unlike $q_1$, the stretch described by $q_3$ shows little correlation with the magnitude at maximum. While $q_3$ has a larger impact on the light curves than $q_2$, the reverse is true for individual colors. At late phases, $q_3$ exhibits a brighter-redder correlation that might be employed to distinguish this vector from $q_1$ when using photometric data only.

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Average light curves and color evolution with phase in synthetic bands, and the predicted effect of a variation of $q_1$ by $\pm 1\sigma$.</figcaption></figure>

The influence of $q_3$ on spectral structures is mainly visible in calcium regions, although variations can also be observed in other regions. Around Ca II IR, between -6 days and maximum, brighter SNe Ia exhibit a deeper absorption trough at high calcium velocity. This high-velocity calcium feature is visible for all $q_3$ values at -12 days, but fades faster for dimmer SNe Ia. Meanwhile the lower velocity counterpart is only visible in $q_3$ for the dimmest supernovae at -12 days and appears later for brighter SN Ia.

# 5. SUGAR performance

## 5.1. Fitting spectra with SUGAR

The measurement of spectral indicators is very sensitive to the signal-to-noise ratio of the spectrum, and at high redshift, measuring the indicators in the observer's reference frame is not possible. However, once the model has been trained, the parameters $q_1$, $q_2$, $q_3$, $A_V$, and $\Delta M_{\text{gray}}$ for a given SN Ia can be directly estimated from its spectral time series. Estimating the parameters $q_1$, $q_2$, $q_3$, $A_V$, and $\Delta M_{\text{gray}}$ is done directly on spectra by minimizing the following $\chi_i^2$:

$$ \chi_i^2 = (\mathbf{M}'_i - \mathbf{A}'\mathbf{x}_i)^T \mathbf{W}\mathbf{M}'_i (\mathbf{M}'_i - \mathbf{A}'\mathbf{x}_i), \quad (29) $$

where $\mathbf{M}'_i$ is the vector containing the entire spectral time series of SN Ia $i$, ordered in the same way as for the SUGAR training:

$$ \mathbf{M}'_i = (\mathbf{y}_{3340\AA} \cdots \mathbf{y}_{\lambda_i} \cdots \mathbf{y}_{8580\AA})_i. \quad (30) $$

where $\mathbf{y}_{\lambda_i}$ corresponds to the light curve at wavelength $\lambda$ for the observed phases. The model is projected onto observation phases using cubic spline interpolation. For notation efficiency, it is thus ordered as:

$$ \mathbf{A}' = (\mathbf{M}'_0, 1, f(\lambda, R_V), \alpha_{1'}, \alpha_{2'}, \alpha_{3'}). \quad (31) $$

The vector $\mathbf{x}_i$ contains the parameters that describe the SN Ia $i$ according to the SUGAR model:

$$ \mathbf{x}_i^T = (1, \Delta M_{\text{gray}} i, A_V i, q_{i,1}, q_{i,2}, q_{i,3}). \quad (32) $$
<figure><img src="image_13.png" /><figcaption>Fig. 12. Average spectral time series and the effect of a variation of $q_2$ by $\pm 1\sigma$. Each phase is separated by a constant magnitude offset.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Average light curves and color evolution with phase in synthetic bands, and the predicted effect of a variation of $q_2$ by $\pm 1\sigma$.</figcaption></figure>

Finally, $W_{M_i'}$ is the weight matrix that comes from the $M_i'$ errors. Estimating $x_i$ amounts to minimizing the $\chi_i^2$ with respect to $x_i$, giving the solution:

$$x_i = (A'^T W_{M_i'} A')^{-1} (A'^T W_{M_i'} M_i') \quad (33)$$

Therefore the covariance on the $x_i$ is given by:

$$\operatorname{cov}(x_i) = (A'^T W_{M_i'} A')^{-1}. \quad (34)$$

<figure><img src="image_15.png" /><figcaption>Fig. 14. Average spectral time series and the effect of a variation of $q_3$ by $\pm 1\sigma$. Each phase is separated by a constant magnitude offset.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Average light curves and color evolution with phase in synthetic bands, and the predicted effect of a variation of $q_3$ by $\pm 1\sigma$.</figcaption></figure>

The term $x_i$ is computed for the 105 SNe Ia that were used to train SUGAR. The values obtained are compared with the factors and $A_V$ measured at maximum light in Fig. 16. Both quantities are highly correlated. However, this correlation decreases with decreasing order of the factor index. This is because the factor uncertainties increase with index due to the noise being greater for higher-order components. The increase in noise therefore reduces the correlation. The extinction term is quite
<figure><img src="image_13.png" /><figcaption>Fig. 12. Average spectral time series and the effect of a variation of $q_2$ by $\pm 1\sigma$. Each phase is separated by a constant magnitude offset.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Average light curves and color evolution with phase in synthetic bands, and the predicted effect of a variation of $q_2$ by $\pm 1\sigma$.</figcaption></figure>

Finally, $W_{M_i}$ is the weight matrix that comes from the $M_i'$ errors. Estimating $x_i$ amounts to minimizing the $\chi_i^2$ with respect to $x_i$, giving the solution:

$$x_i = (A'^T W_{M_i} A')^{-1} (A'^T W_{M_i} M_i') \quad (33)$$

Therefore the covariance on the $x_i$ is given by:

$$\operatorname{cov}(x_i) = (A'^T W_{M_i} A')^{-1}. \quad (34)$$

<figure><img src="image_15.png" /><figcaption>Fig. 14. Average spectral time series and the effect of a variation of $q_3$ by $\pm 1\sigma$. Each phase is separated by a constant magnitude offset.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Average light curves and color evolution with phase in synthetic bands, and the predicted effect of a variation of $q_3$ by $\pm 1\sigma$.</figcaption></figure>

The term $x_i$ is computed for the 105 SNe Ia that were used to train SUGAR. The values obtained are compared with the factors and $A_V$ measured at maximum light in Fig. 16. Both quantities are highly correlated. However, this correlation decreases with decreasing order of the factor index. This is because the factor uncertainties increase with index due to the noise being greater for higher-order components. The increase in noise therefore reduces the correlation. The extinction term is quite
<figure><img src="image_17.png" /><figcaption>Fig. 16. Factors estimated directly from the spectral indicators, as a function of the factors estimated only from the SUGAR model and the $A_V$ estimated at maximum as a function of those estimated with SUGAR. The red lines represent the equation $x = y$.</figcaption></figure>

compatible between its maximum estimate and its estimation with SUGAR.

The term $x_i$ is also computed for the 58 SNe Ia kept for validation. The distribution of the $x_i$ parameters from the validation sample is compared to the $x_i$ parameters from the training in Fig. 17. From this, one can see that the training and validation samples are compatible, as it is confirmed in following sections.

## 5.2. Comparison with SALT2

### 5.2.1. Difference between SUGAR and SALT2 models

Throughout this article, we use the SALT2 model as a benchmark to evaluate the properties of the SUGAR model, because it has become a reference for cosmological analyses. SALT2 describes the SED of SNe Ia in the following functional form (Guy et al. 2007):

$$F(t, \lambda) = X_0[F_0(t, \lambda) + X_1F_1(t, \lambda)] \exp[C C_L(\lambda)], \quad (35)$$

where the observed flux $F$ depends on the individual SN Ia parameters ($X_0, X_1, C$) and on the model average flux $F_0$ (defined up to a multiplicative factor), the model flux variation $F_1$, and an empirical color law $C_L$. In the regime where $X_1F_1$ is small with respect to $F_0$, this model, once expressed in magnitudes, can be identified with the SUGAR model: ($F_0, F_1, C_L$) and ($X_1, C$) take the roles of ($M_0, \alpha_1, f(R_V)$), and ($q_1, A_V$), respectively. $\Delta M_{\text{gray}}$ is the equivalent of a linear combination of ($\log X_0, X_1, C$).

Although the major differences between the models come from the inclusion of two additional intrinsic components in SUGAR, its color law, and the data on which they were trained, there are many other differing characteristics, which are summarized Table 1. The inclusion of two additional intrinsic components is enabled by the level of detail available in our spectrophotometric time series. We therefore expect the SUGAR model to provide a more faithful representation of the spectral details than SALT2. Though analogous in spirit, $X_1$ and $q_1$ are derived employing a quite different paradigm: $F_1$ comes from a PCA analysis performed in flux space, while the $\alpha_i$ are trained in magnitude space and driven to reproduce the prominent spectral indicators. A direct generalization of the SALT2 approach with more components is presented in Saunders et al. (2018). While the empirical SALT2 $C_L$ function models the average color variation, our color curve was shown to be accurately described by a Cardelli et al. (1989) extinction law. The SALT2 color parameter $c$ thus mixes information about reddening by interstellar media (described by $A_V$) with intrinsic SNe Ia properties that were not directly matched with one of the $q_i$ parameters. These differences in the treatment of colors would manifest in long-range spectral variation (because the extinction curve is smooth). In regions where the phase coverage is sufficient, the different phase interpolation methods used are not expected to contribute significantly to differences between models. Otherwise, and for early phases in particular, we expect the differences to be driven by the training samples rather than interpolation technique.

The wavelength and phase coverage of SUGAR are restricted compared to SALT2. This is directly related to the training data: SUGAR is trained on low redshift spectrophotometric time series, while SALT2 is trained on a wide redshift range of photometric data, employing a handful of spectra to help refine the SED details. As SALT2 uses low and high redshift data it offers a better coverage of the UV domain. It also has a larger phase coverage thanks to the rolling cadence strategy, while SNfactory SNe Ia have to be specifically targeted. While SALT2 could be easily trained with SNfactory data added to their sample, the reverse is not true: major adaptation would be needed to allow SUGAR to incorporate photometric data in its training. However, this does not prevent fitting photometric light curves with SUGAR in the same way as with SALT2, provided they fall within the phase and spectral coverage of SUGAR. In contrast, SALT2 was not designed to fit spectrophotometric data.

### 5.2.2. Spectral residuals

The SUGAR model is designed to improve the spectral description of SNe Ia by adding components beyond stretch and color. While the ideal model would perfectly match the data, the level of mismatch can be quantified by the residuals, that is the difference between the left and right term of Eq. (22). The improvement with respect to SALT2 is then studied by comparing the data to the spectral reconstruction predicted by both models. One example of these comparisons is presented Fig. 18, which shows the spectra and residuals for SN-training-4. For this supernova, SUGAR gives a better description than SALT2. The
<figure><img src="image_17.png" /><figcaption>Fig. 16. Factors estimated directly from the spectral indicators, as a function of the factors estimated only from the SUGAR model and the $A_V$ estimated at maximum as a function of those estimated with SUGAR. The red lines represent the equation $x = y$.</figcaption></figure>

compatible between its maximum estimate and its estimation with SUGAR.

The term $x_i$ is also computed for the 58 SNe Ia kept for validation. The distribution of the $x_i$ parameters from the validation sample is compared to the $x_i$ parameters from the training in Fig. 17. From this, one can see that the training and validation samples are compatible, as it is confirmed in following sections.

## 5.2. Comparison with SALT2

### 5.2.1. Difference between SUGAR and SALT2 models

Throughout this article, we use the SALT2 model as a benchmark to evaluate the properties of the SUGAR model, because it has become a reference for cosmological analyses. SALT2 describes the SED of SNe Ia in the following functional form (Guy et al. 2007):

$$F(t, \lambda) = X_0[F_0(t, \lambda) + X_1F_1(t, \lambda)] \exp[C C_L(\lambda)], \quad (35)$$

where the observed flux $F$ depends on the individual SN Ia parameters ($X_0, X_1, C$) and on the model average flux $F_0$ (defined up to a multiplicative factor), the model flux variation $F_1$, and an empirical color law $C_L$. In the regime where $X_1F_1$ is small with respect to $F_0$, this model, once expressed in magnitudes, can be identified with the SUGAR model: ($F_0, F_1, C_L$) and ($X_1, C$) take the roles of ($M_0, \alpha_1, f(R_V)$), and ($q_1, A_V$), respectively. $\Delta M_{\text{gray}}$ is the equivalent of a linear combination of ($\log X_0, X_1, C$).

Although the major differences between the models come from the inclusion of two additional intrinsic components in SUGAR, its color law, and the data on which they were trained, there are many other differing characteristics, which are summarized Table 1. The inclusion of two additional intrinsic components is enabled by the level of detail available in our spectrophotometric time series. We therefore expect the SUGAR model to provide a more faithful representation of the spectral details than SALT2. Though analogous in spirit, $X_1$ and $q_1$ are derived employing a quite different paradigm: $F_1$ comes from a PCA analysis performed in flux space, while the $\alpha_i$ are trained in magnitude space and driven to reproduce the prominent spectral indicators. A direct generalization of the SALT2 approach with more components is presented in Saunders et al. (2018). While the empirical SALT2 $C_L$ function models the average color variation, our color curve was shown to be accurately described by a Cardelli et al. (1989) extinction law. The SALT2 color parameter $c$ thus mixes information about reddening by interstellar media (described by $A_V$) with intrinsic SNe Ia properties that were not directly matched with one of the $q_i$ parameters. These differences in the treatment of colors would manifest in long-range spectral variation (because the extinction curve is smooth). In regions where the phase coverage is sufficient, the different phase interpolation methods used are not expected to contribute significantly to differences between models. Otherwise, and for early phases in particular, we expect the differences to be driven by the training samples rather than interpolation technique.

The wavelength and phase coverage of SUGAR are restricted compared to SALT2. This is directly related to the training data: SUGAR is trained on low redshift spectrophotometric time series, while SALT2 is trained on a wide redshift range of photometric data, employing a handful of spectra to help refine the SED details. As SALT2 uses low and high redshift data it offers a better coverage of the UV domain. It also has a larger phase coverage thanks to the rolling cadence strategy, while SNfactory SNe Ia have to be specifically targeted. While SALT2 could be easily trained with SNfactory data added to their sample, the reverse is not true: major adaptation would be needed to allow SUGAR to incorporate photometric data in its training. However, this does not prevent fitting photometric light curves with SUGAR in the same way as with SALT2, provided they fall within the phase and spectral coverage of SUGAR. In contrast, SALT2 was not designed to fit spectrophotometric data.

### 5.2.2. Spectral residuals

The SUGAR model is designed to improve the spectral description of SNe Ia by adding components beyond stretch and color. While the ideal model would perfectly match the data, the level of mismatch can be quantified by the residuals, that is the difference between the left and right term of Eq. (22). The improvement with respect to SALT2 is then studied by comparing the data to the spectral reconstruction predicted by both models. One example of these comparisons is presented Fig. 18, which shows the spectra and residuals for SN-training-4. For this supernova, SUGAR gives a better description than SALT2. The
<figure><img src="image_18.png" /><figcaption>**Fig. 17.** SUGAR parameters fit directly from spectral data for the 105 SNe Ia used to train SUGAR (blue) and for the 58 SNe Ia from the validation sample (red).</figcaption></figure>

<table><caption>**Table 1.** Main differences between the SALT2.4 model and the SUGAR model.</caption><thead><tr><th>Model properties</th><th>SALT2</th><th>SUGAR</th></tr></thead><tbody><tr><td>Rest frame wavelength coverage</td><td>2000–9200 Å</td><td>3254–8649 Å</td></tr><tr><td>Rest frame phase coverage</td><td>−20 to +50 days</td><td>−12 to +48 days</td></tr><tr><td>Training data</td><td>Spectroscopy and photometry</td><td>Spectrophotometric time series</td></tr><tr><td>Redshift range</td><td>Low and high redshift</td><td>Low redshift</td></tr><tr><td>Interpolation</td><td>Cubic spline</td><td>Gaussian process with squared exponential kernel</td></tr><tr><td>Intrinsic properties from</td><td>Weighted PCA like on SED</td><td>EM-FA on spectral features</td></tr><tr><td>Number of intrinsic component</td><td>1</td><td>3</td></tr><tr><td>Color</td><td>Color law (third order polynomial)</td><td>Cardelli et al. (1989) law</td></tr><tr><td>Error on the model</td><td>Yes</td><td>No</td></tr></tbody></table>

improvement is especially clear in the UV after +17 days and in the infrared at any phase, and particularly in the OI λ7773 and Ca II IR regions. This is not surprising insofar as SALT2 is essentially trained to reproduce the bluer part of the spectrum. The SUGAR description is also better in the regions around lines, even if this description is not totally satisfactory in some areas: Ca II H&K, Si II λ6355, and the bluer part of the spectrum at late phases. In the B-band, the reconstruction of the spectral details is quite similar between SUGAR and SALT2, although SUGAR seems to follow more faithfully the spectral details, as can be seen in the Mg II area.

SN-training-4 is only one example among the 105 SNe Ia on which SUGAR was trained. To have a statistical description of the precision of the models, we derive the dispersion of the residuals obtained for all SNe Ia as a function of wavelength and

phase. In Fig. 19, we present the wRMS of this dispersion: for each wavelength, we compute the wRMS across all phases and weight it by the spectral variance. To show how the accuracy of the model evolves when going from a classic Tripp (1998) relation (when only $q_1$ and $A_V$ are included) to the full model, the wRMS for SUGAR is computed three times with progressive inclusion of $q_1$, $q_2$, and $q_3$. We compute the residual dispersion for both the training and the validation sample. What was observed for SN-training-4 is confirmed in this statistical analysis across all SNe Ia: at any wavelength, SUGAR gives a better description than SALT2. The strongest improvement is in the IR, but is also significant in the UV. All the spectral details are significantly improved, such as the Ca II H&K, Si II λ6355, the OI λ7773, and the Ca II IR. The best performance of the SUGAR model is obtained in the region around 6500 Å, where
<figure><img src="image_18.png" /><figcaption>**Fig. 17.** SUGAR parameters fit directly from spectral data for the 105 SNe Ia used to train SUGAR (blue) and for the 58 SNe Ia from the validation sample (red).</figcaption></figure>

<table><caption>**Table 1.** Main differences between the SALT2.4 model and the SUGAR model.</caption><thead><tr><th>Model properties</th><th>SALT2</th><th>SUGAR</th></tr></thead><tbody><tr><td>Rest frame wavelength coverage</td><td>2000–9200 Å</td><td>3254–8649 Å</td></tr><tr><td>Rest frame phase coverage</td><td>−20 to +50 days</td><td>−12 to +48 days</td></tr><tr><td>Training data</td><td>Spectroscopy and photometry</td><td>Spectrophotometric time series</td></tr><tr><td>Redshift range</td><td>Low and high redshift</td><td>Low redshift</td></tr><tr><td>Interpolation</td><td>Cubic spline</td><td>Gaussian process with squared exponential kernel</td></tr><tr><td>Intrinsic properties from</td><td>Weighted PCA like on SED</td><td>EM-FA on spectral features</td></tr><tr><td>Number of intrinsic component</td><td>1</td><td>3</td></tr><tr><td>Color</td><td>Color law (third order polynomial)</td><td>Cardelli et al. (1989) law</td></tr><tr><td>Error on the model</td><td>Yes</td><td>No</td></tr></tbody></table>

improvement is especially clear in the UV after +17 days and in the infrared at any phase, and particularly in the OI λ7773 and Ca II IR regions. This is not surprising insofar as SALT2 is essentially trained to reproduce the bluer part of the spectrum. The SUGAR description is also better in the regions around lines, even if this description is not totally satisfactory in some areas: Ca II H&K, Si II λ6355, and the bluer part of the spectrum at late phases. In the B-band, the reconstruction of the spectral details is quite similar between SUGAR and SALT2, although SUGAR seems to follow more faithfully the spectral details, as can be seen in the Mg II area.

SN-training-4 is only one example among the 105 SNe Ia on which SUGAR was trained. To have a statistical description of the precision of the models, we derive the dispersion of the residuals obtained for all SNe Ia as a function of wavelength and

phase. In Fig. 19, we present the wRMS of this dispersion: for each wavelength, we compute the wRMS across all phases and weight it by the spectral variance. To show how the accuracy of the model evolves when going from a classic Tripp (1998) relation (when only $q_1$ and $A_V$ are included) to the full model, the wRMS for SUGAR is computed three times with progressive inclusion of $q_1$, $q_2$, and $q_3$. We compute the residual dispersion for both the training and the validation sample. What was observed for SN-training-4 is confirmed in this statistical analysis across all SNe Ia: at any wavelength, SUGAR gives a better description than SALT2. The strongest improvement is in the IR, but is also significant in the UV. All the spectral details are significantly improved, such as the Ca II H&K, Si II λ6355, the OI λ7773, and the Ca II IR. The best performance of the SUGAR model is obtained in the region around 6500 Å, where
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** SUGAR and SALT2 fit of SN-training-4. **Left:** observed spectral series of SN-training-4 (in black), compared to SALT2 (in red) and SUGAR (blue). **Right:** residuals of SALT2 (red) and SUGAR (blue) models to the observations. The black line represents the zero value of the residuals.</figcaption></figure>

the dispersion gets as low as 0.05 mag. It should be kept in mind that the spectral dispersion includes noise variance that may be high in some spectra and is therefore not representative of the dispersion expected in photometry after integration over broad band filters. The phase evolution of the wRMS is presented in Fig. 20: the dispersion in each phase bin is calculated across all wavelengths. This includes wavelengths for which SALT2 is not adapted: this explains relatively high values of wRMS for this model compared to SUGAR.

The successive inclusion of $q_1, q_2$, and $q_3$ can help us understand the origin of SUGAR's improvement over SALT2. As the training samples for both models are different, one may wonder what SALT2 would have given if trained on SNfactory data. This situation can be studied when including only $q_1$: SUGAR$_{|q_1}$ is slightly worse than SALT2 between Si II λ4131 and Si II λ6355 features, and between Si II λ6355 and O I λ7773. These regions cover the wavelength range where SALT2 was trained and optimized, while SUGAR$_{|q_1}$ is only the by-product of the three component model and is not by itself optimal for a single component. Earlier phase are more difficult to reproduce in both models: this could be an effect of higher variability beyond stretch at those phases. For later phases, SUGAR$_{|q_1}$ does not exhibit the same increase in wRMS as SALT2. This is also an effect of the lack of coverage for IR wavelength for SALT2. One could therefore expect SALT2 to behave at least as well as SUGAR$_{|q_1}$ in UV and IR when trained on SNfactory data. As expected, the inclusion of additional components improves the accuracy of SUGAR: $q_2$ improves the description of the spectral features because it is correlated with (Ca II H&K, Si II λ6355, and Ca II IR) while one has to wait for the inclusion of $q_3$ to reach the full SUGAR precision.

The improvement brought by SUGAR is thus enabled by two factors: the spectral coverage and the spectral resolution that allows the inclusion of the new components, $q_2$, which improves the color description, and $q_3$, which helps reduce the dispersion away from maximum light. Moreover, the improvement added by $q_2$ and $q_3$ are also seen in the validation sample, which confirms that these components are not due to overtraining.

### 5.2.3. Relation between SALT2 and SUGAR parameters

In Sect. 3.2 we presented the correlations between the SALT2 parameters and the first five factors found using factor analysis.

After training the full SUGAR model and reconstructing the parameters of individual supernovae as described in Sect. 5.1, we can study the correlations between the SALT2 parameters (X₀, X₁, C) and the SUGAR parameters (ΔM<sub>gray</sub>, q₁, q₂, q₃, A<sub>V</sub>). The direct comparison of X₀ and ΔM<sub>gray</sub> is irrelevant due to the different prescription for lifting the degeneracies of both models. We recognize ΔM<sub>gray</sub> as the Hubble diagram residuals, a quantity that we can also compute from SALT2 parameters (and denote as Δμ<sup>corr</sup>). The results of this comparison between SALT2 and SUGAR parameters are presented in Fig. 21. The lessons from the correlations of SALT2 parameters and the original factors still hold: X₁ is strongly correlated with q₁ and q₃; q₂ has no significant correlation with any of SALT2 parameters, except with Δμ. The small correlation of q₂ and q₃ with Δμ shows that the inclusion of this factor can help improve the standardization. The SALT2 color C is strongly correlated with A<sub>V</sub>, as would be expected since the purpose of these two parameters is to take into account the reddening of SNe Ia. Contrary to the results from Sect. 3.2 with the factors alone, C now exhibits a 4.5σ correlation with q₃. This correlation might come from the slight fluctuation in color with respect to phase visible in Fig. 15 and the inclusion of the validation sample. The slight correlation of C and ΔM<sub>gray</sub> comes from the correlation of ΔM<sub>gray</sub> and A<sub>V</sub>. Indeed, ΔM<sub>gray</sub> and A<sub>V</sub> were determined separately during the training, without fixing potential degeneracy between them. This explains why they are correlated. Finally, Δμ has a strong correlation with ΔM<sub>gray</sub>. This indicates that the inclusion of the additional factors in SUGAR is unable to catch the major source of magnitude variability of SNe Ia, even if it provides a much improved description of spectral details.

## 6. Discussion

### 6.1. Gray dispersion and dispersion matrix fitting

For the model at maximum light, we did not fit a gray term but did include a dispersion matrix while estimating the extinction curve. However, for the full spectral time series, we fit for a gray offset and used the extinction curve determined at maximum light. One might wonder why we did not do a simultaneous fit of the full spectral time series, the gray offset and the extinction curve. Our decision was motivated by the following logic.
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** SUGAR and SALT2 fit of SN-training-4. **Left:** observed spectral series of SN-training-4 (in black), compared to SALT2 (in red) and SUGAR (blue). **Right:** residuals of SALT2 (red) and SUGAR (blue) models to the observations. The black line represents the zero value of the residuals.</figcaption></figure>

the dispersion gets as low as 0.05 mag. It should be kept in mind that the spectral dispersion includes noise variance that may be high in some spectra and is therefore not representative of the dispersion expected in photometry after integration over broad band filters. The phase evolution of the wRMS is presented in Fig. 20: the dispersion in each phase bin is calculated across all wavelengths. This includes wavelengths for which SALT2 is not adapted: this explains relatively high values of wRMS for this model compared to SUGAR.

The successive inclusion of $q_1, q_2$, and $q_3$ can help us understand the origin of SUGAR's improvement over SALT2. As the training samples for both models are different, one may wonder what SALT2 would have given if trained on SNfactory data. This situation can be studied when including only $q_1$: SUGAR$_{|q_1}$ is slightly worse than SALT2 between Si II λ4131 and Si II λ6355 features, and between Si II λ6355 and O I λ7773. These regions cover the wavelength range where SALT2 was trained and optimized, while SUGAR$_{|q_1}$ is only the by-product of the three component model and is not by itself optimal for a single component. Earlier phase are more difficult to reproduce in both models: this could be an effect of higher variability beyond stretch at those phases. For later phases, SUGAR$_{|q_1}$ does not exhibit the same increase in wRMS as SALT2. This is also an effect of the lack of coverage for IR wavelength for SALT2. One could therefore expect SALT2 to behave at least as well as SUGAR$_{|q_1}$ in UV and IR when trained on SNfactory data. As expected, the inclusion of additional components improves the accuracy of SUGAR: $q_2$ improves the description of the spectral features because it is correlated with (Ca II H&K, Si II λ6355, and Ca II IR) while one has to wait for the inclusion of $q_3$ to reach the full SUGAR precision.

The improvement brought by SUGAR is thus enabled by two factors: the spectral coverage and the spectral resolution that allows the inclusion of the new components, $q_2$, which improves the color description, and $q_3$, which helps reduce the dispersion away from maximum light. Moreover, the improvement added by $q_2$ and $q_3$ are also seen in the validation sample, which confirms that these components are not due to overtraining.

### 5.2.3. Relation between SALT2 and SUGAR parameters

In Sect. 3.2 we presented the correlations between the SALT2 parameters and the first five factors found using factor analysis.

After training the full SUGAR model and reconstructing the parameters of individual supernovae as described in Sect. 5.1, we can study the correlations between the SALT2 parameters (X₀, X₁, C) and the SUGAR parameters (ΔM<sub>gray</sub>, q₁, q₂, q₃, A<sub>V</sub>). The direct comparison of X₀ and ΔM<sub>gray</sub> is irrelevant due to the different prescription for lifting the degeneracies of both models. We recognize ΔM<sub>gray</sub> as the Hubble diagram residuals, a quantity that we can also compute from SALT2 parameters (and denote as Δμ<sup>corr</sup>). The results of this comparison between SALT2 and SUGAR parameters are presented in Fig. 21. The lessons from the correlations of SALT2 parameters and the original factors still hold: X₁ is strongly correlated with q₁ and q₃; q₂ has no significant correlation with any of SALT2 parameters, except with Δμ. The small correlation of q₂ and q₃ with Δμ shows that the inclusion of this factor can help improve the standardization. The SALT2 color C is strongly correlated with A<sub>V</sub>, as would be expected since the purpose of these two parameters is to take into account the reddening of SNe Ia. Contrary to the results from Sect. 3.2 with the factors alone, C now exhibits a 4.5σ correlation with q₃. This correlation might come from the slight fluctuation in color with respect to phase visible in Fig. 15 and the inclusion of the validation sample. The slight correlation of C and ΔM<sub>gray</sub> comes from the correlation of ΔM<sub>gray</sub> and A<sub>V</sub>. Indeed, ΔM<sub>gray</sub> and A<sub>V</sub> were determined separately during the training, without fixing potential degeneracy between them. This explains why they are correlated. Finally, Δμ has a strong correlation with ΔM<sub>gray</sub>. This indicates that the inclusion of the additional factors in SUGAR is unable to catch the major source of magnitude variability of SNe Ia, even if it provides a much improved description of spectral details.

## 6. Discussion

### 6.1. Gray dispersion and dispersion matrix fitting

For the model at maximum light, we did not fit a gray term but did include a dispersion matrix while estimating the extinction curve. However, for the full spectral time series, we fit for a gray offset and used the extinction curve determined at maximum light. One might wonder why we did not do a simultaneous fit of the full spectral time series, the gray offset and the extinction curve. Our decision was motivated by the following logic.
<figure><img src="image_20.png" /><figcaption>**Fig. 19.** wRMS of the residuals as a function of the wavelength for the full SUGAR model (blue line), for the SUGAR model only corrected by $q_1$ and $A_V$ (blue dashed dotted line), for the SUGAR model only corrected by $q_1$, $q_2$ and $A_V$ (blue dashed lines), and for the SALT2 model (red line). *Top*: training data. *Bottom*: validation data.</figcaption></figure>

In order to perform a simultaneous fit, one important thing is to have a good estimate of the dispersion matrix, which is essential for the estimation of the extinction curve, as discussed in Scolnic et al. (2014). In addition, this matrix would have the advantage of taking into account temporal correlations of residuals, in addition to the color correlations of the residuals, which would improve the model developed at maximum light.

However, from a numerical point of view, this would break the sparse algebra approximations that significantly speed up the training. Indeed, the sparse algebra (due to Eq. (28)) makes it possible to improve speed efficiency of minimization of Eq. (26) from $O((N_{\lambda}N_{t}N_{\mathrm{param}})^{3})$ to $O((N_{\lambda} \times (N_{t}N_{\mathrm{param}})^{3})$, where $N_{\lambda}$ is the number of wavelength bin, $N_{t}$ the number of bin in terms of time, and $N_{\mathrm{param}}$ is the number of parameters to fit for a given wavelength and epoch. In SUGAR training, $N_{\lambda} \sim 190$, $N_{t} \sim 30$, and $N_{\mathrm{param}} = 5$ (average spectrum + three factors + extinction curve), consequently, not doing sparse algebra significatively slow down the speed of the training algorithm by $O(N_{\lambda}^2)$.

Moreover, there would be problems with numerical stability for the estimation of the dispersion matrix due to its size; in the case of the simultaneous fit: it would be ~20 000 × 20 000. Indeed, it would be necessary to make sure that the dispersion matrix is positive definite, which would involve doing a Singular Value Decomposition that would again slow down the speed

of the training algorithm. In addition, it is evident that the gray offset is degenerate with the dispersion matrix because the dispersion matrix can fully capture a gray dispersion. Therefore we deliberately did not include it explicitly in the fit of the extinction curve within SUGAR in Sect. 3.3. There is however a degeneracy between $A_{\lambda 0}$ and the gray offset, and part of the latter is captured by this parameter. The SUGAR model presented here is already a significant improvement over the SALT2 model and provides insights into understanding SN Ia variability.

## 6.2. Test of adding an additional component

In the Sect. 3.2, we discussed the number of factors needed to describe the final SUGAR SED, and we concluded that this cannot be determined only from the factor description of spectral features at maximum light. Indeed, even if they are strongly related, the main goal is to know the number of components needed to describe the full SED and not the number of components needed to describe the spectral features space at maximum light. One way to check if the choice of three factors used here is optimal is to retrain the SUGAR model with more than three components and observe how the spectral residuals evolve with this change.

In the following, we ran the training of SUGAR twice, each time adding an additional component, that is we reproduced
<figure><img src="image_20.png" /><figcaption>**Fig. 19.** wRMS of the residuals as a function of the wavelength for the full SUGAR model (blue line), for the SUGAR model only corrected by $q_1$ and $A_V$ (blue dashed dotted line), for the SUGAR model only corrected by $q_1$, $q_2$ and $A_V$ (blue dashed lines), and for the SALT2 model (red line). *Top*: training data. *Bottom*: validation data.</figcaption></figure>

In order to perform a simultaneous fit, one important thing is to have a good estimate of the dispersion matrix, which is essential for the estimation of the extinction curve, as discussed in Scolnic et al. (2014). In addition, this matrix would have the advantage of taking into account temporal correlations of residuals, in addition to the color correlations of the residuals, which would improve the model developed at maximum light.

However, from a numerical point of view, this would break the sparse algebra approximations that significantly speed up the training. Indeed, the sparse algebra (due to Eq. (28)) makes it possible to improve speed efficiency of minimization of Eq. (26) from $O((N_{\lambda}N_{t}N_{\mathrm{param}})^{3})$ to $O((N_{\lambda} \times (N_{t}N_{\mathrm{param}})^{3})$, where $N_{\lambda}$ is the number of wavelength bin, $N_{t}$ the number of bin in terms of time, and $N_{\mathrm{param}}$ is the number of parameters to fit for a given wavelength and epoch. In SUGAR training, $N_{\lambda} \sim 190$, $N_{t} \sim 30$, and $N_{\mathrm{param}} = 5$ (average spectrum + three factors + extinction curve), consequently, not doing sparse algebra significatively slow down the speed of the training algorithm by $O(N_{\lambda}^2)$.

Moreover, there would be problems with numerical stability for the estimation of the dispersion matrix due to its size; in the case of the simultaneous fit: it would be ~20 000 × 20 000. Indeed, it would be necessary to make sure that the dispersion matrix is positive definite, which would involve doing a Singular Value Decomposition that would again slow down the speed

of the training algorithm. In addition, it is evident that the gray offset is degenerate with the dispersion matrix because the dispersion matrix can fully capture a gray dispersion. Therefore we deliberately did not include it explicitly in the fit of the extinction curve within SUGAR in Sect. 3.3. There is however a degeneracy between $A_{\lambda 0}$ and the gray offset, and part of the latter is captured by this parameter. The SUGAR model presented here is already a significant improvement over the SALT2 model and provides insights into understanding SN Ia variability.

## 6.2. Test of adding an additional component

In the Sect. 3.2, we discussed the number of factors needed to describe the final SUGAR SED, and we concluded that this cannot be determined only from the factor description of spectral features at maximum light. Indeed, even if they are strongly related, the main goal is to know the number of components needed to describe the full SED and not the number of components needed to describe the spectral features space at maximum light. One way to check if the choice of three factors used here is optimal is to retrain the SUGAR model with more than three components and observe how the spectral residuals evolve with this change.

In the following, we ran the training of SUGAR twice, each time adding an additional component, that is we reproduced
<figure><img src="image_21.png" /><figcaption>**Fig. 20.** wRMS of the residuals as a function of phase for the full SUGAR model (blue line), for the SUGAR model only corrected by $q_1$ and $A_V$ (blue dashed dotted line), for the SUGAR model only corrected by $q_1$, $q_2$ and $A_V$ (blue dashed lines), and for the SALT2 model (red line). Top: training data. Bottom: validation data.</figcaption></figure>

Sec. 4.2 with $q_1 - q_4$ and $q_1 - q_5$ components. In both cases the value of $R_V = 2.6$ found with three factors is fixed in order to focus only on intrinsic parameters. In Fig. 22 we compare the spectral residuals of the SUGAR model to those of SUGAR trained with the additional factors.

As expected, the addition of the two new components from factor analysis does not improve the description of SNe Ia as significantly as the addition of $q_2$ or $q_3$. The factor $q_4$ slightly improves the description within the Ca II H&K area (0.05 mag) and the Si II λ6355 (0.03 mag), but does not improve the SED outside these areas. The factor $q_5$ does not significantly improve the description of the SED. This is confirmation that our choice of using three factors provides a good description of the SED.

## 7. Conclusion

In this paper we have presented a new spectra-temporal empirical model of SNe Ia, named SUGAR. This model significantly improves the spectral description of SNe Ia compared to the current state of the art by going beyond the classical stretch and color parametrization.

In Sect. 2, we presented the SNfactory spectrophotometric dataset that was used to train the model. In Sect. 3, we presented

the intermediary data that were used to train the full SUGAR model. In a first step, we selected a set of 13 spectral indicators near maximum light in B-band that are composed of pseudo-equivalent widths and minima of P-Cygni profiles. Those spectral indicators were chosen to describe the intrinsic part of the SUGAR model because they are easy to define and are independent of host-galaxy dust extinction. Then, we defined a new basis where the 13 spectral indicators are uncorrelated. For this we developed a factor analysis algorithm that is more robust in the presence of errors than PCA algorithms. Three factors seems to be effective enough to describe the spectral indicators space at maximum light. The first factor describes the coherent variation of the pseudo-equivalent widths, mainly of the silicon and calcium lines. Like those lines, this factor is strongly correlated with stretch. The second factor is mainly correlated with the velocities and shows a very weak link with pseudo-equivalent widths, except for the Ca II H&K and S II W. The third factor shows a slight correlation with the stretch parameter. Once the factors have been defined, we established a model of SED at maximum light based on these three first factors from spectral features, using the same underlying method as in Chotard et al. (2011). This allows us to separate the intrinsic behavior from color variation due to dust. Finally, we find that the color curve obtained is compatible with a Cardelli et al. (1989) extinction curve with $R_V$ of 2.6. We then developed an interpolation method using
<figure><img src="image_21.png" /><figcaption>**Fig. 20.** wRMS of the residuals as a function of phase for the full SUGAR model (blue line), for the SUGAR model only corrected by $q_1$ and $A_V$ (blue dashed dotted line), for the SUGAR model only corrected by $q_1$, $q_2$ and $A_V$ (blue dashed lines), and for the SALT2 model (red line). Top: training data. Bottom: validation data.</figcaption></figure>

Sec. 4.2 with $q_1 - q_4$ and $q_1 - q_5$ components. In both cases the value of $R_V = 2.6$ found with three factors is fixed in order to focus only on intrinsic parameters. In Fig. 22 we compare the spectral residuals of the SUGAR model to those of SUGAR trained with the additional factors.

As expected, the addition of the two new components from factor analysis does not improve the description of SNe Ia as significantly as the addition of $q_2$ or $q_3$. The factor $q_4$ slightly improves the description within the Ca II H&K area (0.05 mag) and the Si II λ6355 (0.03 mag), but does not improve the SED outside these areas. The factor $q_5$ does not significantly improve the description of the SED. This is confirmation that our choice of using three factors provides a good description of the SED.

## 7. Conclusion

In this paper we have presented a new spectra-temporal empirical model of SNe Ia, named SUGAR. This model significantly improves the spectral description of SNe Ia compared to the current state of the art by going beyond the classical stretch and color parametrization.

In Sect. 2, we presented the SNfactory spectrophotometric dataset that was used to train the model. In Sect. 3, we presented

the intermediary data that were used to train the full SUGAR model. In a first step, we selected a set of 13 spectral indicators near maximum light in B-band that are composed of pseudo-equivalent widths and minima of P-Cygni profiles. Those spectral indicators were chosen to describe the intrinsic part of the SUGAR model because they are easy to define and are independent of host-galaxy dust extinction. Then, we defined a new basis where the 13 spectral indicators are uncorrelated. For this we developed a factor analysis algorithm that is more robust in the presence of errors than PCA algorithms. Three factors seems to be effective enough to describe the spectral indicators space at maximum light. The first factor describes the coherent variation of the pseudo-equivalent widths, mainly of the silicon and calcium lines. Like those lines, this factor is strongly correlated with stretch. The second factor is mainly correlated with the velocities and shows a very weak link with pseudo-equivalent widths, except for the Ca II H&K and S II W. The third factor shows a slight correlation with the stretch parameter. Once the factors have been defined, we established a model of SED at maximum light based on these three first factors from spectral features, using the same underlying method as in Chotard et al. (2011). This allows us to separate the intrinsic behavior from color variation due to dust. Finally, we find that the color curve obtained is compatible with a Cardelli et al. (1989) extinction curve with $R_V$ of 2.6. We then developed an interpolation method using
<figure><img src="image_22.png" /><figcaption>**Fig. 21.** SALT2 parameters estimated in photometry as a function of the SUGAR parameters estimated in spectroscopy (on training and validation sample). The color code estimates the correlation between the parameters SALT2 and SUGAR.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 22.** wRMS of the residuals as a function of the wavelength for the full SUGAR model (blue line), for the SUGAR model with the addition of the factor $q_4$ (blue dashed dotted line), and for the SUGAR model with the addition of the factor $q_4$ and $q_5$ (blue dashed lines).</figcaption></figure>

Gaussian process in order to train the SUGAR model on a uniform and fixed time grid. The correlation length obtained varies between 5 and 12 days depending on spectral regions, which justifies the use of a 3-day time step for the grid. The uncertainty from the Gaussian processes is generally underestimated by a factor of 1.3 and more investigation is needed to understand why this occurs.

The training process for the SUGAR model is described in Sect. 4, and an interpretation of each new component is provided. Both factors $q_1$ and $q_3$ resemble a stretch effect, but $q_3$ has less impact around maximum light and in color space than $q_1$. The effects of $q_2$ are strongest in the areas of spectral features,

and mainly evident in the infrared as compared to broad band photometry.

After calculating the model, we showed that, instead of going through the calculation of spectral indicators, we can work directly with the spectral time series to recover the three factors and extinction parameter. By studying model residuals as a function of wavelength, it is shown that SUGAR improves the spectral description 0.1–0.4 mag with respect to SALT2. This is valid for both training and validation data sets, which confirms that there was no overtraining resulting from the addition of $q_2$ and $q_3$. This shows that three parameters, defined at a given phase (i.e. maximum light) have predictive power at other phases.
<figure><img src="image_22.png" /><figcaption>**Fig. 21.** SALT2 parameters estimated in photometry as a function of the SUGAR parameters estimated in spectroscopy (on training and validation sample). The color code estimates the correlation between the parameters SALT2 and SUGAR.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 22.** wRMS of the residuals as a function of the wavelength for the full SUGAR model (blue line), for the SUGAR model with the addition of the factor $q_4$ (blue dashed dotted line), and for the SUGAR model with the addition of the factor $q_4$ and $q_5$ (blue dashed lines).</figcaption></figure>

Gaussian process in order to train the SUGAR model on a uniform and fixed time grid. The correlation length obtained varies between 5 and 12 days depending on spectral regions, which justifies the use of a 3-day time step for the grid. The uncertainty from the Gaussian processes is generally underestimated by a factor of 1.3 and more investigation is needed to understand why this occurs.

The training process for the SUGAR model is described in Sect. 4, and an interpretation of each new component is provided. Both factors $q_1$ and $q_3$ resemble a stretch effect, but $q_3$ has less impact around maximum light and in color space than $q_1$. The effects of $q_2$ are strongest in the areas of spectral features,

and mainly evident in the infrared as compared to broad band photometry.

After calculating the model, we showed that, instead of going through the calculation of spectral indicators, we can work directly with the spectral time series to recover the three factors and extinction parameter. By studying model residuals as a function of wavelength, it is shown that SUGAR improves the spectral description 0.1–0.4 mag with respect to SALT2. This is valid for both training and validation data sets, which confirms that there was no overtraining resulting from the addition of $q_2$ and $q_3$. This shows that three parameters, defined at a given phase (i.e. maximum light) have predictive power at other phases.