A MORE GENERAL MODEL FOR THE INTRINSIC SCATTER IN TYPE Ia SUPERNOVA DISTANCE MODULI

Type Ia supernovae (SNe Ia) have been extensively used as "standard candles" to measure astronomical distances. However, they are often called "standardizable candles" to emphasize the fact that the light-curve shape and color can be used to correct the apparent magnitudes, resulting in a more accurate measurement of the distance. The correlation between light-curve shape and magnitude was demonstrated by Phillips et al. (1993), who parameterized the light-curve shape with the parameter Δm₁₅, the decay in B-band magnitude 15 days after the peak. Different parameterizations of the light-curve shape parameter have been used: Δ for MLCS2k2 (Jha et al. 2007), x₁ for Spectral Adaptive Lightcurve Template, version 2 (SALT2; Guy et al. 2007, 2010), and s_B for SiFTO (Conley et al. 2008), but the parameter is essentially determined empirically to produce the tightest correlation between light-curve shape and distance. This work will use the SuperNova ANAlysis (SNANA) software package (Kessler et al. 2009b) implementation of the SALT2 model that uses x₁ as the light-curve shape parameter.

SN Ia color depends on the light-curve shape parameter, but SNe also exhibit variations in color for fixed light-curve shape. The source of color variation is not well understood, but it is widely suspected that variations in the explosion process, host galaxy dust, and the circumstellar environment all play a role. The difference in measured color and the mean color for a given x₁ is denoted by the parameter c. Regardless of the process that produces different colors, it may be expected that the apparent B-band magnitude (m_B) will show some dependence on color, and that an appropriate color correction should improve the precision of SN as standard candles. These considerations have led to the widespread practice of fitting SN light curves with a one parameter family of templates plus a color law to describe the color.

The determination of color is straightforward given multi-band photometric measurements and a color law, and the result is traditionally given as E(B − V). However, the amount of correction to apply to a given SN has been less clear. The correction has often been expressed as an extinction correction, which is related to the color through the parameter R_V = A_V/E(B − V), where A_V is the extinction in V-band. Many efforts to determine this correction have relied on the parameterization of Cardelli, Clayton, and Mathis (Cardelli et al. 1989), which we will call the CCM model, where the parameter R_V changes both the shape of the extinction curve as a function of wavelength and the amount of extinction for a given color excess E(B − V). The simplest approach (Jha et al. 2007) is to fix R_V = 3.1, the average value determined from the studies of extinction in the Milky Way (Sneden et al. 1978). However, many studies of SN light-curve data have indicated that smaller values of R_V provide a better description of the SN data. In the CCM model, the single parameter R_V can be measured by estimating A_V or by comparing the relative extinction in different filters (by considering E(B − V) − E(V − R), for example).

The SALT2 model parameterizes the SN rest-frame spectrum S(λ, p) as a function of wavelength (λ) and phase (p) according to

$$\begin{equation} S(\lambda,p)=x_0 \, [S_0(\lambda,p)+x_1 S_1(\lambda,p)] \, e^{[-c{\rm CL}(\lambda)]}, \end{equation} \tag{ 1 }$$

where the fitted parameters x₀, x₁, and c are the overall scale, the light-curve shape parameter, and the color, respectively. The functions S₀ and S₁ are fixed by a training sample of SN light curves and spectra. The SALT2 light-curve fit does not rely on the CCM parameterization, but uses a color law CL(λ) determined empirically from SN data. The SALT2 parameters are determined by fitting the measured light-curve data for each SN to the synthetic magnitudes calculated from the redshifted spectrum given by Equation (1).

The light-curve parameters x₀, x₁, and c specify a model SN spectrum, but an additional step is needed to use the light-curve parameters to calculate a standard SN magnitude. One can calculate the apparent magnitude m_B, where m_B is the apparent B-band synthesized from the rest-frame spectrum given by Equation (1), and a standardized SN brightness m^s_B is computed with parameters α and β:

$$\begin{equation} m_B^s=m_B + \alpha x_1 - \beta c. \end{equation} \tag{ 2 }$$

Guy et al. (2007) describe a technique that simultaneously determines the cosmological parameters and the coefficients (α and β) that minimize the residuals on the Hubble diagram. The SALT2 approach to correction for color is formally equivalent to the description in terms of R_V and E(B − V), except for the interpretation that the value of β may result from some effect other than extinction by host galaxy dust. Regardless of the interpretation, β, which gives the magnitude correction for B-band, and R_V, which gives the magnitude correction for V-band, may be compared via the relation β = R_B = R_V + 1.

In this work, we present a new method for fitting the parameters α and β from the data. The method has been implemented in the program SALT2mu, which has been added to the SNANA software package. The motivations for a new program are as follows.

• 1.  
  The fit for α and β is decoupled from the cosmology. If the mean values of x₁ and c are a function of redshift (a magnitude-limited survey will have such selection effects), the average residual depends on both the cosmology and the α and β parameters. We prefer to determine the α and β solely by minimizing the rms of the distribution and to let the mean of the distribution determine the cosmology.
• 2.  
  A new description of intrinsic scatter is used as discussed below.
• 3.  
  The SALT2mu output can be used interchangeably with different cosmology fitters since SALT2mu produces distance moduli separately from the cosmology fit. This feature allows a direct comparison of SN distance moduli without reference to a particular cosmology.
• 4.  
  The sample used to determine α and β can be different from the sample used to determine the cosmology.

A significant additional benefit to writing a new program is that SALT2mu reads files from the SALT2 light-curve fits and outputs files that can be read by any of the SNANA cosmology fitters. The fitting model resides in a single function making it easy to modify the χ² function by adding new parameters, for example.

