Type Ia Supernova Distances at Redshift >1.5 from the Hubble Space Telescope Multi-cycle Treasury Programs: The Early Expansion Rate

SNe Ia measure distances most directly; roughly speaking, each SN provides an independent measurement of the luminosity distance to its redshift. For a flat universe, we have

$$\begin{eqnarray}&&{d}_{L}(z)=\displaystyle \frac{c}{{H}_{0}}(1+z){\int }_{0}^{z}\displaystyle \frac{{dz}^{\prime} }{E(z^{\prime} )};\end{eqnarray} \tag{ 1 }$$

therefore the (inverse) Hubble parameter, the derivative of the comoving distance, must be inferred indirectly when starting from raw SN Ia data.

A variety of interrelated methods have been used for this purpose. Some analyses have focused on model-independent reconstruction of an analytical E(z) function or of other dynamical quantities like the deceleration parameter q(z) (Sahni & Starobinsky 2006; Shafieloo et al. 2006; Shafieloo 2007; Ishida & de Souza 2011). Such reconstructions are useful for understanding where the data are most constraining, and they can indicate whether the functional form for H(z) naturally preferred by the data is consistent with that of a physical model like ΛCDM. On the other hand, it is not possible, or at least not straightforward, to subsequently incorporate the reconstructions in a likelihood function, or otherwise in a statistical analysis, in order to constrain cosmological parameters.

Other methods focus on obtaining direct measurements of E(z) at several redshifts by smoothing and/or weighting the individual SNe and differentiating the distance-redshift relation (Tegmark 2002; Daly & Djorgovski 2003, 2004; Daly et al. 2008). One proposed method (Wang & Tegmark 2005), which has been employed in some subsequent analyses (Riess et al. 2007; Avgoustidis et al. 2009; Mortsell & Clarkson 2009), seeks direct, independent estimates of E(z) in redshift bins by first converting SN distance moduli into their corresponding comoving distances r_i, then transforming these r_i into noisy, but locally unbiased, estimates of $E{(z)}^{-1}$ between neighboring SNe. A specific weighted average then yields a minimum-variance estimate of $E{(z)}^{-1}$ over a wider redshift bin. We have verified numerically²⁶ that this procedure is actually equivalent to the familiar weighted least-squares fit of a line to the r_i versus z_i data over the same wide redshift bin, where the slope corresponds to $E{(z)}^{-1}$. Both the least-squares estimator and that of Wang & Tegmark (2005) have been shown to be unbiased and have minimum variance, assuming SN redshifts are known exactly and E(z) is constant over the redshift bin, so it is not surprising that these estimators coincide.

While such an approach is attractive in that it directly transforms the SN distances into independent measurements of E(z) at different redshifts, it has notable problems that make it unsuitable in practice. The first step requires converting SN distance moduli into comoving distances, and one must therefore assume a value for the intercept of the Hubble diagram, which is unknown a priori. As this quantity is partially degenerate with E(z), particularly the lowest-redshift measurement, fixing the intercept to some best-fit value would artificially remove a degree of freedom from the fit, resulting in underestimated uncertainties. One could instead interpret the estimates as estimates of ${AE}{(z)}^{-1}$, where A is an arbitrary constant. In this case, though, properly extracting cosmological information from the E(z) measurements would require fully marginalizing over A in a fit to multiple measurements of ${AE}{(z)}^{-1}$.

Furthermore, an E(z) estimate using this method reflects some average of E(z) over the redshift bin, not necessarily the value at the bin's center. Unless E(z) is constant over the redshift bin, this will lead to a bias, and since only a handful of E(z) values can be constrained robustly with current data, one might expect the bias to be significant. Indeed, by simulating instances of our SN data (see Section 3.2), we have verified that biases in such E(z) estimates are typically a large fraction (∼0.5) of their uncertainty, making the measurements unsuitable for later cosmological inference.

We now describe a somewhat different approach for determining E(z) from SN Ia data. We explain how it avoids the problems discussed in Section 3.1 and provides more meaningful and robust E(z) measurements. We will assume that the true, underlying E(z) function is a continuous, smooth function of redshift, which is certainly the case for most physical and empirical models studied in the literature.

