SALT3: An Improved Type Ia Supernova Model for Measuring Cosmic Distances

The specification of the model given above is degenerate; for example, the scale of the flux surfaces may be changed by reducing M₀ and M₁ and increasing x₀ by the same factor. To remove degeneracies, we apply further model definitions, which are used as constraints during the training process. These definitions are arbitrary, and have been chosen to define a "fiducial" supernova whose SED is the M₀ flux surface at the mean of the observed distributions of light-curve parameters. The definitions are:

• 1.  
  The rest-frame synthetic B-band flux of the M₀ component at peak is fixed such that ${m}_{B}^{\mathrm{peak}}=10.5$ when x₀ = 1.
• 2.  
  The rest-frame synthetic B-band flux of the M₁ component at peak is defined to be zero.
• 3.  
  The distribution of the light-curve parameter x₁ in the training sample is defined to have zero mean.
• 4.  
  The distribution of x₁ has a standard deviation of unity.
• 5.  
  The distribution of the light-curve parameter c in the training sample is defined to have zero mean.
• 6.  
  The color law is defined such that CL(4300 AA) = 0 and CL(5430 Å) = − 1, corresponding to central wavelengths for B and V bandpasses.
• 7.  
  The distributions of x₁ and c have no correlation in the training sample.

The location and scale of the x₁ distribution and the location of the c distribution inferred for any cosmological sample are thus relative to the demographics of the training sample for a particular model. However, the scale of the c distribution has been normalized by definition 6. Our model definitions and SALT2 differ only in the last definition, intended to separate the phase-independent color from the phase-dependent M₁ component. The SALT2 training code does not constrain the correlation between x₁ and c, and instead fixes the V-band flux (and B-band flux following definition 2) of the M₁ component to be 0 at peak brightness. This restricts the SALT2 training code's ability to determine the best parameterization of the SN Ia population by requiring that SALT2 x₁ has no effect on observed B − V color at peak. From a phenomenological perspective, it is easier to make inferences about the latent populations when the stretch and color parameters are uncorrelated, as correlated parameters imply that there is redundant information in the two parameters. Removing the correlation between the parameters x₁ and c also has a physically intuitive meaning; if the color is sourced by dust extinction, it should not depend on the processes associated with light-curve stretch.

The phenomenological c parameter is often referred to as a "color." However, c is not defined by the value of a particular photometric color at a particular phase, but rather by a fixed effect on the B–V color at all phases. ¹⁴ The choice of B–V color as the normalization does not restrict our color model to any particular wavelength dependence. As the color law is defined over the whole wavelength range of the model, the c parameter can be inferred using data from any wavelengths, subject to the model uncertainties derived from the data during the training process. These definitions are sufficient to define the SALT3 flux model without degeneracy, and are enforced by taking them as constraints on the fitted model parameters { m ₀ , m ₁ , cl } during the model training process.

As Rubin (2020) and Rose et al. (2020) suggest that SN Ia SEDs are determined by 3–5 SED parameters and a color, a model with a single SED parameter (x₁ in our case) will not capture the full diversity of the population. We refer to the variation unexplained by our model as "in-sample variance." This variation can be contrasted with the uncertainty in the model parameters due to training with a finite sample of SNe with finite signal-to-noise photometry, which we refer to as "out-of-sample variance." For a photometric observation in a filter with central wavelength λ_c and phase p, in-sample variance is addressed by a "model variance" composed of two terms. The first is a diagonal uncertainty defined as

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{f}^{2}(p,{\lambda }_{c}) & = & {\left[{x}_{0}\exp (c\cdot {CL}({\lambda }_{c}))\displaystyle \int T\left(\displaystyle \frac{\lambda }{1+{z}_{\mathrm{Hel}}}\right)\lambda d\lambda \right]}^{2}\\ & & \times \left[{\sigma }_{{M}_{0}}^{2}(p,{\lambda }_{c})+{x}_{1}^{2}{\sigma }_{{M}_{1}}^{2}(p,{\lambda }_{c})\right.\\ & & \left.+2{x}_{1}{C}_{{M}_{0}{M}_{1}}(p,{\lambda }_{c}){\sigma }_{{M}_{0}}(p,{\lambda }_{c}){\sigma }_{{M}_{1}}(p,{\lambda }_{c})\right],\end{array}\end{eqnarray} \tag{ 3 }$$

where ${\sigma }_{{M}_{0}}(p,{\lambda }_{c};{\sigma }_{{{\rm{M}}}_{0}})$ and ${\sigma }_{{M}_{1}}(p,{\lambda }_{c};{\sigma }_{{{\rm{M}}}_{1}})$ represent variability in the flux surfaces and ${C}_{{M}_{0}{M}_{1}}(p,{\lambda }_{c};{{\boldsymbol{C}}}_{{{\boldsymbol{M}}}_{{\bf{0}}},{{\boldsymbol{M}}}_{{\bf{1}}}})$ is the correlation between flux surfaces. These variance terms are described by zeroth-order B-splines (equivalent to binning the data by phase and wavelength) with eight basis functions in wavelength and 12 basis functions in phase. As detailed in Section 2.5, we use a maximum likelihood estimator to determine the B-spline coefficients ${{\boldsymbol{\sigma }}}_{{{\boldsymbol{M}}}_{{\bf{0}}}},{{\boldsymbol{\sigma }}}_{{{\boldsymbol{M}}}_{{\bf{1}}}},{{\boldsymbol{C}}}_{{{\boldsymbol{M}}}_{{\bf{0}}},{{\boldsymbol{M}}}_{{\bf{1}}}}$. This is distinct from the approach of SALT2, which took the in-sample variance to have the same form as that of the out-of-sample variance, scaling the latter (evaluated by leave-one-out tests) using a smooth function or "error snake" to fix the ${\chi }_{\nu }^{2}$ of the photometry of the training sample to unity. Because the in-sample variance is determined by variability of the underlying population of SNe light curves while the out-of-sample variance is determined by the distribution of available training data, we consider our approach to be a better account of the remaining variance of the SN Ia population.

