# Addressing type Ia supernova color variability with a linear spectral
template

Cássia S. Nascimento
Instituto de Física, Universidade Federal do Rio de Janeiro, 21941-972, Rio de
Janeiro, RJ, Brazil
<EMAIL_ADDRESS>
& João Paulo C. França
Centro Brasileiro de Pesquisas Físicas, Rua Xavier Sigaud st. 150, 22290-180,
Rio de Janeiro, Brazil
Instituto de Física, Universidade Federal do Rio de Janeiro, 21941-972, Rio de
Janeiro, RJ, Brazil
<EMAIL_ADDRESS>
& Ribamar R. R. Reis
Instituto de Física, Universidade Federal do Rio de Janeiro, 21941-972, Rio de
Janeiro, RJ, Brazil
Observatório do Valongo, Universidade Federal do Rio de Janeiro, 20080-090,
Rio de Janeiro, RJ, Brazil
<EMAIL_ADDRESS>

###### Abstract

We investigate the potential of a pure linear expansion for the rest-frame
flux of a type Ia supernova light curve fitter based on the well known
Spectral Adaptive Light Curve Template 2 (SALT2). We generate the expansion
components by performing Principal Component Analysis (PCA) and Factor
Analysis (FA) onto a representative training set. Then, we derive a Tripp-like
expression for the distance modulus and fit the $\Lambda$CDM cosmological
model on the Pantheon sample. The constraining power of the model, dubbed Pure
Expansion Template for Supernovae (PETS), and SALT2 is evaluated and we found
compatible results for $\Omega_{m0}$ and $\Omega_{\Lambda 0}$ within $68\%$
uncertainty between the two models, with PETS’ fit parameters exhibiting non
negligible linear correlations with SALT2’ parameters. We find non negligible
correlations between PETS’s fit parameters and the supernovae host galaxies
masses, while the Hubble Diagram residues show no correlation with fit
parameters, redshift or host galaxy mass. The model nuisance parameters,
$\alpha$ and $\beta$, are slighted correlated and we find evidence for
redshift evolution for $\beta$. The intrinsic scatter,
$\sigma_{\textrm{int}}$, shows a subtle redshift evolution that should be
further investigated increasing the number of high redshift supernovae in the
cosmology sample.

_Keywords_ supernovae: general $\cdot$ dust, extinction $\cdot$ methods: data
analysis $\cdot$ methods: statistical $\cdot$ cosmological parameters

##  1 Introduction

Measuring the acceleration of the universe and its properties has been a
direct consequence of modeling Type Ia supernovae light curves ([47, 42]).
From the relationship between magnitude and decay rate of the light curves,
[43] confirmed the possibility of using Type Ia supernovae (SNe), as standard
candles in cosmology. Since then there has been a constant movement of
revisitation and improvements regarding the training and fitting of Type Ia
supernova models ([4, 39, 44]). This constant improvement is more recently
powered by the need to reduce the systematic errors of the empirical
standardization process as the statistics improved with the dataset growth.
For the next generation surveys such as Rubin (LSST) ([28]) and Nancy Grace
Telescope (WFIRST) ([63]), we can expect to detect around 300.000 supernovae
up to a redshift of $z=3$ ([54]). This will correspond to an increase in the
available Ia SNe sample of 2 orders of magnitude, resulting in a substantial
reduce in statistical errors when measuring cosmological parameters.

As a consequence of the increasing datasets, modeling techniques have been
tailored, in particular those that make use of Spectral Energy Distributions
(SEDs) such as SALT2 ([19, 20]) and SNEMO ([58]). However, despite the well-
established use of SALT2, [4] found after the standardization a remaining
dispersion in magnitude of $m=0.15$ resulting from systematic errors
associated with, among other origins, dust extinction in the host galaxy and
the standardization method itself. Aiming to better understand the impact of
color variation description in Ia SNe rest-frame flux and inspired by the
SALT2 light curve fitting model, in this paper we extensively explore the
description of a generic Type Ia supernova spectral surface by a pure third-
order phase-spectroscopic expansion model, initially proposed by [58].

Following SALT2 and SNEMO we do not separate the sources of color variation.
Not separating the intrinsic and dust extinction reddening contributions can
obscure one of the most important known correlations of type Ia supernovae,
i.e. intrinsic color correlation with luminosity. However, a complete
separation of both effects requires accurate dust distribution estimates for
every host galaxy in the sample. Recent efforts were made using Ia SNe spectra
to separate intrinsic and extrinsic color variations, [9, 57], both finding an
total-to-selective extinction ratio consistent with the Milk-way value.
Nevertheless, from the best of our knowledge, there is no widely accepted
procedure to completely disentangle these contributions using only photometric
data. We also chose not to deredden our data, avoiding the assumption both
intrinsic and dust contributions can be described by tuning an extinction law,
as is usually done for most fitters. Instead, we apply the decomposition
methods to the reddened data after the SEDs reconstruction process and
investigate the model performance when considering $\Lambda$CDM cosmology with
our derived expression for the distance modulus. This description can be
naturally extended to both UV and infrared regions by the expected increasing
training sample.

Our model was trained over the Nearby Supernova Factory Data Release 9 (DR9)
([58]) data sample. On the spectra we applied a forward-backward filter to
reduce systematic errors and then started the reconstruction of the SEDs. We
perform 1-d Gaussian Process Regressions with gaussian priors over the kernel
parameters and a mean behaviour information from [26] to ensure a high
resolution SED reconstruction. Next, we derive a distance modulus from the
observed Ia SN rest-frame flux. Inspired by SUGAR ([37]) and SNEMO ([58])
models, the spectral surfaces were decomposed using decomposition methods such
as Principal Component Analysis and Factor Analysis. In the end, we compared
the results in the context of light curve fitting and cosmological fits,
investigating, among other behaviours, the evolution of the Hubble Diagram
residues with respect to SNe redshifts.

In section 2, we present the pure expansion model and its peculiarities such
as the training grid and limitations. In section 3, we illustrate the pre-
training process of filtering the spectra and reconstructing the SEDs. In
section 4, we present the process of decomposing the spectral surfaces via
Principal Component Analysis (PCA) and Factor Analysis (FA), and analyse the
light curve fits of a validation sample. In section 5, we derive our distance
modulus equation and discuss the Pantheon fitting results. In sections 6 and
7, we discuss the cosmological results for both PCA and FA approaches when
considering $\Lambda$CDM cosmology. Finally, in section 8, we present the
final remarks.

##  2 The pure expansion spectral model

Based on the homogeneity hypothesis for type Ia supernovae, we analyse in
details the light curve fitting and cosmology of a pure linear expansion
model, a variation of the empirical flux modeling by SALT2 (Spectral Adaptive
Light-Curve Template 2), [19, 20], initially proposed by [58] and further
investigated here. The underlying model describes the rest frame flux of a
given type Ia supernova, $\phi(p,\lambda;{\bf x})$, as a pure expansion,

$\phi(p,\lambda;{\bf
x}):=x_{0}[M_{0}(p,\lambda)+x_{1}M_{1}(p,\lambda)+x_{2}M_{2}(p,\lambda)+...].$
(1)

Here $p$ and $\lambda$ are, respectively, phase (i.e. number of days since
maximum light in B-band) and wavelength. The parameter $x_{0}$ controls the
rest frame flux amplitude while the set of surfaces $M_{i}(p,\lambda)$
consists of SEDs that, together with the remaining free parameters, $x_{1}$
and $x_{2}$, should accommodate the observed dispersion for these objects.
From now on, we will refer to the pure expansion model as PETS (Pure Expansion
Template for Supernovae).

This modeling arises from selecting a descriptive training sample and assuming
a general type Ia supernova can be described as a linear combination of them.
Next, applying a transformation to obtain an ordered basis and computing the
amount of variability explained by the first few components enable us to
reduce the dimensionality of the linear expansion while still describing a
great amount of information about the initial set of objects.

The SALT2 model proposes an analogous expression, keeping only the first two
terms, and multiplying this expansion by an exponential of a free parameter,
$c$, times a function of wavelength, a color law $CL(\lambda)$, the latter
also being determined by the training process. The SALT2 model exploit the two
main empirical correlations seen in these objects, luminosity correlation with
light curve shape and luminosity correlation with color, both perceived by
[43]. Their $M_{0}$ surface is identified as the mean surface of the training
sample while $M_{1}$ allows for additional variation to the mean surface,
being interpreted as a stretch-like feature, while the exponential term
accounts for color variations, often without mention to the nature of this
variation.

The first step when obtaining a SALT2-like spectral time series model with an
exponential term is to remove the reddening effect from the spectra and light
curves inputted in the model training. This process relies on the assumption
that dust and intrinsic variations affect the spectra in the same fashion, and
by tuning the $R_{V}$ parameter from an extinction law such as CCM89, [7], or
F99, [13], it is possible to compensate for both effects. The phase
independent reddening is then entirely described by this exponential term that
deforms the time series templates along the wavelength axis. In this
description $c$ contains a environmental dependent fraction of color that is
not correlated with luminosity. In this case, two supernovae with the same
exact parameters are not guaranteed to have the exact same luminosity.

Following SALT2 and SNEMO, we chose not to distinguish reddening due to dust
and intrinsic variations when constructing the model. Instead, we argue that
if our sample is representative it should comprise a good range of supernova
luminosities with different dust extinction contributions, making it still
possible to standardize them even with non negligible dust extinction. On top
of that, the distance modulus here derived removes the necessity of
correlating the environmental dependent fraction of $c$ with luminosity
variations.

Using equation (1) we can obtain the photometric flux model (m) in observer
frame, $F_{Y}^{(m)}(p(1+z))$, see [33], by integrating the spectroscopic flux,
$\phi(p,\lambda;{\bf x})$, through a given bandpass $Y$ with observer frame
transmission function $T_{Y}(\lambda)$:

$F_{Y}^{(m)}(p(1+z))=(1+z)\int_{0}^{\infty}d\lambda[\lambda\phi(p,\lambda)T_{Y}(\lambda(1+z))]\,,$
(2)

where $z$ is the redshift. It is important to reinforce our parameters can not
be immediately recognized as stretch or color. For now, the fit parameters are
non-physical parameters that corrects the type Ia supernova magnitude
dispersion grounded on these object shared similarities.

When considering the linear model portrayed in equation (1) without a previous
derredening, we remove the necessity of approximating the unknown intrinsic
color variation behaviour by the effect of dust reddening. In addition, when
performing this integration for a flux with an exponential term, one may
consider to approximate this term, whether by a first order Taylor expansion
or by choosing $\lambda=\lambda_{eff}$. The former approximation starts to
break with increasing color parameter at the UV region, reaching 30% relative
discrepancy at 3000Å for $c=0.3$, while the latter approximation departs from
the original expression in both directions centered on $\lambda=\lambda_{eff}$
for each filter, reaching 50% of relative discrepancy at 3000Å for $c=0.3$.
Hence the pure linear model avoid yet this other approximation when evaluating
equation (2).

Figure 1: 2-d histogram of SNFactory DR9 spectra used for model training. The
white dashed lines represent the boundaries of our model, selecting the region
with higher data coverage.

Figure 2: SALT2 $M_{0}$ template, in z axis we have flux in units of
erg/s/cm2/Å. The area inside the white dashed lines represent the region
covered by our PETS models.

Our training sample, which we will discuss in further detail in Section 3 is
the Supernova Nearby Factory Data Release 9 from [58]. In Fig. 2 we show the
2-d histogram of these training sample spectra and the model boundaries chosen
to ensure good data coverage.

Fig. 2 represents the density plot of the SALT2 $M_{0}$ spectral surface with
the white dashed lines marking the region where our model is defined. We can
see at low wavelengths around the date of maximum light a non-negligible flux
outside the limits of our model. Even though this is an important region where
type Ia supernovae show greater variability, it is also responsible for less
precise results in every light curve fitting method and is often omitted when
fitting for cosmology. Increasing the number of spectra in this region will
allow for an immediate extension of the pure linear expansion model without
any approximation that breaks in the UV region.