In our approach, we parametrize E(z) by its value at several specific redshifts and employ a basic interpolation scheme to define the complete E(z) function, which can then be numerically integrated to compute the luminosity distance and compare to the data. This allows us to constrain the E(z) parameters using the full SN data set in its raw form, as one would in a standard dark energy analysis. This way, any nuisance parameters associated with the SN data, notably the distance scale or Hubble diagram intercept, can be properly marginalized over in the fit.

While the total number of E(z) values to constrain is somewhat arbitrary, there are several considerations. Choosing too many E(z) parameters results in weaker constraints and posterior distributions that are less likely to be Gaussian. Choosing too few E(z) values increases the chance that the estimates will be biased, as the interpolating function will deviate from the functional form of the underlying cosmology. The specific redshifts, while also somewhat arbitrary, should reflect the redshift range and distribution of the SNe. Since the E(z) measurements, especially those at neighboring redshifts, will naturally be somewhat correlated, choosing too small a separation in redshift between a given pair will lead to undesirably large pairwise correlations in the estimates. We choose redshifts such that the resulting E(z) correlation coefficients are ≲0.5 to avoid such redundancy.

Overall, we find that employing a shape-preserving piecewise-cubic Hermite interpolating polynomial (implemented as pchip in MATLAB; see Kahaner et al. 1988) to interpolate (and extrapolate) the E(z) function works particularly well, though other interpolation schemes (various splines, simple linear interpolation) are also generally suitable. For any specified E(z) (any fiducial cosmology), it is straightforward to determine whether the E(z) estimates resulting from the interpolation and fitting procedure are unbiased. To check this, we repeatedly simulate instances of our SN Ia data; that is, we keep the same SN redshifts and covariance matrix as the real SN data, but repeatedly sample the distance moduli from a multivariate Gaussian centered on the fiducial cosmology. Of course, unbiased constraints for the fiducial cosmology do not guarantee unbiased results for other cosmologies. In principle, one could perform this check for each specific model of interest; however, there is reason to worry only when a model predicts E(z) to vary rapidly or have features too narrow to be captured by the widely spaced E(z) parameters. For the highest-redshift SNe Ia, a modest amount (∼25%) of the integral of $E{(z)}^{-1}$ must be evaluated via extrapolation beyond the last redshift anchor of the E(z) function. However, our simulations indicate that this does not bias this highest-redshift measurement of E(z). Indeed, we have verified that all of the E(z) measurements are biased by ≲10% of their individual statistical uncertainties.

In essence, our procedure trades the ability to make direct, independent measurements of E(z) at redshifts that are somewhat uncertain (and not randomly so) for the ability to obtain precise, unbiased E(z) estimates at specific redshifts. Only this latter type of estimate allows for accurate subsequent cosmological inference with the E(z).

We now constrain E(z) for the Pantheon compilation of 1040 SNe Ia, which we will supplement with the high-redshift CANDELS and CLASH SNe. The Pantheon compilation (Scolnic et al. 2017) includes data from multiple surveys (CfA(1–4), CSP, SDSS, SNLS, Pan-STARRS1, HST) calibrated for a joint cosmological analysis. Below we summarize the key aspects of the Pantheon analysis, and we refer the reader to Scolnic et al. (2017) for additional details and a complete discussion.

The Pantheon analysis presents the full set of spectroscopically confirmed SNe Ia from the Pan-STARRS1 (PS1) Medium Deep Survey, building on the earlier analysis of the first 1.5 yr of PS1 (Rest et al. 2014; Scolnic et al. 2014). It relies on the Supercal cross-calibration procedure presented by Scolnic et al. (2015), which uses the relative consistency of the Pan-STARRS1 photometry over 3π steradians of the sky to tie together the photometric systems of the individual surveys. The Pantheon analysis also incorporates the BBC methodology of Kessler & Scolnic (2017) (see also Scolnic & Kessler 2016), which corrects for distance biases dependent on the light-curve properties of the SNe and the surveys from which they are selected.

The Pantheon analysis employs the SALT2 light-curve fitter (Guy et al. 2007; Betoule et al. 2014), which determines an overall normalization of the log-flux (m_B), a shape parameter (x₁), and a color (c) for each SN light curve, along with associated uncertainties. We standardize the SNe by modeling an individual SN Ia distance modulus as