The second term of the in-sample variance is a covariant "color scatter" that allows light curves of the same supernova in different bands to be coherently offset relative to one another. This model is similar to chromatic models of intrinsic scatter like that of Guy et al. (2010) and the diagonal terms of the covariance matrix from Chotard et al. (2011). The color scatter model describes all photometric observations made using the same bandpass as correlated, with no correlation between measurements in different bandpasses. This is quantified by a relative covariance k²(λ_c ; a ) between two fluxes, which is parameterized by the exponential of a fourth-order polynomial in wavelength, where the polynomial coefficients a are fitted parameters. Thus, the model covariance matrix for two photometric observations i and j in broadband filters X_i ,X_j of a given supernova with model fluxes ${\vec{f}}_{\mathrm{Model}}$ is

$$\begin{eqnarray}&&k(\lambda ;{\boldsymbol{a}})=\exp (\sum _{i=0}^{4}{{\boldsymbol{a}}}_{i}{\lambda }^{i})\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\begin{array}{l}{({{\rm{\Sigma }}}_{\mathrm{Model}})}_{{ij}}={\delta }_{{ij}}{\sigma }_{f}^{2}({p}_{i},{\lambda }_{c(i)})\\ \ \ +\ \left\{\begin{array}{cl}{k}^{2}([{\lambda }_{c}]){\left({\vec{f}}_{\mathrm{Model}}\right)}_{i}{\left({\vec{f}}_{\mathrm{Model}}\right)}_{j} & {X}_{i}={X}_{j}\\ 0 & \mathrm{otherwise}\end{array}\right..\end{array}\end{eqnarray} \tag{ 5 }$$

Although the SALT3 model is intended for use as a photometric light-curve model, spectral data are included in the training process in order to better constrain the shape of spectral features. However, given that spectral data have larger calibration uncertainties compared to broadband fluxes, we follow the SALT2 training code in "recalibrating" spectral data, modulating the model by a smooth function to match the continuum of the observed spectrum. We modify the spectral flux equation (Equation (1)) for use with spectra during the training by removing the color term and replacing it with a recalibration term of similar form,

$$\begin{eqnarray}&&{F}_{\mathrm{spec}}(p,\lambda )={x}_{0}[{M}_{0}(p,\lambda )+{x}_{1}{M}_{1}(p,\lambda )]\exp (\sum ^{{n}_{\mathrm{RC}}}{y}_{i}{\lambda }^{i}/i!),\end{eqnarray} \tag{ 6 }$$

where the recalibration term is $\exp (\sum ^{{n}_{\mathrm{RC}}}{y}_{i}{\lambda }^{i}/i!)$, and the spectral recalibration nuisance parameters y_i are fitted during the training procedure. When spectra are perfectly calibrated, the fitted recalibration term will reproduce the effect of the color parameter on the spectrum; by removing the color law entirely from the spectral flux equation, we remove any impact of miscalibrated or contaminated spectra on the color law. Any broadband or low-frequency miscalibration of the training spectra will thus have no impact on the final model. High-frequency miscalibration of the spectra or contamination from host-galaxy light may propagate biases into the final model. The quantity n_RC controls how many recalibration parameters are allowed for each spectrum, and is determined based on the wavelength extent of the spectrum and the number of filter bands available for that SN (additional filter bands better constrain the recalibration term, allowing for more free parameters). The model variances are defined similarly to the photometry, without the contribution of the color scatter:

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{F}^{2}(p,\lambda ) & = & {\left[{x}_{0}\exp (\sum _{i}^{n}{y}_{i}{\lambda }^{i}/i!))\right]}^{2}\\ & & \times [{\sigma }_{{M}_{0}}^{2}(p,\lambda )+{x}_{1}^{2}{\sigma }_{{M}_{1}}^{2}(p,\lambda )\\ & & +2{x}_{1}{C}_{{M}_{0}{M}_{1}}(p,\lambda ){\sigma }_{{M}_{0}}(p,\lambda ){\sigma }_{{M}_{1}}(p,\lambda )].\end{array}\end{eqnarray} \tag{ 7 }$$

In regions of phase-wavelength space that are poorly constrained by spectra, the M₀ and M₁ components can acquire artifacts such as high-frequency ringing, or deconvolution noise. To reduce these artifacts, we use a "regularization" procedure that penalizes large derivatives in the model where there is an absence of spectroscopic data, as parameterized by a binned spectral density function N_eff(p, λ), with every spectrum assigned equal weight. We implement three kinds of regularization: phase gradient, wave gradient, and dyadic, each of which is similar to the regularization terms of Mosher et al. (2014). These are applied to the two flux surfaces M₀(p, λ) and M₁(p, λ). For a flux surface S(p, λ), the regularization terms are defined as

• 1.  
  Phase gradient regularization penalizes large derivatives with respect to phase in less-constrained regions

  $$\begin{eqnarray}&&{\chi }_{\mathrm{Phase}}^{2}[S(p,\lambda )]={A}_{\mathrm{Phase}}\sum _{i}^{2{N}_{\mathrm{Phase}}}\sum _{j}^{2{N}_{\lambda }}\displaystyle \frac{{\left({\left.\tfrac{\partial S}{\partial p}\right|}_{\genfrac{}{}{0em}{}{p={p}_{i}}{\lambda ={\lambda }_{j}}}\right)}^{2}}{{N}_{\mathrm{eff}}({p}_{i},{\lambda }_{j})}.\end{eqnarray} \tag{ 8 }$$
• 2.  
  Wave gradient regularization similarly penalizes large derivatives with respect to wavelength in less-constrained regions

  $$\begin{eqnarray}&&{\chi }_{\mathrm{Phase}}^{2}[S(p,\lambda )]={A}_{\mathrm{Wave}}\sum _{i}^{2{N}_{\mathrm{Phase}}}\sum _{j}^{2{N}_{\lambda }}\displaystyle \frac{{\left({\left.\tfrac{\partial S}{\partial \lambda }\right|}_{\genfrac{}{}{0em}{}{p={p}_{i}}{\lambda ={\lambda }_{j}}}\right)}^{2}}{{N}_{\mathrm{eff}}({p}_{i},{\lambda }_{j})}.\end{eqnarray} \tag{ 9 }$$
