1 Introduction

In the past two decades, a ‘concordance \(\varLambda \)CDM model’ of cosmology has emerged, in which data are interpreted assuming that the Universe is isotropic and homogeneous, and General Relativity (GR) is the true theory of gravity. In this framework, spacetime is modelled by the maximally-symmetric Friedman–Lemaître–Robertson–Walker class of solutions to the Einstein field equations, obtained by requiring the stress-energy tensor to be exactly isotropic and homogeneous on all scales. In the best-fit \(\varLambda \)CDM model [1] the energy density of the Universe today consists of 4.9% baryonic matter, 26.8% dark (gravitating) matter and 68.3% ‘dark energy’ which has come to dominate the Universe at redshift \(z < 1\) and is causing its expansion rate to accelerate—in its most economical form Einstein’s Cosmological Constant \(\varLambda \).

The Cosmological Constant can also be interpreted as the energy density of fluctuations of the vacuum in the quantum field theories that describe the energy-momentum tensor. Thus the effective \(\varLambda \) appearing in the Friedmann-Lemaître equations is the sum of the arbitrary classical contribution (from geometry), and the quantum contribution (from matter) which cannot be rigorously calculated, even its sign. That these unrelated terms should conspire together so as to make \(\varLambda \) tiny enough to allow the Universe to have lasted for over 10 billion years—rather than gone into perpetual inflation or simply recollapsed at a very early age—is the as yet unsolved Cosmological Constant problem [2].Footnote 1

As is well-documented [5, 6], \(\varLambda \) has had a chequered history with claims for its presence which were subsequently found to be unjustified because of systematic uncertainties (in particular the evolution of ‘standard candles’, ‘standard rods’, galaxy number counts etc). Its recent reincarnation was initiated by indirect arguments [7] arising of a desire to get diverse cosmological datasets to be concordant with the FLRW model, rather than any compelling physical argument why \(\varLambda \) should have a value of \({{\mathcal {O}}}(H_0^2)\), where \(H_0 \equiv 100h\) km s\(^{-1}\)Mpc\(^{-1} \sim 10^{-42}\)GeV is the present expansion rate and the only dimensionful parameter in the model.Footnote 2 It was not until observations of 93 Type Ia supernova in the late 1990s claimed direct evidence for an accelerating Universe [9, 10], that \(\varLambda \) was taken seriously and \(\varLambda \)CDM came to be canonised as the ‘standard model’ of cosmology. It was therefore surprising when some years later a principled statistical analysis of a bigger sample of 740 uniformly analysed SNe Ia in the Joint Lightcurve Analysis (JLA) catalogue [11] demonstrated that the evidence for acceleration is rather marginal i.e. \(<3\sigma \) [12]. In response, it was argued [13] that doing the analysis differently, and in particular taking into account other cosmological data such as on BAO and CMB anisotropies analysed in the FLRW framework, firmly reinstates the acceleration of the universe. To assist the reader in assessing these discussions and to place the history in context, we briefly review supernova cosmology in Sect. 3.

The observed Universe is not quite isotropic. It is also manifestly inhomogeneous. The interpretation of data from the real Universe in terms of the FLRW model has proceeded under the weaker assumptions of statistical isotropy and homogeneity, the modern version of the Cosmological Principle. Since the CMB actually exhibits a large dipole anisotropy, this is assumed to be due to our motion with respect to the cosmic rest frame in which the Universe looks FLRW, and data can be analysed according to the Friedman–Lemaître equations. Hence cosmological data are corrected for this motion using a special relativistic boost. However as we show in Sect. 2, this assumption is no longer tenable. Several independent data sets now argue against the existence of such a frame in the real Universe. At low redshift (\(z < 0.1\)), all local matter appears to have a coherent ‘bulk flow’ aligned with the direction of the CMB dipole. At higher redshift (\(z > 1\)) the matter dipole is far larger than is expected from the kinematic interpretation of the CMB dipole.

