© ESO, 2012

1. Introduction

Observations of type I supernovae (SNe Ia) directly established in the late 1990s that the universe is accelerating (Riess et al. 1998; Perlmutter et al. 1999). This most remarkable result is possibly only surpassed by the discovery in the 1930s that the universe expands. The SNe Ia are still the backbone for the prevailing ΛCDM concordance model, which is also corroborated by the combination of other probes, such as cosmic microwave background anisotropies, baryon acoustic oscillations and galaxy clustering (Komatsu et al. 2011; Eisenstein et al. 2005; Percival et al. 2010; Vikhlinin et al. 2009; Mantz et al. 2010).

Already in the landmark papers from the two original surveying groups, a legitimate concern was expressed about not properly accounted for systematic effects. At that time, however, there were few (on the order of tens) observed SNe Ia in any of the available samples, so the uncertainties were statistically limited; the main task was gathering new, larger, and uniform datasets. For the recent compilations and surveys (Wood-Vasey et al. 2007; Kowalski et al. 2008; Hicken et al. 2009; Conley et al. 2011; Kessler et al. 2009; Amanullah et al. 2010), and more so for even future ones (Pan-STARRS; DES; LSST), we are no longer sample-limited; the urgency is again on the physics of the SNe Ia phenomenon and all the aspects that impact it (for review see, for instance, Kirshner 2010; Howell 2011; Goobar & Leibundgut 2011). Despite the major importance of investigating the systematics, here we focus on a distinct aspect, namely the statistically consistent determination of the parameters when the variance is to be estimated as well.

In contrast to the traditional χ² method in most of the literature, we here discuss a different approach to parameter fitting with unknown variance, which is based on the likelihood method. We recall that, if we aim to estimate both the covariance and the mean of a Gaussian process, the ordinary (uncorrected) χ² approach (or any iterative recipe thereof, for that matter) cannot be straightforwardly applied, because we might miss a nontrivial term in the objective function to be extremized. This is necessary when the covariance itself depends on free parameters of the underlying model. Recently, similar criticisms to the traditional χ² method were presented by Kim (2011, see also Vishwakarma & Narlikar 2010). Kim (2011) used a likelihood approach to simulated data for the distance moduli, and he focused on determining the intrinsic dispersion, without reference to a particular light-curve fitter. We reconsider the problem in the context of the MLCS2K2 and SALT2 light-curve fitters with real data. We point out that, in principle, with SALT2, additional difficulties may arise owing to the astrophysical parameters. We compare the results of the traditional χ² treatment to those of the likelihood method. We restrict ourselves to SNe Ia datasets only, without taking into account other important sources of information (such as CMB, BAOs or clusters), to make our point pristine, by avoiding the masking due to any other such probes.

The structure of the paper is as follows. In Sect. 2, we briefly recall the general procedure to obtain the final estimated parameters from the raw data. In Sect. 3, we describe the two most widely used “fitter pipelines” (MLCS2k2 and SALT2), which use the traditional χ² approach. In Sect. 4, we critically review the aforementioned traditional approach and present the likelihood-based method. In Sect. 5, the main numerical results are shown for both fitters for data from the first-year SDSS-II Supernova Survey together with other supernova data sets (ESSENCE, SNLS, HST, and a compilation of nearby type Ia supernovae), which constitute the sample combination “e” as described in Kessler et al. (2009). Finally, in Sect. 6, we end with some discussions and conclusions.

2. Light-curve fitting

As mentioned in Sect. 1, SNe Ia are standardizable candles; this means that there is a phenomenological recipe according to which the raw light curves, after being subjected to a transformation by an N_st-parameter function, furnish a new set of so-called standardized light curves; in other words the dispersion in the new magnitudes is considerably lower than in the original input set. This phenomenological recipe is not unique at all: there are several light-curve fitters in the literature (Jha et al. 2007; Wang et al. 2003, 2005; Guy et al. 2005, 2007; Conley et al. 2008; Burns et al. 2011). Here we will exploit some features of the two most common ones: the multicolor light-curve shape (MLCS2k2) (Jha et al. 2007) one and the spectral adaptive light-curve template (SALT2) (Guy et al. 2007).

In a nutshell:

• The MLCS2k2 fitting model (Jhaet al. 2007) describes the variationamong SNe Ia light curves with a single parameter (Δ). The MLCS2k2 theoretical magnitude, observed in an arbitrary filter Y, at an epoch γ, is given by (Kessler et al. 2009) ${\mathit{m}}_{\mathit{Y,\gamma }}^{\mathrm{th}}\mathrm{=}{\mathit{M}}_{{\mathit{Y}}^{\mathrm{\prime }}\mathit{,\gamma }}\mathrm{+}{\mathit{p}}_{{\mathit{Y}}^{\mathrm{\prime }}\mathit{,\gamma }}\mathrm{\Delta }\mathrm{+}{\mathit{q}}_{{\mathit{Y}}^{\mathrm{\prime }}\mathit{,\gamma }}{\mathrm{\Delta }}^{\mathrm{2}}\mathrm{+}{\mathit{K}}_{{\mathit{Y}}^{\mathrm{\prime }}\mathit{Y,\gamma }}\mathrm{+}\mathit{\mu }\mathrm{+}{\mathit{X}}_{{\mathit{Y}}^{\mathrm{\prime }}\mathit{,\gamma }}^{\mathrm{host}}\mathrm{+}{\mathit{X}}_{\mathit{Y}}^{\mathrm{MW}}\mathit{,}$(1)where Y′ ∈  { U,B,V,R,I }  is one of the supernova rest-frame filters for which the model is defined, Δ is the MLCS2k2 shape-luminosity parameter that accounts for the correlation between peak luminosity and the shape/duration of the light curve, ${\mathit{X}}_{{\mathit{Y}}^{\mathrm{\prime }}\mathit{,\gamma }}^{\mathrm{host}}$ is the host-galaxy extinction, ${\mathit{X}}_{\mathit{Y}}^{\mathrm{MW}}$ is the Milky Way extinction, K_Y′Y,γ is the K-correction between rest-frame and observer-frame filters, and μ is the distance modulus. The coefficients M_Y′,γ, p_Y′,γ, and q_Y′,γ are model vectors that have been evaluated using nearly 100 well observed low-redshift SNe Ia as a training set. Above, γ = 0 labels quantities at the B-band peak magnitude epoch. Fitting the model to the magnitudes of each SN Ia, usually fixing R_V, gives μ, Δ, A_V and t₀, the B-band peak magnitude epoch; for more details, see (Kessler et al. 2009).