$$\begin{eqnarray}&&\mu ={m}_{B}-M+\alpha \,{x}_{1}-\beta \,c+{{\rm{\Delta }}}_{M}+{{\rm{\Delta }}}_{B}.\end{eqnarray} \tag{ 2 }$$

The ${{\rm{\Delta }}}_{M}$ term is an additional correction for the empirical host-mass step, where SNe in high-stellar-mass host galaxies (${\mathrm{log}}_{10}({M}_{* }/{{\rm{M}}}_{\odot })\gtrsim 10$) are ∼0.05 mag brighter on average, after light-curve standardization.²⁷ The ${{\rm{\Delta }}}_{B}$ term represents the distance bias correction. Note that M, α, β, and the amplitude of the mass step (included in the ${{\rm{\Delta }}}_{M}$ term) are all nuisance parameters that must be determined by a fit to the data. In our analysis, only M (effectively, the Hubble diagram offset) is fit along with the cosmological parameters E(z). The other parameters are well determined independently of cosmology in the Pantheon analysis. The inferred values are $\alpha \approx 0.15\mbox{--}0.16$ and $\beta \approx 3.0\mbox{--}3.7$, where the results vary depending on the intrinsic scatter model.²⁸ Finally, note that the distance modulus as predicted by the cosmological model is given by

$$\begin{eqnarray}&&\mu (z)=5{\mathrm{log}}_{10}\left[\displaystyle \frac{{d}_{L}(z,{\boldsymbol{p}})}{1\,\mathrm{Mpc}}\right]+25,\end{eqnarray} \tag{ 3 }$$

where d_L is the luminosity distance, which is a function of redshift and also depends on the set of cosmological parameters ${\boldsymbol{p}}$.

The statistical uncertainties of SN distance moduli are modeled, in the standard way, as a combination of observational measurement uncertainty, intrinsic scatter, and additional scatter due to gravitational lensing, peculiar velocities, and redshift measurement uncertainty.²⁹ The inferred value for the intrinsic scatter is ${\sigma }_{\mathrm{int}}\approx 0.1$, although, like α and β, the value depends on the intrinsic scatter model. After bulk-flow corrections are applied to the low-redshift SNe, we add a peculiar-velocity scatter of ${\sigma }_{v}=250$ km s⁻¹. We assume a value ${\sigma }_{\mathrm{lens}}=0.055z$ for the lensing scatter (Jönsson et al. 2010). Note that the distribution of the shift in observed magnitude due to lensing is non-Gaussian (e.g., Jönsson et al. 2006), with a tail of strongly magnified SNe; however, by examination of foreground structures we have verified that none of our CANDELS or CLASH SNe are likely to fall in this tail, making the lensing scatter contribution to the distance uncertainty a good approximation. Note that there is also statistical uncertainty in the host-mass correction and the distance bias correction.

The Pantheon analysis also includes a rigorous analysis of systematic errors, adding terms to the covariance matrix of SN distances to account for uncertainties in photometric calibration (including terms for individual survey calibration, the Supercal cross-calibration procedure, and the SALT2 model itself), the intrinsic scatter model, survey selection functions, Milky Way dust extinction, β evolution, the host mass step and its evolution, and peculiar velocity coherent flow corrections.

Standard data-quality cuts were applied to remove SNe that are not expected to follow the empirical standardization relations. Specifically, we keep only SNe with $| {x}_{1}| \lt 3$, ${\sigma }_{{x}_{1}}\lt 1$, $| c| \lt 0.3$, a light-curve fit with ${\chi }^{2}/{N}_{\mathrm{dof}}\lt 3$, and an uncertainty in the time of peak brightness of less than 2 days. Similar cuts have been used in most recent SN Ia cosmological analyses (e.g., Betoule et al. 2014; Rest et al. 2014; Riess et al. 2016). These cuts eliminate 3 of the silver and gold MCT SNe (CLH11Tra, GND13Gar, GND13Jay; see Table 4). Finally, a $4\sigma$ outlier rejection from the best-fit Hubble diagram is applied and removes GND13Cam, leaving 9 HST MCT SNe Ia in the joint analysis.³⁰ Note that here we do include EGS13Rut, which is on the edge of the ${\sigma }_{{x}_{1}}$ cut but has typical light-curve fit parameters. Although the final MCT addition of 9 SNe represents $\lt 1$% of the combined sample, the unusually high redshifts (7 with $z\gt 1.5$) provide unique leverage on E(z) at z = 1.5.

