COSMOLOGY WITH PHOTOMETRIC SURVEYS OF TYPE Ia SUPERNOVAE

The apparent observed flux from a SN at z can be written as the convolution of the spectral energy distribution (SED) and the transmission function of the telescope,

$$\begin{equation} f_{\rm obs} = \int T_X(\lambda _{\rm obs}){\rm SED}\big(\lambda _{\rm obs}, t_{\rm obs}, s, E_{B-V}^{\rm host}, E_{B-V}^{MW}\big)d\lambda _{\rm obs}, \end{equation} \tag{ 1 }$$

where T_X(λ_obs) is the filter response for band X, λ_obs is the observed wavelength, t_obs is the observation date, s is the stretch factor, and E^host_B−V and E^MW_B−V are the color excess for the host galaxy and the Milky Way, respectively.

The transmission functions used in our analysis for 5-bands of LSST are shown in Figure 1 (LSST filter draft 2005).⁹ We also plot rest-frame SEDs for Type Ia SNe at z = 0.2, 0.5, and 0.8. These SED templates are from Nugent et al. (2002) and they cover the spectral wavelength from 1000 to 25000 Å in rest-frame days from −20 to 70 with respect to the B-band maximum light day. The SNe flux for the epochs before −20 are set to be zero.

Zoom In Zoom Out Reset image size

There are now several techniques to parameterize the SNIa light curves, such as the 15 day decline after the B-band maximum light Δm₁₅ (Phillips 1993) and the multicolor light curve shape method (MLCS and the update version MLCS2k2; Riess et al. 1996; Jha et al. 2007). In this paper, we calibrate the SNIa light curve with the timescale and stretch factor relation following the works of Perlmutter et al. (1997, 1999). By stretching and compressing the time axis around the rest frame B-band maximum light day, this method can fit the observed light curve very well using the light curve template (Goldhaber et al. 2001). Then this SED can be re-scaled by the apparent, unextincted B-band peak magnitude

$$\begin{equation} m_B = M_B+5{\rm log_{10}}d_L(z,{\bf \theta })+25-\alpha (s-1)+\Delta m, \end{equation} \tag{ 2 }$$

where M_B is the B-band absolute peak magnitude, d_L is the luminosity distance which is a function of the redshift z and a broad set of cosmological parameters denoted by θ, and α is the coefficient of the relation between s and m_B. Here, we take M_B = −19.3, α = 1.5 (Knop et al. 2003; Astier et al. 2006), and the set of cosmological parameters θ with Ω_m0 = 0.27, Ω_Λ0 = 0.73, and h₀ = 0.71 (Komatsu et al. 2009), where Ω_m0 and Ω_Λ0 have the usual meaning with the present-day matter and dark energy density parameters and h₀ is the dimensionless Hubble constant. Besides, we also consider the dispersion Δm of the rest -frame B-band peak magnitude after the calibration of the stretch factor. The B-band filter we use is from the Johnson–Morgan system (Bessell 1990, 2005). Also, the timescale of the SED is calibrated by the stretch factor s, which is assumed to be available from −15 to 35 around the B-band maximum luminosity day (Astier et al. 2006). We note that there is also an intrinsic color scatter, σ^B−V_int (standard deviation), that should be taken into account when producing mock light curves. Based on prior work, we find this uncertainty to be small with a value of ∼0.05 mag (Phillips et al. 1999; Jha et al. 2007). Hence, it would not affect our results much, so that for simplicity we don not consider it here.

