IS THE EXPANSION OF THE UNIVERSE ACCELERATING? ALL SIGNS POINT TO YES

The discovery of the accelerating universe by two teams (Riess et al. 1998; Perlmutter et al. 1999) in the late 1990s was one of the major breakthroughs in cosmology. Using Type Ia supernovae (SNe Ia) as standard candles, both teams independently determined that high-redshift SNe were fainter than expected in a matter-dominated universe, implying the need for a cosmological constant (or more generally, dark energy) to accelerate the expansion of the universe, increasing cosmological distance as a function of redshift.

SNe Ia are not perfect standard candles, however. Work leading up to the discovery (Phillips 1993; Hamuy et al. 1996; Riess et al. 1996; Perlmutter et al. 1997) demonstrated the need for empirical standardization relations. Peak absolute magnitudes correlate with the width of the light curve (broader light-curve SNe are more luminous) and the color of the supernova (redder SNe are less luminous). In the years since, other empirical standardization relations have been noted, including one related to host galaxy stellar mass (Kelly et al. 2010; Sullivan et al. 2010); perhaps driven significantly by the local star formation rate; Rigault et al. 2013).

The Sloan Digital Sky Survey (SDSS) and SuperNova Legacy Survey (SNLS) SN teams have completed the Joint Light-curve Analysis (JLA; Betoule et al. 2014). This analysis incorporates a thorough recalibration of both surveys (Betoule et al. 2013), and the full set of spectroscopically confirmed SDSS SNe Ia (Sako et al. 2014); it represents the most up-to-date large SNe Ia compilation.⁴ A recent claim (Trøst Nielsen et al. 2016, hereafter N16) was made that this data set provides only "marginal evidence" for acceleration. We examine the statistical model N16 used to make this claim, and find it deficient for the task. In particular, a simple (and well-justified) update of the model to better account for changes in the observed SN light curve parameter distributions with redshift significantly increases the statistical strength of the acceleration evidence.

In the case of JLA, the standardization relations employed are light curve width (x₁ in the framework of SALT2, Guy et al. 2007), color (c), and host galaxy stellar mass. The dependent variable is taken to be the rest-frame B-band magnitude (m_B). The light curve parameters are determined by comparing a rest-frame spectral energy distribution model to the observer-frame photometry; similarly, the host stellar mass is estimated from broadband photometry. The cosmological results rely on the ability of the statistical framework to fit the standardization relations (in JLA, these are taken to be linear in x₁ and c, and a step function in host mass), yet the uncertainties (a general term that we take to include unexplained dispersion around the model) in the dependent variable (m_B) and independent variables (x₁, c, host mass) are of similar size. The JLA analysis itself used a frequentist line-fitting procedure with only modest biases in its regime of applicability (Mosher et al. 2014).

In contrast, the statistical model from N16 uses a Bayesian Hierarchical Model (c.f., Gull 1989). In the N16 model, the latent ("true") parameters for each SN are modeled with nuisance parameters, which are marginalized over to obtain an inference of the global parameters. The distribution of the latent parameters must be adequately described by the prior. For example, flat priors on the latent variables cause a bias in the fit (Gull 1989). Making the parameters of the prior ("hyperparameters") part of the model avoids this bias (this multi-level nature gives rise to the name "Hierarchical").

The key shortcoming of the N16 model is that it assumes redshift-independent distributions for x₁ and c. As shown in Figure 1, the observed distributions (plotted points) are far from redshift-independent. Two effects visible in the data—selection effects and the correlation between older host galaxies and narrower light-curve (lower x₁) SNe (Hamuy et al. 1995)—result in a more luminous distribution of SNe as an increasing function of redshift. The selection effects are particularly evident in color, where only bluer SNe (more negative c) are above the completeness limit for the high-redshift end of each ground-based sample. By incorrectly treating these distributions as redshift-independent, N16 biased their latent x₁ and c toward the global mean, effectively removing some of the standardization (Conley et al. 2007; Wood-Vasey et al. 2007; Karpenka 2015).⁵ The JLA sample was corrected for selection bias in Betoule et al. (2014), but only in the sense that SNe that are selected to be more luminous after standardization are adjusted to be less luminous. The bias correction cannot compensate for a deficient standardization, as provided by a constant-in-redshift model of the distributions.