• The SALT2 fitter (Guy et al. 2007) uses a two-dimensional surface in time and wavelength that describes the temporal evolution of the rest-frame spectral energy distribution (SED) or specific flux for SNe Ia. The original model was trained on a set of combined light-curves and 303 spectra, not only from (very) nearby but also medium- and high-redshift SNe Ia. In SALT2, the rest-frame specific flux at wavelength λ and phase (time) p (p = 0 at B-band maximum) is modeled by $\mathit{\phi }\mathrm{\left(}\mathit{p,\lambda }\mathrm{;}{\mathit{x}}_{\mathrm{0}}\mathit{,}{\mathit{x}}_{\mathrm{1}}\mathit{,c}\mathrm{\right)}\mathrm{=}{\mathit{x}}_{\mathrm{0}}\mathrm{\left[}{\mathit{M}}_{\mathrm{0}}\mathrm{\right(}\mathit{t,\lambda }\mathrm{)}\mathrm{+}{\mathit{x}}_{\mathrm{1}}{\mathit{M}}_{\mathrm{1}}\mathrm{(}\mathit{t,\lambda }\mathrm{\left)}\mathrm{\right]}\exp \mathrm{\left[}\mathit{cC}\mathrm{\right(}\mathit{\lambda }\mathrm{\left)}\mathrm{\right]}$(2)and does depend through the parameters x₀,x₁, and c on the particular type Ia supernova. M₀(t,λ), M₁(t,λ), and C(λ) are determined from the training process described in Guy et al. (2007). M₀(t,λ), M₁(t,λ) are the zeroth and the first moments of the distribution of training sample SEDs. One might consider adding moments of higher order to Eq. (2). To compare with photometric SNe Ia data, the observer-frame flux in passband Y is calculated as ${\mathit{F}}^{\mathit{Y}}\mathrm{\left(}\mathit{p}\mathrm{\right(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{\left)}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathrm{\int }\mathrm{d}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{\left[}{\mathit{\lambda }}^{\mathrm{\prime }}\mathit{\phi }\mathrm{\right(}\mathit{p,}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{\left)}{\mathit{T}}^{\mathit{Y}}\mathrm{\right(}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}{\mathrm{\left)}}^{\mathrm{\right]}}\mathit{,}$(3)where T^Y(λ) defines the transmission curve of the observer-frame filter Y, and z is the redshift. Each SN Ia light curve is fitted separately using Eqs. (2) and (3) to determine the parameters x₀, x₁, and c with corresponding errors. However, the SALT2 light-curve fit does not yield an independent distance modulus estimate for each SN Ia. As we will see in the next section, the distance moduli are estimated as part of a global parameter fit to an ensemble of SN Ia light curves in which cosmological parameters and global SN Ia properties are also obtained.

In the next section we will discuss how to obtain constraints on cosmological parameters using MLCS2k2 and SALT2 output quantities as our data.

3. The traditional χ² approach

The prevailing SNe Ia cosmological analysis is based on the χ² function: ${\mathit{\chi }}^{\mathrm{2}}\mathrm{:}\mathrm{=}{{X}}^{\mathrm{T}}{\mathbf{\Sigma }}^{-1}{X}\mathit{,}$(4)where X: = (μ − μ_th(z,Θ)), μ is the set of distance moduli derived from the light-curve fitting procedure for each SNe Ia event, at redshifts given by z, μ_th(z,Θ) is the theoretical prediction for them, given in terms of a vector Θ of parameters and Σ is the covariance matrix of the events.

3.1. The χ² approach from SALT2 output

The SALT2 light-curve fitter gives three quantities, with corresponding errors, to be used in the analysis of cosmology: ${\mathit{m}}_{\mathit{B}}^{\mathrm{\ast }}\mathrm{:}\mathrm{=}\mathrm{-}\mathrm{2.5}\log \left[{\mathit{x}}_{\mathrm{0}}\mathrm{\int }\mathrm{d}{\mathit{\lambda }}^{\mathrm{\prime }}{\mathit{M}}_{\mathrm{0}}\mathrm{(}\mathit{p}\mathrm{=}\mathrm{0}\mathit{,}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{)}{\mathit{T}}^{\mathit{B}}\mathrm{\left(}{\mathit{\lambda }}^{\mathrm{\prime }}\mathrm{\right)}\right]\mathit{,}$(5)to be interpreted as the peak rest-frame magnitude in the B band, x₁, a parameter related to the stretch of the light-curve and c, related to the color of the supernova alongside the redshift z of the supernova. The distance modulus is modeled, in this context, as a function of (${\mathit{m}}_{\mathit{B}}^{\mathrm{\ast }}$, x₁, c) and two new parameters, δ: = (α,β), plus the peak absolute magnitude, in B band, M_B. Defining the corrected magnitude as ${\mathit{m}}_{\mathit{B}}^{\mathrm{corr}}\mathrm{\left(}{\delta }\mathrm{\right)}\mathrm{:}\mathrm{=}{\mathit{m}}_{\mathit{B}}^{\mathrm{\ast }}\mathrm{+}\mathit{\alpha }{\mathit{x}}_{\mathrm{1}}\mathrm{-}\mathit{\beta c,}$(6)we can write $\mathit{\mu }\mathrm{\left(}{\delta }\mathit{,}{\mathit{M}}_{\mathit{B}}\mathrm{\right)}\mathrm{=}{\mathit{m}}_{\mathit{B}}^{\mathrm{corr}}\mathrm{\left(}{\delta }\mathrm{\right)}\mathrm{-}{\mathit{M}}_{\mathit{B}}\mathit{.}$(7)Assuming that all SN Ia events are independent, one can rewrite Eq. (4) as ${\mathit{\chi }}_{\mathrm{SALT}\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}{\theta }\mathit{,}{\delta }\mathit{,}\mathrm{ℳ}\mathrm{\right(}{\mathit{M}}_{\mathit{B}}\mathit{,h}\mathrm{\left)}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{\mathrm{\left[}{\mathit{\mu }}_{\mathit{i}}\mathrm{\right(}{\delta }\mathit{,}{\mathit{M}}_{\mathit{B}}\mathrm{)}\mathrm{-}{\mathit{\mu }}_{\mathrm{th}}\mathrm{(}{\mathit{z}}_{\mathit{i}}\mathrm{;}{\theta }\mathit{,h}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}\mathrm{\left(}{\delta }\mathrm{\right)}\mathrm{+}{\mathit{\sigma }}_{\mathrm{int}}^{\mathrm{2}}}\mathit{,}$(8)where N is the number of SNe Ia in the sample, θ denotes the cosmological parameters other than h, with the present value of the Hubble parameter given by H₀ = 100   hkm   s^-1   Mpc^-1. The theoretical distance modulus is given by ${\mathit{\mu }}_{\mathrm{th}}\mathrm{\left(}\mathit{z}\mathrm{;}{\theta }\mathit{,h}\mathrm{\right)}\mathrm{=}\mathrm{5}\log \mathrm{\left[}{\mathrm{𝒟}}_{\mathrm{L}}\mathrm{\right(}\mathit{z}\mathrm{;}{\theta }\mathrm{\left)}\mathrm{\right]}\mathrm{+}{\mathit{\mu }}_{\mathrm{0}}\mathrm{\left(}\mathit{h}\mathrm{\right)}\mathit{,}$(9)with ${\mathit{\mu }}_{\mathrm{0}}\mathrm{\left(}\mathit{h}\mathrm{\right)}\mathrm{:}\mathrm{=}\mathrm{5}\log \left(\frac{{\mathrm{10}}^{\mathrm{3}}\mathit{c}\mathit{/}\mathrm{\left(}\mathrm{km}\hspace{0.17em}\mathrm{s}-1\mathrm{\right)}}{\mathit{h}}\right)\mathrm{\simeq }\mathrm{42.38}\mathrm{-}\mathrm{5}\log \mathit{h}\mathit{.}$(10)The dimensionless luminosity distance (in units of the Hubble distance today), , for comoving observers in a Robertson-Walker universe is given by $\begin{array}{ccc}& & {\mathrm{𝒟}}_{\mathrm{L}}\mathrm{\left(}\mathit{z}\mathrm{;}{\theta }\mathrm{\right)}\mathrm{=}\\ & & \{\begin{array}{c}\\ & \\ & \\ & \end{array}\end{array}$(11)where Ω_k0 is the “curvature density parameter”, such that it is proportional to the three-curvature, and $\mathit{E}\mathrm{\left(}\mathit{z}\mathrm{;}{\theta }\mathrm{\right)}\mathrm{:}\mathrm{=}\mathit{H}\mathrm{\left(}\mathit{z}\mathrm{;}{\theta }\mathit{,h}\mathrm{\right)}\mathit{/}{\mathit{H}}_{\mathrm{0}}$(12)is the dimensionless Hubble parameter. As indicated in Eq. (8), the χ² function depends on the parameters M_B and h only through their combination $\mathrm{ℳ}\mathrm{\left(}{\mathit{M}}_{\mathit{B}}\mathit{,h}\mathrm{\right)}\mathrm{:}\mathrm{=}{\mathit{M}}_{\mathit{B}}\mathrm{+}{\mathit{\mu }}_{\mathrm{0}}\mathrm{\left(}\mathit{h}\mathrm{\right)}\mathit{.}$(13)It may thus be thought of as effectively directly dependent on only the parameters θ, δ and ℳ.

A floating dispersion term, σ_int, is added in quadrature to the distance modulus dispersion, which is given by $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}\mathrm{\left(}{\delta }\mathrm{\right)}& \mathrm{=}& {\mathit{\sigma }}_{{\mathit{m}}_{\mathit{B}}^{\mathrm{\ast }}\mathit{,i}}^{\mathrm{2}}\mathrm{+}{\mathit{\alpha }}^{\mathrm{2}}{\mathit{\sigma }}_{{\mathit{x}}_{\mathrm{1}}\mathit{,i}}^{\mathrm{2}}\mathrm{+}{\mathit{\beta }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{c,i}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{\alpha }{\mathit{\sigma }}_{{\mathit{m}}_{\mathit{B}}^{\mathrm{\ast }}\mathit{,}{\mathit{x}}_{\mathrm{1}}\mathit{,i}}\mathrm{-}\mathrm{2}\mathit{\beta }{\mathit{\sigma }}_{{\mathit{m}}_{\mathit{B}}^{\mathrm{\ast }}\mathit{,c,i}}\\ & & \end{array}$(14)where ${\mathit{\sigma }}_{\mathit{\mu ,z,i}}^{\mathrm{2}}$ is the contribution to the distance modulus dispersion owing to redshift uncertainties from peculiar velocities and also from the measurement itself. Following Kessler et al. (2009), we will model this contribution for simplicity with the distance-redshift relation for an empty universe, which gives ${\mathit{\sigma }}_{\mathit{\mu ,z,i}}\mathrm{=}{\mathit{\sigma }}_{\mathit{z,i}}\left(\frac{\mathrm{5}}{\ln \mathrm{10}}\right)\frac{\mathrm{1}\mathrm{+}{\mathit{z}}_{\mathit{i}}}{{\mathit{z}}_{\mathit{i}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{z}}_{\mathit{i}}\mathit{/}\mathrm{2}\mathrm{)}}\mathit{,}$(15)with ${\mathit{\sigma }}_{\mathit{z,i}}^{\mathrm{2}}\mathrm{=}{\mathit{\sigma }}_{\mathrm{spec}\mathit{,i}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{pec}}^{\mathrm{2}}\mathit{,}$(16)where σ_spec,i is the redshift measurement error, and σ_pec = 0.0012 is the redshift uncertainty due to peculiar velocity.

As advocated by some groups (Astier et al. 2006), minimizing Eq. (8) gives rise to a bias towards increasing values of α and β. To circumvent this feature, an iterative method is employed (Astier et al. 2006).

In this iterative method, the χ² presented in Eq. (8) is replaced by ${\mathit{\chi }}_{\mathrm{SALT}\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}{\theta }\mathit{,}{\delta }\mathit{,}\mathrm{ℳ}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{\mathrm{\left[}{\mathit{\mu }}_{\mathit{i}}\mathrm{\right(}{\delta }\mathit{,}{\mathit{M}}_{\mathit{B}}\mathrm{)}\mathrm{-}{\mathit{\mu }}_{\mathrm{th}}\mathrm{(}{\mathit{z}}_{\mathit{i}}\mathrm{;}{\theta }\mathit{,h}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}\mathrm{\left(}{\eta }\mathrm{\right)}\mathrm{+}{\mathit{\sigma }}_{\mathrm{int}}^{\mathrm{2}}}\mathrm{·}$(17)In this expression, η is not a parameter of the ${\mathit{\chi }}_{\mathrm{SALT}\mathrm{2}}^{\mathrm{2}}$. To obtain the best-fit values for the parameters, η is given initial values and the optimization is performed on θ,δ and ℳ. After this step, η is updated with the best-fit value of δ and the optimization is performed again. The process continues until a convergence is achieved, which means that η does not change under the required precision.

During this process σ_int is not considered as a free parameter to be optimized, being determined rather by the following procedure: start with a guess value (usually σ_int = 0.15). Perform the iterative procedure described above. The value of σ_int is then obtained by fine-tuning it so that the reduced χ² equals unity (with all the other parameters fixed on their best-fit values). The iterative procedure is repeated once more with this new value and the final best-fit values are obtained. It is important to note that the value of σ_int affects both the best fit and the confidence levels of the parameters, since it changes the weight given to each supernova in the χ² (cf. Eq. (17)).

3.2. The χ² approach from MLCS2k2 output

The MLCS2k2 light-curve fitter gives a direct estimation of the distance modulus. In this context, the analogue of Eq. (8) is ${\mathit{\chi }}_{\mathrm{MLCS}\mathrm{2}\mathrm{k}\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}{\theta }\mathit{,h}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}}^{\mathit{N}}\frac{\mathrm{[}{\mathit{\mu }}_{\mathit{i}}\mathrm{-}{\mathit{\mu }}_{\mathrm{th}}\mathrm{(}{\mathit{z}}_{\mathit{i}}\mathrm{;}{\theta }\mathit{,h}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{int}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathit{\mu ,z,i}}^{\mathrm{2}}}\mathit{,}$(18)where σ_i is the distance modulus dispersion as given by MLCS2k2.

The procedure to obtain σ_int is similar to the one described in the previous subsection, however, here we use only a subsample of nearby SNe Ia and not the full one, as for the SALT2 analysis. After setting up the value of σ_int, we minimize the ${\mathit{\chi }}_{\mathrm{MLCS}\mathrm{2}\mathrm{k}\mathrm{2}}^{\mathrm{2}}$ using the full SNe Ia sample to obtain the best-fit values for θ and h.

4. The proposed likelihood approach

Considering the SNe Ia light-curve fitting parameters as Gaussian distributed random variables, we propose to take as starting point the likelihood $\mathit{L}\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathit{N}}\det \hspace{0.17em}\mathbf{\Sigma }}}\exp \mathrm{(}\mathrm{-}{{X}}^{\mathrm{T}}{\mathbf{\Sigma }}^{-1}{X}\mathit{/}\mathrm{2}\mathrm{)}\mathit{,}$(19)which is related to the χ² in Eq. (4) by $\mathrm{ℒ:}\mathrm{=}\mathrm{-}\mathrm{2}\ln \mathit{L}\mathrm{=}{\mathit{\chi }}^{\mathrm{2}}\mathrm{+}\ln \det \hspace{0.17em}\mathbf{\Sigma }\mathrm{+}\mathit{N}\ln \mathrm{\left(}\mathrm{2}\mathit{\pi }\mathrm{\right)}\mathit{.}$(20)Equations (19) and (20) are the single basis upon which our whole statistical procedure lies. When the full covariance of the problem is known, minimizing χ² is completely equivalent to minimizing ℒ. However, this is not the case for current SNe Ia observations and neglecting the last but one term in Eq. (20) would, in principle, lead to a biased result. Our proposal is to minimize the following functions for each case discussed in the previous section $\begin{array}{ccc}{\mathrm{ℒ}}_{\mathrm{SALT}\mathrm{2}}\mathrm{\left(}{\theta }\mathit{,}{\delta }\mathit{,}\mathrm{ℳ}\mathit{,}{\mathit{\sigma }}_{\mathrm{int}}\mathrm{\right)}& \mathrm{=}& {\mathit{\chi }}_{\mathrm{SALT}\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}{\theta }\mathit{,}{\delta }\mathit{,}\mathrm{ℳ}\mathit{,}{\mathit{\sigma }}_{\mathrm{int}}\mathrm{\right)}\\ & & \mathrm{+}\sum _{\mathit{i}}^{\mathit{N}}\ln \mathrm{\left(}{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}\mathrm{\right(}{\delta }\mathrm{)}\mathrm{+}{\mathit{\sigma }}_{\mathrm{int}}^{\mathrm{2}}\mathrm{)}\end{array}$(21)and $\begin{array}{ccc}{\mathrm{ℒ}}_{\mathrm{MLCS}\mathrm{2}\mathrm{k}\mathrm{2}}\mathrm{\left(}{\theta }\mathit{,h,}{\mathit{\sigma }}_{\mathrm{int}}\mathrm{\right)}& \mathrm{=}& {\mathit{\chi }}_{\mathrm{MLCS}\mathrm{2}\mathrm{k}\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}{\theta }\mathit{,h,}{\mathit{\sigma }}_{\mathrm{int}}\mathrm{\right)}\\ & & \mathrm{+}\sum _{\mathit{i}}^{\mathit{N}}\ln \mathrm{(}{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{int}}^{\mathrm{2}}\mathrm{)}\mathit{,}\end{array}$(22)where we neglected parameter-independent terms. ${\mathit{\chi }}_{\mathrm{SALT}\mathrm{2}}^{\mathrm{2}}$ (θ,δ,ℳ,σ_int) and ${\mathit{\chi }}_{\mathrm{MLCS}\mathrm{2}\mathrm{k}\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}{\theta }\mathit{,h,}{\mathit{\sigma }}_{\mathrm{int}}\mathrm{\right)}$ are given, respectively, by Eqs. (8) and (18) now considering σ_intalso as a free parameter. With this procedure, we can directly obtain unbiased probability distributions functions for all parameters, including σ_int and δ.