Following the methodology and discussion in Section 3.2, we parametrize $E{(z)}^{-1}$ by its value at six redshifts (chosen to best summarize the sample) and therefore have six free parameters to constrain. It is important to remember that the Hubble diagram offset is a free parameter as well, though we analytically marginalize over this offset, with a flat prior, in the likelihood. We assume a flat universe (${{\rm{\Omega }}}_{k}=0$) throughout, so the E(z) measurements are cosmological-model-dependent in this sense. To obtain the constraints, we sample the likelihood using a custom Markov chain Monte Carlo (MCMC) code employing the basic Metropolis–Hastings algorithm. We impose flat, hard-bound priors on the $E{(z)}^{-1}$ parameters wide enough that extending the bounds does not affect the resulting constraints. The final MCMC chains were inspected to verify convergence.

The resulting marginalized posterior likelihoods for $E{(z)}^{-1}$ are Gaussian to a good approximation, and the constraints are given in Table 6. In Figure 1, we convert the measurements of $E{(z)}^{-1}$ into E(z) measurements by reprocessing the MCMC chains and then compare the results with and without the MCT SNe. It is not surprising that the MCT SNe subsantially improve the measurement of E(z) at z = 1.5. They permit a ∼20% measurement of $E(z=1.5)$, roughly a factor of three improvement over the result without the MCT SNe. While the CANDELS and CLASH SNe mostly affect the measurement at z = 1.5, they also improve and shift some lower-redshift measurements, which are somewhat correlated (≈8% and 4% improvements at z = 0.9 and 0.55, respectively). By eye, the set of E(z) measurements may appear somewhat discrepant with the fiducial ΛCDM model, but the overall ${\chi }^{2}$, which includes the moderate correlations, is 6.7 for the 6 degrees of freedom.

Zoom In Zoom Out Reset image size

Table 6.  Pantheon + MCT SN Ia Measurements of E(z)

z	$E{(z)}^{-1}$ a	Correlation Matrix						E(z)	Distance Residual ${\rm{\Delta }}\mu /(0.01\,\mathrm{mag}$)b
0.07	1.007 ± 0.024	1.00						0.994 ± 0.023	−0.13 ± 0.99
0.2	0.898 ± 0.016	0.40	1.00					1.113 ± 0.020	−0.23 ± 1.26
0.35	0.893 ± 0.029	0.52	−0.13	1.00				1.122 ± 0.037	+0.23 ± 1.32
0.55	0.732 ± 0.033	0.35	0.35	−0.18	1.00			1.369 ± 0.063	+0.11 ± 1.97
0.9	0.652 ± 0.051	0.02	−0.08	0.19	−0.41	1.00		1.54 ± 0.12	+1.15 ± 2.85
1.5	0.337 ± 0.078	0.00	−0.06	−0.05	0.16	−0.21	1.00	${2.69}_{-0.52}^{+0.86}$	−3.42 ± 6.78

Notes.

^aMean and standard deviation of the marginalized likelihood, approximately Gaussian in all cases. ^bEffective distance moduli relative to those of a fiducial ΛCDM cosmology (${{\rm{\Omega }}}_{m}=0.3$), as determined by an interpolated fit to the residuals using the same redshift control points as the E(z) analysis.

Download table as:  ASCIITypeset image

In Figure 2, we scale E(z) by ${(1+z)}^{-1}$ to illustrate the constraints on the time derivative of the scale factor $\dot{a}(z)$, relative to its present value, for the same data shown in Figure 1. In this space, it is clear that the low-redshift and high-redshift E(z) measurements together provide evidence for both recent acceleration and earlier deceleration epochs, as predicted by standard cosmological models. In addition to the fiducial ΛCDM model, we show dynamical models with fixed deceleration parameter q₀. The $\dot{a}(z)$ values track the ${q}_{0}=-0.5$ model at $z\lesssim 0.5$ (where the low-z behavior matches that of a ΛCDM model with ${{\rm{\Omega }}}_{m}\approx 0.3$) but show deceleration with respect to that curve at higher redshifts. The coasting cosmology (${q}_{0}=0$), pure acceleration cosmology (${q}_{0}=-0.5$), and pure deceleration cosmology (${q}_{0}=0.5$, equivalent to a flat CDM model with ${{\rm{\Omega }}}_{m}=1$) are strongly disfavored with ${\rm{\Delta }}{\chi }^{2}=78.9$, ${\rm{\Delta }}{\chi }^{2}=40.5$, and ${\rm{\Delta }}{\chi }^{2}=357.9$, respectively, for 6 degrees of freedom. The measurement at z = 1.5 alone, while consistent with the other models, disfavors ${q}_{0}=-0.5$ with ${\rm{\Delta }}{\chi }^{2}=14.2$.

