Cosmological inference from within the peculiar local universe

The existence of 'peculiar' velocities due to the formation of cosmic structure marks a point of discord between the real Universe and the usually assumed Friedmann-Lema\'{i}tre-Robertson-Walker metric which accomodates only the smooth Hubble expansion on large scales. In the standard $\Lambda$CDM model framework, Type Ia supernovae data are routinely"corrected"for the peculiar velocities of both the observer and the supernova host galaxies relative to the cosmic rest frame, in order to infer evidence for acceleration of the expansion rate from their Hubble diagram. However observations indicate a strong, coherent local bulk flow that continues outward without decaying out to a redshift $z \gtrsim 0.1$, contrary to the $\Lambda$CDM expectation. By querying the halo catalogue of the Dark Sky Hubble-volume N-body simulation, we find that an observer placed in an unusual environment like our local Universe should see correlations between supernovae in the JLA catalogue that are 2-8 times stronger than seen by a typical or Copernican observer. This accounts for our finding that peculiar velocity corrections have a large impact on the value of the Cosmological Constant inferred from supernova data. We also demonstrate that local Universe-like observers will infer a downward biased value of the clustering parameter $S_8$ from comparing the density and velocity fields. More realistic modelling of the peculiar local Universe is thus essential for correctly interpreting cosmological data.