5. Results

In this section we compare the results obtained from the χ² and the likelihood approaches, as described in Sects. 3 and 4, using real data from the first year SDSS-II Supernova Survey together with other supernova data sets (ESSENCE, SNLS, HST, and a compilation of nearby supernovae) as described in Kessler et al. (2009).

Fig. 1 Contours of 68% and 95% confidence level in the plane Ωm0 × ΩΛ0 for the ΛCDM model. We marginalized over all other parameters with flat prior. The solid (dashed) lines are the results for the likelihood (χ2) approach. The black circle (gray square) is the best fit value for the likelihood (χ2). Left panel: results for the MLCS2k2 version of the SDSS compilation. Right panel: results for the SALT2 version of the SDSS compilation.

Fig. 2 Contours of 68% and 95% confidence level in the plane Ωm0 × w for the FwCDM model. We marginalized over all other parameters with flat prior. The solid (dashed) lines are the results for the likelihood (χ2) approach. The black circle (gray square) is the best fit value for the likelihood (χ2). Left panel: results for the MLCS2k2 version of the SDSS compilation. Right panel: results for the SALT2 version of the SDSS compilation.

Fig. 3 Contours of 68% and 95% confidence level in the plane β × α. We marginalized over all other parameters with flat prior. The solid (dashed) lines are the results for the likelihood (χ2) approach. The black circle (gray square) is the best fit value for the likelihood (χ2). Left panel: results for the ΛCDM model with the SDSS compilation. Right panel: results for the FwCDM model with the SDSS compilation.