During the transit, the SNe light will be partly absorbed by the dust of the host galaxy. We employ the reddening law of Cardelli et al. (1989) with R_V = 3.1, from infrared to far-ultraviolet (0.3 μm⁻¹ ⩽ x ⩽ 10 μm⁻¹, where x = 1/λ). For the optical to near-ultraviolet wavelength range (1.1 μm⁻¹ ⩽ x ⩽ 3.3 μm⁻¹), we use an updated version for extinction given by O'Donnell (1994). The latter uses the same analytical form for extinction as Cardelli et al. (1989) but with values of the fitting parameters revised slightly from the previous version. The level of extinction we assume here is consistent with the one measured recently by Menard et al. (2009) corresponding to large angular scales based on galaxy–quasi-stellar object cross-correlation in Sloan Digital Sky Survey. Since the SNe are at a different redshift than the observer, the spectrum is redshifted for both wavelength and the phase, i.e., λ' = λ(1 + z) and t' = t(1 + z). We also apply an extinction associated with dust in the Milky Way (Burstein & Heiles 1982; Schlegel et al. 1998), and assume that we have a perfect measurement of $E_{B-V}^{\rm {MW}}$. This assumption has no effect on our final conclusions. While the extinction of the Milky Way has different values for different sky regions, we do not have any information on the exact field selection of future SNe surveys from ground. Thus, we do not account for sky variation of extinction and simplify by just taking an average value with $E_{B-V}^{\rm {MW}} \approx 0.03$ and R_V = 3.1 for the extinction law. Finally, the spectrum is integrated with the LSST filters to get the mock light curve sampling in each of the five LSST filters.

Since the mean redshift of the SNe detections with LSST main survey is expected to be about 0.5 and the deeper, but smaller, survey can potentially detect SNe out to ∼1, we choose the redshift range from 0.01 to 1.1 when making mock SNe samples.

When creating large samples, we assume the flat $\rm \Lambda CDM$ model with Ω_m0 = 0.27 and h₀ = 0.71. Then, the redshift z, stretch factor s, the extinction of the host galaxy E^host_B−V, and the magnitude dispersion Δm are generated from the Gaussian distribution with truncated tails as follows: 0.01 ⩽ z ⩽ 1.1 with $\bar{z}=0.5$ and σ_z = 0.4, 0.6 ⩽ s ⩽ 1.4 with $\bar{s}=1$ and σ_s = 0.3, −0.1 ⩽ E^host_B−V ⩽ 0.3 with $\overline{E_{B-V}^{\rm host}}=0.0$ and σ_E = 0.2 and −0.3 ⩽ Δm ⩽ 0.3 with $\overline{\Delta m}=0.0$ and σ_Δm = 0.17. This extra dispersion acts as an extra source of noise in our mock data (Sullivan et al. 2006; Hamuy et al. 1995, 1996; Phillips et al. 1999; Guy et al. 2005).

The photometric error we use for LSST comes from Ivezic et al. (2008) and takes the form of

$$\begin{equation} \sigma _{\rm phot}^2 = \sigma _{\rm sys}^2 + \sigma _{\rm zero}^2 + \sigma _{\rm rand}^2, \end{equation} \tag{ 3 }$$

where σ_sys is the systematic photometric error which is designed to be very small (<0.005 mag), σ_zero is the absolute photometric error that we set to be σ_zero = 0.02 mag (Astier et al. 2006; Sullivan et al. 2006). We note that, in practice, there is only one zero-point realization in any given experiment that is applied to all SNe. This would result in a non-diagonal covariance matrix for the distance modulus (Kim & Miquel 2006). However, since the inclusion or non-inclusion of this covariance does not change the principles of our methodology, for simplicity, we ignore this correlation here. We do suggest that it must be considered in an analysis of real data.

In Equation (3), σ²_rand is the random photometric error for point sources given by

$$\begin{equation} \sigma _{\rm rand}^2=(0.04-\gamma)x+\gamma x^2. \end{equation} \tag{ 4 }$$

Here, γ is a parameter related to the sky brightness and readout noise, among others, and $x=10^{0.4(m-m_5)}$, where m is the magnitude and m₅ is 5σ depth for a detection of a point source in each of LSST bands. The m₅ is a function of the sky brightness, the seeing, the exposure time, atmospheric extinction, the air mass, and the overall throughput of the instrument. All of the values of these parameter can be found in Table 2 of Ivezic et al. (2008).

We randomly generate the first observational day from −20 to 35 rest-frame days to ensure that we always have enough data to establish the stretch factor. We next randomly select the data point to occur every 3 or 4 days based on the cadence of the LSST (Ivezic et al. 2008). Finally, about 1000 mock SNIa flux data are generated for each of the five filters.

In Figure 2, we show the examples of the mock light curves in g, r, and i bands at different redshifts. The mock flux is created from the Gaussian distribution with the mean on the light curve. Here, we set the first observe day t^obs₀ = −10 in the observer frame.