Zoom In Zoom Out Reset image size

As an illustration of the power of the E(z) measurements in constraining (spatially flat) cosmologies, we compare constraints on common dark energy parameterizations in Figure 3. Remarkably, the constraints are nearly identical whether the parameters are constrained with the SN Ia data directly or with the E(z) measurements in Table 6. It may not be too surprising that E(z) captures the constraining power of the SNe for simple one-or-two-parameter models. One would not expect the same for fits with many degrees of freedom (e.g., more complicated dark energy models); in practice, however, current and near-future SN Ia data can only meaningfully constrain 2–3 expansion parameters anyway. For models that assume a flat universe and predict fairly smooth, featureless H(z), the E(z) constraints will be an efficient summary of the present SN Ia data.

Zoom In Zoom Out Reset image size

The use of SNe Ia as standardizable candles across redshift relies on the understanding that their uncommonly homogeneous luminosities and colors follow from their nature as carbon–oxygen white dwarfs close to the Chandrasekhar mass. While uncertainty persists regarding how these degenerate stars approach that mass limit, either by accretion from a nondegenerate companion or through the tidal disruption followed by accretion of a degenerate companion, there has long been agreement about this model based on the well-understood physics of degenerate stars (Hoyle & Fowler 1960; Arnett 1969; Colgate & McKee 1969). The thermonuclear detonation of a Chandrasekhar-mass carbon–oxygen white dwarf yields a mass of radioactive nickel whose energy output matches that of a SN Ia (Arnett et al. 1985) and whose modeled nucleosynthesis matches its spectral elements (Nomoto et al. 1984). More recently, prediscovery observations of SN 2011fe, a prototypical SN Ia in M101, demonstrated that the progenitor did not exceed a radius of 2% solar, fully consistent with the expected white dwarf (Li et al. 2011; Nugent et al. 2011; Bloom et al. 2012). Yet the difficulty and low likelihood of ever directly observing a white dwarf system before it becomes a SN Ia leaves enough uncertainty and model freedom to support the consideration of redshift evolution of the standardized SN Ia luminosity.

From SN Ia observations spanning a wide range of redshifts and sampling the epochs when cosmic expansion accelerates and decelerates, it is possible to distinguish such evolution from the uncertain properties of dark energy (Riess & Livio 2006). As an illustration of the power of SNe Ia at $z\gt 1$ to separate evolution from cosmology, we briefly reconsider the analysis of Tutusaus et al. (2017), which shows that power-law cosmology, where the scale factor evolves as $a(t)\propto {t}^{n}$ for some exponent n, is an equally good fit to SN Ia data (primarily at $z\lt 1$) as the ΛCDM model (with ${{\rm{\Omega }}}_{m}$ free) when the standardized luminosity is also allowed to vary with redshift according to some simplistic, empirical models of SN Ia evolution. Although such models are not astrophysically motivated, they may be useful for exploring the separation of other SN distance-dependent effects (e.g., gray extinction) from cosmological parameters.

Here, as an illustration, we consider Model B (${\rm{\Delta }}M(z)=\epsilon {z}^{\delta }$) from Tutusaus et al. (2017) with fixed $\delta =0.3$. We separately fit both ΛCDM and power-law cosmology to our combined (Pantheon + MCT) data; in each case, we fit for the Hubble diagram intercept, a cosmological parameter (${{\rm{\Omega }}}_{m}$ or n), and the amplitude of the assumed intrinsic luminosity evolution. We compare these fits in Figure 4. Fitting only the SNe at $z\lt 1$, a power law with n = 1.1 is a slightly better fit to the SN Ia data than ΛCDM. Indeed, when analyzing the JLA compilation, which features only ∼5 SNe at $z\gt 1$, Tutusaus et al. (2017) claims a mild preference for the power law (note that their analysis also included BAO and H(z) information).

