Marginal evidence for cosmic acceleration from Type Ia supernovae

The ‘standard’ model of cosmology is founded on the basis that the expansion rate of the universe is accelerating at present — as was inferred originally from the Hubble diagram of Type Ia supernovae. There exists now a much bigger database of supernovae so we can perform rigorous statistical tests to check whether these ‘standardisable candles’ indeed indicate cosmic acceleration. Taking account of the empirical procedure by which corrections are made to their absolute magnitudes to allow for the varying shape of the light curve and extinction by dust, we find, rather surprisingly, that the data are still quite consistent with a constant rate of expansion.

In the late 1990's, studies of Type Ia supernovae (SN Ia) showed that the expansion rate of the universe appears to be accelerating as if dominated by a cosmological constant [1][2][3] . Since then supernova cosmology has developed rapidly as an important probe of 'dark energy' . Empirical corrections are made to reduce the scatter in the observed magnitudes by exploiting the observed (anti) correlation between the peak luminosity and the light curve width and the colour 4,5 . Other such correlations have since been found e.g. with the host galaxy mass 6 and metallicity 7 . Cosmological parameters are then fitted, along with the parameters determining the light curves, by simple χ 2 minimisation 1,[8][9][10][11] . This method has a number of pitfalls as has been emphasised earlier 12,13 .
With ever increasing precision and size of SN Ia datasets, it is important to also improve the statistical analysis of the data. To accomodate model comparison, previous work [14][15][16] has introduced likelihood maximisation. In this work we present an improved maximum likelihood analysis, finding rather different results.

Supernova Cosmology
There are several approaches to making SN Ia 'standardiseable candles' . The different philosophies lead to mildly different results but the overall picture seems consistent 17 . In this paper we adopt the widely used approach of 'Spectral Adaptive Lightcurve Template 2′ (SALT2) 18,19 wherein the SN Ia are standardised by fitting their light curve to an empirical template, and the parameters of this fit are used in the cosmological analysis. (A more comprehensive statistical model of light curves spanning optical through near-infrared data has subsequently been constructed in a hierarchical Bayesian framework 20 ). Every SN Ia is assigned three parameters, one being ⁎ m B , the apparent magnitude at maximum (in the rest frame 'B-band'), while the other two describe the light curve shape and colour corrections: x 1 and c. The distance modulus is then taken to be: where M is the absolute magnitude, and α and β are assumed to be constants for all SN Ia. These global constants are fitted along with the cosmological parameters. The physical mechanism(s) which give rise to the correlations that underlie these corrections remain uncertain 21,22 . The SN Ia distance modulus is then compared to the expectation in the standard Λ CDM cosmological model:   , the probability density of the true parameters writes: , the observed X , and the estimated experimental covariance matrix Σ d (including both statistical and systematic errors), the probability density of the data given some set of true parameters is: To combine the exponentials we introduce the vector which can be integrated analytically to obtain: This is the likelihood (equation (3)) for the simple model of equation (4), and the quantity which we maximise in order to derive confidence limits. The 10 parameters we fit are . We stress that it is necessary to consider all of these together and Ω m and Ω Λ have no special status in this regard. The advantage of our method is that we get a goodness-of-fit statistic in the likelihood which can be used to compare models or judge whether a particular model is a good fit. Note that the model is not just the cosmology, but includes modelling the distributions of x 1 and c.
With this MLE, we can construct a confidence region in the 10-dimensional parameter space by defining its boundary as one of constant . So long as we do not cross a boundary in parameter space, this volume will asymptotically have the coverage probability cov 0 is the pdf of a chi-squared random variable with ν degrees of freedom, and max  is the maximum likelihood.
To eliminate the so-called 'nuisance parameters' , we set similar bounds on the profile likelihood. Writing the interesting parameters as θ and nuisance parameters as φ, the profile likelihood is defined as Comparison to other methods. It is illuminating to relate our work to previously used methods in SN Ia analyses. One method 14 maximises a likelihood, which is written in the case of uncorrelated magnitudes as  but this cannot be used to compare models, since it is tuned to be 1 per degree of freedom for the Λ CDM model by adjusting an arbitrary error σ int added to each data point. This has been criticised 12,13 , nevertheless the method continues to be widely used and the results presented without emphasising that it is intended only for parameter estimation for the assumed (Λ CDM) model, rather than determining if this is indeed the best model.