The inclusion of a more general description of the intrinsic scatter of the fitted SALT2 parameters is a major feature, as is the determination of α and β in a way that is independent of the cosmology. The details of the incorporation of the intrinsic scatter are given in the mathematical description (Section 2) below, and most of the paper is devoted to the issues associated with the more general description of the intrinsic scatter. The program is then applied to the Sloan Digital Sky Survey supernova (SDSS-II SN) data, which is used to illustrate some of the features of the program and issues with its use.

The SDSS-II SN data (Frieman et al. 2008) were obtained during 2004–2007 as part of extension of the original SDSS (York et al. 2000). The SDSS telescope (Gunn et al. 2006) and imaging camera (Gunn et al. 1998) produce photometric measurements in each of the five SDSS filters (Fukugita et al. 1996) spanning the range of 350–1000 nm. The most useful filters for observing SDSS SN, however, are g, r, and i because the SNe are too faint in u and z to be well measured by the SDSS SN survey except at low redshifts. Further details of the survey, data processing, and selections are given in Section 3.

The standardized 10 pc SN magnitude is determined by the expression

$$\begin{equation} M_0^{(B)}=m_B-\mu (z)+\alpha (z) x_1 - \beta (z) c, \end{equation} \tag{ 3 }$$

where μ(z) is the distance modulus and M^(B)₀ is the standardized SN magnitude for B-band. We have generalized the unknown parameters α and β to functions of redshift to accommodate possible effects of evolution of SNe or their environment. These functions are not known from first principles, but must be determined from the data. We use a nearly equivalent formulation

$$\begin{equation} M_0(z)=m_x-\mu (z)+\alpha (z) x_1 - \beta (z) c, \end{equation} \tag{ 4 }$$

where

$$\begin{equation} m_x = -2.5\log _{10}(x_0). \end{equation} \tag{ 5 }$$

The definition shown in Equation (5) has the advantage that the magnitude m_x does not depend on the shape of the B-band filter. It differs from the B-band standard magnitude by about 10 mag. The functions M₀(z), α(z), and β(z) are to be found by the minimization of the expression for χ²

$$\begin{equation} \chi ^2 = \displaystyle \sum _{n=1}^N\frac{\left[m_{xn}-\mu (z_n)+\alpha (z_n)x_{1n}-\beta (z_n)c_{n}-M_0(z_n)\right]^2}{\sigma _n^2+\sigma _{\rm {int}}^2} \end{equation} \tag{ 6 }$$

where m_xn, x_1n, and c_n are the fitted SALT2 parameters for the nth SN and N is the number of SNe. The quantity σ_n is the error on the quantity in the numerator due to measurement uncertainty for the nth SN. The quantity σ_int is discussed below and its form is displayed in Equation (7). The function M₀(z) is, by construction, degenerate with cosmological information contained in the distance modulus μ(z). The introduction of a nuisance function M₀(z) allows our fit to minimize the scatter in the Hubble diagram while being insensitive to the cosmology. In order to reduce this function to a few parameters we assume that M₀(z) is an arbitrary constant in each of several redshift bins centered at z = z_b, and we use the exact expression for the distance modulus μ(z) for some assumed cosmology. Effectively, we are using an assumed cosmology to extrapolate the residual to the redshift at the center of the bin and minimizing the Hubble diagram scatter at those reference points. Since our method introduces several new parameters, the statistical accuracy of the fit is reduced by the additional parameters. And since our method relies on the minimization of the scatter of distance moduli in a single bin, it requires a bin population of at least two SNe—preferably many more. These considerations are of no consequence for the SDSS sample, but may be of some concern in applying our method to samples with sparse data or unconventional cosmologies. We show in Section 4 that the fit results are insensitive to the choice of cosmological parameters if the number of redshift bins is sufficiently large.

Although this formalism allows fitting an arbitrary redshift dependence¹¹ for α and β, we will assume that α and β are constants, independent of redshift, in this paper. Recently, a number of papers (Kelly et al. 2010; Lampeitl et al. 2010; Sullivan et al. 2010) have found that the value of M₀ shows a dependence on host galaxy properties. Other papers have found correlations between line velocities in SN spectra (Wang et al. 2009; Foley & Kasen 2011) and their standardized magnitudes and other correlations with spectral features (Nordin et al. 2011). Those aspects will not be explored in this paper, but it is straightforward to generalize Equation (4) by adding a new term γy_n to the numerator of Equation (6), where y_n is some measured property of the nth SN and γ is a parameter to be determined by the minimization. In addition, the measurement uncertainty and intrinsic scatter in y_n should be incorporated into σ_n and σ_int, respectively, following the formalism described below.

When the SALT2 light-curve parameters are used according to Equation (4) to construct an SN Hubble diagram, the dispersion exceeds that expected from the light-curve fit parameter errors. A simple way to parameterize the Hubble diagram scatter is to assert that the distribution of light-curve parameters can be described by an arbitrary population in the plane described by Equation (4) plus a Gaussian spread that is described by a covariance matrix, which we will refer to as the intrinsic scatter. The quantity σ_int in the denominator of Equation (6) describes the expected error in the numerator due to the intrinsic scatter. We further assume that the covariance matrix is constant, independent of the SN light-curve parameters. These considerations result in an expression for the intrinsic scatter

$$\begin{eqnarray} \sigma _{\rm {int}}^2 &=& \Sigma _{00} + \alpha (z)^2\Sigma _{11}+\beta (z)^2 \Sigma _{{\rm cc}} +2\alpha (z)\Sigma _{01}\nonumber\\ &&-2\beta (z)\Sigma _{0c}-2\alpha (z)\beta (z)\Sigma _{1c}, \end{eqnarray} \tag{ 7 }$$