The peculiar Raychaudhury equation [16] has been employed to show [17] that the acceleration of the expansion of space as inferred by an observer embedded within such a bulk flow should show a scale-dependent dipolar modulation. We have tested this hypothesis [18](hereafter: CMRS19) and find that the inferred acceleration of the expansion rate is indeed directional and its dipole component is 50 times larger than and statistically far more significant (at \(3.9 \sigma \)) than the isotropic component which is non-zero at only \(1.4 \sigma \), suggesting that the inferred acceleration cannot be due to \(\varLambda \). We also showed (Figure 2 of CMRS19 [18]) that SNe Ia data as disseminated publicly and thus employed by many authors in fitting for concordance cosmology, not only acknowledge the existence of this bulk flow, but ‘correct’ the observations towards an assumed cosmic rest frame in an unrealistic manner using unreliable data sets [19]. The methodology in CMRS19 [18] has been subjected to further criticism [20] (RH19) which we have rebutted seriatim [21].

General Relativity (GR), as expressed in the Einstein Field Equations and the Geodesic Equation, is a geometric description of gravity in which spacetime is modelled as a pseudo-Riemannian manifold warped locally by the presence of matter and energy. The real Universe contains structure over a vast range of scales. While any volume element in the FLRW model expands at the same rate as all others, as described by a single scale factor a, the expansion of the real Universe in GR is an average effect arising from the coarse-graining of small scales and is thus locally different everywhere [23]. This is properly described by the Raychaudhury equation—which reduces to the Friedmann-Lemaître equation as a limiting case [24].Footnote 3 Due to these cardinal differences between the real Universe and its simplified description in the FLRW model, even as the ‘precision cosmology’ programme has in the past two decades focused on measuring the \(\varLambda \)CDM cosmological parameters as precisely as possible, there is an ongoing theoretical debate as to whether this really makes physical sense [26, 27]. We emphasize that the ‘fitting problem’ in cosmology [28, 29] is at the heart of the problem and dedicate Sect. 4 to a discussion of how this relates to the corrections employed in supernova cosmology. We conclude in Sect. 5 that the evidence for the acceleration of the Hubble expansion rate from supernovae remains marginal.

2 The universe is anisotropic

The ‘Cosmological Principle’ began as an assumption [30] but is now said to be supported by the smoothness of the cosmic microwave background (CMB), which has temperature fluctuations of only about 1 part in 100,000 on small angular scales. These high multipoles in the CMB angular power spectrum are attributed to Gaussian density fluctuations created in the early universe with a nearly scale-invariant spectrum which have grown through gravitational instability to create the large-scale structure in the present Universe.Footnote 4 The observed dipole anisotropy of the CMB is however \(\sim 100\) times larger and is believed to be due to our motion with respect to the ‘CMB frame’ in which it would look isotropic [32, 33], so is called the kinematic dipole. This motion was attributed to the gravitational effect of the inhomogeneous distribution of matter on local scales, dubbed the “Great Attractor” [35]. In subsequent studies however this ‘bulk flow’ has been seen to be significantly faster and extend deeper than the expectation in the \(\varLambda \)CDM model [36,37,38,39]. In fact it extends up to and beyond the Shapley supercluster, i.e. out to \(> 300\) Mpc [14, 15], which is well beyond the scale of \(\sim 100\) Mpc on which the universe is supposedly sensibly homogeneous according to galaxy counts in large-scale surveys [40, 41]. A consistency check would be to measure the concommitant effects on higher multipoles of the CMB anisotropy [42], however, the significance of the detection of these is \(<3\sigma \) even in the precise measurements by the Planck satellite [43], which allows up to 40% of the observed dipole to be due to effects other than our motion [31].

Yet another probe is to examine the directional behaviour of the X-ray luminosity-temperature correlation of galaxy clusters which have been mapped out to several hundred Mpc. Examination of the ROSAT catalogue of several hundred clusters too reveals an anisotropy in roughly the same direction at high significance \(>5\sigma \) [53, 54].