##  3 The data set and pre-processing steps

Aiming to train our template surfaces ${\bf M}=(M_{0},M_{1},M_{2})$, we chose
the Data Release 9 from Nearby Supernova Factory111Data Release 9 available in
https://snfactory.lbl.gov/, [1]. The sample used here for model training and
validation consists of 2466 spectra from 171 spectroscopically confirmed type
Ia supernovae with redshift ranging from $0.01$ to $0.08$. The spectra are
already shifted to a common rest-frame ($z=0$) and corrected for Milky Way
dust extinction following [59, 7] (see [8, 52] for more details) and the
observer-frame B-band date of maximum light from a SALT2 analysis is already
subtracted. For more information about these procedures see [58].

This data set consists of high-quality selected supernovae that have at least
five spectra, of those at least one is prior to the date of maximum light, and
at least four are between 10 and 35 days. All dates are measured relative to
the observer-frame B-band maximum light estimate, that defines the rest-frame
phase as $p=(t-t_{B,max})/(1+z)$. This higher number of spectroscopic data
than is usually available allows us to construct a SED for each supernova, in
a similar fashion as seen in [58]. These SEDs are surfaces of specific flux in
units of erg/s/$\textnormal{cm}^{2}$/Å (multiplied by an arbitrary factor) as
a function of phase, in units of days, and rest-frame wavelength, in units of
Å. To extract the components that best explain the variability in our data
set, we begin by interpolating our SEDs and evaluating them in a regular grid.

The majority of the spectral data from our training sample range from 3300Å to
8500Å in wavelength and from -15 to 50 days in phase. Based on the 2-d
histogram seen in Fig. 2, we chose the model boundaries and mesh grid size as
a regular grid from -10 to 50 days, with 1-day bins and from 3400Å to 8400Å,
with 10Å bins. This choice ensures good data coverage while preventing the
loss of bandpasses when fitting photometric data. Our model then includes the
optical region and a small portion of near-infrared, where type Ia supernovae
luminosities usually show less dispersion.

###  3.1 Pre Processing: Filtering the spectra

In order to avoid nonphysical fluctuations of the specific flux we applied a
forward-backward digital filter Filtfilt, from the Scientific Python library
SciPy, [70]. Following [37] we decided to apply a cutoff of $100$ Å, which
represents the lowest half-width absorbing line in type Ia supernovae typical
spectra.

We derive the uncertainty on the flux filtered spectra as the sum in
quadrature of the contributions from measurement error and a filtering one.
The latter is necessary to take into account the effect seen in some high
oscillatory portions of the data, where the filtering can lead to values
considerably far from the original data. This contribution is equal to one-
third of the difference between the real and filtered data. This procedure
will be important since, when applying the Gaussian Process Regression (GPR),
oscillations around the real values in conjunction with underestimation of
uncertainties can wrongfully lead to a nonphysical highly oscillatory
regression.

###  3.2 Pre Processing: Applying Gaussian Process Regression

After the filtering process we interpolate each spectra and its uncertainties
and evaluate them onto a regular wavelength grid. Then, for each supernova we
position their regularized data in a three-dimensional space of specific flux
per wavelength per phase and to completely reconstruct each supernova SED we
perform a Gaussian Process Regression (GPR) on each plane of constant
wavelength.

Gaussian processes (GP) are supervised learning methods used in regression and
classification problems. For a more detailed discussion see [45]. Here we seek
one different function that explains the data for each wavelength plane. Each
of these functions will be depicted by a one-dimensional GP with mean
$m(\textbf{p})=\mathbb{E}[f(\textbf{p})]$ and covariance function
$K(\textbf{p},\textbf{p}^{\prime})=\mathbb{E}[(f(\textbf{p})-m(\textbf{p}))(f(\textbf{p}^{\prime})-m(\textbf{p}^{\prime}))]$,

$f(\textbf{p})\sim\mathcal{GP}(m(\textbf{p}),K(\textbf{p},\textbf{p}^{\prime})).$
(3)

The input values in the regression, $y_{i}$, are called the target values and
they are, apart from a gaussian noise $\epsilon_{i}$, equal to the function we
wish to model via a GP,

$y_{i}=f(p_{i})+\epsilon_{i},\\\
\epsilon_{i}\sim\mathcal{N}(0,\sigma_{i}^{2}),$ (4)

in our case we perform a heteroscedastic regression since each input point has
a different noise value.

The outline is to choose a kernel (i.e. a covariance function) that generates
a family of curves with the desired characteristics and retain, after several
samplings, only the ones that can describe the target values. Averaging over
these curves provides the predictive function and uncertainties estimates.
This can be accomplished in a bayesian modeling by describing a joint
distribution for the target values $y_{i}$ and for the function $f$, evaluated
at the test points (i.e. new points where we wish to calculate the function),
$f(p_{*})=:f_{*}$, and later conditioning the joint distribution on the
observations.

We start from the assumption that the observed target values, $y_{i}$, and the
function $f$ evaluated at the test points, $f_{*}$, can be described by the
same multivariate Gaussian distribution. This allows us to construct the joint
distribution,

$\begin{bmatrix}\textbf{y}\\\
\textbf{f}_{*}\end{bmatrix}\sim\mathcal{N}\left(\textbf{m}(P),\begin{bmatrix}K(P,P)+\sigma^{2}&K(P,P_{*})\\\
K(P_{*},P)&K(P_{*},P_{*})\end{bmatrix}\right),$ (5)

where the multivariate distribution covariance matrix, $\boldsymbol{\Sigma}$,
was depicted in terms of the kernel function evaluated at each pair of points,
whether from the train points vector, $P$, or test points vector, $P_{*}$.
Lastly, $\sigma^{2}$ is a diagonal matrix carrying the target values
variances. This distribution is completely defined when specifying a mean
vector, m(P), and covariance matrix, $\boldsymbol{\Sigma}$.

We can then obtain the conditional distribution that provides the predictions
for the function evaluated at the test points given the test points itself,
the training points, and the corresponding observed target values,

$\textbf{f}_{*}|P,\textbf{y},P_{*}\sim\mathcal{N}(\bar{\textbf{f}}_{*},\textnormal{cov}(\textbf{f}_{*})),$
(6)

where the predictive mean is

$\bar{\textbf{f}}_{*}=\textbf{m}(P_{*})+K(P_{*},P)[K(P,P)+\sigma^{2}]^{-1}(\textbf{y}-\textbf{m}(P)),$
(7)

and the predictive covariance is

$\textnormal{cov}(\textbf{f}_{*})=K(P_{*},P_{*})-K(P_{*},P)[K(P,P)+\sigma^{2}]^{-1}K(P,P_{*}).$
(8)

Depending on the data and the prior knowledge about the function behavior we
can choose different kernels. A widely used one is the Radial Basis Function
kernel, RBF, which describes the covariance between two points as an
exponential decline with their squared distances and a set of two parameters
that tune the variance and the correlation range.

We employ the Matérn Kernel for our GPRs, which includes an additional
parameter, $\nu$, to regulate the smoothness of the function. The lower this
parameter, the less smooth the function becomes, as compared to the RBF
Kernel. In the limit that $\nu\rightarrow\infty$, we recover RBF. The Matérn
kernel function is expressed as

$k(p_{i},p_{j})=\frac{\sigma^{2}}{\Gamma(\nu)2^{\nu-1}}\left[\frac{\sqrt{2\nu}(p_{i}-p_{j})^{2}}{\Delta
l}\right]^{\nu}K_{\nu}\left[\frac{\sqrt{2\nu}}{\Delta
l}(p_{i}-p_{j})^{2}\right],$ (9)

where $\sigma^{2}$ is the variance and $\Delta l$ is the length scale
parameter which controls the correlation range. $K_{\nu}$ is the modified
Bessel functions and $\Gamma(\nu)$ is the Gamma function for a given $\nu$.
Here we are specially interested in $\nu=5/2$ to keep a reasonable smoothness.

We apply this GPR formalism to our data via the Python library [17]. A very
important step when performing the GPR is using the mean information,
$\textbf{m}(P)$, otherwise, the spectra without pre-maximum data lead to a
wrong behavior at lower phases, as can be seen in Fig. 3 (in this context the
maximum is specific to each monochromatic light-curve and does not coincide
with B-band date of maximum light). This issue affects several supernovae and
directly influences the SEDs reconstruction mainly at the region of low
wavelengths and before maximum B-band light, profoundly impacting the
forthcoming feature extraction. To avoid this incorrect regression we adopted
the type Ia supernova template from [26] as the aforementioned mean.

Figure 3: Gaussian process regression for the constant wavelength plane
$\lambda=3450$Å. The black error bars are the filtered spectra data. In orange
and blue we have GPR performed over this set of points with and without the
mean information, respectively. The solid black curve represents the template
used as mean in the former case.

The original data was named by [58] as either Train or Test supernova followed
by an identification number. Thus, the Train and Test flags do not represent
our choice of training or validation objects, it just identifies a specific
supernova. From the initial sample of 171 supernovas, six were excluded, two
due to problems in GPR and four due to poor quality in their reconstructed
SEDs. Regarding the latter, they had too few spectra taken after maximum, and
even with template information they showed a large amount of nonphysical
negative flux in this unconstrained region, they are Train_SN93, Test_SN15,
Test_SN26 and Train_SN96. Train_SN30 was also excluded since it showed several
variations from the mean which are not observed in any other object in this
sample.

For each GPR we perform a Maximum Likelihood Estimation (MLE) to tune the
kernel parameters. A parameter inference showed many GPRs had difficulty in
constraining the length scale parameter, $\Delta l$, often reaching as
extremely low values as $10^{-3}$ days, when it should be in the range of $9$
days to $30$ days since it represents the horizontal data scale. To compensate
for this issue we added a gaussian prior on both kernel parameters when
performing the MLE, successfully avoiding these low values that lead to
overfitted regressions.

The length scale gaussian prior changes as wavelength grows. In general, the
data shows a peak before $t_{max,B}$ at lower wavelengths with a horizontal
scale of about 15 days, value which was chosen as the initial prior. At around
5500Å these monochromatic light curves get wider, the peak moves smoothly and
hence we readjust the prior. Growing even further in wavelength another peak
appears, changing the horizontal scale again. Essentially, at every point
where we observed a specific behavior common to most SNe, we programmed
changes to the prior. During this entire process, the monochromatic light
curves are normalized to one, avoiding readjusts to the variance gaussian
prior. This stage of the training process is carried out manually.

Another interesting behavior occurs at higher wavelengths, where the MLE
favors lower length scales because of what seems to be a random fluctuation of
data not accompanied by an explanatory uncertainty. This apparent
underestimation of uncertainties leads to many undesirable regressions. Then,
for every GPR at planes of $\lambda\geq 5000$Å , a small fixed white kernel of
variance $5\times 10^{-4}$ was added to the Matérn kernel, allowing the
predictive curve to deviate slightly more from the target values. This white
kernel can not be added to every GPR, otherwise, as we are using a
personalized mean function, if the data is already really well explained by
it, the prediction turns out to be exactly equal to the mean, since all small
variations on the data will be solely explained by noise. For lower
wavelengths we made fewer exceptions, only when directly detecting this same
behavior of a wiggly prediction (low $\Delta l$), a new regression was
performed adding the same fixed white kernel.

##  4 Model Training

###  4.1 Dimensionality reduction of supernovae SEDs with Principal Component
Analysis

From the set of reconstructed SEDs, we select a random subsample of 147 to
compose the training set and save the remaining objects for a further
validation step. Based on the homogeneity hypothesis of type Ia supernovae,
one may argue most of the information on this sample should be correlated and
partially redundant. Thus, following [19, 35, 56, 21], we investigate
initially an orthonormal transformation to an ordered lower dimension basis
via a Principal Component decomposition, initially proposed by [58] and here
more thoroughly investigated.

The Principal Component Analysis (PCA), [25, 40], aims to find the
uncorrelated directions that successively maximize the explained variance of
the original data. This process translates into diagonalizing the sample
covariance matrix, ordering the eigenvectors according to the highest
eigenvalues, and projecting the original data onto the transpose eigenvectors,
leading us to the ordered Principal Components (PCs). More detailed
discussions can be found in [30] and [3].