Zoom In Zoom Out Reset image size

The best-fit value for each light curve parameter usually can be found using the nonlinear least-squares fitting technique (e.g., Sullivan et al. 2006). Here, considering the number of the parameters, the efficiency, and the accuracy, we would like to employ the MCMC technique to perform the fitting process. This method does not require to assume the Gaussian distribution for the likelihood, and it is easy to perform the marginalization over other parameters when quoting error for one parameter. Most importantly, it is very efficient for the multiparameter fitting (Neil 1993; Lewis & Bridle 2002; MacKay 2003; Doran & Mueller 2004; Gong & Chen 2007; Trotta 2008).

Our purpose is to estimate the posterior probability P(θ|D) for the parameter set θ given the observational data set D. Based on the Bayes theorem

$$\begin{equation} P(\theta |{\bf D}) = \frac{\mathcal {L}({\bf D}|\theta)P(\theta)}{P({\bf D})}, \end{equation} \tag{ 6 }$$

where $\mathcal {L}({\bf D}|\theta)\sim e^{-\chi ^2/2}$ is the likelihood which denotes the probability to get D given the parameters θ, P(θ) is the prior probability for θ, and P(D) is the normalization factor which would not affect our analysis here.

The Metropolis–Hastings algorithm is applied in our MCMC technique to decide if a new point should be accepted by an acceptance probability:

$$\begin{eqnarray} {\bf a}({\theta _{n+1} | \theta _n}) & = & \min \Bigg \lbrace \frac{P(\theta _{n+1} | {\bf D})\; {\bf q}(\theta _n | \theta _{n+1})}{P(\theta _n|{\bf D})\; {\bf q}(\theta _{n+1}|\theta _n)}\, 1\Bigg \rbrace \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} & = & \min \Bigg \lbrace \frac{{\mathcal {L}}({\bf D}|\theta _{n+1})\; {\bf q}(\theta _n|\theta _{n+1})}{{\mathcal {L}}({\bf D}|\theta _n)\; {\bf q}(\theta _{n+1}|\theta _n)}\, 1\Bigg \rbrace, \end{eqnarray} \tag{ 8 }$$

where q(θ_n+1|θ_n) is the proposal density to propose a new point θ_new given a current point θ_n in the chain. Here, we assume uniform prior probabilities for the parameters which are canceled in Equation (8). If a = 1, the new point θ_new is accepted; otherwise, the new point is accepted with probability a. This process is repeated until a new point is accepted, and then we set θ_n+1 = θ_new. Also, we set a uniform Gaussian-distributed proposal density for every point, so that it is independent of the position on the chain, i.e., q(θ_n+1|θ_n) = q(θ_n|θ_n+1), we then have

$$\begin{equation} {\bf a}({\theta _{n+1}|\theta _n}) = \min \Bigg \lbrace \frac{{\mathcal {L}}({\bf D}|\theta _{n+1})}{{\mathcal {L}}({\bf D}|\theta _n)}\, 1\Bigg \rbrace. \end{equation} \tag{ 9 }$$

Since the proposal density determines the step size of the MCMC process, it is closely related to the convergence and mixing of the chain. Here, we adopt the adaptive step size Gaussian sampler given by Doran & Mueller (2004). The criterion of the convergence we use was described in Gelman & Rubin (1992), and after convergence we freeze the step size (Doran & Mueller 2004).

The ranges of the parameters when we run the MCMC are set as follows: z ∈ (0, 2), m_B ∈ (10, 30), s ∈ (0.5, 1.5), E^host_B−V ∈ (−0.5, 0.5), Δm ∈ (−0.5, 0.5), and t^rest₀ ∈ (−20, 40). For each mock SNIa, we take about 10, 000 chain points to illustrate the probability distribution of the parameters after the burn-in and thinning process.