Fig. 4 Distributions for σint for the ΛCDM model from the likelihood approach. We marginalized over all other parameters with flat prior. The horizontal lines depict the 68%, 95% and 99% confidence intervals. The dashed line is the result for the χ2 approach. Left panel: results for the MLCS2k2 version of the SDSS compilation. Right panel: results for the SALT2 version of the SDSS compilation.

To perform the comparison, we considered the following cosmological models:

• ΛCDM, in which we can write the Friedmann equation in terms of the parameters θ = (Ω_m0,Ω_k0) as ${\mathit{E}}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{;}{\theta }\mathrm{\right)}\mathrm{=}{\mathrm{\Omega }}_{\mathit{m}\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{3}}\mathrm{+}{\mathrm{\Omega }}_{\mathit{k}\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{\Omega }}_{\mathit{m}\mathrm{0}}\mathrm{-}{\mathrm{\Omega }}_{\mathit{k}\mathrm{0}}\mathrm{)}\mathit{.}$(23)

• FwCDM, described by θ = (Ω_m0,w) and ${\mathit{E}}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{;}{\theta }\mathrm{\right)}\mathrm{=}{\mathrm{\Omega }}_{\mathit{m}\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{3}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{\Omega }}_{\mathit{m}\mathrm{0}}\mathrm{-}{\mathrm{\Omega }}_{\mathit{k}\mathrm{0}}\mathrm{\left)}\mathrm{\right(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{3}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{w}\mathrm{)}}\mathit{.}$(24)