This process assumes the training set is representative of the general pair of
type Ia supernova plus host galaxy and thus form a suitable basis for a linear
description. Additionally, it assumes that there is a lower dimension
hyperplane capable of describing the original data with a few first Principal
Components. It is important to note PCA does not assume the existence of
hidden variables and the Principal Components which are linear combinations of
the original basis are not affected by the choice of the new basis dimension
and neither are directly associated with any physical interpretation.

To find accurate PCs using Singular Value Decomposition (SVD) the mean of each
feature, in our case each supernova, is required to be zero. Therefore, we
remove the mean prior to applying PCA and then, we transform the original
vector of means to the new basis, returning a contribution for each PC. This
process is necessary otherwise the first surface that we will further verify
that resembles an average type Ia supernova flux would show non physical
negative values. A further treatment we choose not to apply is scaling each
feature such that the standard deviation is equal to one. This is common when
dealing with features of different scales, in which cases PCA could prioritize
some over others. However, the scaling would treat the most common objects in
equity with some less common, diverse objects with lower relative rest-frame
fluxes seen in our sample.

We perform this dimensionality reduction through the python package Scikit-
learn, [41]. The input data is a $30561\times 164$ matrix with the flattened
SEDs placed as columns. Fig. 4 shows the cumulative explained variance as a
function of the number of components, $N_{c}$. For $N_{c}=3$ we explain
$98.2\%$ of the sample variance while for $N_{c}=10$ we explain about
$99.2\%$, but at the cost of including seven new components, which later
translates to seven new light-curve fit parameters. This highly explained
variance by the first component alone reassures the homogeneity of the type Ia
supernovae composing our sample.

Figure 4: Cumulative explained variance in terms of number of components, both
for Principal Component Analysis and Factor Analysis. As expected, first order
templates holds a higher variance. Together, $M_{0}$, $M_{1}$ and $M_{2}$
dominates the variability, exceeding/reaching $98\%$ explained variability in
PCA/FA description. PCA as expected explains better the sample variability
than FA, since the former focus on maximizing this quantity while FA focus on
the off-diagonal elements of the data covariance matrix.

In Fig. 6 we see the projection of the original data onto the first three
components of the ordered basis for different planes of constant wavelength.
The top image of Fig. 6 shows our $M_{0}(p,\lambda)$ element defined in
equation 1, this surface has only positive values and resembles the typical
behavior of a normal Ia SN SED, as occurs for most spectral template fitting
models. The following $M_{1}(p,\lambda)$ and $M_{2}(p,\lambda)$ elements, also
depicted in the same figure, add the finer details when describing a
particular SED. Their main variations concentrate close to the day of maximum
light. And while $M_{1}(p,\lambda)$ shows a more symmetric behavior around
zero flux, $M_{2}(p,\lambda)$ shows higher amplitudes at specific wavelengths.

Figure 5: The original components projected onto the first three PCA
components. From top to bottom we show $M_{0}$, $M_{1}$ and $M_{2}$. This
components covers more than $98\%$ of the total explained variance.

Figure 6: The original data projected onto the first three FA components in
the absence of rotation

.

We verify that the major explained variability surface is less dependent on
the training set in contrast to the ones that hold less explained variability.
Indeed, $M_{0}(p,\lambda)$ is the most similar among the other spectral
fitting methods, as SALT2, [19], SNEMO2, SNEMO7, and SNEMO 15, [58], while the
remaining surfaces are more sample dependent.

The first surface, $M_{0}(p,\lambda)$, is often associated with the effective
mean SED of the training sample. One could then expect the second and third
elements to show some similarities with the next higher-order statistical
moments, i.e standard deviation and skewness. Nonetheless, Fig. 6 shows
$M_{1}(p,\lambda)$ is not positive definite, and substituting the components
for these statistical moments showed a worse performance.

###  4.2 Dimensionality reduction of supernovae SEDs with Factor Analysis

Another method that was implemented for type Ia supernova feature extraction
by [58] and [37] is Factor Analysis (FA). FA relies on an important underlying
assumption that there are latent (i.e. unobserved) variables which can explain
our original data, often with a reduced set of components.

As we discussed previously, the PCs are constructed as linear combinations of
the original variables, without any mention of their physical interpretations
nor with an explicit model. Nonetheless, for Factor Analysis, the observed
quantities, x, are the ones assumed to be linear combinations of the latent
variables, f, except for an error term, e, (i.e.
$\textbf{x}=\boldsymbol{\Lambda}\textbf{f}+\textbf{e}$), with the
$\boldsymbol{\Lambda}$ matrix carrying the coefficients. This generative model
also requires a number of assumptions over the distributions and correlations
of the data, the latent variables, also known as common factors, and the error
terms, also known as the specific factors. These three quantities are assumed
to have null expected values and the common and specific factors are assumed
to be uncorrelated, with themselves and with each other.

Without loss of generality the common and specific factors are considered to
be described by gaussian distributions,
$\textbf{f}\sim\mathcal{N}(0,\textbf{I})$ and
$\textbf{e}\sim\mathcal{N}(0,\boldsymbol{\Psi})$, respectively. Then the
conditional distribution of the observed variables is given by

$p(\textbf{x}|\textbf{f})=\mathcal{N}(\textbf{x}|\boldsymbol{\Lambda}\textbf{f}+\boldsymbol{\mu},\boldsymbol{\Psi}).$
(10)

Here $\boldsymbol{\mu}$ is an offset and $\boldsymbol{\Psi}$ is a diagonal
matrix with different entries, characterizing a heteroscedastic noise. This
decomposition consists in an iterative method using maximum likelihood
estimate with an SVD approach and it was also implemented using the python
package Scikit-learn, [41]. It is interesting to note that there is a
probabilistic description of PCA that is closely related to FA, Probabilistic
Principal Component Analysis, known as PPCA. This also linear-Gaussian
framework differs from FA only by defining an homoscedastic, isotropic noise,
where $\boldsymbol{\Psi}=\sigma^{2}\textbf{I}$. In the limit
$\sigma^{2}\rightarrow 0$ it recovers the PCA results, hence the name. For a
more detailed discussion on these methods see [5] and [68].

Overall, as argued by [30], PCA and FA can be understood as explaining
different aspects of the sample covariance matrix. While PCA concentrates in
explaining the diagonal elements through maximizing explained variance, FA
concentrates on the off-diagonal elements by forcing the specific factors to
be uncorrelated and leaving the off-diagonal elements to be fully explained by
the common factors.

In addition, the solution for f is not unique, arbitrary rotations can be
performed over this matrix and some rotations can facilitate recognizing the
role of the new variables. We can perform orthogonal and oblique rotations,
the latter allowing for correlated factors. These rotations will be further
explored in a future work, for now we provide the results for FA in absence of
rotations.

Since FA searches for latent variables, it is common to analyse the
correlations between the variables in the new basis with the variables in the
old one, and recognize which older variables are more relevant, gaining some
interpretation to these new projected values.

As in PCA, for accurate results it is necessary to center each feature before
applying the algorithm. The mean values are stored, transformed and added back
to each new basis component, after the model fitting.

In Fig. 6 we represent the surfaces that define a basis for describing a
general type Ia supernova in FA context. The input was the same subset used
for PCA feature extraction consisting of 147 SEDs constructed with SNFactory
spectra and Gaussian Processes. As previously seen for our PCA surfaces, for
SNEMO and SALT2, the first surface resembles the structure of an average Ia SN
SED. The remaining surfaces then will be responsible for correcting finer
details that differentiate each supernova. As the majority of the emitted flux
is concentrated closer to the day of maximum light and around $3000$Å to
$6000$Å, it is expected that the corrections also concentrate on the same
region. The surfaces aside from the first one show portions of negative flux
and it is interesting to note the second one shows a symmetry with respect to
the horizontal axis, with the correction changing sign when growing in
wavelength. And last, the third surface shows a peak around this same region
where the second surface changes sign.

To have the same number of free parameters as the most widely used fitters, we
assume the existence of three hidden variables. In Fig. 4 we see a comparison
in explained variance between the two approaches. As expected the PCA
performance is slightly better since this algorithm focuses on maximizing the
explained variance. The model depicted in this figure, for comparison only, is
FA with 10 hidden variables. FA with three components will also explain
roughly $98\%$ of this sample variability.

###  4.3 Validation fits

The remaining 17 Ia supernovae that did not take part into the feature
extraction form a test set for model validation. For this step we used the
Python package for Supernova Cosmology SNCosmo, [2]. With SNCosmo we were able
to create a spectral time series as a function of arbitrary parameters and fit
the light curves by generating synthetic photometry from the SEDs of the
remaining Ia supernovae. To draw maximum information from this data we chose
the set of filters used by The Carnegie Supernova Project (CSP), [36]. Our
photometric system consists on filters CSP-g, CSP-r and CSP-V9844 with AB
magnitudes. The first two filters cover a range from about 4000Å to 7000Å but
with a lower coverage around 5500Å, this same region is mainly contained in
CSP-V9844. These three filters extract an important amount of information from
our original SEDs and capture a good portion of the optical region.

The considered models were created truncating the rest-frame flux expansion,
equation (1), whether after the first three terms or after the first ten
terms. The fits quality was measured through chi squared per degrees of
freedom. This quality fit ruler is not ideal since it strongly depends on the
choices made to estimate the measurement errors in this synthetic photometry.
However, as the validation fits only intends to compare the fitting ability of
the methods and establish that the PCA/FA process produce reasonable model
components, we maintain this parameter as a ruler for now.

Whether using PCA or FA the ten components versions performed better. When
adopting PCA, 12 of the 17 validation fits had really good agreement between
the three and ten components versions, the latter performing something similar
to a fine tuning. The other five fits showed a greater improvement going from
the PCA version to the FA version. For three components PCA, the chi squared
per degrees of freedom was usually no lower than ten and sometimes reaching up
to one extra order of magnitude, while for the ten components PCA in many
occasions it was lower than one, probably indicating overfitting. Most of the
light curves detailed behavior was captured by this ten components model but
at the cost of adding seven new parameters.

It is important to note three components PCA has explained variability of
$98.2\%$, but this number refers to fitting the SEDs and not the flux
integrated along a given set of filters. In other words, we should not expect
the ability of explaining a SED as a linear combination of the model
components to be similar in performance when dealing with correlated
photometric data. We argue the three component versions struggles with the
fine details due to this highly correlated synthetic data, since there is
great overlap in filter transmissions.

Regarding the FA versions, the three components model performs better than its
PCA analogous and this can be visualized in Fig. 7. The transition between
three to ten FA components models seems to just tune the peaks height.
Possibles rotation considered for FA, as varimax and quartimax, will only
provide different parameters values due to projections on rotated axis, the
fit quality itself does not change. Overall, this test set supernova portrayed
in Fig. 7 sums up the main behavior observed when comparing FA and PCA fits,
the former generally offers a lower $\chi^{2}/ndof$ and apparently is better
at reproducing finer details. Nevertheless, both methods produce satisfactory
rest-frame model components and we proceed to incorporate both methods with
three components in our cosmology analyses and observe which performs better
when fitting real photometry data and how their different approaches will
affect the cosmology results.

Figure 7: Light curve fitting of the representative test set supernova
Train_SN58, as named by [58]. On the left figure we have three components and
ten components PCA fits, in solid and dashed lines, respectively. On the right
figure we have the three component FA. The $\chi^{2}/ndof$ for PCA versions
are respectively, 5.38 and 0.47, and 5.22 for the FA version.

##  5 Cosmological parameters constraints

The next step is to apply the pure linear expansion model to a larger sample
named Pantheon, introduced by [62], which covers a greater redshift range. In
Section 5.3 we will provide further information about the sample and the
fitting process. For now, we describe how to construct and compare the Hubble
diagrams for PETS and SALT2 models.

The theoretical distance modulus is obtained from the luminosity distance in
units of Hubble distance today through