Zoom In Zoom Out Reset image size

2.1. Redshift-independent Distributions

As a starting point, we perform an analysis similar to that in N16, using Hamiltonian Monte Carlo to sample from the posterior (we describe the details in the Appendix). We assume a cosmological model with cold matter and a cosmological constant (ΛCDM). We make four measurements: a ΛCDM universe allowed to have spatial curvature (i.e., ${{\rm{\Omega }}}_{m}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}\ne 1$), and a flat ΛCDM universe (the assumption of flatness is discussed more in the discussion section), each with both sets of model assumptions. We compute the deceleration parameter q₀, (${q}_{0}\equiv {-\tfrac{\ddot{a}}{a{H}^{2}}|}_{t={t}_{0}}$, equal to ${{\rm{\Omega }}}_{m}/2-{{\rm{\Omega }}}_{{\rm{\Lambda }}}$ for a ΛCDM cosmology). We evaluate the statistical significance of acceleration (${q}_{0}\lt 0$) by comparing the 50th percentile of the posterior with the difference of the 50th percentile and the 84th percentile (taken as an estimate of 1σ in the $+{q}_{0}$ direction), and then rounding to 0.1σ. The statistical significance of the acceleration is $3.1\sigma$ with no constraint on curvature, and $8.7\sigma$ assuming a flat universe (see Figures 2, 3, left panels). The difference between this estimate and one derived from explicitly measuring the fraction of samples with ${q}_{0}\gt 0$ is modest (66 posterior samples out of 60,000 have ${q}_{0}\gt 0$).

Zoom In Zoom Out Reset image size

2.2. Redshift-dependent Distributions

Next, we introduce a simple model of the observed distributions as a function of redshift. We allow each source of SN discovery (Nearby, Sloan Digital Sky Survey, SuperNova Legacy Survey, Hubble Space Telescope) to have a linear variation in the mean with redshift (for the Hubble Space Telescope SNe, we use only a constant mean in redshift, as this sample is too small to constrain any variation). This model is shown in Figure 1; the variation with redshift is highly statistically significant. We also try a more flexible model in redshift (Rubin et al. 2015), and it makes only a small difference (the only requirement on the model is to be at least as flexible in redshift as the cosmological model under consideration; Rubin et al. 2015). The statistical significance of the acceleration increases to $4.2\sigma$ and $11.2\sigma$ assuming a flat universe (see Figures 2, 3, right panels). Again, the difference between this estimate and one derived from explicitly measuring the fraction of samples with ${q}_{0}\gt 0$ is modest (only one posterior sample out of 60,000 has ${q}_{0}\gt 0$).

Zoom In Zoom Out Reset image size

2.3. Other Cosmological Models

For a result that relies only on kinematics, we also compute q₀ constraints using the Visser (2004) series expansion of luminosity distance as a function of redshift. We take the first three terms (including q₀ and ${j}_{0}\equiv {\tfrac{\dddot{a}}{a{H}^{3}}|}_{t={t}_{0}}$). The q₀ constraints are illustrated on the third row of Figure 3. Even with flat priors on q₀ and j₀ (allowing the kinematics to venture into regions of parameter space that may be hard to realize dynamically), we find 3.7σ evidence for acceleration (2.8σ with the N16 model). SN data alone can be used to derive constraints on the joint posterior of ${{\rm{\Omega }}}_{m}$ and the dark energy equation of state parameter ($w={P}_{\mathrm{DE}}/{\rho }_{\mathrm{DE}}$), assuming a flat universe (Garnavich et al. 1998; Perlmutter et al. 1999); we compute constraints for this model as well. For simplicity, we take a flat prior on both ${{\rm{\Omega }}}_{m}$ and w, and find strong evidence for ${q}_{0}\lt 0$, as shown in the bottom row of Figure 3. The constraints on q₀ are non-Gaussian, but in both the N16 model and ours, no samples (out of 60,000) have ${q}_{0}\gt 0$.