We chose the ΛCDM and FwCDM models to directly compare the best-fit and the 68% and 95% confidence contours for the parameters δ, σ_int and θ, for both SALT2 and MLCS2k2 data. The best-fit values were obtained with the MIGRAD minimization of the Minuit (James & Roos 1975) implementation in ROOT (Antcheva et al. 2009) and the probability distributions were obtained with Monte Carlo Markov Chains (MCMC). We considered as the probability distributions, in the context of χ² approach, the following functions: $\begin{array}{ccc}{\mathrm{𝒫}}_{\mathrm{SALT}\mathrm{2}}\mathrm{\left(}{\theta }\mathit{,}{\delta }\mathit{,}\mathrm{ℳ}\mathrm{\right)}& \mathrm{=}& {\mathit{N}}_{\mathrm{SALT}\mathrm{2}}\exp \mathrm{[}\mathrm{-}{\mathit{\chi }}_{\mathrm{SALT}\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}{\theta }\mathit{,}{\delta }\mathit{,}\mathrm{ℳ}\mathit{,}{\mathit{\sigma }}_{\mathrm{int}}\mathrm{\right)}\mathit{/}{\mathrm{2}}^{\mathrm{]}}\mathit{,}\\ {\mathrm{𝒫}}_{\mathrm{MLCS}\mathrm{2}\mathrm{k}\mathrm{2}}\mathrm{\left(}{\theta }\mathit{,h}\mathrm{\right)}& \mathrm{=}& {\mathit{N}}_{\mathrm{MLCS}\mathrm{2}\mathrm{k}\mathrm{2}}\exp \mathrm{[}\mathrm{-}{\mathit{\chi }}_{\mathrm{MLCS}\mathrm{2}\mathrm{k}\mathrm{2}}^{\mathrm{2}}\mathrm{\left(}{\theta }\mathit{,h,}{\mathit{\sigma }}_{\mathrm{int}}\mathrm{\right)}\mathit{/}{\mathrm{2}}^{\mathrm{]}}\mathit{,}\end{array}$where the normalization factors N_SALT2 and N_MLCS2k2 are independent of the parameters to be estimated. Note that for the traditional χ² method σ_int is fixed so the probability distribution does not depend on it.