$\mu_{th}(z;\boldsymbol{\theta},h)=5\log_{10}(\mathcal{D}_{L}(z;\boldsymbol{\theta}))+\mu_{0}(h),$
(11)

with
$\mu_{0}(h)=5\log_{10}(1000c/(\textnormal{km}/\textnormal{s}))-5\log_{10}h$.
Here $h$ is a dimensionless quantity parametrizing the Hubble constant,
$H_{0}=100h$ km s-1 Mpc-1 and $\boldsymbol{\theta}$ stands for the remaining
cosmological parameters for a specific cosmological model. The luminosity
distance in units of Hubble distance today is defined as

$\mathcal{D}_{L}(z;\boldsymbol{\theta}):=\begin{cases}\displaystyle\frac{1+z}{\sqrt{\Omega_{k0}}}\sinh\left(\sqrt{\Omega_{k0}}\int_{0}^{z}\frac{dz^{\prime}}{E(z^{\prime};\boldsymbol{\theta})}\right)&\textnormal{if}\;\Omega_{k0}>0,\\\
(1+z)\displaystyle\int_{0}^{z^{\prime}}\frac{dz^{\prime}}{E(z^{\prime};\boldsymbol{\theta})}&\textnormal{if}\;\Omega_{k0}=0,\\\
\displaystyle\frac{1+z}{\sqrt{-\Omega_{k0}}}\sin\left(\sqrt{-\Omega_{k0}}\int_{0}^{z}\frac{dz^{\prime}}{E(z^{\prime};\boldsymbol{\theta})}\right)&\textnormal{if}\;\Omega_{k0}<0,\\\
\end{cases}$ (12)

where $E(z;\boldsymbol{\theta})$ is the Hubble parameter in units of the
Hubble constant and $\Omega_{k0}$ is the dimensionless curvature parameter
measured today.

The cosmological model considered is $\Lambda$CDM, assuming a universe with
curvature, cold dark matter and cosmological constant, which give us the
following dimensionless Hubble parameter

$E_{\Lambda\textnormal{CDM}}(z;\boldsymbol{\theta})=\sqrt{\Omega_{m0}(1+z)^{3}+\Omega_{\Lambda
0}+\Omega_{k_{0}}(1+z)^{2}}.$ (13)

Recalling that $\Omega_{k_{0}}=1-\Omega_{m_{0}}-\Omega_{\Lambda_{0}}$. Here
$\boldsymbol{\theta}=(\Omega_{m0},\Omega_{\Lambda 0})$, representing the
cosmological parameters we are interested in constraining and $\Omega_{r0}$,
the dimensionless radiation density parameter measured today, is neglected
since it is about 5 orders of magnitudes smaller than the cosmological
parameters $\boldsymbol{\theta}$.

###  5.1 SALT2 distance modulus

The SALT2 empirical distance modulus, that carries information from the
derived light-curve fitting parameters, is given by

$\mu_{\textrm{SALT2}}=m_{B}^{corr}-M_{B}=m_{B}^{*}-M_{B}+\alpha x_{1}-\beta
c,$ (14)

where $m_{B}^{*}$ is the rest-frame B-band magnitude at maximum light, and
according to [39] this quantity can be obtained from the light curve parameter
$x_{0}$ using $m_{B}^{*}=-2.5\log_{10}x_{0}+10.635$. $M_{B}$ is the type Ia
supernovae absolute magnitude, the parameters $x_{1}$ and $c$ corrects the
supernova apparent magnitude to account for the observed dispersion. Lastly,
$\alpha$ and $\beta$ are nuisance parameters that need to be included in the
parameter inference.

The statistical uncertainty matrix for this distance modulus can be portrayed
as

$\sigma^{2}_{\mu_{\textrm{SALT2}}}=\sigma^{2}_{m_{B}^{*}}+\alpha^{2}\sigma_{x_{1}}^{2}+\beta^{2}\sigma_{c}^{2}+\alpha\sigma_{m_{B}^{*},x_{1}}-\beta\sigma_{m_{B}^{*},c}-\alpha\beta\sigma_{x_{1},c}+\sigma^{2}_{\mu,z},$
(15)

where $\sigma^{2}_{m_{B}^{*}}$, $\sigma_{x_{1}}^{2}$, $\sigma_{c}^{2}$,
$\sigma_{m_{B}^{*},x_{1}}$, $\sigma_{m_{B}^{*},c}$ and $\sigma_{x_{1},c}$ are
elements of the covariance matrix for the light curve parameters. The last
term propagates to the distance modulus an error due to the redshift
measurement and due to peculiar velocities contamination. Following [33], we
adopt the distance-redshift relation for an empty universe leading to

$\sigma_{\mu,z}=\sigma_{z}\left(\frac{5}{\log 10}\right)\frac{1+z}{z(1+z/2)},$
(16)

where $\sigma^{2}_{z}=\sigma^{2}_{\textrm{spec}}+\sigma^{2}_{\textrm{pec}}$.
The term $\sigma^{2}_{\textrm{spec}}$ represents the measurement uncertainty
and $\sigma^{2}_{\textrm{pec}}$ is the contribution due to peculiar velocity
uncertainties, estimated as $0.0012$ by [34].

###  5.2 PETS distance modulus

The distance modulus can be written as the difference between apparent and
absolute bolometric magnitudes, $\mu(z)=m(z)-M$. The above description for
SALT2’s distance modulus relies on the empirical observed correlations between
stretch and color with absolute magnitude in B-band. Recognizing the
parameters $x_{1}$ and $c$ proposed in the empirical flux modelling as stretch
and color, respectively, allows to insert those values directly in the
distance modulus as a correction accounting for the absolute magnitude
dispersion, under the assumption these two parameters describe completely the
dispersion seen in Ia SNe light curves.

PETS’ $x_{1}$ and $x_{2}$ parameters can not, by construction, be immediately
recognized as either stretch or color, instead, we derive an expression for
the distance modulus from equation (2), that already contains K-corrections.
The rest-frame flux, as defined in equation (1), can be incorporated into the
observed flux expression (equation (1)). Taking the base 10 logarithm of this
quantity and dividing it by a reference value in the AB magnitude system yield
the apparent magnitude measured by an observer in relative velocity with
respect to the targeted Type Ia supernovae. We can manipulate this expression
to recover a distance modulus that resembles that of SALT2.

If we define

$I_{i}(p;z):=\int_{0}^{\infty}d\lambda
M_{i}(p,\lambda)\lambda(1+z)T_{B}(\lambda(1+z)),$ (17)

we can rewrite the observer-frame flux as

$F_{B}(p(1+z))=x_{0}I_{0}(p;z)\left[1+x_{1}\frac{I_{1}(p;z)}{I_{0}(p;z)}+x_{2}\frac{I_{2}(p;z)}{I_{0}(p;z)}\right],$
(18)

and at B-band maximum light we obtain

$F_{B}(0;z)=x_{0}I_{0}(0;z)\left[1+x_{1}\alpha(z)+x_{2}\beta(z)\right],$ (19)

where we defined two nuisance parameters dependent on redshift

$\alpha(z):=\frac{I_{1}(0;z)}{I_{0}(0;z)},\quad\quad\beta(z):=\frac{I_{2}(0;z)}{I_{0}(0;z)}.$
(20)

Using that $\mu(z)=m_{B}(z)-M_{B}-K_{BB}(z)$ we have

$\displaystyle\mu(z)$
$\displaystyle=-2.5\log_{10}\left[\frac{F_{B}(0;z)}{F_{B,\textnormal{ref}}}\right]-M_{B}$
(21)
$\displaystyle=m_{B}(z)-2.5\log_{10}[1+x_{1}\alpha(z)+x_{2}\beta(z)]-M_{B},$

with

$m_{B}(z)=-2.5\log_{10}\left[\frac{x_{0}I_{0}(0;z)}{F_{B,\textnormal{ref}}}\right].$
(22)

Now we can express the distance modulus as

$\mu(z)=m_{B,\textnormal{corr}}(z)-M_{B},$ (23)

with the corrected apparent magnitude

$m_{B,\textnormal{corr}}(z)=m_{B}(z)-2.5\log_{10}[1+x_{1}\alpha(z)+x_{2}\beta(z)],$
(24)

where $x_{1}$ and $x_{2}$ are the same parameters used to account for
dispersion in the rest-frame flux of type Ia supernovae.

Here, we need to approximate $m_{B}$ for the rest-frame counterpart,
$m_{B}^{*}$, since when performing the integral $I_{0}(0;z)$ as the redshift
grows the overlap between $M_{0}(0;z)$ and the redshifted filter transmission
rapidly goes to zero.

In a first approximation we also neglect the redshift dependence of our
nuisance parameters, $\alpha(z)=\alpha$ and $\beta(z)=\beta$. Afterwards we
will allow for evolution with redshift by solving for a fixed cosmology in a
redshift binned Hubble Diagram.

PETS distance modulus can be written as

$\mu_{\textrm{PETS}}(z)=m_{B}^{*}(z)-2.5\log_{10}(1+\alpha x_{1}+\beta
x_{2})-M_{B},$ (25)

where the redshift dependence of $m_{B}^{*}$ is fully contained in $x_{0}$. We
have the following statistical uncertainty contribution

$\displaystyle\sigma^{2}_{\mu_{\textrm{PETS}}}$
$\displaystyle=\sigma_{m_{B}^{*}}^{2}+\alpha^{2}\gamma^{2}\sigma_{x_{1}}^{2}+\beta^{2}\gamma^{2}\sigma_{x_{2}}^{2}-2\alpha\gamma\sigma_{m_{B},x_{1}}-2\beta\gamma\sigma_{m_{B},c}$
(26)
$\displaystyle+2\alpha\beta\gamma^{2}\sigma_{x_{1},x_{2}}+\sigma_{z,\mu}^{2},$

where the $\gamma$ factor depends only on the nuisance parameters,

$\gamma(\alpha,\beta):=\frac{2.5\log_{10}e}{1+\alpha x_{1}+\beta x_{2}}.$ (27)

We do not expect the parameters other than $m_{B}^{*}$ to carry the same
information as their SALT2 counterparts. In our description
$m_{B}^{*}=-2.5\log_{10}x_{0}+$const is the rest-frame apparent magnitude,
i.e. apparent magnitude measured by an observer with no relative velocity with
respect to the target. The constant term depends on the $M_{0}(p,\lambda)$
surface of the underlying model and also on the shape of the chosen reference
filter transmission, then this offset varies depending on the fitting model.

We also consider a correction due to different host galaxies masses, if
$M_{stellar}>10^{10}M_{\odot}$ we replace in the aforementioned equations
$M_{B}$ for $M_{B}+\Delta_{M}$.

The parameter inference is performed through a Monte Carlo Markov Chain (MCMC)
using the Python package emcee, [15]. Only with type Ia supernovae data we can
not simultaneously constrain $M_{B}$ and $h$, then we define
$\mathcal{M}(M_{B},h)=M_{B}+\mu_{0}(h)$. We can write the $\chi^{2}$ for PETS
model, neglecting correlations between different supernovae, as

$\chi^{2}_{\textrm{PETS}}(\theta,\delta,\mathcal{M})=\sum^{N}_{i}\frac{[\mu_{i,\textrm{PETS}}(z_{i},\delta,M_{B})-\mu_{th}(z_{i};\boldsymbol{\theta},h)]^{2}}{\sigma_{i,\textrm{PETS}}^{2}(\delta)+\sigma_{\textrm{int}}^{2}}.$
(28)

An analogous expression is defined for SALT2. The input log-likelihood is

$\log
L=-0.5\left[\chi^{2}_{\textrm{PETS}}(\boldsymbol{\theta},\mathcal{M},\boldsymbol{\delta},\sigma_{\textrm{int}})+\sum_{i}^{N}\log{[\sigma^{2}_{i,\textrm{PETS}}(\boldsymbol{\theta},\boldsymbol{\delta})+\sigma_{\textrm{int}}^{2}]}\right],$
(29)