In Figure 3, we compare the input redshift of each of our 1000 SNeIa in the simulation with the redshift obtained from MCMC fitting of SEDs to the multiwavelength light curves. We find that when z < 0.2 the SNeIa light curves are adequately sampled with enough accuracy to allow good redshift estimates along with other parameters, with uncertainties as small as 0.001. Such an error is comparable with the spectroscopic measurements, and even if the spectroscopic measurements could provide a higher precision on the measurement of the redshift, in any case the unknown bulk flows (Cooray & Caldwell 2006; Zhang & Chen 2008) would produce an error on the redshift at this level. For the medium redshift 0.2 < z < 0.5, the estimated redshift is still useful but the uncertainty is about 0.1. For z>0.5, the apparent magnitude becomes large, and since σ_phot ∼ 10^0.4m, the redshift errors increase quickly with increasing redshift and can reach ∼1. The limitation at high redshift is also due to lack of near-IR photometric coverage and addition of IR bands beyond the z band will improve photometric determinations when z>0.5.

Zoom In Zoom Out Reset image size

The residuals and 1σ errors for the total six fitting parameters in the MCMC analysis are shown in Figure 4. Similar to the redshift, the multiwavelength light curve model fitting leads to parameter accuracies that are remarkably accurate when z < 0.2, except for Δm as it acts as an extra source of noise independent of the redshift. As shown in Figure 5, over the whole redshift range studied out to z of 1.1, the dispersion of the fitting z, m_B, s, E^host_B−V, Δm, and t^rest₀ are mainly less than ±0.3, ±0.4, ±0.2, ±0.2, ±0.4, and ±4, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To establish the overall effect of the uncertainty from photometric redshift for cosmological studies, we also generate 1000 SNeIa with spectroscopic redshifts z_spec with the LSST photometric error σ_phot, i.e., we just fix the redshift and only model fit the other five parameters. The Hubble diagram for the spectroscopic and photometric cases is shown in Figure 6. Only one-sixth of the whole data are shown in each figure.

Zoom In Zoom Out Reset image size

We use the MCMC approach to fit the cosmological parameters from the two Hubble diagrams. Two cosmological scenarios are considered, that first one is $\rm \Lambda CDM$ with non-flat geometry and the second is $\rm wCDM$ with the time-evolved equation of state for the dark energy with w(z) = w₀ + w₁z/(1 + z).

In Figure 7, we show the contour maps of Ω_m0 versus Ω_Λ0 and w₀ versus w₁ with 1σ and 2σ errors. As can be seen, the 1σ contours using z_phot and its error in cosmological parameter fits nearly overlap with the 2σ contours of the case where redshift is known precisely using z_spec. Also, we find little deviation for the directions of the main axis of the contours for the two cases. Thus, for a survey such as those planned for LSST, we effectively find a factor of ∼ 2 degradation in parameter uncertainties when using the SNIa sample with only photometric redshifts compared with one with also spectroscopic redshifts.

Zoom In Zoom Out Reset image size

In terms of the dark energy figure of merit that involves the inverse area of the w₀ versus w₁ ellipse, photometric SNe samples lead to a factor of 4 degradation compared to spectroscopic sample. This difference, however, is likely to be a minor issue: compared with the planned SNeIa surveys that will involve spectroscopic measurements of a few thousand SNe per year, photometric-only surveys such as the one with LSST will produce a sample of a few hundred thousand SNe. Moreover, we have to note that for a real survey the sample is always magnitude limited, so that some z>1.1 objects could contaminate the z < 1.1 sample and lead to a bias in the fitting results. Given that we cannot quantify this bias fraction, we don not include this effect in our analysis. In an upcoming paper, we hope to implement a new technique to account for such biases in large SNe samples.

To estimate the level of contamination, we create mock samples of light curve data for SNIb/c and SNII. The spectral templates for SNIb/c are from Levan et al. (2005), and for SNII we use the templates of SNIIP and SNIIL given by Gilliland et al. (1999) and Baron et al. (2004).¹⁰ Here, for simplicity, we consider primarily the SNIIP and SNIIL. The SNIIn which have "unusual" progenitors (Mobberley 2007) may be an important contamination for the SNIa (Poznanski et al. 2007a). We may discuss these objects in future work. We set the percentage of SNIIP and SNIIL as 50% and 50%, respectively, for the SNII sample. The mock flux data of the SNIb/c and SNII samples are generated with the same procedure as those used to generate the SNIa mock data in Section 2.