• 3.  
  Dyadic regularization encourages the flux surfaces to be separable in phase and wavelength, and is zero when a flux surface takes the form S(p, λ) = g(p) × h(λ)

  $$\begin{eqnarray}&&\begin{array}{l}{\chi }_{\mathrm{Dyadic}}^{2}[S(p,\lambda )]={A}_{\mathrm{Dyadic}}\displaystyle \sum _{i}^{2{N}_{\mathrm{Phase}}}\displaystyle \sum _{j}^{2{N}_{\lambda }}\displaystyle \frac{1}{{N}_{\mathrm{eff}}({p}_{i},{\lambda }_{j})}\\ \ \times \ {\left[{\left.\displaystyle \frac{\partial S}{\partial p}\right|}_{\displaystyle \genfrac{}{}{0em}{}{p={p}_{i}}{\lambda ={\lambda }_{j}}}{\left.\displaystyle \frac{\partial S}{\partial \lambda }\right|}_{\displaystyle \genfrac{}{}{0em}{}{p={p}_{i}}{\lambda ={\lambda }_{j}}}-S({p}_{i},{\lambda }_{j}){\left.\displaystyle \frac{{\partial }^{2}S}{\partial p\partial \lambda }\right|}_{\displaystyle \genfrac{}{}{0em}{}{p={p}_{i}}{\lambda ={\lambda }_{j}}}\right]}^{2}.\end{array}\end{eqnarray} \tag{ 10 }$$

The relative strength of the three regularization terms is determined by the weights {A_Phase, A_Wave, A_Dyadic}, which are inputs to the SALTshaker code. The summations over phase and wavelength points are evenly spaced over the model phase and wavelength ranges, with twice as many points as basis functions along each axis. As regularization can bias the model surfaces by oversmoothing them (Mosher et al. 2014), we tune the weights to ensure that the regularization terms do not contribute significantly to the total χ². Further work to determine how this regularization scheme affects the model and to choose optimal model configurations will require applying our training and analysis framework to simulations (M. Dai et al. 2021, in preparation; Pierel et al. 2021).

Based on the model definitions above, we construct the training χ² for the model as

$$\begin{eqnarray}&&{\chi }_{\mathrm{Total}}^{2}={\chi }_{\mathrm{Phot}}^{2}+{A}_{\mathrm{Spec}}\displaystyle \frac{{N}_{\mathrm{Phot}}}{{N}_{\mathrm{Spec}}}{\chi }_{\mathrm{Spec}}^{2}+{\chi }_{\mathrm{Constraints}}^{2}+{\chi }_{\mathrm{Reg}}^{2},\end{eqnarray} \tag{ 11 }$$

where N_Phot, N_Spec are, respectively, the number of photometric and spectroscopic data points.

The photometric χ² term is

$$\begin{eqnarray}\begin{array}{rcl}{\chi }_{\mathrm{Phot}}^{2} & = & \displaystyle \sum _{n}^{{N}_{\mathrm{SN}}}{\left({\vec{f}}_{\mathrm{obs}}^{(n)}-{\vec{f}}_{\mathrm{Model}}^{(n)}\right)}^{T}{\left({{\rm{\Sigma }}}_{\mathrm{Total}}^{(n)}\right)}^{-1}({\vec{f}}_{\mathrm{obs}}^{(n)}-{\vec{f}}_{\mathrm{Model}}^{(n)})\\ {{\rm{\Sigma }}}_{\mathrm{Total}}^{(n)} & = & \mathrm{diag}\left(\left({\sigma }_{{Phot}}^{(n)}\right)2\right)+{{\rm{\Sigma }}}_{\mathrm{Model}}^{(n)},\end{array}\end{eqnarray} \tag{ 12 }$$

where for the nth supernova, ${\vec{f}}_{\mathrm{obs}}^{(n)}$ is the vector of observed photometric fluxes and ${\vec{f}}_{\mathrm{Model}}^{(n)}$ is the vector of the model fluxes integrated over the photometric bandpass. The covariance ${{\rm{\Sigma }}}_{\mathrm{Total}}^{(n)}$ combines photometric uncertainties and the model uncertainty covariance described in Equation (5). The factor A_Spec is a constant "spectral suppression" term that downweights the contribution of spectra to the training χ² in order to reduce the sensitivity of the training to unknown systematic errors in the spectral data (Guy et al. 2007). We set this term such that the spectral and photometric data have roughly equal contributions to the total training χ². The spectral χ² term is then defined

$$\begin{eqnarray}&&{\chi }_{\mathrm{Spec}}^{2}=\sum _{n}^{{N}_{{SN}}}\sum _{i}^{{N}_{\mathrm{spec}}}\sum _{j}^{{N}_{\mathrm{points}}}\displaystyle \frac{{\left[{F}_{\mathrm{spec}}({p}^{(n,i)},{\lambda }_{j}^{(n,i)})-{\left({\vec{f}}_{\mathrm{obs}}^{(n,i)}\right)}_{j}\right]}^{2}}{{\left({\sigma }_{\mathrm{obs}}\right)}_{j}^{2}+{\sigma }_{\mathrm{Model}}{\left({p}^{(n,i)},{\lambda }_{j}^{(n,i)}\right)}^{2}},\end{eqnarray} \tag{ 13 }$$

where for the nth supernova, ${F}_{\mathrm{spec}}({p}^{(n,i)},{\lambda }_{j}^{(n,i)})$ is the model spectral flux at the jth wavelength bin of the ith spectrum, ${\left({\vec{f}}_{\mathrm{obs}}^{(n,i)}\right)}_{j}$ is the observed spectral flux, ${\left({\sigma }_{\mathrm{obs}}\right)}_{j}$ is the photometric uncertainty, and ${\sigma }_{\mathrm{Model}}({p}^{(n,i)},{\lambda }_{j}^{(n,i)})$ is the model uncertainty in Equation (7). While calculating ${\chi }_{\mathrm{phot}}^{2}$ and ${\chi }_{\mathrm{spec}}^{2}$, we redden the model to account for the Milky Way E(B − V) along the line of sight to each SN using the reddening law of Fitzpatrick (1999), but we do not include this term in the equations above, for simplicity. The ${\chi }_{\mathrm{Constraints}}^{2}$ is composed of penalty terms used to enforce the model definitions described in Section 2.1, and the regularization term is defined as