where $\boldsymbol{\delta}=(\alpha,\beta)$ are the model nuisance parameters.
And $\sigma_{\textrm{int}}$, the intrinsic scatter, is a free parameter that
will store any remaining unmodeled coherent variability.

###  5.3 Fitting Pantheon sample

To perform the cosmological analysis we chose the Pantheon sample from [62].
This sample consists of 1048 spectroscopically confirmed SNe Ia. In the low
redshift region it includes CfA1-CfA4, [48, 29, 22, 23, 24] and CSP, [14, 12,
64]. Populating the intermediate redshift region we have SDDS, [16, 33, 55],
SNLS [11, 65] and PS1 [46, 61]. And in high-z region we have data from HST,
[66, 49, 50, 18, 53, 51]. This sample is already cross-calibrated and the
light curves were retrieved from SNANA, [34].

The spectral time series model was constructed in SNCosmo using as input the
first three components returned by either our Principal Component Analysis or
Factor Analysis. We restrict our fitting to the region -10 days to 40 days and
3400Å to 7000Å so that the training set would be representative in this
region. We correct for extinction by dust in Milky Way given the color excess,
$E(B-V)_{MW}$, following [13]. The fitting process returns for each supernova
a value of $t_{0}$ (date of maximum light in B-band), $x_{0}$, $x_{1}$ and
$x_{2}$, with corresponding uncertainties.

###  5.4 Fitting results for PCA PETS

From the 1048 SNe in the initial sample, 1017 had a successful fit according
to SNCosmo. However, after visual inspection is possible to note some low-
quality fits. In most cases this happens when there is mainly data before or
after the maximum, increasing the difficulty in constraining the date of
maximum and, consequently, in obtaining accurate values for $x_{1}$ and
$x_{2}$. It is worth to mention an apparent underestimation of uncertainties
for the low redshift sample, leading to higher values of $\chi^{2}/ndof$ even
for visually high quality fits. Hence, we decide to apply cuts to the
parameter space.

We encountered additional challenges when trying to determine the date of
maximum for the low redshift sample. To address this issue, we used a flat
prior $\mathcal{U}(-5,5)$ around the date of maximum reported in the SNANA
files. For the remaining objects we conducted fits for $t_{0}$, $x_{0}$,
$x_{1}$, and $x_{2}$ without using any prior knowledge.

Together with this occasional underestimation of uncertainties, our lack of a
complete model covariance makes the $\chi^{2}/ndof$ parameter a bad quality
fit ruler. The best combination found was to restrict the values of $x_{1}$
and $x_{2}$. The farther from the center of these distributions, the more
likely was for the fit to be of low quality.

The cuts proposed for the PCA PETS parameters are shown in Table 1, where we
can also see how many supernovae are accepted in each step. The cuts on
$\chi^{2}/ndof$, $\sigma_{x_{1}}$ and $\sigma_{x_{2}}$ are purposefully mild
only to eliminate the ones that really distinguish from the underlying
distributions. Fig. 8 shows the distributions for each parameter. In blue the
original sample and in orange the reduced sample, after cuts. Most of the
eliminated supernovas do not fail in all three cuts, many fits with high
values of $x_{1}$ and/or $x_{2}$ have low $\chi^{2}/ndof$, leaving us with a
reduced sample of 691 Ia SNe.

Figure 8: Parameters distribution for PCA PETS fits, in blue the distribution prior to the cuts, in orange after the cuts. The dashed line indicates the applied cuts. Table 1: Cuts applied to the light curve fit parameters for Pantheon sample with PCA PETS. Cuts | # pass cut SNe
---|---
fit | 1017
$-3.5\leq x_{1}\leq 3$ | 820
$|x_{2}|\leq 3$ | 801
$\chi^{2}/\textnormal{ndof}\leq 15$ | 964
$\sigma_{x_{1}}\leq 10$ | 1006
$\sigma_{x_{2}}\leq 10$ | 1003

In Table 2 we have the number of SNe selected per survey and in Fig. 9 we have
the histograms of redshift distributions for each survey. As already expected
from [62] analysis, the fit parameters cuts with respect to the LOWZ sample
affects more strongly CSP, we lost almost 2/3 of the original data, presumably
also due to a great underestimation of uncertainties.

Figure 9: Redshift distribution of the supernovae sample used in the cosmological analysis from PCA PETS model. Table 2: Accepted SNe per survey when applying the cuts for Pantheon subsample using PCA PETS. Survey | Total # of SNe | # pass cut SNe
---|---|---
LOWZ | 172 | 116
SDSS | 313 | 224
SNLS | 236 | 152
PS1 | 271 | 184
HST | 25 | 15
Figure 10: Scatter plots of SALT2 and PCA PETS parameters with the Pearson
correlation coefficients for each pair of parameters.

This same subsample is selected from SALT2 fits which were also performed over
the same light curves using SNCosmo built-in "salt2 v. 2.4". In Fig. 10 we
have a scatter plot for SALT2 parameters versus PCA PETS parameters. To
address their relationships we calculate the Pearson Correlation Coefficient.
This coefficient measures the linear correlation between a pair of parameters.
Values closer to unity indicate the measurement of one feature provides a good
estimate of the second. On the other hand, values closer to zero indicate that
measuring one feature no information about the other is gained, regarding
linear relationships.

We see positive correlations between the pairs $x_{1}$ (SALT2), $x_{1}$ (PCA
PETS) and $c$ (SALT2), $x_{2}$ (PCA PETS). And negative correlations between
the pairs $x_{1}$ (SALT2), $x_{2}$ (PCA PETS) and $c$ (SALT2), $x_{1}$ (PCA
PETS). The stronger the correlation, the strongest is the indication of a
similar origin for those features, i.e. the mechanism which drives the $x_{1}$
(SALT2) for higher or lower values also affects $x_{1}$ (PETS) in the same
direction. As we constructed our model based on a similar hypothesis for pairs
of Ia SNe plus host galaxies without light curve dereddening nor color
extinction separation we do not expect a direct association of PCA PETS and
SALT2 parameters. However, this higher Person coefficients indicate the
Brighter-Broader and Brighter-Bluer correlations are still present but in a
different fashion than in usual Tripp-like estimators.

Models that extract Principal Components from a set of Ia SNe undergo an
approximation to reduce dimensionality, this decision balances the cumulative
explained variance and number of components. In this process we necessarily
discard information, what leads to unmodelled behavior. Our model eliminated
other approximations at the cost of losing direct association with physical
features, albeit through scatter plots (Fig. 10), we verify that our
parameters still carry important information from Ia SNe color and stretch
which are directly correlated with luminosity.

Figure 11: Scatter plot of fit parameters for PCA decomposition components
method as a function of $m_{B}^{*}$, a quantity that correlates linearly with
$\log_{10}(z)$. In blue we have the binned scatter with corresponding standard
deviations.

We can also analyse if the PCA PETS parameters show redshift dependence, this
is portrayed in Fig. 11. For intermediate and high redshift there is no
apparent dependence with redshift. For the low-z sample we observe a deviation
for negative values of $x_{1}$ and positive values of $x_{2}$. However, this
apparent dependence is due to the lower region being populated only by one
sample, here called LOWZ, which favor host galaxies masses, measured through
$\log_{10}(M_{stellar}/M_{\odot})$ higher than 10. This is clear from the
subsample plots in Fig. 12 and in Fig. 19 that shows the same preferences in
the parameters dependence with host galaxy mass, hence reflecting a selection
effect.

Figure 12: Host galaxy mass distribution in selected samples for PCA PETS from
each survey, except HST.

###  5.5 Fitting results for FA PETS

We now analyse the fitting results when equation (1) $M_{i}(p,\lambda)$
surfaces are the FA components. Table 3 shows the applied cuts over the
Pantheon sample after light curve fitting along with the number of supernovas
passing each cut. As for the PCA version, the cuts on $\chi^{2}/ndof$,
$\sigma_{x_{1}}$ and $\sigma_{x_{2}}$ were also purposefully mild.

Table 3: Cuts applied to the light curve fit parameters for Pantheon sample with FA PETS. Cuts | # pass cut SNe
---|---
fit | 1018
$-3.5\leq x_{1}\leq 0.25$ | 888
$|x_{2}|\leq 0.35$ | 898
$\chi^{2}/\textnormal{ndof}\leq 15$ | 994
$\sigma_{x_{1}}\leq 1$ | 998
$\sigma_{x_{2}}\leq 1$ | 999

Since the cuts criteria were different for each training approach, we are not
able to directly compare both resulting tables. Nonetheless, we can
investigate their impact on each survey. From Table 2 and Table 4 we should be
able to assert, that for all surveys (LOWZ, SDSS, SNLS, PS1, HST), the FA
training process shows a better performance, giving a relative number of
passed supernovae higher than PCA ($\text{\\# pass cut SNe}/\text{Total \\# of
SNe}$).

Table 4: Accepted SNe per survey when applying the cuts for Pantheon subsample using FA PETS. Survey | Total # of SNe | # pass cut SNe
---|---|---
LOWZ | 172 | 148
SDSS | 315 | 243
SNLS | 234 | 174
PS1 | 272 | 220
HST | 25 | 20

As the FA components seemed to capture more detailed behavior in the
validation fits, it is expected a better performance when adopting this
description. Indeed, when fitting the Pantheon sample, the parameters
distribution showed less dispersion, with the number of outliers significantly
reduced when comparing to PCA fits.

Figure 13: Scatter plots of SALT2 and FA PETS parameters and the Pearson
correlation coefficients for each pair of parameters.

We can also compare the scatter plots for FA PETS and SALT2 pairs of
parameters. Fig. 13 reveals greater Pearson correlations coefficients between
the two models than seen for the PCA approach. The greatest correlations are
seen for the pairs $x_{2}$ (PETS), $x_{1}$ (SALT2) and $x_{1}$ (PETS), $c$
(SALT2) which can imply that in our model as no derredening was performed, the
first order correction is mainly due to color variations instead of light
curve shape variations.

Comparing Fig. 10 and Fig. 13 we see that in PCA the stretch and color
information seem to be more diluted, with each PETS fit parameter showing a
non-negligible Pearson correlation with SALT2’ parameters. For the FA scatter
plot the situation is different, the stretch and color information seem to be
less diluted but appear in a different order of correction.

Following the analysis performed over the PCA fitting results, an evolution
with $m_{B}^{*}$, or equivalently $\log_{10}z$, is shown in Fig. 14. We see no
apparent behavior with redshift apart from $x_{2}$ tendency to positive values
at low redshift, which is consequence of the host galaxy masses distribution
of the LOWZ sample portrayed in Fig. 23. Summing up, $x_{1}$ and $x_{2}$
parameters show correlations with host galaxy mass when dealing with PCA
components while only $x_{1}$ shows this same correlation when considering FA
components. This reinforces the behavior seen previously that PCA dilutes
color and stretch information, while FA components act on different aspects of
Ia SN variability.

Figure 14: Scatter plot of fit parameters for FA decomposition components
method as a function of $m_{B}^{*}$, a quantity that correlates linearly with
$\log_{10}z$. In blue we have the binned scatter with corresponding standard
deviations.

##  6 Cosmology results for PCA PETS

With light curve fitting results for PCA PETS we perform a cosmological
parameter inference using emcee sampling package [15], with the log-likelihood
described in equation (29) and flat prior over all parameters. The
marginalized contours and one-dimensional marginalized distributions for the
cosmological parameters are shown in Fig. 15. The constraining power of our
PCA PETS model is smaller than SALT2’s but the marginalized one-dimensional
distributions show good agreement inside the 95% confidence region. The
reduced $\chi^{2}$ for the best fit cosmology is 0.97 for PCA PETS while it is
1.012 for SALT2, considering the same subsample.

Table 5 shows the MCMC results for $\Lambda$CDM model, here only including
this subsample of type Ia supernova as cosmological probes. The $\Delta_{M}$
parameters are in agreement within the 95% confidence region, but our PCA PETS
model estimates a lower correction in distance due to host galaxy mass
differences.