An independent test would be to check if the reference frame of matter at still greater distances converges to that of the CMB, i.e. the dipole in the distribution of cosmologically distant sources, induced by our motion via special relativistic aberration and Doppler shifting effects, should align both in direction and in amplitude with that inferred from the CMB dipole [44]. However, this test when done with the NVSS catalogue of radio sources yielded an anomalously high result [45]. This was criticised as it was not an all-sky survey, also since the redshifts are not directly measured there is potential for contamination by a ‘clustering dipole’ due to a local inhomogeneity. However the anomaly persists even after such concerns are addressed [46,47,48,49] and is present in many such studies done with different radio catalogues [50, 51]. Recently the CatWISE catalogue of 1.36 million quasars (with median redshift of 1.2) was analysed in this context and the kinematic expectation for the quasar dipole rejected at \(4.9\sigma \)—the highest significance achieved to date [52]. Knowing the redshift distribution of the quasars, the expected clustering dipole according to the \(\varLambda \)CDM model [46] could be explicitly calculated and shown to be negligible.

Figure 1 shows the directions of the anisotropies reported—in the local bulk flow traced out to \(z \sim 0.1\), in X-ray clusters and SNe Ia out to \(z \sim 1\), and in radio sources and quasars at \(z > 1\).

Fig. 1
figure 1

Examples of directional anisotropy reported in studies of the local bulk flow [14, 15, 36, 38, 39], X-ray clusters [53, 54], SNe Ia [21], high redshift radio sources [49, 50] and quasars [52]. These are all close to the CMB dipole direction [43] which is also marked

3 Supernova cosmology

Type Ia Supernovae (SNe Ia) are believed to be thermonuclear explosions in low-mass stars, e.g. triggered when the mass of a white dwarf is driven by the accretion of material from a companion over the maximum that can be supported by electron degeneracy pressure. Since this happens at a critical mass, the Chandrasekhar limit of \(\sim 1.4 M_\odot \), all Type Ia supernovae are taken to have the same intrinsic luminosity, i.e. a ‘standard candle’. In practice the intrinsic magnitudes of nearby SNe Ia (to which distances are known via independent means) exhibit a rather large scatter. However by exploiting the observed linear correlation of the (colour-dependent) luminosity decline rate with the peak magnitude [55], this scatter can be considerably reduced. This makes SNe Ia ‘standardisable’ candles, i.e. the intrinsic magnitude can be inferred with relatively low scatter (0.1–0.2 mag) by measuring the lightcurves in different (colour) bands [56]. Further assuming that the intrinsic properties themselves do not evolve with redshift, observations of SNe Ia can be used to measure the cosmological evolution of the luminosity distance (i.e. of the scale factor) as a function of redshift.

In detail, however, the different empirical techniques for implementing the Phillips corrections [55], viz. the Multi Colour Lightcurve Shape (MLCS) strategy [10], the ‘stretch factor’ corrections [9] and the template fitting or \(\varDelta m_{15}\) method [57, 58], do not agree with each other—see Figure 4 of Ref. [56]. As the sample of SNe Ia has grown, the tension between the methods has in fact increased [59]. The MLCS strategy was to simultaneously infer the Phillips corrections and the cosmological parameters using Bayesian inference. However, a two-step process—the ‘Spectral Adaptive Lightcurve Template’ (SALT)—is now adopted, wherein the shape as well as the colour [60] parameters required for the Phillips corrections are first derived from the lightcurve data, and the cosmological parameters are then extracted in a separate step [61]. The current incarnation of this method is SALT2, employed in analysis of recent public SNe Ia data sets [11, 62], in which every SNe Ia is assigned three parameters, \(m_B^*\), \(x_1\) and c—respectively, the apparent magnitude at maximum (in the rest frame ‘B-band’), the lightcurve shape, and the lightcurve colour correction. This can be used to construct the distance modulus using the Tripp formula [60]:

$$\begin{aligned} \mu _\mathrm {SN} = m_B^* - M^0_B + \alpha x_1 - \beta c , \end{aligned}$$
(1)

where \(M^0_B\) is the absolute magnitude (degenerate with the absolute distance scale i.e. the value of \(H_0\)) while \(\alpha \) and \(\beta \) are parameters which are assumed to be constants for all SNe Ia. (Further parameters can be added, e.g. a ‘mass step correction’ according to the mass of the SNe Ia host galaxy, but this turns out to be irrelevant in the fitting exercise, whereas the stretch and colour corrections above are both important and uncorrelated with each other [12].) This is related to the luminosity distance \(d_\mathrm {L}\) as

$$\begin{aligned} \mu = 25 + 5 \mathrm {log}_{10} (d_\mathrm {L}/\mathrm {Mpc}), \end{aligned}$$
(2)

where \(d_\mathrm {L}\) is a function of the \(\varLambda \)CDM model parameters:

$$\begin{aligned} d_\mathrm {L}= & {} (1 + z) \frac{d_\mathrm {H}}{\sqrt{\varOmega _k}} \mathrm {sin}\left( \sqrt{\varOmega _k} \int _0^{z} \frac{H_0 \mathrm {d}z'}{H(z')}\right) , \mathrm {for \ } \varOmega _k > 0 \nonumber \\= & {} (1 + z) d_\mathrm {H} \int _0^{z} \frac{H_0 \mathrm {d}z'}{H(z')}, \mathrm {for \ } \varOmega _k = 0 \nonumber \\= & {} (1 + z) \frac{d_\mathrm {H}}{\sqrt{\varOmega _k}} \mathrm {sinh}\left( \sqrt{\varOmega _k} \int _0^{z} \frac{H_0 \mathrm {d}z'}{H(z')}\right) , \mathrm {for \ } \varOmega _k < 0 \nonumber \\ \mathrm {where:} d_\mathrm {H}= & {} c/H_0, \quad H_0 \equiv 100h~\mathrm {km}\,\mathrm {s}^{-1}\mathrm {Mpc}^{-1}, \nonumber \\ H= & {} H_0 \sqrt{\varOmega _\mathrm {m} (1 + z)^3 + \varOmega _k (1 + z)^2 + \varOmega _\varLambda } . \end{aligned}$$
(3)

Here H the Hubble parameter (\(H_0\) being its present value), \(d_\mathrm {H}\) is called the ‘Hubble distance’, and \(\varOmega _\mathrm {m} \equiv \rho _\mathrm {m}/(3H_0^2/8\pi G_\mathrm {N}), \varOmega _\varLambda \equiv \varLambda /3H_0^2, \varOmega _k \equiv -kc^2/H^2_0 a^2_0\) are the matter, cosmological constant and curvature densities in units of the critical density which are related in the \(\varLambda \)CDM model by the ‘cosmic sum rule’: \(1=\varOmega _\mathrm {m} + \varOmega _\varLambda + \varOmega _k\).

To extract cosmological parameters, this distance modulus needs to be compared to the theoretical prediction from a model, using an appropriate and principled statistical procedure.

3.1 Standardisable candles?

Using SNe Ia to measure cosmological parameters involves the implicit assumption that their intrinsic properties do not evolve with redshift. This requires that the SN luminosity after standardisation (\( m_B^* + \alpha x_1 -\beta c\)) should be invariant or equivalently that \(M^0_B\) should not evolve with look back time. This assumption has been challenged [63, 64], although this too is disputed [65]. In addition, it has been argued that the distributions of \(x_1\) and c are explicitly sample- and redshift-dependent [13, 20] even though this requires adding 12 new parameters to a 10-parameter model, so is not justified by any information criterion [21]. If that is indeed the case, it is unclear how SNe Ia can still be thought of as standardisable candles. This would require the other terms in Eq.(1) to evolve with redshift at just the right rate such that \(M^0_B\) remains invariant. If the absolute magnitude were to evolve with redshift this would of course trivially undermine the inference of accelerated expansion [66, 67].