$$\begin{eqnarray}\begin{array}{rcl}{\chi }_{\mathrm{Reg}}^{2} & = & \displaystyle \sum _{S(p,\lambda )}^{\{{M}_{0},{M}_{1}\}}{\chi }_{\mathrm{Phase}}^{2}[S(p,\lambda )]+{\chi }_{\mathrm{Wave}}^{2}[S(p,\lambda )]\\ & & +{\chi }_{\mathrm{Dyadic}}^{2}[S(p,\lambda )].\end{array}\end{eqnarray} \tag{ 14 }$$

The SALTshaker code is initialized with the configuration parameters shown in Table 1. SALTshaker then determines best-fit values of the model parameters. These parameters are shown in Table 2. We define flux model parameters as those that control the flux surfaces and color law, uncertainty model parameters as those that control the output model uncertainties, and nuisance parameters as those that describe individual supernovae.

SALTshaker alternates between simultaneously fitting the nuisance parameters and model flux parameters while keeping uncertainties fixed and fitting the model uncertainty parameters while keeping model fluxes fixed. The color scatter is kept fixed with the corresponding covariance term set to k(λ_c ) = 0 during this stage of the training process, as we find that allowing nonzero color scatter results in biased flux surfaces due to the regularization procedure. The flux and nuisance parameters are fit using an iterative Levenberg–Marquardt algorithm. To reduce the number of (computationally expensive) Jacobian evaluations of the model residuals required, we use Schubert's method to perform a rank one update on the Jacobian after each iteration while maintaining its sparsity structure (Schubert 1970; Marwil 1979). However, this technique is unsuitable for determining the model uncertainties.

To fit the model uncertainty parameters defined in Equation (5), we define a log-likelihood $({{ \mathcal L }}_{\mathrm{Phot}.\mathrm{Errs}.}$) given as

$$\begin{eqnarray}&&-2\mathrm{log}({{ \mathcal L }}_{\mathrm{Phot}.\mathrm{Errs}.})={\chi }_{\mathrm{Phot}}^{2}-\sum _{n}^{{N}_{{SN}}}\mathrm{log}(| {{\rm{\Sigma }}}_{\mathrm{Total}}^{(n)}| ).\end{eqnarray} \tag{ 15 }$$

The model uncertainty parameters are chosen to maximize the log-likelihood using the optimizer iMinuit while the flux model and nuisance parameters are kept fixed (James & Roos 1975; Dembinski et al. 2020). After the training has converged with the color scatter fixed such that k(λ) = 0, we use iMinuit to fit the parameters controlling the color scatter; during this final fit, the color law is allowed to vary (having previously been fit as a flux model parameter). The out-of-sample variance is then estimated by inverting the Hessian matrix of the flux parameters obtained from the model fitting process, with the regularization terms suppressed by a factor of 100, and propagating those parameter uncertainties into the flux surfaces. The model uncertainties and out-of-sample variance are added together as the total model uncertainty surface.

The computation time required for the SALT3.JLA training sample with SALTshaker is approximately three minutes per iteration, with ∼25 iterations required for convergence. On the larger SALT3.K21 sample, iterations are significantly longer, with ∼25 minutes required per iteration. We note that 29-Å spectral binning, which was used in the SALT2.JLA training, improves speeds substantially by reducing the amount of data by up to an order of magnitude. We include this as an option in SALTshaker, but the models in this work use the native binning of the input spectra. Finally, the slowest component of the code is the iterative fitting of the error model, which for SALT3.K21 requires approximately 4 hours to perform 80 iterations and reach convergence through iMinuit. Error model iterations are performed every five iterations by default, but could likely be performed less often without adversely affecting the final model. Faster methods of error model fitting will be an important avenue for future improvement and should improve speeds considerably.

An overview of the SALTshaker training procedure is given in Figure 2.

Zoom In Zoom Out Reset image size

To distinguish between models trained on different input data, we define "SALT3.X" as the SALT3 model created with SALTshaker using training sample X. We refer to the samples used to evaluate the performance of a given SALT3 model as "validation" samples, as it will be expedient to evaluate the trained model in some cases on data that were not used in the model training. We use the SNANA light-curve fitting program to fit validation samples with both SALT2.JLA and SALT3.X.

Given fitted SALT parameters m_B , x₁, c, the variances for each parameter ${\sigma }_{x}^{2}$, and the covariances between the parameters σ_x,y , the Tripp estimator for distance modulus is

$$\begin{eqnarray}&&\mu ={m}_{B}+\alpha \cdot {x}_{1}-\beta \cdot c-{ \mathcal M }\end{eqnarray} \tag{ 16 }$$

with distance uncertainties

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{\mu } & = & {\sigma }_{\mathrm{int}}^{2}+{\sigma }_{\mu ,z}^{2}+{\sigma }_{\mathrm{lens}}^{2}+{\sigma }_{{m}_{B}}^{2}+{\left(\alpha {\sigma }_{{x}_{1}}\right)}^{2}+{\left(\beta {\sigma }_{c}\right)}^{2}\\ \ & & +2\alpha \beta {\sigma }_{c,{x}_{1}}+2\alpha {\sigma }_{{m}_{B},{x}_{1}}+2\beta {\sigma }_{{m}_{B}{,}_{c}},\end{array}\end{eqnarray} \tag{ 17 }$$

where σ_int, α, β, and ${ \mathcal M }$ are nuisance parameters, σ_μ,z is computed from a peculiar velocity uncertainty of 250 kms⁻¹, and σ_lens = 0.055z. We use the SALT2mu method (Marriner et al. 2011; Kessler & Scolnic 2017) implemented in SNANA to estimate the nuisance parameters as well as distances in five redshift bins. Allowing for a shift in location and scale of the light-curve parameters, the observed distributions are similar as compared to SALT2.JLA (see Section 5.2.2). We thus expect selection biases to be common between the two models, and therefore we do not use SALT2mu to correct for selection biases.