Zoom In Zoom Out Reset image size

In contrast, when we include our 24 SNe at $z\gt 1$, a nearly coasting (marginally accelerating) power-law cosmology (best fit n = 1.04) together with simplistic SN Ia evolution is no longer as good a fit as the ΛCDM model, with a relative probability of $\exp (-{\rm{\Delta }}{\chi }^{2}/2)\approx 20 \%$. Without invoking evolution (that is, fixing $\epsilon \equiv 0$), the ΛCDM model is a much better fit than the power-law model, with the latter strongly disfavored with ${\rm{\Delta }}{\chi }_{{\rm{\Lambda }}\mathrm{CDM}}^{2}=8.3$, a relative probability of 1.6%, when including the new SNe at $z\gt 1.5$. Meanwhile, assuming ΛCDM and fitting for the evolution amplitude yields a value consistent with zero, $\epsilon =0.08\pm 0.15$, so there is no motivation for including it based on astrophysical or empirical considerations. A more comprehensive investigation of SN Ia evolution and cosmology is underway (D. L. Shafer et al., 2018 in preparation).

We note the addition of the MCT SNe to the Pantheon compilation also further reduces the already-low likelihood of the "empty universe" solution where ${{\rm{\Omega }}}_{m}\approx 0$ and ${{\rm{\Omega }}}_{{\rm{\Lambda }}}\approx 0$ in an open ΛCDM universe, a location Nielsen et al. (2016) claimed to be marginally consistent ($\sim 3\sigma$) with SN data alone using unconventional priors on SN distributions, to the boundary of the 6σ contour.

WFIRST was the top space-based recommendation of the 2010 U.S. astronomy and astrophysics decadal survey. The mission is still in formulation, but current plans specify a 2.4 m primary mirror and include a wide-field instrument for cosmology. The cosmology science objectives, as detailed in the most recent report from the Science Definition Team (Spergel et al. 2015), will be accomplished through a combination of SN Ia, galaxy, and weak-lensing surveys.

The WFIRST SN survey is anticipated to yield a large sample of thousands of SNe, many at $z\gt 1$ with precise distances. These SNe will vastly improve upon the high-redshift E(z) measurements available today, allowing nontrivial and precision tests of the ΛCDM model independent of the BAO and weak-lensing constraints in a redshift range that is currently not well constrained.

Here we wish to forecast realistic constraints on E(z) from WFIRST. Typical forecasts (e.g., for dark energy figures of merit) rely on Fisher matrix formalism, which is exact only for Gaussian posterior distributions and otherwise underestimates parameter uncertainties. For SN Ia forecasts, one typically assumes idealized, or roughly estimated, redshift distributions and makes simple assumptions about the measurement error. Here instead we employ a detailed simulation of one potential observing strategy for the WFIRST SN survey (Hounsell et al. 2017). We then constrain the E(z) parameters using the methodology of Section 3.2 that was employed in Section 3.3 for our current Pantheon + MCT data.

For our illustration, we consider the Imaging All-z strategy described by Hounsell et al. (2017). This particular strategy relies on multi-band imaging for classification and assumes follow-up spectroscopy will provide host-galaxy redshifts. Hounsell et al. (2017) also assumes a large external sample of 800 SNe at $z\lt 0.1$. As the size of future systematic uncertainties is hard to predict, Hounsell et al. (2017) simulates a range of scenarios, and here we opt for all-around optimistic assumptions about future systematic errors (for what this entails, see Hounsell et al. 2017). In this scenario, the contribution of systematic errors is not negligible but is subdominant in the error budget.

In Figure 5, we compare our current Pantheon + MCT constraints on E(z) with simulated constraints from the WFIRST Imaging All-z strategy. We find that we are able to constrain E(z) robustly, albeit with moderate pairwise correlations, at 9 redshifts in the range $0.07\lt z\lt 2.5$. In Table 7, we list the percent errors for E(z) corresponding to our forecast. Note that these results are negligibly changed whether we quote percent errors on E(z) or its inverse. We find that WFIRST allows 8 measurements of E(z) at the 1%–3% level, along with a robust but less precise measurement at $z\approx 2.5$. Notably, this is a constraint on the expansion rate at a redshift higher than any SN Ia has even been observed to date.

Zoom In Zoom Out Reset image size