Analysis of JLA catalogue
We focus on the Joint Lightcurve Analysis (JLA) catalogue 11 . (All data used are available on http://supernovae. in2p3.fr/sdss_snls_jla/ReadMe.html -we use the covmat_v6.) As shown already in Fig. 1, the distributions of the light curve fit parameters x 1 and ĉ are well modelled as gaussians. Maximisation of the likelihood under specific constraints is summarised in Table 1 and the profile likelihood contours in the Ω m − Ω Λ plane are shown in Fig. 2. In Fig. 3 we compare the measured distance modulus, with its expected value in two models: 'Λ CDM' is the best fit (Table 1) accelerating universe, while 'Milne' is an universe expanding with constant velocity. The error bars are the square root of the diagonal elements of Σ l + A T−1 Σ d A −1 so include both experimental uncertainties and intrinsic dispersion. We show also the residuals with respect to the Milne model (which has been raised to take into account the change in M 0 ).
To assess how well our Gaussian model for the latent variables describes the data, we show the 'pull' distribution in Fig. 4. These are defined as the normalised, decorrelated residuals of the data, where U is the upper triangular Cholesky factor of the covariance matrix Σ d + A T Σ l A. Performing a K-S test, comparing the pull distribution to a unit variance gaussian gives a p-value of 0.1389.
To check the validity of our method and approximations, we do a Monte Carlo simulation of experimental outcomes from a model with parameters matching our best fit (see Table 1). Figure 5 shows the distribution of −2 log [ / ] true max   , which is just as is expected.

Discussion
That the SN Ia Hubble diagram appears consistent with an uniform rate of expansion has been noted earlier 16,[23][24][25] . We have confirmed this by a statistically principled analysis, using the JLA catalogue of 740 SN Ia processed by the SALT2 method. We find marginal (i.e.  σ 3 ) evidence for the widely accepted claim that the expansion of the universe is presently accelerating 3 .
The Bayesian equivalent of this method (a "Bayesian Hierarchical Model") has been presented elsewhere 13 and has recently been applied to the same dataset, finding results consistent with ours 26 . We note that a Bayesian consistency test 27 has been applied (albeit using the flawed 'likelihood' (equation 12) and 'constrained χ 2 ' (equation 13) methods) to determine the consistency between the SN Ia data sets acquired with different telescopes 28 . These authors do find inconsistencies in the UNION2 catalogue but none in JLA. This test had been applied earlier to the UNION2.1 compilation finding no contamination, but those authors 29 fixed the light curve fit 'nuisance' parameters, so their result is inconclusive. Including a 'mass step' correction for the host galaxies of SN Ia 11 has little effect.  While our gaussian model (4) is not perfect, it appears to be an adequate first step towards understanding SN Ia standardisation. One might be concerned that various selection effects (e.g. Malmquist bias) affect the data. Such effects may not be amenable to our approximate method and are better addressed in a Bayesian approach 26 . We are concerned here solely with performing the analysis in a statistically sound manner to highlight the different conclusion from previous analyses 11 of the same data.
Whether the expansion rate is accelerating or not is a kinematic test and it is only for ease of comparison with previous results that we have chosen to show the impact of doing the correct statistical analysis in the Λ CDM framework. In particular the 'Milne model' refers here to an equation of state p = − ρ/3 and should not be taken to mean an empty universe. For example the deceleration due to gravity may be countered by bulk viscosity associated with the formation of structure, resulting in expansion at approximately constant velocity even in an universe containing matter but no dark energy 30 . Such a cosmology is not prima facie in conflict with observations of the angular scale of fluctuations in the cosmic microwave background or of baryonic acoustic oscillations, although this does require further investigation. In any case, both of these are geometric rather than dynamical measures and do not provide compelling direct evidence for a cosmological constant -rather its value is inferred from the assumed 'cosmic sum rule': Ω Λ = 1 − Ω m + Ω k . This would be altered if e.g. an additional term due to the 'back reaction' of inhomogeneities is included in the Friedmann equations 31 .
The CODEX experiment on the European Extremely Large Telescope will aim to measure the 'redshift drift' over a 10-15 year period to determine whether the expansion rate is really accelerating 32 .

Methods: Confidence ellipsoids
The confidence ellipsoid is the collection of points , which obey where  is a symmetric matrix and x MLE is the MLE. The enclosed volume is a confidence region with coverage probability corresponding with high precision to the value obtained from Equation (10). The eigenvectors of  are then the principal axes of the ellipsoid, and the eigenvalues are the inverse squares of the lengths of the principal axes. We approximate this matrix with the sample covariance from the MC of section 3 as  = − x x cov( , ) 1 . To make reading the matrix of eigenvectors easier, we round all numbers to 0.1. Thus, we get the following approximate eigenvectors of  , in columns with respective lengths of semi-axes We also list the rounded correlation matrix, We see that the only pronounced correlations are between Ω m , Ω Λ and M 0 . This is also apparent from Table 1.
Code Availability. The code and data used in the analysis are available at: http://dx.doi.org/10.5281/zenodo.34487.