Since the discussion on this issue [63, 65] hinges on how the intrinsic scatter (of unknown origin) of SNe Ia data in the Hubble diagram is treated statistically, a closer examination of the usual method is in order.

3.2 The statistical method

Historically, analyses in supernova cosmology [68,69,70,71,72] have relied on a ‘constrained \(\chi ^2\) statistic’ defined as follows:

$$\begin{aligned} \chi ^2 = \sum _i \frac{(\mu _i - \mu _i^\mathrm {obs})^2}{\sigma _{\mu i}^2} . \end{aligned}$$
(4)

Here \(\mu _i^\mathrm {obs}\) is the distance modulus for the \(i^{th}\) SNe Ia derived from observations by employing the Phillips corrections according to Eq.(1) and \(\mu _i\) is the corresponding model prediction based on its redshift (Eq.2). The variance \(\sigma _{\mu i}^2\) is the sum of different sources of uncertainty added in quadrature.

$$\begin{aligned} \sigma _{\mu i}^2 = (\sigma _{\mu i}^\mathrm {fit})^2 + (\sigma _{\mu i}^{z})^2 + (\sigma _{\mu ,i}^l)^2 + (\sigma _{\mu }^\mathrm {int})^2 , \end{aligned}$$
(5)

where \(\sigma _{\mu i}^\mathrm {fit}\) is the uncertainty from the lightcurve fitting process, \(\sigma _{\mu i}^{z}\) is the uncertainty of the SN redshift measurement, from spectroscopy as well the peculiar velocities of, and within, the host galaxy, and \(\sigma _{\mu ,i}^l\) is the uncertainty due to gravitational lensing. Since each SNe Ia observation corresponds to a datum that is still intrinsically scattered around the expectation from the fiducial cosmological model, an additional term \(\sigma _{\mu }^\mathrm {int}\) is introduced and estimated from the fit by requiring that \(\chi ^2/\mathrm {d.o.f.} = 1\), sometimes on a sample-by-sample basis. This is very ad hoc and statistically unprincipled [73, 74]. Since \(\chi ^2/\mathrm {d.o.f.}\) is the statistic that is usually used to judge goodness of fit, requiring it to be 1 (i.e. a perfect fit) as a condition to estimate yet another parameter (\(\sigma _{\mu }^\mathrm {int}\)) forsakes any ability to judge the goodness-of-fit of the assumed model, and no confidence interval can be provided on the \(\sigma _{\mu }^\mathrm {int}\) thus estimated. In essence a large effective number of trials associated with trying various values of \(\sigma _{\mu }^\mathrm {int}\) in the process of arriving at the final value have been ignored.

Table 1 Summary of approaches in major SNe Ia cosmology analyses. \(N_\mathrm {out}\) is the number of outliers rejected, while CRF is the choice made of Cosmic Rest Frame. The first analyses argued that gravitational lensing is a negligible effect, however, subsequently a z-dependent systematic uncertainty to account for it has been introduced. \(^a\) Explicitly noted as “3 \(\sigma \) outlier” rejections
Fig. 2
figure 2

Sky distribution of the SNe Ia used by the HZT [10] (left) and the SCP [9] (right). Supernovae at \(z < 0.1\) are indicated in red while those at \(z > 0.1\) are in green

Fig. 3
figure 3

Left: The profile likelihood for the 10-parameter FLRW analysis (3), using the JLA [11] colour c and stretch \(x_1\) corrections—see Table 2. Right: The same for the cosmographic analysis (6)—see Table 3

It is motivated by these considerations that the Maximum Likelihood Estimator method [12, 75] was introduced for supernova cosmology. Ref. [74] notes that if \(x_1\) and c evolve with redshift, the likelihood-based methods return biased values of the parameters (while the ‘constrained \(\chi ^2\)’ method continues to be robust), however, this conclusion is arrived at using Monte Carlo simulations which assume the \(\varLambda \)CDM model and is therefore a circular argument. It has been emphasised [76] that systematic uncertainties and selection biases in the data need to be corrected for in a model-independent manner, before fitting to a particular cosmological model.