In Figs. 1 and 2 we show the confidence contours for the parameters θ for ΛCDM and FwCDM models, respectively. For these models, we note that the differences between the best fit and the area of the contours for the χ² and likelihood approaches are more significant when the MLCS2k2 fitter is used. Indeed, for the SALT2 fitter, the differences are not significant (less than 13%) – see also Fig. 4 and discussion below. Whether this is a general feature or depends on the models or dataset used has to be further investigated.

In Fig. 3 we show the confidence contours for the SALT2 parameters δ for both ΛCDM and FwCDM models. We can see that there is no significant difference in the constraints for α. For β we find a bias that is, however, small compared to the 68% confidence interval for this parameter.

In Fig. 4 we show the probability distributions for σ_int, given by the likelihood approach, for MLCS2k2 and SALT2 data. The traditional χ² approach gives only one value for this parameter without uncertainty and we represent it by the dashed vertical line in the figure. We can see that the discrepancy between the value obtained from the χ² approach and the best-fit value obtained from the likelihood approach is greater for the MLCS2k2 data. The results are incompatible at more than 99% confidence interval, which does not happen for the SALT2 data. This can possibly be due to the fact that σ_int is obtained using only a nearby sample in the χ² approach for MLCS2k2, while the whole sample is used in the likelihood approach. This issue deserves additional investigation and will be the subject of future work.

Fig. 5 Evolution of SALT2 parameters α, β, ℳ, and σint with redshift for the SDSS compilation for model ΛCDM. The cosmological parameters were kept fixed at their global best-fit values. The black circles and the gray squares are the results for the likelihood and the χ2 approaches, respectively. The error bars represent the 68% confidence intervals, marginalizing over all other parameters with flat priors.

We also allow for a possible variation of the parameters α, β, ℳ and σ_int with redshift, in the context of the SALT2 data. To perform such analysis the dataset was divided into redshift bins and the cosmological parameters were fixed at the best-fit values obtained from the global fit, then releasing α, β, ℳ and σ_int to be determined in each bin. The results are shown in Fig. 5. We found evidence of evolution for the parameter β, in agreement with Kessler et al. (2009); furthermore we also found evidence of evolution for σ_int, which might support the use of a variable σ_int instead of a constant one.

In addition to the SN Ia sample combination “e” from Kessler et al. (2009), we also used the Union2 compilation (Amanullah et al. 2010) and the SNLS third-year data (Guy et al. 2010). However, these references do not provide the covariances between the SALT2 parameters (${\mathit{m}}_{\mathit{B}}^{\mathrm{\ast }}$,x₁,c). To make a fair comparison with the SDSS data, we repeated the analyses for all three samples, this time neglecting the cross terms in Eq. (14). Without taking into account the covariances, the differences in the areas for the SDSS compilation are essentially the same. The main difference is that, for the FwCDM model, the constraints are tighter for the χ² approach than for the likelihood one. For the Union 2 compilation, the difference in the areas ranged from 13% to 14%, while the SNLS data showed the greater discrepancies, reaching up to  ~ 53%.

6. Conclusion

In this work we considered data from the first-year SDSS-II Supernova Survey together with other supernova data sets (ESSENCE, SNLS, HST, and a compilation of nearby supernovae) and compared two different approaches (traditional χ² and likelihood) to determine constraints from SNe Ia processed by two of the most used light-curve fitters in the literature, the MLCS2k2 and the SALT2. The MLCS2k2 gives an estimation of the distance modulus for each SN Ia, with its corresponding variance, which can be directly compared with the model prediction for this quantity. With SALT2 we can only have an estimate of the distance modulus depending on parameters to be obtained simultaneously with the cosmological ones. Furthermore, in the SNe Ia analysis, it is common to introduce a residual, unknown, contribution σ_int to the variance which, in the traditional approach, is determined by imposing that the reduced χ² be unity, when considering the full sample, in the SALT2 case, or only nearby SNe Ia, in the MLCS2k2 case.