where Σ is the intrinsic covariance matrix of the parameters and Σ₀₀, Σ₁₁, and Σ_cc are the diagonal elements that correspond to the SALT2 parameters m_x, x₁, and c. Our introduction of a covariance matrix to describe the intrinsic scatter differs from the conventional approach where the intrinsic scatter is assumed to be a single number. Equation (7) shows that the conventional approach is equivalent to our model if Σ₀₀ is the only non-zero term, leading to the conclusion that the conventional approach applies the intrinsic scatter to m_x but not to x₁ or c. We introduce the more general approach because the additional terms have not been excluded on either experimental or theoretical grounds and, in fact, there is good reason to believe that there is a significant scatter that can be attributed to the SALT2 color parameter. The measurement error term (σ²_n) in Equation (6) is dominated by photometric errors and is calculated by the SALT2 light-curve fit for each SN and is given by a covariance matrix that is of the same form as the intrinsic scatter covariance displayed in Equation (7). The denominator in Equation (6) is therefore the sum of a constant intrinsic scatter covariance matrix, which we have introduced above, and the measurement covariance, which has traditionally been included. The fact that the measurement error and the intrinsic error are both of the form shown in Equation (7) means that each of the SALT2 parameters has both an intrinsic scatter and a measurement error that are taken into account in computing χ². The formulation of the minimization problem defined by Equations (6) and (7) is an extension of the formula given in Press et al. (2007) for the case of fitting a line to two variables with known errors. A more extensive treatment of the linear regression problem from an astronomical perspective has been given, for example, by Kelly (2007).

After the parameters α and β are determined by the minimization of χ² the distance moduli are calculated according to

$$\begin{equation} \mu _n = m_{xn}+\alpha (z_n) x_{1n}-\beta (z_n) c_{n}-M_0(z_b). \end{equation} \tag{ 8 }$$

These distance moduli can then be used for input to cosmological fits.

We now apply this method to the SDSS four-year sample of 529 spectroscopically confirmed SNe. The SDSS SN survey is a so-called rolling search, where a portion of the sky is repeatedly scanned to discover new SNe and to measure the light curves of ones already discovered. The survey measured an equatorial stripe about 25 wide in declination between a right ascension of 20^h and 04 Full coverage of the stripe can be obtained in two nights, but the average cadence was about four nights because of inclement weather and interference from moonlight. The data were obtained in the Fall of 2005, 2006, and 2007 and a smaller amount of data was taken in a test run in 2004. The survey is sensitive to SN beyond a redshift of 0.4, but beyond a redshift of 0.2 the completeness and the ability to obtain accurate photometry deteriorates. The SDSS camera images were processed by the SDSS imaging software and SN were identified via a frame-subtraction technique. More details and references can be found in Frieman et al. (2008). The candidate selection and spectroscopic identification has been described by Sako et al. (2008). The SN photometry used in this analysis was obtained using the scene modeling technique of Holtzman et al. (2008), where the sequence of SN observations is modeled as a variable point source and a host galaxy background that is constant in time.

We select a sample of well-measured light curves for this analysis based on the criteria used by Kessler et al. (2009a). We require the following.

• 1.  
  Spectroscopically confirmed as an SN Ia using techniques similar to those used for the first year sample (Zheng et al. 2008).
• 2.  
  Redshift in the range 0.02 < z < 0.42.
• 3.  
  A minimum number of measurements without any requirement on the measurement error and a minimum number of well-measured points, where well-measured means an error of 20% or less (i.e., S/N > 5)
• 4.  
  An acceptable fit to the SALT2 model (confidence level greater than 0.001).

These selections result in a total sample of 343 SNe with negligible contamination from sources that are not SNe Ia. The redshift range spans 0.036–0.419, with a mean redshift of 0.22. The SDSS sample thus defined does not have uniform selection criteria as a function of redshift and, therefore, there is a significant change in the population as a function of redshift. In particular, the spectroscopically selected sample of SDSS SNe is significantly biased toward bluer, brighter SNe as the redshift increases; the selection bias is discussed in detail in Kessler et al. (2009a). To the extent that the data are described by the model of Equation (4) with an accurate parameterization of α(z) and β(z) and the SALT2 fit parameters are unbiased, we expect to be insensitive to having a biased population.

An application of SALT2mu to the SDSS data using four redshift bins yields the results shown in Table 1, which shows the fit χ² for 337 degrees of freedom (343 SNe), the fitted values of α and β, and their errors σ_α and σ_β. In Section 4, we show that four bins are adequate to decouple the SDSS data from the determination of the cosmological parameters. The first row in Table 1 shows the fit when the intrinsic scatter is assumed to be zero resulting in χ² = 1305, reproducing the well-known result that the scatter in the Hubble diagram is larger than what is expected from the measurement error alone. Subsequent rows illustrate how the fitted parameters change with different assumptions on the nature of the intrinsic scatter. Each of the fit variables is assumed, in turn, to be totally responsible for the intrinsic scatter with the magnitude of the scatter being adjusted to achieve χ²/dof ≈ 1. In Table 1 we have introduced new variables to indicate the size of the intrinsic scatter: $\sigma _0=\sqrt{\Sigma _{00}}$, $\sigma _1=\sqrt{\Sigma _{11}}$, and $\sigma _c=\sqrt{\Sigma _{{\rm cc}}}$. From Table 1 we see that the different assumptions on the nature of the intrinsic scatter result in rather different estimates of the α and β parameters compared to the reported statistical error. An assumed scatter in m_x (σ₀ ≠ 0) decreases both α and β compared to the case where the intrinsic scatter is taken to be zero. An assumed scatter in x₁ (σ₁ ≠ 0) increases α by 15.7σ while decreasing β by 1.4σ relative to the fit where the intrinsic scatter is assumed to be entirely due to m_x. The opposite effect is seen when the scatter is assumed to be due to c (σ_c ≠ 0): β increases by 7.1σ and α decreases by 3.1σ. Even more variation could be exhibited by a more general intrinsic scatter covariance matrix. However, the variations in α and β are not without limit. The values are all the same sign and approximately the same magnitude. Roughly speaking, we are changing our estimate of how much to correct the observed slope for the flattening caused by the scatter described by the intrinsic covariance matrix (Equation (7)).