Figure 15: Marginalized cosmology results comparison of SALT2 and PCA PETS for
$\Lambda$CDM model. Table 5: Result of marginalized parameters for
$\Lambda$CDM cosmological model for our PCA PETS and FA PETS models as well as
for SALT2 without distance bias correction.

Model | $\Omega_{m0}$ | $\Omega_{\Lambda 0}$ | $\mathcal{M}$ | $\Delta_{M}$ | $\alpha$ | $\beta$ | $\sigma_{\textrm{int}}$
---|---|---|---|---|---|---|---
PCA PETS | 0.34 $\pm$ 0.08 | 0.49 ${}^{+0.13}_{-0.14}$ | 23.93 $\pm$ 0.02 | -0.028 $\pm$ 0.015 | 0.129 $\pm$ 0.005 | 0.124 $\pm$ 0.005 | 0.147 $\pm$ 0.006
SALT2 (PCA subsample) | 0.31 ${}^{+0.06}_{-0.07}$ | 0.54 ${}^{+0.10}_{-0.11}$ | 24.13 $\pm$ 0.01 | -0.043 $\pm$ 0.011 | 0.131 $\pm$ 0.006 | 2.62 $\pm$ 0.07 | 0.094 $\pm$ 0.005
FA PETS | 0.36 $\pm$ 0.07 | 0.48 ${}^{+0.13}_{-0.12}$ | 23.98 $\pm$ 0.02 | -0.030 $\pm$ 0.014 | 1.72 $\pm$ 0.06 | 0.88 $\pm$ 0.05 | 0.145 $\pm$ 0.006
SALT2 (FA subsample) | 0.29 $\pm$ 0.06 | 0.49 $\pm$ 0.10 | 24.14 $\pm$ 0.014 | -0.048 $\pm$ 0.010 | 0.126 $\pm$ 0.006 | 2.65 $\pm$ 0.07 | 0.097 $\pm$ 0.005

It is also worth to mention the slight correlation between the nuisance
parameters seen in PCA PETS that is not observed in SALT2, as shows Fig. 16.
As the principal components are by definition uncorrelated, a pure linear
expansion has uncorrelated components, different from SALT2.

Figure 16: Marginalized contours in $\alpha$ x $\beta$ plane, on the left for
PCA PETS and on the right for SALT2 (PCA subsample).

The intrinsic scatter, $\sigma_{\textrm{int}}$, is a term added in quadrature
to the distance modulus uncertainty of each supernova and is meant to capture
any behavior unaccounted by the model. PCA PETS returns a intrinsic scatter
value 66% higher than SALT2.

Figure 17: Upper panel shows the Hubble Diagram for PCA PETS model. In the
bottom graph is the HD residue as function of redshift.

In the top panel of Fig. 17 we show the Hubble Diagram (HD) for our Pantheon
subsample, corrected for intrinsic dispersion using PCA PETS. The black solid
line is the theoretical distance modulus with parameters fixed at the best fit
displayed in Table 5, and the blue scatter points are drawn from equation (25)
at each supernova redshift, with statistical errors given by the square root
of equation (26) with parameters fixed at their best fits. In the bottom panel
we have the residues for this HD, here defined as
$\mu_{\textrm{PETS}}-\mu_{th}$. A histogram of this residues show a mean of
0.006 mag, indicating no visible rigid translation of the distribution, what
would indicate a bias or an uncounted systematic error. It also shows a
dispersion of 0.19 mag and skewness of 0.21.

In Fig. 18 the HD residues dependence with redshift and fit parameters,
$m_{B}^{*}$, $x_{1}$ and $x_{2}$, are displayed for PCA PETS model. Each gray
scatter point in the background represents a Ia SN from the underlying
subsample, the binning results shown in blue helps to clarify the subgroups
behavior. All parameters show normal dispersion around the solid black line
which represents null residue. For intermediate redshift however, we see a
greater dispersion, that could be related to the higher concentration of data
on this region. Overall, no parameter showed correlation with higher residues.

Figure 18: Hubble Diagram residues for PCA PETS dependence with redshift and
fit parameters.

In Fig. 19 we show the light-curve fit parameters $x_{1}$, $x_{2}$ and the HD
residues dependence with host galaxy mass. The dependence seen here resemble
that observed for SALT2 $x_{1}$ and $c$ parameters reported in [62]. The
binned scatter points show a preference for positive values of $x_{1}$ and
negative values of $x_{2}$ for lower host galaxy masses and the opposite for
higher masses, growing/declining in a linear fashion, while the HD residues
again does not exhibit any dependence.

Figure 19: Fit parameters and HD residue dependence with host galaxy mass for
PCA PETS fits.

Our first description of distance modulus, equation (23), shows redshift
dependence not accommodated by our analysis. Keeping the apparent magnitude as
its rest-frame counterpart but allowing the nuisance parameters and the
intrinsic scatter to evolve with redshift and solving an MCMC simulation for
all cosmological parameters fixed at their best fit values, provides the
results displayed in Fig. 20. The binned results are calculated for groups of
similar number of supernovas, around 40 per group. The nuisance parameters
show a more clear dependence, tending to lower values with increasing
redshift, while the intrinsic scatter show a slower rise with redshift.
However, we need to take into account that the error bars shown here are only
statistical contributions, including systematic terms would reduce this
apparent evolution primarily for the intrinsic scatter. It is important to
note the evolution on our nuisance parameters, must compensate each other when
assuming they have a constant value since the HD residues does not absorb this
redshift evolution.

Figure 20: Evolution of PCA PETS fit parameters with redshifts.

The distance modulus and corresponding uncertainties for SALT2 can be obtained
following the same procedure, using equation (14) and the square root of
equation (15). In Fig. 21 we see no clear dependence with redshift for the
difference $\mu_{SALT2}-\mu_{\textrm{PETS}}$. There is an offset, which is due
to the differences in assumptions when calculating the additive term for
$m_{B}^{*}$, since this value will depend, among others factors, on the shape
of the B-band transmission function. There is also a clear growth in
dispersion for higher redshifts, this region is the main responsible for the
tail in the distribution show on the right, which has skewness of -0.08 and a
standard deviation of 0.15 mag. These results do not include distance bias
corrections, which according to the usual practise requires simulations of
thousand of type Ia Supernovae to reproduce selection effects characteristics
that are specific for each survey. As we chose to compare our results with
SALT2’s sharing the same input Pantheon subsample, we would have to perform a
new set of simulations using BBC (BEAMS with Bias Correction), [32], to
account for the selection bias in our chosen subsample, the same would have to
be performed for each SALT2 subsample used here since PCA and each FA
rotations select a different set of objects, based on cuts applied to each fit
parameters distributions.

Since our main goal on this paper was to evaluate the constraining power of a
pure linear expansion model while observing the behavior of the nuisance
parameters and possible residue correlations in contrast with SALT2 and we do
not have yet a bias correction for our model, we chose to compare both models
in equity, i.e. we do not correct SALT2 results for selection effects. The
correction performed via BBC alters directly the values of $m_{B}^{*}$,
$x_{1}$ and $c$ parameters and amount to a total distance correction of almost
$-0.1$ mag for $z\sim 0.6$. The correction is dependent with redshift and
becomes more visible at $z\sim 0.3$ and further. BBC outputs a distance
modulus for SALT2 since the nuisance parameters are not estimated together
with cosmological parameters, which results in an increase of up to 30%
compared with non-BBC fitting version values. Apart from this effect, the
cosmological parameters are also significantly altered, what we could already
expect from the correction being concentrated mainly on the high redshift
supernovae. After applying BBC, $\Omega_{m0}$ is driven to a slight lower
value while $\Omega_{\Lambda 0}$ increases more significantly.

Figure 21: The left panel shows the dependence with redshift of the difference
$\mu_{\textrm{SALT2}}-\mu_{\textrm{PCA PETS}}$. On the right panel is the
histogram for those differences for all redshifts.

##  7 Cosmology results for FA PETS

We now analyse the cosmological results when considering FA PETS model. As
seen previously, the components retrieved from FA method concentrates in
explaining the off-diagonal elements of the sample covariance matrix and show
less explained variability when fitting the Ia SN SEDs in contrast to PCA
PETS. This slight better performance in explaining variability does not
necessarily translate to a better performance when fitting correlated
photometric data, as also mentioned previously. Indeed, as seen in Section 5.5
the FA PETS fit parameters distributions for the Pantheon sample is better
behaved and we lose less objects when applying the selection cuts.

Figure 22: Marginalized cosmology results comparison of SALT2 and FA PETS for
$\Lambda$CDM model.

Fig. 22 shows the confidence levels in the ($\Omega_{m0}$, $\Omega_{\Lambda
0}$) plane and the given parameters posteriors. In addition to the confidence
plot, the best fit parameters are given in Table 5 for both FA PETS and SALT2
considering the same subsample. Here the reduced $\chi^{2}$ for the best fit
cosmology is 0.96 for FA PETS and 1.013 for SALT2. For FA PETS constraints in
Fig. 22 it is possible to detect a rigid translation in phase space towards
higher $\Omega_{m0}$ and lower $\Omega_{\Lambda 0}$, the same behavior
observed for PCA PETS relative to SALT2.

Figure 23: Fit parameters and HD residue dependence with host galaxy mass for
FA fits.

Fig. 23 shows the dependence of the HD residues and light curve fit parameters
with the host galaxies mass. For $x_{1}$ there is no apparent correlation
while for $x_{2}$ there is a similar behavior previously seen for PCA PETS,
now with a more steep slope. This reinforces the differences between the
components of both methods, with the principal components and common factors
clearly carrying distinct information. Regarding the HD residues, no host
galaxy mass dependence is observed.

Apart from these distinctions, performing the MCMC analysis for redshift
binned nuisance parameters and $\sigma_{\textrm{int}}$ with the cosmological
parameters fixed at their best fit values, we encounter a similar redshift
evolution previously seen for PCA PETS. HD residues do not show dependencies
with redshift nor fit parameters. FA PETS also predicts compatible results of
$\Omega_{m0}$ and $\Omega_{\Lambda 0}$ parameters with SALT2 predictions, with
a higher estimate of coherent intrinsic scatter, $\sigma_{\textnormal{int}}$.
On top of that, a similar correlation of $\alpha\times\beta$ is verified and
no redshift dependent behaviour is observed for $\mu_{\textnormal{FA
PETS}}-\mu_{\textnormal{SALT2}}$.

##  8 Conclusion

In this paper we investigated two different models describing the rest-frame
flux of type Ia supernova as pure linear expansions, called PCA PETS and FA
PETS. The pure linear expansion rest-frame flux was originally proposed by
[58] and here investigated in greater details and applied to $\Lambda$CDM
cosmology. We obtained the model components through Principal Component
Analysis and Factor Analysis after a robust reconstruction of the training
sample SEDs using 1-d Gaussian Process Regression for each plane of constant
wavelength with [26] surface as the mean. This reconstruction provided a SED
resolution of 10Å $\times$ 1 day. Following, we derived directly from the
rest-frame flux definition an expression for the model distance modulus and
fitted the Pantheon sample to constrain the $\Lambda$CDM cosmological
parameters.

The driven idea was to take advantage of Ia SNe similarities and adopt an
alternative description using decomposition methods over a training set of
SEDs. Unlike the usual rest-frame flux description that includes an
exponential term to model two different color-related effects, reddening due
to dust and intrinsic color variations, here a pure expansion description was
applied, removing the need to approximate the color variation effect by an
extinction law. We verified this pure expansion is a reasonable assumption
when constructing a light curve fitter and both PCA and FA methods showed
showed reliable light-curve modeling results. The cosmological constraints for
PETS favours higher values of $\Omega_{m0}$ and slightly lower values of
$\Omega_{\Lambda 0}$, when compared to SALT2. With respect to this matter,
further investigation of PETS model systematics is needed to conclude if this
behaviour is related to the rest-frame flux description or if the distance
bias correction would reconcile the best fits. In this description we lost the
link between the fit parameters with the effects of color and stretch, but we
were able to encounter a Tripp-like formula when bringing our model to
observer-frame with minimal approximations.