Given estimated distance moduli μ and Hubble residuals relative to a nominal ΛCDM cosmology (Δμ) for both models, we define two metrics. First, RMS(Δμ) across the validation sample and the relative distance difference between models defined as

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{Diff}({{\rm{\Delta }}}_{z}\mu ) & = & \left[\mu (0.40\lt z\lt 0.6| \mathrm{SALT}3.X)\right.\\ & & \left.-\mu (0.40\lt z\lt 0.6| \mathrm{SALT}2.\mathrm{JLA})\right]\\ & & -\left[\mu (0.01\lt z\lt 0.2| \mathrm{SALT}3.X)\right.\\ & & \left.-\mu (0.01\lt z\lt 0.2| \mathrm{SALT}2.\mathrm{JLA})\right],\end{array}\end{eqnarray} \tag{ 18 }$$

where μ(Z∣M) indicates a weighted average distance across a redshift range Z given a model M. Additionally, we will show binned Hubble residuals across the redshift range. We discuss differences in the nuisance parameter β, as this has physical implications for dust properties; however, α depends on both the demographics of the training sample and our revised separation of stretch and color, and thus has no useful comparison across models (see Section 2.1). Similarly, we compare synthetic photometry, but these comparisons are most relevant for simulated data when the underlying model is known. We train on multiple simulated and real data samples to assess how our models perform; in Table 3, we summarize the abbreviations used for these training samples.

The JLA data used to train the SALT2.JLA model consist of 420 SNe with light curves, 83 of which include spectroscopy, from a compilation of low-z samples (see Table 4 for details and citations), the Sloan Digital Sky Survey (SDSS; Holtzman et al. 2008; Kessler et al. 2009b; Sako et al. 2018), and the Supernova Legacy Survey (Astier et al. 2006, with spectra from Walker et al. 2011; Balland et al. 2018; and private communication with M. Betoule and C. Balland). These data are summarized in Table 4, along with the new training data discussed in Section 4, and are included in the data release at http://saltshaker.readthedocs.io/.

We use the "Supercal" procedure (Scolnic et al. 2015) to update the photometric calibration (including filter function central wavelengths and photometric zero points) of the JLA training sample. Supercal uses the 3π sky coverage of the PS1 system and its observations of secondary standard stars to precisely determine the offsets between different photometric systems. Scolnic et al. (2015) found that applying these calibration corrections to the sample of SN Ia used to measure the Hubble flow without updating the training calibration could shift w by 0.026. While we do not study the impact of the calibration on the trained model, Taylor et al. (2021) examines this with a full reanalysis of the DES-SNIa cosmology results using a retrained SALT2 model. We also use the Schlafly & Finkbeiner (2011) corrections to the Schlegel et al. (1998) Milky Way dust maps.

To test SALTshaker with a known input model and cosmology, we first simulate the SALT2.JLA training data to produce our "simJLA" training sample. Training on the actual JLA training sample is also a useful test (Section 3.5), but it does not provide a known truth model to validate the training outputs. For this simJLA test, every simulated SN that goes into the training sample has SALT2.JLA as the simulated model, allowing us to test the consistency of input and output models.

We use the SNANA software (Kessler et al. 2009c) to generate Monte Carlo realizations of SN photometry and spectroscopy mimicking the demographics of the JLA training sample of 420 SNe. SNANA simulations are frequently used to simulate a random realization for a sample, in order to explore SN distance biases as a function of redshift. Our goal here is different: we aim to accurately simulate every individual event in the JLA sample, including cadence and signal-to-noise for both photometric and spectroscopic observations. We therefore use the cadence, redshift, and best-fit c and x₁ from each JLA light curve as input to the simulation. Each set of measured SN properties (z, t₀, x₁, c) is used to simulate a rest-frame SED with the SALT2.JLA model. Spectral variations (intrinsic scatter) are added to the rest-frame SED using the method in K13 and the covariance model from the (Chotard et al. 2011, hereafter C11) scatter model. ¹⁵ We do not use the default (Guy et al. 2010, hereafter G10), as the G10 model has nonphysical scatter at extremely red wavelengths that does not match observations. As described in Kessler et al. (2019) the simulation applies cosmological dimming, lensing, peculiar velocity, galactic extinction, and redshifting to produce a redshifted SED at the top of the atmosphere. The filter transmissions and cadence are used to determine measured fluxes and uncertainties. Finally, the noise properties of each measured spectrum are applied to the simulated SED.

This simulation does not exactly reproduce the data set, because of random fluctuations in photometric noise and intrinsic scatter—and also because the underlying SALT2 model formalism is an approximation, as discussed in Pierel et al. (2021). Nonetheless, the simulation is similar to the data and is therefore sufficient for testing SALTshaker, as illustrated in Figure 3 for a representative low-z SN Ia. Comparisons between the parameters of simulations and data after fitting with the SALT2.JLA model are shown in Figure 4, with only slight observed differences in average m_B uncertainty.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We use SALTshaker to train a model using our simulated JLA sample, producing a model we call SALT3.simJLA. In Figure 5, we show the relative difference between the input SALT2.JLA synthetic photometry and the corresponding synthetic photometry recovered from the SALT3.simJLA model. There is significant discrepancy in the ultraviolet, where regularization strongly impacts the recovered model spectra, but at central passband wavelengths 4000 Å < λ < 7000 Å and − 10 days < p < 30 days, the other three light curves (B, V, and R) are consistent with the simulated model at a level better than 1%. To illustrate the impact of changes in the color law on a "typical" light curve in units of magnitudes, we use the quantity σ_c · ΔCL(λ), where the standard deviation of the distribution of the SALT color parameter c is ∼0.1. As can be seen from Figure 6, the color-law difference is >0.05 mag when λ < 3000 Å, the regime where the color law is least constrained.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We find that the RMS of σ_c ΔCL(λ) is 0.01mag between wavelengths 3500 Å < λ < 7000 Å, central filter wavelengths for which the SALT2.JLA model is considered reliable. We conclude that, because few SNe constrain the SED at λ < 3000 Å, the limited set of light-curve parameters c, x₁ in this wavelength region poorly constrain the color law and the spectral components. In Sections 4 and 5, we substantially expand the training data in order to address this issue.