If we are only interested in SN as distance indicators, α and β are nuisance parameters of no particular interest. The average SN brightness at a given redshift will be given by the mean of Equation (4): if the mean values of x₁ and c do not evolve with redshift, then different values of α and β will result in a redshift-independent offset that is not significant in determining the cosmology. The mean values of x₁ and c for the SDSS data are plotted in Figure 1 showing, however, a substantial increase in x₁ and a decrease in c with redshift. These trends are primarily caused by the spectroscopic selection bias toward brighter SNe Ia. The different values listed in Table 1 would result in differences of ∼0.1 in the average values of μ(z), a rather large error for cosmological studies. In practice, the effect can be mitigated by limiting the sample to lower redshifts, where the selection efficiency is higher, correcting for biases, and combining high and low redshift experiments that have different efficiencies. Nonetheless, the possibility of an uncorrected bias is a concern.

Zoom In Zoom Out Reset image size

Our results can be compared to the results of Amanullah et al. (2010) who combined data from several sources to obtain a value of β = 2.51 ± 0.07 with an analysis most similar to our case with σ₀ = 0.14. Guy et al. (2007) found β = 1.77 ± 0.16 for the first year data from the SNLS and Kessler et al. (2009a) found β = 2.65 ± 0.08 for the first year of SDSS data. Amanullah et al. (2010) also quote numbers for their four largest samples: β = 2.38 ± 0.14 (SDSS), 2.73 ± 0.13 (CfA), 2.50 ± 0.17 (ESSENCE), and 1.72 ± 0.18 (SNLS). More recently Conley et al. (2011) found β = 3.18 ± 0.10 for the SNLS three-year sample using a chi-squared minimization technique. The study of Nobili and Goobar (2008) found that the variations in SN colors were most consistent with a CCM parameterization when R_V = 1.75 ± 0.27. However, the Nobili and Goobar (2008) study differs significantly from this analysis in that their value of R_V is based on their measurement of the wavelength dependence of the reddening using the Cardelli parameterization. Kessler et al. (2009a) also investigated SN color and found R_V = 2.18 ± 0.14(stat) ± 0.48(syst), but that result relies on modeling the color changes due to the magnitude-limited sample. All these results are broadly consistent in finding β < 4.1, but they exhibit more scatter than one would expect from the quoted errors.

We will return to the SDSS SN data in Section 5, but it is useful to first examine the performance of SALT2mu on simulated SN data. We first determine the number of data bins that are required to decouple our fits from the cosmology over a wide range of parameters and conclude that four bins is adequate for the SDSS data. We then verify the formalism with a simplified simulation and show that we recover the input values of α and β without bias. We next look at a realistic simulation and show that the fitted α and β dependence on the intrinsic covariance matrix Σ follow the trends seen in the data. We compute the intrinsic scatter covariance matrix from the simulation (which requires knowing the simulation parameters), and show that the input values of α and β are recovered when the correct intrinsic scatter covariance matrix is used.

We generate the simulated SNe by using SNANA to generate SN light-curve data according to the SALT2 model. The SNe are generated with a redshift distribution appropriate for the SDSS survey and with α = 0.18 and β = 3.4. Realistic cadence, measurement errors, and selection effects are applied to the simulated SN, and they are processed with the same light-curve fitting program and selection criteria as the SDSS SN, resulting in SALT2 fit parameters for the simulated SN. In addition to the photometric measurement errors, the simulation applies color smearing to model the intrinsic scatter. Color smearing consists of random magnitude offsets that are generated for each SN and observed filter and added to each light-curve point measured in that filter. The magnitude of the smearing is adjusted to 0.08 in each observed filter to give approximately the same Hubble diagram scatter as seen in the real data.

The first test is a relatively trivial one to verify that the fit results for α and β are, indeed, independent of cosmology. We expect that a large number of parameters M₀(z_b) will decouple the determination of α and β from the cosmological parameters, but wish to minimize the number of nuisance parameters that are introduced. We illustrate the fit dependence on the number of redshift bins in Figure 2 using the simulation of the SDSS data set, which has over 10 times the number of SNe (4417) as the SDSS data. Figure 2 shows the fit results for α and β assuming three different cosmologies: Ω_M = 0.27, $\Omega _\Lambda =0.73$, w = −1 (red), Ω_M = 0.27, $\Omega _\Lambda =0.0$ (green), and Ω_M = 0.27, $\Omega _\Lambda =0.73$, w = 0 (blue). The difference in the fitted parameters is noticeable when only one bin is used, but becomes negligible as the number of bins becomes greater than about five. We conclude that using four redshift bins is adequate for the SDSS data, given the lower statistical precision of the data and the established constraints on the cosmological parameters. All results in this paper are based on dividing the SDSS data (real or simulated) into four equal sized bins spanning the redshift range 0.02 < z < 0.42.

Zoom In Zoom Out Reset image size

While the fit results in Figure 2 show that we have successfully decoupled the fit for α and β from the cosmology, the dependence on the assumptions about the intrinsic scatter matrix remains. For the fit results shown in Figure 2 we took σ₀ = 0.14 and σ₁ = σ_c = 0, and this choice results in fitted values of α = 0.184 and β = 2.86 which differ from those used to generate the simulated SNe, namely, α = 0.18 and β = 3.4. The discrepancy is not surprising since our simulation includes significant color smearing that we have not taken into account. We will determine more accurate parameters for the intrinsic scatter covariance below, but first we check the formalism with a simplified simulation. We use the same simulated SN but replace the simulated values of m_x, x₁, and c with the exact values that were used to generate the light curves but smeared according to the 3 × 3 error matrix including any assumed intrinsic scatter. This procedure tests the self-consistency of the SALT2mu model and eliminates concerns about biases in the SALT2 parameter fitting and inaccuracies in the calculation of the parameter errors. The results of the fits are shown in Table 2 with the same four assumptions about the intrinsic scatter that were used in Table 1. We have also tested the formalism when all the elements of the intrinsic scatter covariance matrix are non-zero; the last line in Table 2 is an example. For each case we recover the input values that are consistent within the statistical error of about 1%. The χ²/dof consistent with 1. We do not observe any variation in the parameter values since the parameters of each SN are generated according to a Gaussian error distribution that is the combination of the measurement error and whatever intrinsic scatter was assumed.