An examination of previous SNe Ia analyses (see Table 1) reveals that the choice of the component of \(\sigma _{\mu i}^{z}\) from peculiar velocities as well as \(\sigma _{\mu i}^{l}\) the uncertainty due to lensing, has varied considerably. As is shown in figure B.1. of CMRS19, while \(\sigma _{\mu i}^\mathrm {int}\) is independent of redshift, the peculiar velocity component of \(\sigma _{\mu i}^{z}\) varies inversely with the redshift of the SNe Ia, while \(\sigma _{\mu i}^{l}\) is proportional to the redshift.

3.3 The accelerated expansion of the Universe

The original analyses by the Supernova Cosmology Project (SCP) [9] and High z Supernova Search Team (HZT) [10] employed, respectively, 60 and 50 SNe Ia, 17 of which were in common. Subsequent tests on larger datasets have shown that the evidence is strongly dependent on the assumption of isotropy [77], while supernova datasets themselves show evidence for a local bulk flow, i.e. anisotropy [14, 15]. CMRS19 [18] found that the inferred acceleration is also highest in this general direction when the analysis is done in the heliocentric frame. In Fig. 2 we show that most of the SNe Ia observed by SCP and HZT were in fact also in the same general direction.

Fig. 4
figure 4

The profile likelihood for the 22-parameter cosmographic analysis of RH19 [20], employing the sample- and redshift-dependent treatment of \(c_0\) and \(x_{1,0}\) [13]—see Table 4

Table 2 Results of a 10-parameter fit to the \(\varLambda \)CDM model (3), using the JLA [11] colour (c) and stretch (\(x_1\)) corrections
Table 3 Results of a 10-parameter fit in the cosmography framework (6), using the JLA [11] stretch (\(x_1\)) and colour (c) corrections

Nielsen et al. [12] showed using a Maximum likelihood estimator that even with the 740 SNe Ia in the much enlarged JLA catalogue [11], uniformly analysed with the SALT2 template, the evidence for acceleration is quite marginal (\(\sim 2.8\sigma \)). This analysis employed the JLA catalogue as it ships, i.e. with peculiar velocity (PV) corrections. Since such corrections introduce an arbitrary discontinuity within the data and modify the low redshift SNe Ia which serve to fix the lever arm of the Hubble diagram [18], the impact of these corrections on the evidence for cosmic acceleration is worth examining. (Note that PV corrections in SNe Ia analyses were first introduced only in 2011 [71].)

Over half of the 740 SNe Ia in the JLA catalogue are relatively local (\(z < 0.1\)) and the furthest is at \(z = 1.4\). Hence the luminosity distance \(d_\mathrm {L}\) can be quite accurately, to within 7% of Eq. (3), expanded as a Taylor series in terms of the Hubble parameter \(H_0\), the deceleration parameter \(q_0 \equiv -\ddot{a} a/{\dot{a}}^2\), and the jerk \(j_0 \equiv \dot{\ddot{a}}/a H^3\). This is ‘cosmography’, i.e. independent of assumptions about the content of the universe. (Specifically in the \(\varLambda \)CDM model, \(q_0 \equiv \varOmega _\mathrm {m}/2 - \varOmega _\varLambda \).) Modified to explicitly show the dependence on the measured heliocentric redshifts, this writes (RH19):

$$\begin{aligned}&d_\mathrm {L} (z, z_\mathrm {hel})\nonumber \\&\quad = \frac{cz}{H_0} \left[ 1 + \frac{1}{2}(1 - q_0)z - \frac{1}{6} (1-q_0 - 3q^2_0 + j_0 - \varOmega _k) z^2 \right] \nonumber \\&\quad \times \frac{1+z_\mathrm {hel}}{1+z} , \end{aligned}$$
(6)