Introduction
The flat ΛCDM 'standard model' of cosmology, which has a dominant fraction of its energy density Ω Λ ∼ 0.7 in the form of a cosmological constant and a fraction Ω m ∼ 0.3 in matter (of which ∼ 85% is cold dark matter and ∼ 15% baryons), is said to be a "good approximation to reality" [1].It is nevertheless experiencing a crisis due to the significant tension between the value of the present expansion rate H 0 (≡ 100h km s −1 Mpc, h ≃ 0.7) determined using the 'cosmic distance ladder' anchored in the local universe, and the value inferred from the Cosmic Microwave Background (CMB) in the ΛCDM model framework [2].Another tension is between the growth parameter S 8 = σ 8 (Ω m /0.3) 0.5 determined from observations of weak gravitational lensing, and from CMB data, where σ 8 is the variance of mass fluctuations in a top-hat sphere of radius 8h −1 Mpc [3].
In the ΛCDM model, data are interpreted using the Friedmann-Lemaítre equations, obtained from general relativity assuming the Cosmological Principle.In its modern form this assumes statistical isotropy and homogeneity in the distribution of matter and radiation in the 'cosmic rest frame' (CRF) in which the CMB dipole, assumed to be of kinematical origin, is presumed to vanish.We have recently shown, however, that the distribution of distant matter as traced by quasars and radio galaxies is not isotropic in the CRF [4,5], thus challenging this foundational assumption of the Friedmann-Lemaítre-Robertson-Walker metric.We have also shown that the acceleration of the Hubble expansion rate inferred from Type Ia supernovae (SNe Ia) is anisotropic in the heliocentric frame, and so cannot be interpreted as due to Λ [6].It is likely an artefact [7] because of our being 'tilted' observers embedded in a coherent bulk flow, which gives rise to the prominent dipole arXiv:2003.10420v5[astro-ph.CO] 8 May 2024 anisotropy in the CMB.Moreover, the local bulk flow extends out significantly further than is expected in the standard ΛCDM model, and no convergence is seen to the CRF out as far as ∼ 200 h −1 Mpc [8].In this paper, we focus on the impact of this anomaly on the S 8 tension and, more importantly, on the estimation of Ω Λ from SNe Ia data.
The relativistic viewpoint [9] is that such deviations from the Hubble flow should be thought of as variations in the expansion velocity field of the universe, rather than as 'peculiar velocities' with respect to a uniformly expanding space of the 'background cosmology', which is assumed a priori to be described by the Friedmann-Lemaítre equations.However, it is the latter approach that has become standard in cosmology, in particular for measuring the Hubble expansion rate.SNe Ia data [10], for example, are analysed in the framework of concordance cosmology (e.g., Ref. [11]) by making special relativistic corrections based on models of the local peculiar velocity field.For example, in order to obtain the 'cosmological' redshift from the measured value, corrections are made for (non-Hubble) velocities using Equation (1).These are typically a few hundreds of km s −1 , so would appear to only be relevant at low redshift z ≲ 0.1.However, they can affect the analysis at higher redshift too, since evidence for accelerated expansion is a dimming of the high-redshift SNe Ia in relation to the low-redshift SNe [12,13].Moreover, the local velocity field is quite noisy, reflecting our rather inhomogeneous neighbourhood, and it is not clear exactly where the separation should be made between the nearby and distant universe.Various empirical methods for accounting for peculiar velocities have thus been proposed [14]; the above problem has often been circumvented by simply excluding from cosmological fits all the SNe Ia at low redshifts.However this severely deprecates the sample statistics (since about a quarter of all known SNe Ia are in fact very local) and moreover it is somewhat arbitrary, e.g., cuts have been applied at both z = 0.01 and at z = 0.025 [15].Another option is to allow for an uncorrelated, and somewhat arbitrary, dispersion in the velocities of SNe Ia, e.g., Perlmutter et al. [12] took the redshift uncertainty due to peculiar velocities to be cσ z = 300 km s −1 , while Riess et al. [13] used cσ z = 200 km s −1 .
The alternative is to correct for the peculiar velocities [16]: for this purpose the IRAS PSCZ catalogue [17], the SMAC catalogue of clusters [14] and the 2M++ catalogue [18] have all been used to infer the peculiar velocity field from the underlying density field.However, these catalogues are rather limited and can be biased; moreover, linear Newtonian perturbation theory is used, hence the extracted velocities are model-dependent.Furthermore, these analyses assume convergence to the CRF at ≳ 100h −1 Mpc, as is expected in the framework of the standard ΛCDM model, even though this is contradicted by many independent observations [8,14,[19][20][21][22][23][24].It is in any case inappropriate to use standard ΛCDM to make such corrections, since the model is itself a subject of the test being carried out.
Specifically, in analyses of the SDSS-II/SNLS3 Joint Lightcurve Analysis (JLA) catalogue of 740 SNe Ia [25], and the subsequent Pantheon catalogue of 1048 SNe Ia (which includes 279 SNe Ia from Pan-STARRS1) [26], the low-redshift SNe Ia have been retained in the cosmological fits by thus "correcting" the individual redshifts and magnitudes of the SNe for the local 'bulk flow' inferred from density field surveys out to z ∼ 0.04 [14] and z ∼ 0.067 [18].However, as noted by [6], in both analyses SNe Ia immediately outside the survey volume of the peculiar velocity field were arbitrarily assumed to be at rest with respect to the cosmic rest frame, despite the fact that the same surveys detected a bulk flow extending beyond the survey volume of 372 ± 127 km s −1 and 159 ± 23 km s −1 , respectively.Moreover the JLA and Pantheon analyses adopted different values of 150 km s −1 and 250 km s −1 , respectively, for the dispersion cσ z of the bulk flow velocity.
Moreover the peculiar velocity corrections applied to both JLA and Pantheon contain significant errors and inconsistencies [6,27].Since the covariance matrices for peculiar velocity corrections have not been provided separately, the impact of these errors on cosmological parameter estimation is hard to quantify.It should also be of concern that the the residual bulk flows of the peculiar velocity surveys align approximately with the directions of maximum hemispherical asymmetry in the sky coverage of the two catalogues.
On the theoretical side, the correlated fluctuations of SNe Ia magnitudes due to peculiar velocities and the impact on cosmological parameter estimation of making such corrections have been extensively studied [16,[28][29][30].However, all these studies assumed that the peculiar velocity statistics are those expected around a typical (aka Copernican) observer in a ΛCDM universe.Such an observer should in fact not observe a bulk flow exceeding ∼200 km s −1 beyond 100 h −1 Mpc (z ≃ 0.033), independently of the form of the matter power spectrum [31], so clearly this assumption is in tension with reality.In this work, we use a cosmological N-body simulations to examine the validity of the peculiar velocity covariances proposed by Hui & Greene [29] in the light of additional information about our local universe.Essentially, we account for the cosmic variance due to observer location, which in our specific case is especially rare for a ΛCDM universe-as we quantify below.
The bulk flow observations suggest that we are not typical observers in a ΛCDM universe [32][33][34].We discuss here the correlated fluctuations of SNe Ia magnitudes and redshifts due to the peculiar velocities and bulk flows in and around 'Local Universe (LU)-like' environments in the z = 0 halo catalogue of the DarkSky ΛCDM simulations [35].We find that previous theoretical predictions [29] for randomly selected typical observers have underestimated the actual covariances for observers like ourselves by a factor of 2-8.
Next we discuss (Section 2) the peculiar velocity corrections employed in JLA [25] and show that these are both arbitrary and incomplete (Section 2.1).We compare the magnitude of the velocities used for the corrections in JLA against those obtained from the Cosmicflows-3 (CF-3) compilation [36] and demonstrate that the JLA values are underestimated by ∼48% on average.We also review the various relevant sources of uncertainties and dispersions that go into the JLA cosmological fits.We then explore (Section 2.2) various methods to fit for the extent of the bulk flow in the LU and present our likelihood analysis (Section 2.3).This is followed by a discussion of related work (Section 4).An Appendix presents the standard methodology of cosmology with SNe Ia (Section A), and the JLA catalogue (Appendix B).
We find that for any consistent treatment of the peculiar velocities (including ignoring them altogether), the JLA dataset favours Ω Λ ≲ 0.45 and is consistent with a non-accelerating universe at ∼2σ.Larger values of Ω Λ which have been found in other analyses of the JLA catalogue [37,38] are in fact due to the incomplete peculiar velocity "corrections" applied.We demonstrate that, consistent with the recent finding from the Cosmicflows-4 survey [8], the JLA data favour a fast (>250 km s −1 ) bulk flow extending out to >200h −1 Mpc, which is quite unexpected in the standard ΛCDM model.
Subsequently (Section 3), we examine the S 8 parameter inferred by randomly selected Copernican observers as well as constrained 'Local-Universe like' observers by comparing the peculiar velocities around them with those expected from the density field.While the variance in S 8 as inferred by Copernican observers is already larger than the disagreement between a recent measurement [39] and the ΛCDM fiducial value from Planck, LU-like observers see both an additional downward bias on S 8 and a larger cosmic variance.

Selecting Local Universe Like Environments
Making the usual assumption that the CMB dipole is due to our motion with respect to the CRF (also known as the 'CMB frame') in which the universe looks isotropic, so the luminosity distance d L is related to the redshift z as in Equation (A2), the redshift of a supernova in the heliocentric frame z hel (obtained by correcting the actually measured redshift for the Earth's motion around the Sun) is related to its redshift z in the CMB frame (sometimes labelled z CMB ) as [28]: where z ⊙ is the redshift induced by our motion with respect to the CMB and z SN is the redshift due to the peculiar motion of supernova host galaxy in the CMB frame.The luminosity distance is similarly corrected as: to obtain d L as a function of z Equation (A2) for the standard ΛCDM model (see Appendix A).The covariance of SNe Ia magnitudes due to peculiar velocities is then given by [28][29][30]: where Here D i is the linear structure growth factor at the redshift of the ith SNe, j ′ l is the derivative of the lth spherical Bessel function and P l is the Legendre polynomial of order l.Note that according to this expression, the covariance in magnitudes between two SNe depends only on their relative angular separation (which comes in through P l ) and is independent of their absolute directions.
N-body simulations can be used to estimate ξ ij Equation (3) for different assumed observers.Figure 1 compares S ij evaluated using the ΛCDM expectation for ξ ij with that read off the z = 0 snapshot halo catalogue of Dark Sky , a Hubble volume, trillionparticle simulation [35], for two very different classes of observers.For the 'Copernican observer' in Figure 1 (left), the halo containing the observer and its orientation are selected at random-such an observer sees the universe as isotropic and homogeneous.However, for the constrained 'LU-like' observer' (right), only halos satisfying the following criteria are considered (the first three being the same as in Ref. [32]): (i) The observer halo has a Milky Way (MW)-like mass, in the range 2.2 × 10 11 < M 200 < 1.4 × 10 12 M ⊙ [40] for the halo mass contained within 200 kpc.(ii) The bulk velocity in a sphere of radius R = 3.125h −1 Mpc centred on the observer is Mpc from the observer.(iv) The angle between the bulk flow of (ii) and the direction to the Virgo-like halo of (iii) is (44.5 ± 5) • (v) The bulk velocity in a sphere of R = 180h −1 Mpc centered on the observer is 260 ± 100 km s −1 [19].(vi) The angle between the bulk flow of (v) and the direction to the Virgo-like halo of (iii) is (69.9 ± 7.5) • .(vii) The angle between the bulk flows of (ii) and (v) is (35.6 ± 7.5) • .
Note that such observers are very rare-we had to examine 8,721,498 halos in order to find 1000 that satisfied these criteria.Hence the above probability of 1000/8721498 = 0.0115% is an upper limit on the probability of finding our Local Universe in ΛCDM cosmology.After an observer satisfying the above criteria is found, the entire system is rotated so that the direction of the bulk flow of criterion (ii) and the direction to the Virgo-like halo of criterion (iii) correspond to the real observed directions.The criterion on the bulk flow direction is exact, while the criterion on the direction to the Virgolike halo is imposed only on the azimuthal angle in a coordinate system in which the z-axis points towards the bulk flow direction.Criterion (iv) then suffices to orient the system.Note that the angular tolerances in (iv), (vi), and (vii) are set to be less stringent than actual observational constraints in order to limit the required computation time.Subsequently, the halos closest to each JLA supernova with z < 0.1 (of which there are 152) are identified in the DarkSky simulation as a proxy for the supernova host galaxy, and their velocities are queried.From these velocities, ξ ij can be calculated for each pair in the set of 152 supernovae.For the typical observer of Figure 1 (left), none of the steps regarding directional orientation discussed above are considered and observers are simply picked at random.As seen in Figure 1 (right), a realistic LU-like observer sees on average correlations between the supernovae of a JLA-like catalogue that are 2-8 times stronger than does a typical observer.This illustrates that the theoretical covariances of Hui & Greene [29], as given in Equation ( 4), is valid only for idealised observers who see neither a local bulk flow nor a preferred orientation in the sky.

Peculiar Velocity Corrections in JLA
It had been noted in Ref. [6] that the peculiar velocity 'corrections' applied to the SNe Ia redshifts and magnitudes in the JLA catalogue (see Appendix B) are neither consistent nor complete.SNe Ia immediately beyond z ∼ 0.06 were taken to be stationary with respect to the CMB and assumed to only have an uncorrelated velocity dispersion cσ z = 150 km s −1 in the cosmological fits, even though observations of clusters indicate a bulk velocity of 372 ± 127 km s −1 due to sources beyond 200h −1 Mpc [14].Unlike the intrinsic dispersion σ M 0 which is assumed to be redshift-independent, the dispersion in the magnitudes as a result of the velocity dispersion is 5σ z /(zln10) i.e., the magnitudes of lower redshift supernovae are selectively more dispersed.As seen in Figure 2, the typical bulk flow in a ΛCDM universe (see, e.g., Ref. [41]) continues to much larger distances, with the velocity decreasing gradually.In some environments, the bulk velocity may even increase beyond a certain scale (as seen in the CosmicFlow-4 survey [8]), although the overall trend has to be a decreasing one if the Universe is to become homogeneous when averaged on large scales.
Cosmicflows-3 (CF-3) [36] is a compilation of the peculiar velocities of 17,669 nearby galaxies, using various independent distance estimators such as the Tully-Fisher relationship.In Figure 3, we compare the velocities that have been used to correct the JLA redshifts with those from CF-3.The galaxy in the CF-3 dataset corresponding to a JLA supernova is identified by cross-matching with a tolerance of 0.01 • , using the tool k3match.Out of 119 JLA SNe Ia at z CMB < 0.06, 112 have CF-3 counterparts within 0.01 • .It is seen from the regression line [42] in Figure 3 that peculiar velocities have been systematically underestimated by 48% in the corrections applied in the JLA analysis [25], compared to the actual measurements by CF-3.shows the best-fit orthogonal distance regression which has a slope of 1.61, i.e., the JLA velocities have been underestimated on average by 48%.(Note that the outlier (SN1992bh) has a peculiar velocity of ∼ 1000 km s −1 according to JLA, but zero according to CF-3.)

Fitting for a Bulk Flow
We consider two illustrative profiles for the bulk flow velocity.A linear ∼1/r fall-off is expected if we are, e.g., being pulled by a 'Great Attractor': where Q ′ is a (dimensionless) scale parameter.We ensure that the velocity never goes negative by setting it to zero above d L = P/Q ′ (see Figure 4).Our unusual bulk flow may however be due to a different physical cause, so we also consider an an exponentially falling form ⟨v⟩ where Q is the scale of the flow.The free parameters P and Q or Q ′ in our modelling of the bulk flow can be determined along with the 10 other parameters used to fit SNe Ia data.The red and green lines show, respectively, the linear and exponential fits to the bulk flow using JLA data -rows 9 and 10 of Table 1 (the dashed lines indicate the extrapolated fits).Note that the fits are not done to the other selected measurements shown.
The effect of a bulk flow term is to modify the distance modulus of Equation (A1) such that: for the exponentially falling bulk flow, and likewise for the linearly falling bulk flow with the expression from Equation (5) now substituted in the square brackets.

The Likelihood Analysis
We can now rewrite z SN and z in Equation (1) as functions of z hel , P and Q for the exponential (Equation ( 6)) or Q ′ for the linear (Equation ( 5)) bulk flow models, respectively.A MLE [37] is then used with the two additional parameters P and Q (or Q ′ ) for the bulk flow, in addition to the usual light curve fitting parameters in the SALT2 template [25]: α, x 1,0 , σ x 1,0 , β, c 0 , σ c 0 , M 0 , σ M 0 .Along with the two ΛCDM model parameters Ω m and Ω Λ we then have 12 parameters in total.As shown in Table 1, the following fits are performed (including an additional dispersion of cσ z = 150 km s −1 as recommended by Ref. [25]): 1.
The same 10-parameter fit as in Ref. [37], using only the z CMB values provided by JLA.

2.
The 10-parameter fit using Ref. [28]-and the JLA provided z hel and z CMB values.

3.
The 10-parameter fit as in 1, using only z hel values (JLA provided) as was done in all SNe Ia analyses until 2011.4.
The 10-parameter fit as in 1, using JLA provided z hel values, after subtracting out bias corrections to m * B .

5.
Exponentially falling bulk flow: 12-parameter fit (including the P and Q parameters of Equation ( 6), using only JLA provided z hel values.No peculiar velocity corrections are applied.6.
Linearly falling bulk flow: 12-parameter fit (including the P and Q parameters of Equation ( 5) using only JLA provided z hel values.No peculiar velocity corrections are applied.7.
JLA-corrected redshifts + Exponential bulk flow: 12-parameter fit: SNe Ia with peculiar velocity corrections applied by JLA, are treated as in (ii) above, while an exponentially falling bulk flow is fitted to the remaining SNe.

8.
JLA-corrected redshifts + Linear bulk flow: As in 7, but with the linear parametrisation of the bulk flow.9.
CF-3 data and the exponential bulk flow fit: 12-parameter fit using Equation ( 2) with the -derived values of z hel and z CMB (see Section 2.1) used for the low z SNe Ia to which the velocity correction can be applied.For the remaining objects, we use the JLA z hel values, and an Exponential bulk flow is fitted using Equation ( 6) as described above.10.CF-3 data and the linear bulk flow fit: 12-parameter fit using Equation ( 5) with the -derived values of z hel and z CMB (see Section 2.1) used for the low z SNe Ia to which the velocity correction can be applied.For the remaining objects, we use the JLA z hel values, and a linear bulk flow is fitted using Equation ( 5).
In all these fits, the direction of the bulk flow is fixed to be in the CMB dipole direction as most previous analyses have shown large dipoles at intermediate redshifts converging to this direction [19,22,24,34,43].In Table 1, we also show the fit results after imposing the additional constraint of 'No acceleration' for a ΛCDM universe i.e.: q 0 ≡ Ω Λ /2 − Ω m = 0.For the last two fits, we also show the effect of imposing the constraint of zero curvature ('Flat'): i.e., Ω Λ + Ω m = 1.
The bulk flow fit is ⟨v⟩ = 478e −d L /458 Mpc km s −1 for the exponential decay form ( 6), and ⟨v⟩ = [246 − 0.175(d L /Mpc)]km s −1 for the linearly falling form (5). Including the bulk flow always improves the quality of the fit as can be seen from the smaller values of −2 log L max .This justifies adding the two parameters characterising it.In all the above fits apart from the 'No acceleration' ones, the best-fit bulk flow extends beyond 200h −1 Mpc at 250 km s −1 .Figure 4 shows our results along with selected recent observations.
Using the CF-3 data and the linear bulk flow fit, as well as other fits of similar quality, the difference in the goodness of fit of the best model (with the lowest value of −2 log L max ) with respect to the corresponding 'No acceleration' fit is now significantly smaller compared to previous studies.Figure 5 demonstrates the degeneracy between the derived value of Ω Λ and the local bulk flow, illustrating that the latter is an essential nuisance parameter to be added to cosmological fits when analysing SNe Ia.Allowing for the bulk flow in the fit demonstrates that the evidence for acceleration using SNe Ia data alone is even weaker than was found previously [37].Table 1.Best-fit parameters and results for the fits described in Section 2.3 using the Maximum Likelihood Estimator [53].Including the bulk flow improves the quality of the fit and decreases the significance of accelerated expansion, as seen from the decrease of −2 log L max ."No acceleration" corresponds to As in Ref. [ The results in Table 1 may be summarised as follows: • Of all the fits, the only ones favouring Ω Λ > 0.5 are just those that include the incorrect and incomplete peculiar velocity 'corrections' of JLA [25].• Fit 4, which has no peculiar velocity corrections at all, as in the cosmic acceleration discovery papers [12] and [13], prefers Ω Λ = 0.396 with <2σ evidence for acceleration.

•
While previous work has suggested that bulk flows should not bias Ω Λ , it in fact drops by ∼30% if we undo the peculiar velocity 'corrections' of JLA and instead use the kinematic data from .This illustrates the huge impact of considering a realistic LU-like observer such as ourselves, rather than the randomly located observer assumed in all previous analyses [16,[28][29][30].In particular this contradicts what is stated in Table 11 of Ref. [25].
The discovery papers [12,13] assumed the uncertainty due to peculiar velocities to be cσ z = 300 km s −1 and 200 km s −1 respectively, but neither made SN-by-SN corrections.The JLA [25] and Pantheon [26] analyses employ incorrect peculiar velocity 'corrections', and adopt arbitrary redshift uncertainties of cσ z = 150 km s −1 and 250 km s −1 , respectively.

Extracting S 8
The peculiar velocity field v(r) is defined as where f (Ω m ) is the logarithmic growth rate of fluctuations (≃Ω 0.55 m in the ΛCDM model) and δ(r) is the matter density contrast field: The density in the above equation is that of gravitating matter, which is dominated by unobservable dark matter in the standard paradigm.It is usually assumed that observed luminous objects trace out the underlying matter density contrast with only linear bias where b g and δ g are the bias and the density fluctuation field of the tracers, respectively.The predicted peculiar velocity field can now be compared to the 'observed' (in the case of the present study, simulated) to estimate the term β t for the tracer t [39,44].Thus, β t can be obtained by fitting a large number of measured tracer velocities against the predictions from the density field.For straightforward comparison with results from weak lensing, we convert this, as in Ref. [39], to S 8 ≡ σ 8 (Ω m /0.3) 0.5 , using the input value of Ω m = 0.2952 used for the Dark Sky simulations.
In practice, the observational tracers used to map the velocity field need not be the same as are used to map the density field; in fact they usually are not [39].However, since our aim is to study the effect of the local environment on velocity-density correlations, we use the same tracers for both, viz. the halos (we use halos rather than particles since the z = 0 halo catalogue of Dark Sky is computationally more tractable than the z = 0 particle snapshot).In order to keep our study as similar as possible to Ref. [44] we use a Gaussian kernel with a smoothing length of 5 h −1 Mpc to smooth out the density fluctuation field.
It is evident from Figure 6 that there is a downward bias in measurements of S 8 by Local Universe-like observers such as ourselves and the cosmic variance is also higher.This may well account for the tension between the value derived from Planck data, and that obtained by comparing the reconstructed velocity field from the 2M++ galaxy redshift compilation to supernovae, Fundamental Plane and Tully-Fisher distances [18,39,45].

Discussion
Other authors [30,46] have come to different conclusions regarding the impact of the bulk flow on cosmological parameter determination, so we comment on why this is so.In general relativity, space-time evolves according to the Einstein field equations, but for simplicity and tractability, structure formation is studied by linearising these around the maximally symmetric FLRW solution.By further making restricted 'gauge' choices [47], cosmological reality is simulated using Newtonian N-body simulations wherein there is a background space that expands homogeneously and isotropically, governed by just one scale factor, while density perturbations evolve around this background according to Newtonian gravity.However, solutions to the linearised field equations can only be linearisations of the solutions to the fully non-linear equations [48,49]; hence, this approach hides known physical phenomena.Indeed, a relativistic universe exhibits locally inhomogeneous expansion beyond that evident in linear perturbation theory around a maximally symmetric background [50].It was therefore proposed to consider the expansion of the universe at late times as an average effect, arising out of the coarse-graining of physics at smaller scales [51].Peculiar velocities are thus at the interface between the real universe and its idealisation modelled by first-order perturbation of the fluid equations or in a Newtonian N-body simulation.Note that in the real universe, peculiar velocities are differences in the expansion rate of the universe at different space-time points, while in an N-body simulation they arise by construction from Newtonian gravity acting on top of a hypothetical uniformly expanding space.Refs.[29,30] use the restricted longitudinal or conformal-Newtonian 'gauge' [47] to derive the covariance (Equation 3) for a typical observer.But, as Figure 1 shows, this would be relevant for cosmology only if each SNe Ia were being observed from a different, randomly sampled, host galaxy.In practice we observe the real universe from only one unique vantage point.Nevertheless, Ref. [30] adds this covariance as a "guaranteed theoretical signal" to the uncertainty budget of the JLA data, thus weakening the preference for a bulk flow to ≲2σ (as seen in their Figure 4).
Ref. [46] claims further that any bias in the inference of dark energy parameters due to the effect of peculiar velocities can be determined a priori via simulations.This misses the point, however, that ΛCDM is a model and N-body simulations contain only as much physics as have been coded into them, i.e., neither capture the real universe.In fact Ref. [46] acknowledges that ignoring the velocity covariance altogether would lead to larger effects due to peculiar velocities-just as we have established here.

Conclusions
To summarise, we are not typical (Copernican) observers-we are embedded in a fast and deep bulk flow and this has significant impact on the covariances used in supernova cosmology. 1We have studied the effect on SNe Ia redshifts in the JLA catalogue of the peculiar velocities of their host galaxies.Using direct measurements of these from Cosmicflows-3, we find that the effect of peculiar velocities for low redshift SNe Ia has been underestimated by 48%.We show that the usual procedure of adding a constant velocity dispersion of a few hundred km s −1 to account for peculiar velocities at high redshift does not take into account the correlated flow of the galaxies.By analysing the DarkSky simulation [35], we demonstrate that 'Local universe-like' observers like ourselves see a 2-8 times stronger correlation between the SNe Ia than a randomly located observer does.The JLA analysis [25] corrected the data assuming the CMB dipole to be entirely kinematic in origin and that convergence to the CMB rest frame occurs abruptly at redshift z ∼ 0.06.Since neither assumption is fully supported by observations, we have adopted a general model of the bulk flow which introduces two additional parameters in the analysis.We do not adopt the ΛCDM model a priori, nor do we make assumptions about the origin of the CMB dipole.This provides an independent estimate of the bulk flow and we find that it persists out to distances beyond 200 h −1 Mpc, with a speed of ∼250 km s −1 .Our maximum likelihood analysis then shows that the accelerated expansion of the universe cannot be inferred as a statistically significant result from the SNe Ia data alone.We use the publicly available SALT2 light curve fits carried out by the JLA collaboration [25], but rather than use their 'constrained' χ 2 statistic which is unprincipled, we employ the Maximum Likelihood Estimator (MLE) due to Nielsen et al. [37].Note that the distribution of 'pulls' using this MLE is gaussian, while the best fit and its uncertainties explicitly satisfy Wilk's theorem [37].Our approach is frequentist but equivalent to the 'Bayesian Hierarchical Model' [38,54,55].It has been used in independent analyses of SNe Ia [56,57].
Ref. [38] advocated that the shape and colour parameters be allowed to depend on both the SNe Ia sample and the redshift; however, this introduces 12 additional parameters (to the 10 used above) and thus violates the Bayesian information criterion [6,58].Ref. [59] noted that if x 1 and c evolve with redshift, the likelihood-based methods return biased values of the parameters (while the 'constrained χ 2 ' method continues to be robust); however, this conclusion is arrived at using Monte Carlo simulations which assume the ΛCDM model and is therefore a circular argument.It has been emphasised by [56] 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.For further discussion of these and other related issues, see Ref. [60].

Notes 1
There are other corrections too such as for gravitational lensing, which become more important than the effect of peculiar velocities at redshift z > 0.15-see Figure B.1 of Ref. [6].

Figure 1 .
Figure 1.The theoretically expected covariance S ij (Equation (3)) plotted against the value found in N-body simulations-in regions around typical observers (left) and constrained 'Local Universe-like' observers (right).Each point is an average over 1000 observers.

Figure 2 .
Figure 2. The bulk flow velocity profiles around 10 random 'Local Universe-like' observers satisfying the criteria in Section 2. Note that the velocity profile around an individual observer need not decrease monotonically, even though the ensemble average in the ΛCDM model (dark blue curve) does so.Indeed, the Cosmicflows-4 data [8] indicates a rising velocity with depth.The shaded blue region is the ±1σ band around the mean value.

Figure 3 .
Figure 3.The line-of-sight velocity of SNe Ia inferred from their z hel and z CMB values quoted by JLA, plotted versus the line-of-sight component of the velocity of the group the object belongs to in the dataset (⟨V CMB ⟩ − gp).The blue horizontal bars are the diagonal errors in the JLA cosmology fit (statistical plus systematic), while the blue vertical bars indicate the random error of 250 km s −1 in the measurement.The green dashed line indicates when the two are equal, while the red dashed lineshows the best-fit orthogonal distance regression which has a slope of 1.61, i.e., the JLA velocities have been underestimated on average by 48%.(Note that the outlier (SN1992bh) has a peculiar velocity of ∼ 1000 km s −1 according to JLA, but zero according to CF-3.)

Figure 4 .
Figure 4.The profile of bulk flow expected in a ΛCDM universe using a top-hat window function is shown as the solid blue line, while the shaded region shows the ±1σ range.Several discrepant measurements (with ±1σ uncertainties) using various surveys are shown for comparison[8,[19][20][21]24].The red and green lines show, respectively, the linear and exponential fits to the bulk flow using JLA data -rows 9 and 10 of Table1(the dashed lines indicate the extrapolated fits).Note that the fits are not done to the other selected measurements shown.

Figure 5 .
Figure 5. 1, 2 and 3 σ contours corresponding to the fit in row 9 of Table 1 wherein peculiar velocities from Cosmicflows-3 are used for SN-by-SN corrections and the flow is allowed to continue beyond the survey volume with an exponential fall-off (Equation 6).The velocity of the bulk flow at a top-hat smoothing scale of radius 200 Mpc is shown in the right histogram of the posterior.The top histogram shows the extracted value of Ω Λ , which is seen to be degenerate with this bulk flow velocity.

Figure 6 .
Figure 6.Distribution of the S 8 parameter extracted by comparing the 'observed' (in this case, simulated) peculiar velocity with the prediction from the density contrast field, for 1000 observers selected at random (Copernican) as well as selected to be in Local Universe-like (LU-like) environments.The median values are shown in black.For comparison, the measurements [18,39,45] made using the same method are shown, as well as the value derived from the CMB measurements by Planck.

Figure A2 .
Figure A2.The redshift distribution of the 4 samples that make up the SDSS-II/SNLS3 Joint Lightcurve Analysis catalogue.