Using the parameters of the SALT2 fits to the simulated data provides a more complicated test since the color-smearing model is different than applying an intrinsic scatter to the simulated SALT2 parameters used to generate the light curves. We apply SALT2mu to the simulated SN using the same procedure used for the SDSS data, and display the results in Table 3. The first line has the intrinsic scatter set to zero, and the succeeding lines have the intrinsic scatter that results in χ²/dof ≈ 1. The general trends shown in Table 3 are very similar to those seen in Table 1: assuming an intrinsic scatter in x₁ and c increases the fitted value of α and β, respectively. The magnitudes of the intrinsic scatter required to achieve χ²/dof ≈ 1 are similar to those required by the SDSS data. But we can also compare the results with the "true" values for the simulated data. While the assumption that the intrinsic scatter is entirely due to color errors (last line in Table 3) produces values of α and β that are close to those used to generate the simulation (α differs by 2.6σ), it is not a very accurate representation of the intrinsic scattering covariance matrix that describes the simulated SN as shown below.

In principle, we know the smearing formulae applied in the simulation and can calculate the intrinsic scatter covariance matrix Σ. But the calculation of the intrinsic smearing matrix is complicated because the error model used in the simulation smears individual light-curve points in a way that is not easily expressed in terms of an intrinsic smearing matrix in m_x, x₁, and c. However, it is relatively straightforward to calculate the error a posteriori. The full error matrix is determined by the variance of the simulated measurements of m_x, x₁, and c relative to the true values (which, of course, are known for the simulated data). The light-curve fitter calculates covariance matrix of measurement errors, which, when added to the intrinsic scatter covariance matrix, should yield the total error matrix. The intrinsic scatter covariance matrix can therefore be calculated for the simulated data by subtracting the observed variance from the prediction for the measurement error. The total error covariance matrix calculated from the variance of the reconstructed SALT2 parameters from the simulated ones is

$$\begin{equation} \begin{array}{lll}\sigma _0= 0.0781, & \xi _{01}=0.253, & \xi _{0c}=0.612, \\ \sigma _1=0.351, & \xi _{1c}=0.183, &\hbox{and} \\ \sigma _c=0.0634, \end{array} \end{equation} \tag{ 9 }$$

where the various ξ are the correlation coefficients (ξ₀₁ = Σ₀₁/(σ₀σ₁), ξ_0c = Σ_0c/(σ₀σ_c), and ξ_1c = Σ_1c/(σ₁σ_c)). Using the same notation to represent the elements of the measurement error matrix we find

$$\begin{equation} \begin{array}{lll}\sigma _0= 0.0365, & \xi _{01}=0.581, & \xi _{0c}=0.694, \\ \sigma _1=0.347, & \xi _{1c}=0.307, &\hbox{and}\\ \sigma _c=0.0286. \end{array} \end{equation} \tag{ 10 }$$

Subtracting Equation (10) from Equation (9) gives an excess error, which we attribute to the intrinsic scatter. The result is

$$\begin{equation} \begin{array}{lll}\sigma _0= 0.0690, & \xi _{01}=-0.119, & \xi _{0c}=0.590, \\ \sigma _1=0.0506, & \xi _{1c}=0.362, &\hbox{and}\\ \sigma _c=0.0566. \end{array} \end{equation} \tag{ 11 }$$

We are now in position to apply the intrinsic scatter from Equation (11) above to the fit from SALT2mu. The result is shown in the first row in Table 4. The result is consistent with the generated values (α = 0.18 and β = 3.4). The small discrepancy in β (2.6σ) not a serious concern because it is quite small compared to the accuracy of the SDSS data set. However, we suspect that uncertainties in the calculation of the measurement errors, including treatment of nonlinearities and a consistent handling of model errors, may contribute a systematic bias at this level of precision. We have performed additional fits with other intrinsic scatter matrices. The second row is a fit with all the non-diagonal elements of the intrinsic scatter matrix set to zero. The fitted parameters are only slightly changed although there is a decrease in χ². The next row shows a fit that has been selected to give a larger change in χ² by tuning the correlations for maximum effect. The last row in Table 4 displays a fit, which may be counterintuitive: the error in x₁ is decreased, but the χ² increases. These results illustrate how complicated the fit behavior can be for the five free parameters that live in the intrinsic scatter matrix.

We do not have a rigorous approach to estimating the intrinsic scatter matrix for actual SN data. Since we do not know the underlying, unsmeared distribution or the actual form of the smear- ing, the observed distribution is difficult to interpret. The scatter in the Hubble diagram sets the overall scale, but more information is needed to determine the individual elements of the intrinsic scatter matrix.

We can derive some information from the color distribution with the help of some plausible assumptions, which were first used by Jha et al. (2007). We assume that the color distribution arises from random color smearing that we parameterize as a Gaussian of unknown width (σ_c) and host galaxy extinction that we characterize as an exponential distribution of unknown slope (τ). We do not expect to be sensitive to the exact shape assumed for galaxy extinction: the important point is that it is one-sided (extinction can only make objects redder, not bluer) and that the most probable value of extinction is zero. In addition, we do not know the position of the edge of the extinction distribution, so we must introduce an additional parameter (c₀) to be determined from the data. Assuming that the intrinsic scatter and host galaxy extinction are independent, the convolution of the two processes results in the following distribution for c:

$$\begin{equation} \rho (c)=A \int _{c-c_o}^\infty \rm exp\left(-\frac{c-c_o-c^\prime }{\tau }\right) \rm exp\left(-\frac{{c^\prime }^2}{2 \sigma _c^2}\right) d c^\prime. \end{equation} \tag{ 12 }$$