By comparing the results obtained from the traditional χ² approach with those of the likelihood method, we showed that for current data and chosen cosmological models, there is a small difference in the best-fit values and confidence contours (~30% difference in area) (cf. Figs. 1 and 2) if the MLCS2k2 fitter is adopted. For SALT2 the difference is less significant ( ≲ 13% difference in areas). In both cases the likelihood approach gives more restrictive constraints. We can understand these results by observing, from Fig. 4, that the estimated value of σ_int in the traditional approach, is higher than the peak of the σ_int distribution obtained with the likelihood method. Indeed, for MLCS2k2 the σ_int value obtained with the traditional χ² approach is outside the 99% confidence level of the distribution obtained with the likelihood method. For SALT2 it is on the order of 68% confidence level and this might explain why for SALT2 the differences between these two approaches are less significant. We also remark that the covariance between the SALT2 parameters (${\mathit{m}}_{\mathit{B}}^{\mathrm{\ast }}\mathit{,}{\mathit{x}}_{\mathrm{1}}\mathit{,c}$) has an important role in the above result since we obtained, for the 68% confidence level contour, an area ratio of 1.13, when considering the covariance, and of 0.95, when neglecting the covariance.

We also studied the possible evolution of the SALT2 parameters α, β and ℳ with the redshift, splitting the samples in redshift bins and performing the fit separately for each one. In this case, we found evidence of evolution in the parameter β and also in σ_int.

While this paper was in preparation two articles appeared that addressed the issue of using χ² with unknown variances. The first (Kim 2011) focused on the question of determining σ_int by applying the likelihood approach to simulated data for the distance modulus, in which case there is no need to specify the light-curve fitter. In contrast, we used real data here and compared the constraints arising from the two most common light-curve fitters. The second paper (March et al. 2011) employed SALT2 and proposed a more sophisticated, Bayesian analysis instead of the likelihood method proposed here and also by Kim (2011). Furthermore, we did not find in our work, with the SDSS first year compilation, any kind of catastrophic biases for parameter estimation when adopting the likelihood approach, a possibility suggested in March et al. (2011).

In summary, we used current data to compare the traditional χ² and the likelihood approaches to determine best fit and confidence regions from SNe Ia. We argued that when the variance is not completely known, minimizing the traditional χ² is not formally equivalent to maximizing the likelihood function, since the normalization of the likelihood, assumed to be Gaussian, is also a function of parameters to be determined. We conclude by suggesting to adopt the likelihood framework instead of the traditional χ² one, since it is more straightforward, numerically more efficient, and self-consistent.

Acknowledgments

We thank Masao Sako for helpful discussions and suggestions. B.L.L. and I.W. also thank CNPq, Brazil, and S. E. J. also thanks ICTP, Italy, for support.

References

All Figures

Fig. 1 Contours of 68% and 95% confidence level in the plane Ωm0 × ΩΛ0 for the ΛCDM model. We marginalized over all other parameters with flat prior. The solid (dashed) lines are the results for the likelihood (χ2) approach. The black circle (gray square) is the best fit value for the likelihood (χ2). Left panel: results for the MLCS2k2 version of the SDSS compilation. Right panel: results for the SALT2 version of the SDSS compilation.

In the text

Fig. 2 Contours of 68% and 95% confidence level in the plane Ωm0 × w for the FwCDM model. We marginalized over all other parameters with flat prior. The solid (dashed) lines are the results for the likelihood (χ2) approach. The black circle (gray square) is the best fit value for the likelihood (χ2). Left panel: results for the MLCS2k2 version of the SDSS compilation. Right panel: results for the SALT2 version of the SDSS compilation.

In the text

Fig. 3 Contours of 68% and 95% confidence level in the plane β × α. We marginalized over all other parameters with flat prior. The solid (dashed) lines are the results for the likelihood (χ2) approach. The black circle (gray square) is the best fit value for the likelihood (χ2). Left panel: results for the ΛCDM model with the SDSS compilation. Right panel: results for the FwCDM model with the SDSS compilation.

In the text

	Fig. 4 Distributions for σint for the ΛCDM model from the likelihood approach. We marginalized over all other parameters with flat prior. The horizontal lines depict the 68%, 95% and 99% confidence intervals. The dashed line is the result for the χ2 approach. Left panel: results for the MLCS2k2 version of the SDSS compilation. Right panel: results for the SALT2 version of the SDSS compilation.
In the text

Fig. 5 Evolution of SALT2 parameters α, β, ℳ, and σint with redshift for the SDSS compilation for model ΛCDM. The cosmological parameters were kept fixed at their global best-fit values. The black circles and the gray squares are the results for the likelihood and the χ2 approaches, respectively. The error bars represent the 68% confidence intervals, marginalizing over all other parameters with flat priors.

In the text