Our analysis over PCA and FA PETS cosmological results showed inferences of
cosmological parameters within 68% agreement with the most widely used fitter
SALT2 results, with our models predicting more conservative statistical
uncertainties over the cosmological parameters. The Hubble Diagram residues
showed no evolution neither with redshift nor with host galaxy mass for both
PCA and FA models. And the residues w.r.t. SALT2 distance modulus predictions
show no dependence with redshift. In future works we plan on extending the
model coverage to include UV region and better understand the model
covariance. We also plan on including higher redshift Ia SNe data, when
available, to re-access the nuisance parameters evolution with redshift and
apply the effects of distance bias selection to our analysis.

Overall, both PCA and FA PETS models performance will heavily depend on how
representative the training set is. These methods undergo an approximation
when keeping only the first two terms of the rest-frame flux expansion, as the
others empirical descriptions of Ia SN rest-frame fluxes. The main differences
reside in what kind of information is being neglected. In order to avoid
correlating dust extinction reddening with luminosity, we decided not to
separate the reddening effect from our SEDs, but at the same time this raised
questions whether important stretch and color information is being neglected.
We argue that the fitting and constraining power of both PCA PETS and FA PETS
are good indicators that this information is not being completely neglected
but would instead be more diluted over the model components.

It is also important to access some aspects of the training procedure.
SNFactory DR9 provides good quality spectra at a number of different epochs
for 171 spectroscopically confirmed nearby Ia SNe and our resulting models are
trained over 164 SEDs reconstructed with this spectra. Even with this high
quality data, the SEDs reconstruction using gaussian process regression will
many times recover wrong early phase behavior if we do not consider a
supernova template. On top of that, underestimated uncertainties can also lead
to poor regressions. These wrong predictions translate to bad reconstructed
SEDs that can result in nonphysical behavior recovered by the principal
components or by the common factors. To obtain a closer description to the
reality and correctly recognize the aspects most common to the underlying
group of Ia SNe, the training procedure depends not only on the amount of good
spectra quality but also on the SED reconstruction procedure.

The high explained variance in the first component on both methods can
indicate a lack of diversity in our training sample. Increasing the sample
will help increase the sample variability completion. We should also be
attentive to the difference in $\Delta_{M}$ constraints, our FA PETS and PCA
PETS obtained lower values than the SALT2 corresponding analysis. This result
can be due to the different color treatments but for now the difference is not
significant to make a firm statement. It is important to note we only
considered a coherent intrinsic scatter which entered as a free parameter in
the cosmological analysis. Other intrinsic scatter should be explored, giving
further insight whether the neglected information was wavelength-dependent. In
this scenario the best intrinsic scatter modeling may not be wavelength
dependent for a model that does not correlates extinction and luminosity
variations.

Regarding $\alpha$ and $\beta$ evolution with redshift seen in Fig. 20, as
these nuisance parameters do not describe the correlations strengths, the
evolution with redshift does not offer extra information around a possible
type Ia SN evolution with redshift. Using oblique FA PETS will allow for a
deeper analysis relating evolution. It is important to note this behaviour
does not seem to impact the HD residues and inclusion of systematics
uncertainties is likely to diminish the evolution seen for
$\sigma_{\textrm{int}}$ also in Fig. 20.

Mapping the systematics is essential to the current cosmological precision era
and therefore supports improvements in current empirical models. FA PETS can
be further explored if we consider common factor rotations. These rotations do
not improve the fitting quality but can provide a different insight about the
hidden variables. It is possible that a specific oblique rotation, that
correlates the factors, lead to a set of hidden variables describing directly
some intrinsic variation effects, like color and stretch, allowing us to
recover physical interpretation for the model components. It is important to
note that even though these two quantities are known to correlate with
luminosity, a complete description of the explosion mechanism behind Ia SNe is
not completely understood, so other hidden variables not included by most
fitters, may play an important role in these empirical descriptions.

## Acknowledgements

We would like to express our gratitude to Martin Makler for his generosity and
support throughout this collaboration. We would also like to thank Kyle
Barbary and Kyle Boone for their helpful contributions towards resolving
SNCOSMO GitHub issues, which significantly eased the fitting analysis process.

This work has also been supported by the Brazilian funding agencies CAPES and
CNPq. CSN thanks the support of CNPq for the PhD scholarship no.
155994/2019-0. JPCF thanks the Brazilian funding agencies CAPES for the MS
scholarship no. 88887.336370/2019-00 and CNPq for the PhD scholarship no.
140210/2021-0. RRRR acknowledges CNPq (grant no. 309868/2021-1).

## Data Availability

The type Ia supernova spectra used in this paper to construct the PETS models
are from Data Release 9 from Nearby Supernova Factory publicly available in
https://snfactory.lbl.gov/. The light curves used for cosmology are from
SNANA, available for download in https://zenodo.org/record/4015325. Our codes
for training, light curve fitting and cosmology will be available after
publication at
https://github.com/CassiaNascimento/PETS_model_for_SN_Ia_LC_fitting.

## SOFTWARE CITATIONS

This work uses the following software packages:

  * •

Astropy ([67, 10]).

  * •

Emcee ([15]).

  * •

GetDist ([38]).

  * •

Matplotlib ([27]).

  * •

Quadpy ([60]).

  * •

Scipy ([70]).

  * •

SNANA ([34]).

  * •

SNCosmo ([2]).

  * •

Scikit-learn ([6, 41]).

  * •

Pillow ([31]).

  * •

Python ([69]).

## References

  * [1] Greg Aldering et al.  “Overview of the Nearby Supernova Factory”  In _Survey and Other Telescope Technologies and Discoveries_ , 4836  SPIE, 2002, pp. 61