In order to minimize the contribution of measurement errors to the width of the Gaussian component, we include only SN where the estimated measurement error for c is less than 0.04, reducing our sample of 343 SNe to 247. Fitting Equation (12) to the data, we find σ_c = 0.061 ± 0.007, τ = 0.076 ± 0.010, and c₀ = −0.073 ± 0.009. The data and resulting fit are shown in Figure 3. As can be seen, the functional form of Equation (12) provides a good description of the data. The fitted value for τ is substantially smaller than the value τ = 0.334 ± 0.088 quoted by Kessler et al. (2009a), but our result is an effective τ for the spectroscopically selected sample with the additional requirement of a well-measured color and is therefore biased against large extinctions. Measurement errors are calculated by the SALT2 light-curve fit and contribute, on average, $\sqrt{\langle \sigma _c^2 \rangle }=0.026$ to the width of the Gaussian component, resulting in an estimate of 0.055 ± 0.007 for the color smearing. We do not expect our values of c₀ and σ_c to be significantly affected by the selection bias, and we do not make a correction for the bias. This result is comparable to the amount of color smearing (σ_c = 0.05) that was required in Table 1 to explain, by itself, the scatter in the Hubble diagram.

Zoom In Zoom Out Reset image size

Unfortunately, we do not have any useful a priori model for the x₁ distribution. We can, however, make the somewhat trivial observation that the intrinsic scatter cannot be larger than the width of the distribution. It is known that passive and star-forming galaxies have different distributions in x₁ and, if we assume that the intrinsic scatter is the same for both, we can obtain an upper limit for the intrinsic scatter from the distribution that excludes passive galaxies. The distributions in x₁ were shown by Lampeitl et al. (2010) for the SDSS data, and we have used that analysis for the 296 SNe this paper has in common with that analysis. The resulting distribution is shown in Figure 4 along with a Gaussian of width 1.0. The Gaussian is not a good fit to the tails of the distribution, but it is a good match to the core and we conclude that an upper limit for the intrinsic scatter in x₁ is 1.0, a limit that is greater than the amount required to explain the scatter in the Hubble diagram.

Zoom In Zoom Out Reset image size

The simulated data exhibit a significant correlation between m_x and c, and it is interesting to ask whether that correlation should be expected in the data. In the simulation, the correlation arises because the SALT2 model defines B-band to be the reference point for the color and because the measured data are largely redward of B-band. The assumed independent variation of each measured band results in the correlation between peak magnitude and color. While the model of color smearing is arbitrarily constructed to reproduce the scatter in the Hubble diagram, there is some evidence to support independent fluctuations in each of the filters. As pointed out in Kessler et al. (2010), if the redshift is fit from the light-curve data, the dispersion of the results around the spectroscopically determined redshifts is larger than can be accounted for by photometric errors alone, but the color-smearing model of SNANA gives approximately the right dispersion. Based on these considerations, we conclude that the correlation between m_x and c is likely to be positive (as seen in the simulations) for the SDSS data.

While we have been able to make some observations on the elements of the intrinsic scatter matrix, our picture is incomplete and rests on assumptions of uncertain validity. An alternative approach is to marginalize over all possible intrinsic scatter matrices. We implement this marginalization by running SALT2mu repeatedly with randomly generated intrinsic scatter matrices. We start with a vector randomly generated on the unit sphere subject to the restriction that all components of the vector are positive. We then multiply the three components of the vector by 0.2, 1.0, and 0.07, respectively, and set σ₀, σ₁, and σ_c to the scaled vector components. The constants are chosen to be about 50% larger than the values (see Table 1) that were required to explain the intrinsic scatter. The normalized correlations are chosen randomly between −1 and +1.

A total of 10,000 intrinsic scatter matrices were generated, but the procedure does not guarantee that the matrices are positive. A total of 6380 of the matrices were found by SALT2mu to be positive definite and were used to determine values for α and β. Each intrinsic scatter matrix is scaled by an arbitrary factor to achieve a χ²/dof = 1. The red (dotted) curve in Figure 5 shows the results of the fits for (a) α and (b) β. The median value is α = 0.150^+0.083_{− 0.021} and β = 2.91^+0.30_{− 0.35}, where the error is the range that includes 68% of the results. Our intrinsic scatter covariance matrix defines an intrinsic error in three dimensions, but the chi-squared is only sensitive to the error in one direction. The intrinsic errors in the other two directions just broaden the distribution of points in the plane defined by Equation (8). If the intrinsic errors are large enough to explain, by themselves, the observed population, the slopes in Equation (8) become undefined and produce the outliers that are observed in Figure 5. The data show a secondary peak near β = 2.5, but this should not be interpreted as evidence for two populations with different slopes. The shape of the plot in Figure 5, of course, depends on the distribution of covariance matrices, which we chose to be uniform. The lower values of β result primarily from intrinsic scatter matrices with a small σ_c. The fit becomes insensitive to σ_c as σ_c → 0 and the peak at β = 2.5 is an artifact of our choice of distribution for the intrinsic scatter matrices and the insensitivity of the fit parameter β in the vicinity of σ_c = 0.

Zoom In Zoom Out Reset image size