Next, we compare distances from the two models. To avoid the statistics of the training sample limiting the precision of our validation, we simulate a "large" JLA simulation for validation, with summary statistics for the Hubble diagram fits shown in Table 5, row 1. Instead of simulating exact x₁ and c values for an "apples-to-apples" comparison with the SALT2.JLA training sample, we generate ∼3000 SNe mimicking a combination of CfA3 (Hicken et al. 2009a), SDSS, and SNLS by using x₁ and c distributions from Scolnic & Kessler (2016). We show the Hubble residuals of this simulated cosmology sample in Figure 7. Given that the Scolnic & Kessler (2016) distributions of c and x₁ are assumed uncorrelated and the training sample has demographics similar to those of the JLA sample, we can expect to recover the same α, β parameters as were measured using SALT2, which will not be true of the real data. We observe consistent measurements of β and RMS (Δμ), with σ_int higher by 0.005 mag due to our treatment of uncertainties. We also observe a slightly higher value of α, which is attributed to the x₁/M₁ degeneracy. The distance difference between the models (Equation (18)) is consistent (Diff(Δ_z μ) = 7 ± 11 mmag) despite the U-band discrepancies in the light curves seen in Figures 5 and 6. As the validation sample here has demographics similar to those of the training sample, the same sparsity of data that allows the observed ultraviolet divergences causes those same divergences to have little effect on the inferred distances. We conclude that, for a SALT model to be reliable at these wavelengths, the density of data in this region must be increased to match the density of data where we recover the input model to ∼1%, i.e., at least a factor of two increase (see Section 4, where we discuss the density of photometric data across the JLA wavelength range).

Zoom In Zoom Out Reset image size

Next, we run SALTshaker on the real JLA training sample; we refer to this trained model as SALT3.JLA. We follow the validation procedure of Section 3.1 using the JLA training set for our validation sample. As shown in row 2 of Table 5, we find a similar α, a slightly lower β, and consistent RMS (Δμ). We find that Diff(Δ_z μ) = 32±28 mmag is consistent with zero. The value of σ_int is slightly higher for SALT3.JLA, attributable to greatly decreased model uncertainties. Given the similar RMS scatter, this σ_int difference does not result in reduced distance precision. The distance moduli are consistent at the 1σ level between SALT2.JLA and the SALT3.JLA, and the model surfaces are consistent within the uncertainties. We show the Hubble residuals of this model in Figure 8 (along with our extended SALT3 models discussed in Section 5).

Zoom In Zoom Out Reset image size

The Carnegie Supernova Project (Stritzinger et al. 2011) observed supernovae discovered by both amateur observers and professional surveys, with the goal of producing a stably calibrated and homogeneous set of light curves throughout the optical and NIR. Photometric observations were conducted using the Henrietta Swope Telescope and du Pont Telescope in optical (ugriBV) and NIR (YJHK_s ) wavelengths.

We make use of optical observations from the second data release of the CSP (Stritzinger et al. 2011) as well as optical spectroscopic data published in Folatelli et al. (2013). Stritzinger et al. (2011) reported photometric observations of 50 objects, with median redshift of z ∼ 0.02, of which 26 were included in the Pantheon analysis. A third data release with revised photometric calibration and more observations has been made public in Krisciunas et al. (2017); however, there is not an available cross-calibration of this data relative to other SN Ia surveys. Our compilation includes 13 SNe included in the Pantheon analysis, which additionally pass our more stringent requirements for density of photometric data. Future work plans to make use the NIR photometric observations of the CSP, as well as other available NIR data, to extend our model into the NIR (J. D. R. Pierel et al. 2021, in preparation).

The CfA4 sample consists of 94 SNe observed by the Center for Astrophysics Supernova Group (Hicken et al. 2012). These objects were discovered by a mixture of amateur and professional supernova searches. Photometric observations were made in the uUBVri bands using the F. L. Whipple Observatory and the KeplerCam instrument. Depending on photometric band, there are an average of ∼10 points in each light curve.

A total of 41 SNe from the CfA4 survey were included in the Pantheon analysis, and 30 pass our more stringent cuts for light-curve quality and are included in our compilation of training data.

The Foundation Supernova Survey (Foley et al. 2018) followed SNe using the Pan-STARRS1 telescope, and measured well-calibrated SN light curves in griz filters with a five-day cadence near maximum light and an approximately eight-day cadence beginning at +10 days after maximum light. To achieve reduced selection effects compared to previous surveys that targeted bright, preselected low-z galaxies, Foundation primarily followed SNe discovered by untargeted surveys such as the All-Sky Automated Survey for Supernovae (Shappee et al. 2014), the Asteroid Terrestrial-impact Last Alert System (Tonry et al. 2018), Gaia (Gaia Collaboration et al. 2016), and the Pan-STARRS Survey for Transients (Huber et al. 2015).

The Foundation first data release in Foley et al. (2018) contains 225 SNe Ia, 180 of which are cosmologically useful, and 153 of which pass our own light-curve quality cuts. These data have been used to measure cosmological parameters in Jones et al. (2019) and the correlation of host galaxy properties with SN distances in Jones et al. (2018). The iz band coverage of Foundation is particularly critical to creating a SALT3 model that is trained to redder wavelengths than enabled by the JLA data alone. We include spectra for 114 Foundation SNe from the Dettman et al. (2021) data release.

The Pan-STARRS medium deep survey covered 70 square degrees of sky over four years, discovering approximately 5200 SNe (Jones et al. 2017; Villar et al. 2020) and spectroscopically classifying ∼10% of these at a median redshift of ∼0.35 (Rest et al. 2014). The Pantheon analysis, which combined these data with JLA, includes 279 PS1-observed SNe Ia with an average of approximately 6 observations per 10 days in griz (Scolnic et al. 2018). We use the 266 of these SNe Ia that pass our light-curve quality cuts in our training data.