bibrangedash 72  International Society for OpticsPhotonics  DOI:
10.1117/12.458107

  * [2] Kyle Barbary et al.  “SNCosmo” If you use SNCosmo, please cite it using these metadata.  In _Zenodo_ , 2022  DOI: 10.5281/zenodo.7117347
  * [3] David Barber  “Bayesian reasoning and machine learning”  Cambridge University Press, 2012
  * [4] Betoule, M. et al.  “Improved cosmological constraints from a joint analysis of the SDSS-II and SNLS supernova samples”  In _A &A_ 568, 2014, pp. A22  DOI: 10.1051/0004-6361/201423413
  * [5] Christopher M. Bishop  “Pattern Recognition and Machine Learning”, 1613-9011  Springer, 2006
  * [6] Lars Buitinck et al.  “API design for machine learning software: experiences from the scikit-learn project”  In _ECML PKDD Workshop: Languages for Data Mining and Machine Learning_ , 2013, pp. 108–122
  * [7] Jason A. Cardelli, Geoffrey C. Clayton and John S. Mathis  “The relationship between infrared, optical, and ultraviolet extinction”  In _The Astrophysical Journal_ 345, 1989  DOI: 10.1086/167900
  * [8] M. Childress et al.  “HOST GALAXIES OF TYPE Ia SUPERNOVAE FROM THE NEARBY SUPERNOVA FACTORY”  In _The Astrophysical Journal_ 770.2 The American Astronomical Society, 2013, pp. 107  DOI: 10.1088/0004-637X/770/2/107
  * [9] Chotard, N. et al.  “The reddening law of type Ia supernovae: separating intrinsic variability from dust using equivalent widths”  In _Astronomy and Astrophysics_ 529, 2011, pp. L4  DOI: 10.1051/0004-6361/201116723
  * [10] The Astropy Collaboration et al.  “The Astropy Project: Building an Open-science Project and Status of the v2.0 Core Package*”  In _The Astronomical Journal_ 156.3 The American Astronomical Society, 2018, pp. 123  DOI: 10.3847/1538-3881/aabc4f
  * [11] A. Conley et al.  “SUPERNOVA CONSTRAINTS AND SYSTEMATIC UNCERTAINTIES FROM THE FIRST THREE YEARS OF THE SUPERNOVA LEGACY SURVEY*”  In _The Astrophysical Journal Supplement Series_ 192.1 The American Astronomical Society, 2010, pp. 1  DOI: 10.1088/0067-0049/192/1/1
  * [12] Carlos Contreras et al.  “THE CARNEGIE SUPERNOVA PROJECT: FIRST PHOTOMETRY DATA RELEASE OF LOW-REDSHIFT TYPE Ia SUPERNOVAE*”  In _The Astronomical Journal_ 139.2 The American Astronomical Society, 2010, pp. 519  DOI: 10.1088/0004-6256/139/2/519
  * [13] Edward L. Fitzpatrick  “Correcting for the Effects of Interstellar Extinction”  In _Publications of the Astronomical Society of the Pacific_ 111, 1999  DOI: 10.1086/316293
  * [14] Gastón Folatelli et al.  “THE CARNEGIE SUPERNOVA PROJECT: ANALYSIS OF THE FIRST SAMPLE OF LOW-REDSHIFT TYPE-Ia SUPERNOVAE*”  In _The Astronomical Journal_ 139.1 The American Astronomical Society, 2009, pp. 120  DOI: 10.1088/0004-6256/139/1/120
  * [15] Daniel Foreman-Mackey, David W. Hogg, Dustin Lang and Jonathan Goodman  “emcee: The MCMC Hammer”  In _Publications of the Astronomical Society of the Pacific_ 125.925, 2013, pp. 306  DOI: 10.1086/670067
  * [16] Joshua A. Frieman et al.  “THE SLOAN DIGITAL SKY SURVEY-II SUPERNOVA SURVEY: TECHNICAL SUMMARY”  In _The Astronomical Journal_ 135.1 The American Astronomical Society, 2007, pp. 338  DOI: 10.1088/0004-6256/135/1/338
  * [17] GPy  “GPy: A Gaussian process framework in python”, http://github.com/SheffieldML/GPy, 2012
  * [18] O. Graur et al.  “TYPE-Ia SUPERNOVA RATES TO REDSHIFT 2.4 FROM CLASH: THE CLUSTER LENSING AND SUPERNOVA SURVEY WITH HUBBLE”  In _The Astrophysical Journal_ 783.1 The American Astronomical Society, 2014, pp. 28  DOI: 10.1088/0004-637X/783/1/28
  * [19] Guy, J. et al.  “SALT2: using distant supernovae to improve the use of type Ia supernovae as distance indicators ***”  In _A &A_ 466.1, 2007, pp. 11–21  DOI: 10.1051/0004-6361:20066930
  * [20] Guy, J. et al.  “The Supernova Legacy Survey 3-year sample: Type Ia supernovae photometric distances and cosmological constraints”  In _A &A_ 523, 2010, pp. A7  DOI: 10.1051/0004-6361/201014468
  * [21] Shiyuan He, Lifan Wang and Jianhua Z. Huang  “Characterization of Type Ia Supernova Light Curves Using Principal Component Analysis of Sparse Functional Data”  In _The Astrophysical Journal_ 857.2 The American Astronomical Society, 2018, pp. 110  DOI: 10.3847/1538-4357/aab0a8
  * [22] Malcolm Hicken et al.  “CfA3: 185 TYPE Ia SUPERNOVA LIGHT CURVES FROM THE CfA”  In _The Astrophysical Journal_ 700.1 The American Astronomical Society, 2009, pp. 331  DOI: 10.1088/0004-637X/700/1/331
  * [23] Malcolm Hicken et al.  “IMPROVED DARK ENERGY CONSTRAINTS FROM 100 NEW CfA SUPERNOVA TYPE Ia LIGHT CURVES”  In _The Astrophysical Journal_ 700.2 The American Astronomical Society, 2009, pp. 1097  DOI: 10.1088/0004-637X/700/2/1097
  * [24] Malcolm Hicken et al.  “CfA4: LIGHT CURVES FOR 94 TYPE Ia SUPERNOVAE”  In _The Astrophysical Journal Supplement Series_ 200.2 The American Astronomical Society, 2012, pp. 12  DOI: 10.1088/0067-0049/200/2/12
  * [25] Harold Hotelling  “Analysis of a complex of statistical variables into principal components”  In _Journal of Educational Psychology_ 24.6 Warwick & York, 1933, pp. 417  DOI: 10.1037/h0071325
  * [26] E. Y. Hsiao et al.  “K-Corrections and Spectral Templates of Type Ia Supernovae”  In _The Astrophysical Journal_ 663.2, 2007, pp. 1187  DOI: 10.1086/518232
  * [27] J. D. Hunter  “Matplotlib: A 2D graphics environment”  In _Computing in Science & Engineering_ 9.3 IEEE COMPUTER SOC, 2007, pp. 90–95  DOI: 10.1109/MCSE.2007.55
  * [28] Ivezić et al.  “LSST: From Science Drivers to Reference Design and Anticipated Data Products”  In _The Astrophysical Journal_ 873.2 The American Astronomical Society, 2019, pp. 111  DOI: 10.3847/1538-4357/ab042c
  * [29] Saurabh Jha et al.  “UBVRI Light Curves of 44 Type Ia Supernovae”  In _The Astronomical Journal_ 131.1, 2006, pp. 527  DOI: 10.1086/497989
  * [30] I. T. Jolliffe  “Principal Component Analysis”  Springer, 2002
  * [31] Hugo Kemenade et al.  “python-pillow/Pillow: 9.0.1”  In _Zenodo_ , 2022  DOI: 10.5281/zenodo.5953590
  * [32] R. Kessler and D. Scolnic  “Correcting Type Ia Supernova Distances for Selection Biases and Contamination in Photometrically Identified Samples”  In _The Astrophysical Journal_ 836.1 The American Astronomical Society, 2017, pp. 56  DOI: 10.3847/1538-4357/836/1/56
  * [33] Richard Kessler et al.  “FIRST-YEAR SLOAN DIGITAL SKY SURVEY-II SUPERNOVA RESULTS: HUBBLE DIAGRAM AND COSMOLOGICAL PARAMETERS”  In _The Astrophysical Journal Supplement Series_ 185.1 The American Astronomical Society, 2009, pp. 32  DOI: 10.1088/0067-0049/185/1/32
  * [34] Richard Kessler et al.  “SNANA: A Public Software Package for Supernova Analysis”  In _Publications of the Astronomical Society of the Pacific_ 121.883 University of Chicago Press, 2009, pp. 1028  DOI: 10.1086/605984
  * [35] A. G. Kim et al.  “Standardizing Type Ia Supernova Absolute Magnitudes Using Gaussian Process Data Regression”  In _The Astrophysical Journal_ 766.2, 2013, pp. 84  DOI: 10.1088/0004-637X/766/2/84
  * [36] Kevin Krisciunas et al.  “The Carnegie Supernova Project. I. Third Photometry Data Release of Low-redshift Type Ia Supernovae and Other White Dwarf Explosions”  In _The Astronomical Journal_ 154.5 The American Astronomical Society, 2017, pp. 211  DOI: 10.3847/1538-3881/aa8df0
  * [37] Léget, P.-F. et al.  “SUGAR: An improved empirical model of Type Ia supernovae based on spectral features”  In _A &A_ 636, 2020, pp. A46  DOI: 10.1051/0004-6361/201834954
  * [38] Antony Lewis  “GetDist: a Python package for analysing Monte Carlo samples”  In _arXiv e-prints_ , 2019, pp. arXiv:1910.13970  DOI: 10.48550/arXiv.1910.13970
  * [39] J. Mosher et al.  “COSMOLOGICAL PARAMETER UNCERTAINTIES FROM SALT-II TYPE IA SUPERNOVA LIGHT CURVE MODELS”  In _The Astrophysical Journal_ 793.1 The American Astronomical Society, 2014, pp. 16  DOI: 10.1088/0004-637X/793/1/16
  * [40] Karl Pearson  “LIII. On lines and planes of closest fit to systems of points in space”  In _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ 2.11 Taylor & Francis, 1901, pp. 559–572  DOI: 10.1080/14786440109462720
  * [41] F. Pedregosa et al.  “Scikit-learn: Machine Learning in Python”  In _Journal of Machine Learning Research_ 12, 2011, pp. 2825–2830
  * [42] S. Perlmutter et al.  “Measurements of $\Omega$ and $\Lambda$ from 42 High-Redshift Supernovae”  In _The Astrophysical Journal_ 517.2, 1999, pp. 565  DOI: 10.1086/307221
  * [43] Mark M Phillips  “The absolute magnitudes of Type IA supernovae”  In _The Astrophysical Journal_ 413, 1993, pp. L105–L108  DOI: 10.1086/186970
  * [44] J. D. R. Pierel et al.  “Extending Supernova Spectral Templates for Next-generation Space Telescope Observations”  In _Publications of the Astronomical Society of the Pacific_ 130.993 The Astronomical Society of the Pacific, 2018, pp. 114504  DOI: 10.1088/1538-3873/aadb7a
  * [45] Carl Edward Rasmussen and Christopher K. I. Williams  “Gaussian Processes for Machine Learning”  The MIT Press, 2005  DOI: https://doi.org/10.7551/mitpress/3206.001.0001
  * [46] A. Rest et al.  “COSMOLOGICAL CONSTRAINTS FROM MEASUREMENTS OF TYPE Ia SUPERNOVAE DISCOVERED DURING THE FIRST 1.5 yr OF THE Pan-STARRS1 SURVEY”  In _The Astrophysical Journal_ 795.1 The American Astronomical Society, 2014, pp. 44  DOI: 10.1088/0004-637X/795/1/44
  * [47] Adam G. Riess et al.  “Observational Evidence from Supernovae for an Accelerating Universe and a Cosmological Constant”  In _The Astronomical Journal_ 116.3, 1998, pp. 1009  DOI: 10.1086/300499
  * [48] Adam G. Riess et al.  “BVRI Light Curves for 22 Type Ia Supernovae”  In _The Astronomical Journal_ 117.2, 1999, pp. 707  DOI: 10.1086/300738
  * [49] Adam G. Riess et al.  “Type Ia Supernova Discoveries at z $>$ 1 from the Hubble Space Telescope: Evidence for Past Deceleration and Constraints on Dark Energy Evolution*”  In _The Astrophysical Journal_ 607.2, 2004, pp. 665  DOI: 10.1086/383612
  * [50] Adam G. Riess et al.  “New Hubble Space Telescope Discoveries of Type Ia Supernovae at z $\geq$ 1: Narrowing Constraints on the Early Behavior of Dark Energy*”  In _The Astrophysical Journal_ 659.1, 2007, pp. 98  DOI: 10.1086/510378
  * [51] Adam G. Riess et al.  “New Parallaxes of Galactic Cepheids from Spatially Scanning the Hubble Space Telescope: Implications for the Hubble Constant”  In _The Astrophysical Journal_ 855.2 The American Astronomical Society, 2018, pp. 136  DOI: 10.3847/1538-4357/aaadb7
  * [52] Rigault, M. et al.  “Strong dependence of Type Ia supernova standardization on the local specific star formation rate”  In _A &A_ 644, 2020, pp. A176  DOI: 10.1051/0004-6361/201730404
  * [53] Steven A. Rodney et al.  “TYPE Ia SUPERNOVA RATE MEASUREMENTS TO REDSHIFT 2.5 FROM CANDELS: SEARCHING FOR PROMPT EXPLOSIONS IN THE EARLY UNIVERSE”  In _The Astronomical Journal_ 148.1 The American Astronomical Society, 2014, pp. 13  DOI: 10.1088/0004-6256/148/1/13
  * [54] B. M. Rose et al.  “Synergies between Vera C. Rubin Observatory, Nancy Grace Roman Space Telescope, and Euclid Mission: Constraining Dark Energy with Type Ia Supernovae”  In _arXiv e-prints_ , 2021, pp. arXiv:2104.01199  DOI: 10.48550/arXiv.2104.01199
  * [55] Masao Sako et al.  “The Data Release of the Sloan Digital Sky Survey-II Supernova Survey”  In _Publications of the Astronomical Society of the Pacific_ 130.988 The Astronomical Society of the Pacific, 2018, pp. 064002  DOI: 10.1088/1538-3873/aab4e0
  * [56] Michele Sasdelli et al.  “A metric space for Type Ia supernova spectra”  In _Monthly Notices of the Royal Astronomical Society_ 447.2, 2014, pp. 1247–1266  DOI: 10.1093/mnras/stu2416
  * [57] Michele Sasdelli et al.  “Breaking the colour-reddening degeneracy in Type Ia supernovae”  In _Monthly Notices of the Royal Astronomical Society_ 460.1, 2016, pp. 373–382  DOI: 10.1093/mnras/stw900
  * [58] C. Saunders et al.  “SNEMO: Improved Empirical Models for Type Ia Supernovae”  In _The Astrophysical Journal_ 869.2 The American Astronomical Society, 2018, pp. 167  DOI: 10.3847/1538-4357/aaec7e
  * [59] David J. Schlegel, Douglas P. Finkbeiner and Marc Davis  “Maps of Dust Infrared Emission for Use in Estimation of Reddening and Cosmic Microwave Background Radiation Foregrounds”  In _The Astrophysical Journal_ 500.2, 1998, pp. 525  DOI: 10.1086/305772
  * [60] Nico Schlömer et al.  “nschloe/quadpy: None”  In _Zenodo_ , 2021  DOI: 10.5281/zenodo.5218525
  * [61] D. Scolnic et al.  “SYSTEMATIC UNCERTAINTIES ASSOCIATED WITH THE COSMOLOGICAL ANALYSIS OF THE FIRST PAN-STARRS1 TYPE Ia SUPERNOVA SAMPLE”  In _The Astrophysical Journal_ 795.1 The American Astronomical Society, 2014, pp. 45  DOI: 10.1088/0004-637X/795/1/45
  * [62] D. M. Scolnic et al.  “The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample”  In _The Astrophysical Journal_ 859.2 The American Astronomical Society, 2018, pp. 101  DOI: 10.3847/1538-4357/aab9bb
  * [63] D Spergel et al.  “Wide-field infrarred survey telescope-astrophysics focused telescope assets WFIRST-AFTA 2015 report”  In _arXiv preprint arXiv:1503.03757_ , 2015
  * [64] Maximilian D. Stritzinger et al.  “THE CARNEGIE SUPERNOVA PROJECT: SECOND PHOTOMETRY DATA RELEASE OF LOW-REDSHIFT TYPE Ia SUPERNOVAE”  In _The Astronomical Journal_ 142.5 The American Astronomical Society, 2011, pp. 156  DOI: 10.1088/0004-6256/142/5/156
  * [65] M. Sullivan et al.  “SNLS3: CONSTRAINTS ON DARK ENERGY COMBINING THE SUPERNOVA LEGACY SURVEY THREE-YEAR DATA WITH OTHER PROBES”  In _The Astrophysical Journal_ 737.2 The American Astronomical Society, 2011, pp. 102  DOI: 10.1088/0004-637X/737/2/102
  * [66] N. Suzuki et al.  “THE HUBBLE SPACE TELESCOPE CLUSTER SUPERNOVA SURVEY. V. IMPROVING THE DARK-ENERGY CONSTRAINTS ABOVE z $>$ 1 AND BUILDING AN EARLY-TYPE-HOSTED SUPERNOVA SAMPLE*”  In _The Astrophysical Journal_ 746.1 The American Astronomical Society, 2012, pp. 85  DOI: 10.1088/0004-637X/746/1/85
  * [67] The Astropy Collaboration et al.  “Astropy: A community Python package for astronomy”  In _A &A_ 558, 2013, pp. A33  DOI: 10.1051/0004-6361/201322068
  * [68] Michael E Tipping and Christopher M Bishop  “Mixtures of probabilistic principal component analyzers”  In _Neural computation_ 11.2 MIT Press One, 1999, pp. 443–482  DOI: 10.1162/089976699300016728
  * [69] Guido Van Rossum and Fred L. Drake  “Python 3 Reference Manual”  Scotts Valley, CA: CreateSpace, 2009
  * [70] Pauli Virtanen et al.  “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python”  In _Nature Methods_ 17, 2020, pp. 261–272  DOI: 10.1038/s41592-019-0686-2