where z can be the measured heliocentric redshift, boosted to the CMB frame, or boosted to the CMB frame with further PV ‘corrections’ applied.

We redo the MLE analysis [12] in four different ways and present the results in Figs. 3 and 4, and in Tables 2 and 3, respectively—both for the standard \(\varLambda \)CDM model (3) and the cosmographic Taylor expansion (6). For each case we also show the fit quality when \(q_0\) is held at zero (“No accn.”).

  1. 1.

    \(z_\mathrm {CMB}\) with PV corr.: This reproduces the earlier analysis of Ref. [12]. The CMB frame redshifts are used, with further corrections made for the peculiar velocities of the SNe Ia w.r.t. the CMB frame, and the peculiar velocity covariance matrix is included.

  2. 2.

    \(z_\mathrm {CMB}\) without PV corr.: Now CMB frame redshifts are used without correcting for the flow of the SNe Ia w.r.t. this frame and the peculiar velocity component of the covariance matrix is excluded. Note that transforming from heliocentric to CMB frame redshifts requires assuming that the CMB dipole is kinematic in origin.

  3. 3.

    \(z_\mathrm {hel}\): Heliocentric redshifts are used, no corrections are employed and the peculiar velocity component of the covariance matrix is excluded.

  4. 4.

    \(z_\mathrm {LG}\): Finally, redshifts in the Local Group frame are used and the peculiar velocity component of the covariance matrix is excluded. This is similar to what was done by the SCP [9].

Our results in Tables 2 and 3 illustrate that the peculiar velocity corrections serve to bias the data towards higher acceleration (more negative \(q_0\)). Using observables corrected to the Local Group, as employed by SCP [9], the change in \(2\mathrm {log}\mathcal {L}\) between the best-fit model and the one with zero acceleration is only 4.0, indicating that the preference for acceleration is \(<1.6\sigma \).

In Table 4 we adopt the redshift-dependent treatment of the stretch and colour corrections advocated by RH19 [13, 20] (\(x_{1,0}\) and \(c_0\) being their present values) and use redshifts transformed to the CMB frame. Nevertheless, the evidence for acceleration remains \(<3\sigma \) unless further ‘corrections’ are made for the peculiar velocities of the SNe Ia with respect to the CMB frame. As shown by CMRS19 [18], these corrections asssume that convergence to the CMB frame is achieved by \(\sim 150\) Mpc even though the data say otherwise [14, 15, 36, 38, 39].

Table 4 Selected results from the 22-parameter fit in the cosmography framework (6) using sample- and redshift-dependent stretch (\(x_{1}\)) and colour (c) corrections [13, 20]

It is unrealistic to make such corrections for peculiar velocities which leave the SNe Ia immediately outside the flow volume uncorrected [18], in particular such a procedure induces a directional bias (which in the JLA uncertainty budget is simply assigned an uncorrelated variance of \(c\sigma _z = 150\) km \(s^{-1}\)). An isotropic acceleration can be extracted from the data only by ‘correcting’ over half of all observed supernovae to the CMB frame—to which convergence has never been demonstrated.

In fact, the peculiar velocity corrections affect the lever arm of the Hubble diagram in a non-obvious manner. The subsequent Pantheon compilation [62] initially included peculiar velocity corrections far beyond the extent of any actual survey or even flow model.Footnote 5 When this issue was pointed out and the bug was fixed, both the magnitudes and heliocentric redshifts of the corresponding supernovae were noted to be discrepant [19].Footnote 6 We cannot analyse data from subsequent surveys such as the Carnegie Supernova Project [79] and the Dark Energy Survey [80] since neither dataset is publicly available in an usable form, unlike the JLA catalogue which provides full details of the lightcurve fitting parameters and all individual covariances.

4 Peculiar velocity corrections and the ‘fitting problem’