Since the core-collapse SNe are intrinsically fainter than the SNeIa and have no magnitude–phase relation, we take the absolute peak magnitude from Richardson et al. (2002) and the Gaussian distribution with $\bar{z}=0.4$ for SNIb/c and SNII, and then set s = 1 when mimicking their B-band peak magnitude. Also, given that the core-collapse SNe are usually found in star-forming regions, they are expected to suffer more extinction from the host galaxy. We set −0.2 < E^host_B−V < 0.4 with $\overline{E_{B-V}^{\rm host}}=0.0$ and σ_E = 0.3. Once simulated, we continue to use the SNIa SED light curves to fit them just as we did in Section 3.

The distribution of the difference of χ²_min (i.e., relative χ²) for the SNIa, SNIb/c, and SNII is shown in Figure 8. We find that when fitted with Ia SED light curves, the χ² values for SNeII are so large that they are easily distinguished from the SNIa even with photometric data alone. However, for SNeIb/c the χ² peak overlaps with that of the SNeIa, so they are the cause of primary contamination to the total sample.

Zoom In Zoom Out Reset image size

To obtain a less-contaminated sample, as a first cut we note that the χ² distribution of the SNeIb/c has a long tail which can extend to tens of thousands, and some of the SNeIb/c can be removed by an overall restriction on the χ² values in the fitting to SNIa light curve template. If the selection is restricted to χ²_rel < 20, keeping all real Type Ia's, this results in a removal of about 40% of the SNeIb/c's. Since the ratio of the rate of SNIa to SNIb/c out to z ∼ 1 is about 10 to 7 (Calura & Matteucci 2006; de Plaa et al. 2007; Sato et al. 2007; Poznanski et al. 2007b; Eldridge et al. 2008; Smartt et al. 2009; Georgy et al. 2009), we expect about 250 SNeIb/c to remain and contaminate a sample that contains 1000 SNeIa selected photometrically. Note that here we just use a simplified assumption with the ratio of the rate of SNIa to SNIb/c to be redshift independent. In the next section, we discuss a statistical method to further reduce the contamination of Ib/c during model fits to the total sample.

We employ the Bayesian estimation method proposed by Press (1996) and Kunz et al. (2007) to further reduce the contamination from SNIb/c. We note that our proposed statistical method cannot distinguish each SNIb/c from a Type Ia individually, but statistically it reduces the overall contamination and the associated bias in cosmological parameters. As we illustrate here, the same method also allows us to jointly estimate the fraction of the SNIa (or Ib/c's) within the whole SNe sample used for cosmology.

We take the case that the observational sample of supposedly Type Ia's D contains a mixture of true SNeIa data D₁ and SNeIb/c's D₀ which mimics Ia. We define a vector v of length the total number of SNe N with the value of v_i taking either 1 or 0 if D_i is or not a SNIa. We also define a quantity p to account for the total fraction of the true SNeIa in the total SNe data set D. Using v and p, the posterior probability can be written as

$$\begin{eqnarray} P(\theta |{\bf D}) &=& \sum _{{p,\bf v}}P(\theta, {\bf v}, p|{\bf D}) \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} &\propto & \sum _{p,{\bf v}}\mathcal {L}({\bf D}|\theta,{\bf v},p)P(\theta, {\bf v}, p) \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} &\propto & \sum _{p,{\bf v}}\mathcal {L}({\bf D}|\theta,{\bf v},p)P(p)P(\theta |p)P({\bf v}|\theta, p). \end{eqnarray} \tag{ 12 }$$

In Equation (10), the sum over p will be the integration if the value of p is continuous, and the sum of v goes through all 2^N possible values of v. Equation (11) is derived from the Bayes theorem, and $\mathcal {L}({\bf D}|\theta, {\bf v},p)$ is the likelihood. The P(θ|p) in Equation (12) can be reduced to P(θ) since there is no reason to believe the parameter p affect the cosmological evolution of the universe (it is not a cosmological parameter).

Thus, we can simplify to

$$\begin{equation} P(\theta |{\bf D}) \propto P(\theta)\sum _p P(p)\sum _{\bf v} \mathcal {L}({\bf D}|\theta,{\bf v},p)P({\bf v}|\theta,p). \end{equation} \tag{ 13 }$$

For any value of p, P(v|θ, p) is 0 when v_i involving the ith datum is a SNIb/c. When normalized, P(v_i = 1|θ, p) = p and P(v_i = 0|θ, p) = 1 − p. Therefore, we find