The Dark Energy Survey (DES) three-year spectroscopically classified SN sample contains 207 SNe Ia at a median redshift of 0.36 (Abbott et al. 2019) (of which 206 pass our light-curve quality cuts). These SNe were discovered by imaging eight 2.7 deg² "shallow" fields (depth ≈ 23.5 mag) and two 2.7 deg² "deep" fields (depth ≈ 24.5 mag) approximately once per week (Smith et al. 2020b). Transients were discovered using a difference-imaging pipeline (Kessler et al. 2015), and final photometry was performed with a "scene modeling" pipeline described in Brout et al. (2019b). See Abbott et al. (2019) and Brout et al. (2019a) for additional details regarding the DES SN Ia data and analyses. These data have a maximum redshift of ∼0.85 and complement SNLS in probing rest-frame near-UV wavelengths with well-calibrated (subpercent) photometric data (Burke et al. 2018).

To ensure that SNe are suitable for inclusion in a training sample, we first require that every SN in the compilation is a spectroscopically classified, Branch-normal SN Ia (see Branch et al. 1993; we also include 1991T-like SNe Ia, following the original SALT2 training procedure of Guy et al. 2007). For the SN light curves, we make the following selection requirements (cuts):

• 1.  
  At least four epochs at phases between −10 < p < 35 days, where p is the rest-frame phase relative to time of B-band peak.
• 2.  
  At least one measurement after peak brightness (5 < p < 20), to constrain the shape.
• 3.  
  At least one measurement in each of at least two filters at −8 < p < 20, to constrain the color of the SN.
• 4.  
  At least one measurement prior to peak brightness (−10 < p < −1), to ensure a well-measured time of maximum light. This is the only cut that was not included in the original SALT2 training ¹⁶ .

For the spectra, we include the original SALT2.JLA training spectra in addition to spectra taken from the kaepora database and spectra taken as part of the Foundation Supernova Survey. A number of these spectra were taken by the Foundation team, but most have been published on the Transient Name Server ¹⁷ as classification spectra.

To ensure minimal host-galaxy contamination in the high-z SN spectra used for model training, Guy et al. (2007) fit the spectra with a combined model including the predicted SN spectrum, the spectral recalibration parameters, and a galaxy model (elliptical, S0, Sa, Sb, and Sc templates). They removed spectra for which there was evidence for host-galaxy contamination at the 68% confidence level.

For Foundation and the additional low-z SN spectra included here, low-z SNe are much brighter relative to their host galaxies than at the redshifts probed by SNLS and SDSS. We clip host galaxy lines, mask regions with uncorrected telluric features, and remove excessively noisy or poorly calibrated regions of each spectrum, but do not attempt to subtract a host galaxy continuum. We remove a handful of spectra with poor quality from visual inspection. After cuts, there are 114 spectra from Foundation and 693 from Kaepora, a subset of which are shown in Section 5. All but five Foundation SNe have redshifts measured from host galaxy features. The complete training data are available at https://saltshaker.readthedocs.io.

Before presenting our model trained on all data in Section 5.2, we check the performance on separate training and validation samples. We take the "K21valid" compilation as our cosmology sample, using our "SALT3.K21train" model as described in Section 3.1, and show the resulting Hubble residuals in Figure 8 and in row 3 of Table 5. We find that nuisance parameters are similar between SALT2.JLA and the SALT3.K21train model, with σ_int slightly higher by 0.01 mag but with a consistent RMS(Δμ) (the SALT3 RMS is negligibly smaller). Diff(Δ_z μ) is consistent with zero at the mmag level.

Finally, we train our best SALT3 model, which we call SALT3.K21, using all of the data described in previous sections as a training sample. The sample includes 1083 SNe, a factor of 2.5 more SNe than the JLA training sample, and 1207 spectra, a factor of three increase in the number of spectra. Synthetic light curves from SALT3.K21 and SALT2.JLA are compared in Figure 11, and the model uncertainties are compared in Figure 12. We see good consistency between SALT2.JLA and SALT3.K21, with modest differences in the u-band and some additional differences in redder bands; at both wavelength ranges, we have substantially increased the training sample (see Figure 10). Similarly, as shown in Figure 13, the color law is consistent with that of SALT2.JLA to within 1% across the entire wavelength range. As has previously been noted, the inferred color law does not resemble observed reddening laws determined from Milky Way extinction measures such as those of Fitzpatrick (1999) and Cardelli et al. (1989) in the ultraviolet and red optical wavelengths, as might be expected if SN Ia color variations are caused exclusively by host-galaxy dust similar to that of the Milky Way.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We compare Hubble residuals of the SALT3.K21 and SALT2.JLA models in Figure 8. Individual standardized distances are consistent to 0.05 mag between the two models, and the effects on the Hubble diagram are found in row 4 of Table 5. Diff(Δ_z μ) is consistent with zero at 3 ± 14 mmag. Finally, row 4 of Table 5 shows nuisance parameters and Hubble diagram metrics, demonstrating that SALTshaker produces a new SALT3 model with slightly lower total dispersion, consistent σ_int, and consistent distances.

We note that the β parameter is lower by 0.17 in the SALT3 model, likely due to a reconsidered separation of color and stretch. The nuisance parameter β has been physically interpreted by comparison to R_B = R_V + 1, the selective-to-absolute extinction ratio of a dust reddening law, finding that it is lower than expected from the typical Milky Way R_V = 3.1. These comparisons have been and continue to be informative, as they demonstrate that either the dust reddening laws of SN Ia host galaxies are distinct from that of the Milky Way, or that circumstellar dust may play some role in the reddening of SN Ia (see also Jha et al. 2007; Chotard et al. 2011, Burns et al. 2014; Amanullah et al. 2015; Brout et al. 2021). As we believe that the previous SALT model had confused the parameter c with some intrinsic color sourced by the parameter x₁ under the assumption that x₁ could have no effect on peak color of SN Ia, a different value of β is expected. As the β we observe is still inconsistent with expected R_B values from the mean reddening of the Milky Way, the aforementioned conclusions are not substantially affected.