We can observe only one Universe from our unique vantage point. As pointed out in Ref. [29] any anisotropies we observe are direct measures of the non-FLRW nature of the real Universe. Fitting an idealised FLRW model to the real Universe necessarily involves a choice of ‘corresponding 2-spheres’ which then fixes the cosmic rest frame (CRF) in the real Universe with respect to which all non-Hubble velocities are defined.

The cosmological redshift (\(z_\mathrm {c}\)) in the cosmic rest frame is modified by peculiar velocities of the source and the observer as [29]:

$$\begin{aligned} 1+z = (1 + z_\mathrm {O})(1 + z_\mathrm {c})(1 + z_\mathrm {s}) \end{aligned}$$
(7)

where z is the measured redshift in some reference frame and \(z_\mathrm {O}\) (\(z_\mathrm {s}\)) is the redshift due to the peculiar velocity of the observer (source) with respect to this frame. (Until 2008, the approximate expression \(z \sim z_\mathrm {c} + z_\mathrm {O} + \cdots \) was used instead of the correct one above [34].) Ref. [29] emphasised that in addition to a best-fit FLRW model given available data, a statement of the goodness-of-fit is also necessary to determine the extent of the non-FLRW nature of the real Universe. The ‘constrained \(\chi ^2\)’ method is however by construction unsuited to judge goodness-of-fit, while the choice of CRF has also evolved (see Table 1). Thus the procedures adopted in supernova cosmology are at odds with the fundamental challenge of the ‘fitting problem’ in cosmology.

5 Conclusions

We have shown [18] that the acceleration of the Hubble expansion rate inferred from SNe Ia magnitudes and redshifts as measured (in the heliocentric frame) is described by a dipole anisotropy. To infer an isotropic component of the acceleration i.e. a monopole (such as can be attributed to \(\varLambda \)), it is necessary to boost to the (possibly mythical) CMB frame, and ‘correct’ the redshifts of the low-z supernovae further for their motion w.r.t. the CMB frame [20]. This is done using bulk flow models [36, 38] which already assume that the universe is well-described by the standard \(\varLambda \)CDM model. Moreover to boost the significance of the (monopole) acceleration above \(\sim 3\sigma \) the supernova lightcurve parameters need to be empirically modelled a posteriori as being both sample- and redshift-dependent [13]—violating both basic principles of unbiased hypothesis testing and the Bayesian Information Criterion [81]. Only by doing so can the significance of acceleration be raised to \(4.2 \sigma \) [13, 20].

It is often said that there is independent evidence for acceleration from observations of CMB anisotropies and large-scale structure. This is not the case, in particular, all low-redshift probes such as BAO, cosmic chronometers H(z) and the growth rate of structure \(\sigma _8(z)\), are also consistent with a non-accelerating universe [66]. The CMB does not have direct sensitivity to \(\varLambda \) and can only determine \(\varOmega _k\) and \(\varOmega _\mathrm {m}\).Footnote 7 It is only via the cosmic sum rule that \(\varOmega _\varLambda \) is then inferred ( \(=1-\varOmega _k-\varOmega _\mathrm {m}\)) to be non-zero. However, this sum rule is just the Friedmann-Lemaître equation which is valid only for the FLRW model. We wish to emphasise that it is the widespread (and at times convoluted) use of the FLRW model in cosmological data analysis that has created a strong bias towards the \(\varLambda \)CDM Universe.

Whether dark energy is just a manifestation of inhomogeneities is an open question [27]. This can impact on analysis of cosmological data in several different ways—since in the ‘cosmic web’ that we inhabit, our location (e.g. whether in a low or high-density region) is important, as is the direction in which we look (e.g. along filaments or voids), and whether we are Copernican or ‘tilted’ observers in a bulk flow. The spatial scales being probed by observations are important as the inhomogeneous Universe looks different when averaged on different scales. Ultimately we would like a scale-dependent description of the Universe (akin to the renormalisation group) but this is a long way from the (perturbed) FLRW model which has served so far as the workhorse of cosmology and led us to infer the existence of unphysical dark energy.