$$\begin{equation} P(\theta |{\bf D})\propto P(\theta)\sum _p P(p)\sum _{\bf v} \left[\prod _{v_i=1}\mathcal {L}_i^1p\prod _{v_i=0}\mathcal {L}_i^0(1-p)\right]. \end{equation} \tag{ 14 }$$

Here $\mathcal {L}_i^1$ is the likelihood that the ith SN is a Ia and this is taken to be

$$\begin{equation} \mathcal {L}_i^1 = \frac{1}{\sqrt{2\pi }\sigma _i}e^{-\chi _i^2/2}, \end{equation} \tag{ 15 }$$

where σ_i is the error and χ²_i = (μ^obs_i − μ^th_i)²/σ²_i, and the μ^obs_i and μ^th_i are the observational and theoretical distance moduli, respectively.

The $\mathcal {L}_i^{0}$ is the likelihood for the SNIb/c samples though we don not know a priori the exact distribution. We use two parameters b and σ₀ belonging to the parameter set θ to describe P⁰_i as

$$\begin{equation} \mathcal {L}_i^0 = \frac{1}{\sqrt{2\pi }\sigma _0}e^{-\chi _0^2/2}, \end{equation} \tag{ 16 }$$

where χ²₀ = (μ^obs_i − μ^th_i − b)²/σ²₀.

We can simplify Equation (14) further by noting that the 2^N summation term can be written as the product of N terms. We finally get

$$\begin{equation} P(\theta |{\bf D})\propto P(\theta)\sum _p P(p)\prod _N \big[\mathcal {L}_i^1 p + \mathcal {L}_i^0(1-p)\big]. \end{equation} \tag{ 17 }$$

The sum over p is easily performed with MCMC runs, and we assume P(p) is a uniform distribution. When analyzing our mock samples we take the ranges for p, b, and σ₀ of p ∈ (0.5, 1), b ∈ (−20, 20), and 1/σ₀ ∈ (0, 1000).

Extending the discussion in Section 4.1, we add 250 SNIb/c data with χ²_rel < 20 to the 1000 SNIa data, and extract cosmological constraints on the time-evolving equation of state of dark energy with an analysis which implements the Bayesian estimation method described above, in addition to a method where all data are analyzed with a MCMC run without making an attempt to account for Ib/c contamination to the total sample.

We show the contour maps of w₀ versus w₁ in Figure 9. The red solid and blue dotted contours are the fitting results with and without the Bayesian estimation, respectively. As can be seen in Figure 9, the constraint on w₀ and w₁ is completely wrong when we ignore the Ib/c contamination in the total sample and just do direct fitting to the Hubble diagram. This is caused by the large differences of the distribution between the SNIa and the SNIb/c data and MCMC chains are easily trapped in a wrong likelihood value. The χ² value for the overall best fit in this case is also very large reaching as high as 2000.

Zoom In Zoom Out Reset image size

When we implement the Bayesian estimation method, the result is improved significantly. Although there is a difference between the best-fit and the actual (fiducial) value, the fiducial value of the cosmological parameter set (w₀, w₁) = (−1, 0) lies safely within the 2σ contour around the best fit. Also, as discussed, we also jointly estimate the fraction of SNeIa in the total data set. We plot the likelihood for P(p) in Figure 10. The fraction of the SNIa in this particular mock data should be 80% while the fitting leads to the result of p = 0.75^+0.05_−0.06 with errors at 1σ. While there still remains a bias associated with the contaminating Type Ib/c's, we have reduced this bias to the level of a few percent. Also, if the distributions of b and σ₀ are better measured in the future, this method would get better typing result.

Zoom In Zoom Out Reset image size

We believe that the Bayesian estimation method provides a useful statistical tool to reduce contamination and to evaluate the fraction of the contaminating SNe in the LSST photometric SNIa survey. Of course, to identify whether each individual SN is a core-collapse one or a Ia, the method discussed above is inadequate, but the Bayesian statistical analysis can be a valuable guide for further advanced study (Poznanski et al. 2002, 2007a; Gal-Yam et al. 2004; Johnson & Crotts 2006; Kuznetsova & Connolly 2007).