We next project the constraints from each of our four models to the q₀–$[{j}_{0}-{{\rm{\Omega }}}_{k}]$ plane, shown in Figure 4 (both j₀ and ${{\rm{\Omega }}}_{k}$ contribute linearly at the same order in luminosity distance, so we cannot distinguish them in this plane). The constraints for 2D models (models other than flat ΛCDM) are similar. However, the high q₀/low $[{j}_{0}-{{\rm{\Omega }}}_{k}]$ region is disfavored by the dynamical models, as ${{\rm{\Omega }}}_{m}$–w cannot reach ${j}_{0}\lt -1/8$, and even an empty universe ("Milne Model" Milne 1935) only has ${j}_{0}-{{\rm{\Omega }}}_{k}=-1$.

Zoom In Zoom Out Reset image size

Our results (flat universe ${{\rm{\Omega }}}_{m}={0.298}_{-0.031}^{+0.033}$) are similar to the frequentist JLA analysis (flat universe ${{\rm{\Omega }}}_{m}=0.295\pm 0.034$). This is unsurprising; frequentist and Bayesian analyses will converge to exactly the same results under a set of assumptions not far from those made here (Rubin et al. 2015). We also note that more advanced analyses can better take into account statistical properties of the data (modeling selection effects, nonlinear standardization relations, a redshift-dependent host-mass relation, outliers, and a model of unexplained dispersion incorporating x₁ and c) (Rubin et al. 2015). N16 did not include the host-mass standardization; excluding this only has a small impact on our results.⁶ However, we focus our attention on the N16 model of the x₁ and c distributions, as it is this model that drives the difference from the JLA analysis.

While constraints derived from SNe Ia alone require a ∼30% flatness constraint to push the supernova measurement of acceleration above 5σ, current experiments have constrained curvature to much better precision than 1% (Planck Collaboration et al. 2016). Even constraints on ${{\rm{\Omega }}}_{m}$ (e.g., galaxy clusters, Allen et al. 2011), which imply ${{\rm{\Omega }}}_{m}\gt 0.2$, cut off the tail of the SN-only posterior extending down to a Milne universe and ${q}_{0}\gt 0$, allowing the acceleration to again reach 5σ confidence.⁷ With the combination of current experiments (SNe Ia, baryon acoustic oscillations, cosmic microwave background, and the Hubble constant), ${{\rm{\Omega }}}_{{\rm{\Lambda }}}$ is constrained to be 0.6911 ± 0.0062 (Planck Collaboration et al. 2016). In order to claim that the evidence for acceleration is "marginal," it is necessary to fully reject all measurements of the curvature of the universe, the basic constraints on the matter density of the universe, and other cosmological data sets.

Even without external constraints, this work demonstrates that a more accurate model for the supernova analysis greatly increases the significance of acceleration. We conclude that the analysis in N16 is both incorrect in its method and unreasonable in its assumptions, leading the authors to question a result that is quite secure when addressed properly.

We appreciate the feedback we received from Greg Aldering, Peter Nugent, Saurabh Jha, Saul Perlmutter, Alex Kim, Peter Garnavich, and Mike Hobson. Support was provided by the Director, Office of Science, Office of High Energy Physics, of the U.S. Department of Energy under contract No. DE-AC02-05CH11231 and NASA ROSES-14 WFIRST Preparatory Science program 14-WPS14-0050.

We sample from the posterior using Stan (Carpenter et al. 2016) through PyStan (https://pystan.readthedocs.io). Following Trøst Nielsen et al. (2016), we assume flat priors on all parameters, but require ${{\rm{\Omega }}}_{m}\gt 0$. Our chains are 2500 samples each (after warmup), and show excellent convergence (the diagnostic of Gelman & Rubin 1992 is smaller than 1.01). We run 24 chains for the results with curvature, and 8 chains for the flat universe results.

As in Rubin et al. (2015), in order to speed up sampling, we decompose the light-curve fit covariance matrix into its eigenvectors and sample over the projections onto these (this results in a covariance matrix that has no correlations between SNe). In addition to a fully Bayesian analysis, the only other change we make from N16 is to include the host-mass standardization, as was done in JLA. We make our code available at https://zenodo.org/badge/latestdoi/71841571 (Rubin 2016).