The blue (solid) curve is the result of restricting the intrinsic scatter matrix to a range favored by the arguments in this section. Specifically we require the intrinsic scatter matrix color element to be within three standard deviations of the value determined by the fit to Equation (12), namely, 0.040 < σ_c < 0.082. We also require σ₁ < 1 since the intrinsic scatter cannot exceed the measured width of the distribution, and ξ_0c > 0 since we expect a positive correlation when the color smearing includes a component that is uncorrelated between the observed filter bands. While the exact values of these cuts are somewhat arbitrary, we have chosen a wide range to account for model uncertainties. Only 1650 of the 6380 intrinsic scatter matrices satisfy these three criteria, but even with the restrictions a wide range of values are found. However, the median values and the ranges that contain 68% of the data are α = 0.135^+0.031_{− 0.013} and β = 3.19^+0.10_{− 0.22}. In addition to the range of values that result from different assumptions about the intrinsic scatter matrix, the statistical errors reported by SALT2mu are σ_α = 0.011 and σ_β = 0.10. Adding these errors in quadrature, the final result is α = 0.135^+0.033_{− 0.017} and β = 3.26^+0.14_{− 0.24}. For each entry in Figure 5, we compute the statistical probability that the entry is consistent with β ⩾ 4.1 using the computed statistical error and assuming a normal distribution. The uncertainty in the intrinsic scatter matrix is computed by averaging the statistical probability for all the entries in Figure 5. From this calculation, we conclude that the probability that our result is consistent with β ⩾ 4.1 is 2%. The probability is higher than would be naively estimated from the value β = 3.26^+0.14_{− 0.24} because our calculation includes the effect of the non-Gaussian tail of the distribution shown in Figure 5(b).

The parameters α and β are usually assumed to be constants, independent of redshift, but the possibility of a redshift dependence has been considered previously by, for example, Kessler et al. (2009a). We might observe a redshift dependence of the fits for α and β for a number of reasons.

• 1.  
  The properties of the SN, themselves, may evolve with redshift.
• 2.  
  The model of Equation (4) may not be an exact model and the resulting fit parameters may be in tension between different SN populations. If the population is a function of redshift, the amount of tension may vary, resulting in different results for α and β.
• 3.  
  Our SN selection has biases, particularly at higher redshift where we tend to select the brighter SN.

While we could directly fit for α(z) and β(z), the likelihood that selection effects will play a significant role suggests that a better approach might be to split the SDSS sample into four redshift bins. The results of the fits are shown in Figure 6. We have used the same ensemble of intrinsic scatter matrices with σ₁ < 1, 0.040 < σ_c < 0.082, and ξ_0c > 0 that was used to construct Figure 5. The errors shown are the combined statistical and systematic error that arise from the uncertainty in the intrinsic scatter matrix. The results shown for α in Figure 6(a) are consistent with the sample average although there is an indication that the value of α may be rising as a function of redshift. Figure 6(b) shows a value of β that is 2.1σ lower in the lowest redshift bin, 2.2σ higher in the range 0.22 < z < 0.32, and 3.7σ lower in the highest redshift bin than the result for the complete sample. The accuracy in the determination of β as a function of redshift varies not only because of the statistics but also because highly reddened SNe are too dim to meet our selection criteria at high redshift. The number of SNe in each redshift bin and the number of highly reddened (c > 0.3) SNe are 36 and 4 (0.02 < z < 0.12), 149 and 2 (0.12 < z < 0.22), 116 and 0 (0.22 < z < 0.32), and 41 and 0 (0.32 < z < 0.42), respectively. The reduced range of c at the higher redshifts reduces the accuracy of the fit and makes the fit considerably more sensitive to the amount of intrinsic scatter in the data. The highest redshift sample consists of 41 SNe with −0.319 < c < 0.126 that show very little correlation between the distance modulus and c, resulting in a low value of β. The range of observed values of c is most likely primarily due to measurement error and intrinsic scatter resulting in the weak correlation between distance modulus and observed c. In principle, our formalism should recover the correct value of β when the correct measurement errors and intrinsic scattering are included, but selection effects have significantly distorted the natural spread. In addition to the restricted color range, the absolute magnitudes of SNe with 0.32 < z < 0.42 are 0.069 ± 0.042 mag brighter on average, assuming standard cosmological parameters ($\Omega _\Lambda =0.7$ and w = −1). As a consequence, our fit to the last redshift bin should not be considered to be a good representation of what would be obtained with an unbiased sample. A more robust result would be possible with an improved simulation of the selection effects and augmenting the spectroscopically confirmed sample with those identified photometrically (Sako et al. 2011) to obtain a more complete sample. Selection effects should be much less important for the midrange (0.22 < z < 0.32), where β is larger than the sample average by 2.2σ. While the lack of highly reddened SNe (c > 0.3) degrades our ability to determine β in this redshift region, we can see if their absence biases β to higher values by excluding SNe with c > 0.3 in the low redshift bins. The β values for the two lowest redshift regions change by less than 0.1 when the highly reddened SNe are removed from the sample.

Zoom In Zoom Out Reset image size

The histograms in Figure 5 show the same trend that was observed in Table 1: when the intrinsic scatter in color is large, the best-fit value of β will increase while α decreases slightly. The fact that our value for β is the highest value of those reported in Section 3 depends on our incorporation of a large intrinsic scatter in SN color. We have not made any correction for the selection bias of our sample based on the results of the simulation, which is subject to the same selection criteria, and shows no significant bias in the determination of α and β.

The error on our result is dominated by an uncertainty in the intrinsic scatter matrix, but there are other potential problems that have not been taken into account. Certainly, our assumption that the intrinsic scatter is a constant, independent of any SN properties, is suspect given our lack of knowledge of the processes involved. In addition, our model assumes that the data follow the relationship of Equation (4) with a fixed value of β. However, we know that the value of R_V varies significantly in our own Galaxy (Schlafly et al. 2010), and we should expect to see at least some variation in SN host galaxies. Variations in β produce effects similar to the intrinsic scatter in color, but the effect can be magnified in a magnitude-limited sample like SDSS where there are few highly reddened SNe, which have higher weight in determining β. Variations in β could explain the small differences seen between the values of β in Figure 6 in the two lowest redshift bins.