5.2.1. Model Uncertainties and Hubble Scatter of SALT3.K21

The uncertainties in Figure 12 show that both the color scatter and light-curve uncertainties in redder bands are much lower in SALT3.K21 than SALT2.JLA. The color scatter, representing an intrinsic variance of the SN Ia population unexplained by the model, is not in general expected to improve as the size of a training sample grows. Other light-curve uncertainties visualized in the right panels of Figure 12 include contributions from both the aforementioned "in-sample" variations, which are beyond the ability of a two-parameter model to address, as well as the "out-of-sample" uncertainties in the determination of the model, which scale with the amount of training data. However, the SALT2.JLA color scatter was itself unconstrained by the JLA data past mean passband wavelengths ∼8000 Å, and we find that, while the color scatter is significant at these wavelengths, it is much smaller than implied by the SALT2.JLA model extrapolated into this region. Our additional data constrain this region up to ∼8500 Å, and redder data will be required to see how this effect carries into the NIR. We also note lower color scatter by ∼1% in the blue, which we attribute to improved relative calibration of the training sample.

Fitting the light curves using the SALT3.K21 model, we compare the uncertainties in distances and light-curve parameters to those found using SALT2.JLA. In Figure 14, we compare performance across the redshift range. We have chosen to show binned values of $\alpha {\sigma }_{{x}_{1}}$ and β σ_c , as these are the more physically relevant quantities; any simple rescaling of the x₁, c parameters would leave these quantities unaffected. Our model shows reduced Hubble scatter and distance uncertainties over the SALT2.JLA model at nearly every redshift, with the least improvement at moderate redshifts z ∼ 0.2, where S/N is high and the SALT2.JLA model is already performing well. Light-curve parameters from SALT3.K21train have smaller uncertainties across the redshift range, with the exception of x₁ uncertainties. Our largest improvements are found at low redshift, where the improved red wavelength coverage of our model allows the use of additional light curve bands in fitting these SNe, providing stronger constraints on light-curve parameters. As shown in panel (d) of Figure 14, the color parameter c in particular is better constrained by the use of redder data, as it provides a longer lever-arm with which to infer the underlying parameter. That these redder bands may have other underlying processes not well-represented by a single color parameter is accounted for in our construction by the color scatter determined from the data, but even in z-band these uncertainties do not exceed ∼0.04 mag. We conclude that photometric observations at these wavelengths, previously discarded in SALT2 cosmology fits, continue to provide cosmologically useful information about the SALT3 light-curve parameters.

Zoom In Zoom Out Reset image size

There is a slight improvement in RMS Hubble residuals, largely attributable to decreased color uncertainties in the low-redshift sample. Breaking this down further, Figure 15 shows the Hubble scatter binned by the number of additional filters used in the SALT3.K21 light-curve fit as compared to the SALT2.JLA light-curve fit. Where two additional filters are available, RMS (Δμ) improves by ∼10%. We conclude improvement in the SALT3 model is most noticeable when it allows us to fit SNe with existing light curves in filters out of the SALT2.JLA wavelength range. At high redshift, SALT3's improved constraints on the NUV model reduce Hubble scatter by ∼0.01 mag, although there is typically not sufficiently red data to take full advantage of the extended wavelength range at these redshifts. Similarly, when breaking down results by survey, as we show in Figure 15, the greatest improvement is in the low-z CFA3 and Foundation samples, where we can make use of I- and iz-band observations (respectively) previously unused in SALT2.JLA-based analyses at low redshift. We are able to reduce the RMS (Δμ) of the Foundation sample from 0.144 ± 0.001 mag to 0.125 mag, an improvement of 15%.

Zoom In Zoom Out Reset image size

5.2.2. Comparison of SALT3.K21 and SALT2.JLA Light-curve Parameters

The SALTshaker training procedure results in different distributions of light-curve parameters compared to SALT2.JLA. Some differences are due to changes in how the SALT3 model is defined, while others are from the demographics of the training sample. For example, a x₁/M₁ degeneracy is broken by setting a constraint on the standard deviation ${\sigma }_{{x}_{1}}=1$, both in our procedure and in the original training procedure; therefore, including additional data with higher stretch increases the scale of the M₁ component.

In Figure 16, we compare the distributions of the x₁ and c parameters from both models. We note a slight rotation, shift, and scale change of the color/stretch distribution relative to SALT2.JLA. These linear transformations cancel in the Tripp standardization of distance, absorbed into the nuisance parameters $\alpha ,\beta ,{ \mathcal M }$. The offset in color in particular is from including bluer high-redshift SNe from PS1 and DES in the training sample. The rotation is due to the distinct procedures for separating stretch and color between the SALT2 training code and SALTshaker. SALT3.K21 x₁ and c values can be approximated from the SALT2.JLA values by the transformations

$$\begin{eqnarray}&&{x}_{1}^{(\mathrm{SALT}3)}\approx 1.028{x}_{1}^{(\mathrm{SALT}2)}+0.138{c}^{(\mathrm{SALT}2)}+0.005\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{c}^{(\mathrm{SALT}3)}\approx 0.002{x}_{1}^{(\mathrm{SALT}2)}+0.985{c}^{(\mathrm{SALT}2)}+0.013.\end{eqnarray} \tag{ 20 }$$

Zoom In Zoom Out Reset image size

To test the similarity of the shapes of the distributions, which can potentially affect selection bias corrections for cosmology, we use the two-dimensional Kolmogorov–Smirnov—like test statistic of Press & Teukolsky (1988). As the means and scales of the distributions lack physical information, we modify the procedure by linearly transforming both sets of light-curve parameters to have the same mean, standard deviation, and x₁-c correlation before calculating the statistic. We find the test statistic for the SALT2.JLA and SALT3.K21 x₁ and c parameters, then bootstrap-resample the data and calculate the test statistic with the resampled data. We derive a p-value = 0.895, and conclude that we cannot distinguish the shapes of the underlying distributions.