We have seen that neglecting intrinsic scatter in the SALT2 c parameter results in a value of β that is biased low. The amount of bias depends on the amount of intrinsic scatter relative to the range in the variable c. Since the range in c decreases as the redshift increases, neglecting or underestimating the amount of intrinsic scatter in c would be expected to result in a decreasing value of β as a function of redshift. The fits to the SDSS data as a function of redshift include intrinsic scatter in color so this effect should not be the explanation for the SDSS data—unless the scatter was underestimated. However, an overestimate of the intrinsic scatter in m_x leads to similar biases and is more likely explanation for the lower value of β that is observed in the highest redshift bin.

A more problematic type of variation in β could arise from the color smearing. We have argued that color variations of SNe must have some component of host galaxy extinction and also some intrinsic color variation, which could be due to variations in the explosion process or local environment or both. We have modeled the intrinsic color variation as color smearing with an effective value β = 0, but it is possible that each mechanism that contributes to the observed color requires a different correction to the distance modulus. Our model can only produce an effective value of β that could be a function of color, redshift, host galaxy type, or other parameters and could, therefore, differ depending on the characteristics of the particular SN sample.

Although we have not included correlations with galaxy types or spectral line features, the introduction of additional dependent variables into Equation (4) is straightforward. In fact, having a more complete description of the SN explosion will not only provide a smaller dispersion for the numerator in Equation (6), it will also decrease the size of the intrinsic scatter that is required and decrease the importance of the uncertainty as to its form.

Recently, a number of papers (Kim 2011; Lago et al. 2011; March et al. 2011) have explored the technique of determining the intrinsic scatter from the data. The approach of March et al. (2011) is similar to ours in that it included our more general form for the intrinsic error matrix but assumed that the off-diagonal terms involving M_B were zero. We have not considered including the intrinsic scattering error as a fit parameter partly for simplicity and partly because the underlying physics of the intrinsic scatter is still uncertain. However, the more rigorous statistical approach described by these authors is a possible direction for future work.

We argue from the SDSS SN photometric data that there is likely some intrinsic scatter in SN color that is responsible for some of the scatter in the Hubble diagram. A similar conclusion is reached by Foley & Kasen (2011), based on much more detailed information involving differences in the line velocities measured from the SN spectra, which are modeled as arising from differences in viewing angle. While this type of detailed spectral information may not be available for all future SN surveys, an understanding of the underlying physical processes will be important for the most accurate treatment of the data.

We have described a new formalism to fit the parameters α and β that are used to determine the standard magnitudes of SNe Ia. The determination of these parameters is accomplished in a way that is independent of cosmology by introducing a few nuisance parameters to accommodate any reasonable cosmology. The formalism also introduces a more general form for the intrinsic scatter which introduces a large uncertainty into the fit if the parameters describing the intrinsic scattering are not known accurately. We have shown by simulation that the mathematical formalism is self-consistent, but have pointed out ways in which the model may fail to be a complete description of SNe.

We find that the SDSS data, when fit to the form of Equation (4), are described by α = 0.135^+0.033_{− 0.017} and β = 3.19^+0.14_{− 0.24} by marginalizing over the uncertainty in the intrinsic scatter matrix. Our result relies on our conclusion that SNe are subject to a substantial color smearing in addition to reddening from host galaxy dust, as indicated by our fit to the color distribution where we determined that σ_c = 0.055 ± 0.007. The uncertainty in the parameters of the intrinsic scatter matrix results in a much larger error than is obtained if the parameters are assumed to be known. However, even with the larger value of β and the larger error, we find that SDSS data differ at the 98% confidence level from β = 4.1, the value expected for extinction by dust in the Milky Way.

The SDSS-II SN survey was managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions were the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, Cambridge University, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPA), the Max-Planck-Institute for Astrophysics (MPiA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington.

This work is based in part on observations made at the following telescopes. The Hobby–Eberly Telescope (HET) is a joint project of the University of Texas at Austin, the Pennsylvania State University, Stanford University, Ludwig-Maximillians-Universität München, and Georg-August-Universität Göttingen. The HET is named in honor of its principal benefactors, William P. Hobby and Robert E. Eberly. The Marcario Low-Resolution Spectrograph is named for Mike Marcario of High Lonesome Optics, who fabricated several optical elements for the instrument but died before its completion; it is a joint project of the HET partnership and the Instituto de Astronomía de la Universidad Nacional Autónoma de México. The Apache Point Observatory 3.5 m telescope is owned and operated by the Astrophysical Research Consortium. We thank the observatory director, Suzanne Hawley, and former site manager, Bruce Gillespie, for their support to this project. The Subaru Telescope is operated by the National Astronomical Observatory of Japan. The William Herschel Telescope (WHT) is operated by the Isaac Newton Group, the Nordic Optical Telescope (NOT) is operated jointly by Denmark, Finland, Iceland, Norway, and Sweden, and the Telescopio Nazionale Galileo (TNG) is operated by the Fundación Galileo Galilei of the Italian INAF (Istituto Nazionale di Astrofisica) all on the island of La Palma in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. Observations at the ESO New Technology Telescope at La Silla Observatory were made under programme IDs 77.A-0437, 78.A-0325, and 79.A-0715. Kitt Peak National Observatory, National Optical Astronomy Observatories (NOAO), is operated by the Association of Universities for Research in Astronomy, Inc. (AURA) under cooperative agreement with the NSF. The South African Large Telescope (SALT) of the South African Astronomical Observatory is operated by a partnership between the National Research Foundation of South Africa, Nicolaus Copernicus Astronomical Center of the Polish Academy of Sciences, the Hobby–Eberly Telescope Board, Rutgers University, Georg-August-Universität Göttingen, University of Wisconsin-Madison, University of Canterbury, University of North Carolina-Chapel Hill, Dartmouth College, Carnegie Mellon University, and the United Kingdom SALT consortium. The WIYN Observatory is a joint facility of the University of Wisconsin-Madison, Indiana University, Yale University, and NOAO. The W. M. Keck Observatory is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.

This work was supported in part by the U. S. Department of Energy under contract number DE-AC0276SF00515.