Effects of nonlinear inhomogeneity on the cosmic expansion with numerical relativity

Effects of nonlinear inhomogeneity on the cosmic expansion with numerical relativity


We construct a three-dimensional, fully relativistic numerical model of a universe filled with an inhomogeneous pressureless fluid, starting from initial data that represent a perturbation of the Einstein-de Sitter model. We then measure the departure of the average expansion rate with respect to this homogeneous and isotropic reference model, comparing local quantities to the predictions of linear perturbation theory. We find that collapsing perturbations reach the turnaround point much earlier than expected from the reference spherical top-hat collapse model and that the local deviation of the expansion rate from the homogeneous one can be as high as at an underdensity, for an initial density contrast of . We then study, for the first time, the exact behavior of the backreaction term . We find that, for small values of the initial perturbations, this term exhibits a scaling, and that it is negative with a linearly growing absolute value for larger perturbation amplitudes, thereby contributing to an overall deceleration of the expansion. Its magnitude, on the other hand, remains very small even for relatively large perturbations.

04.25.dg, 04.20.Ex, 98.80.Jk

Cosmology as a physical theory of the Universe was born soon after the formulation of general relativity one hundred years ago Ellis (2015), yet the extent to which relativistic nonlinearity may affect structure formation remains largely unexplored. With the increasing volume of cosmological data and their precision, more sophisticated modelling is required, and thus it is becoming timely to quantify these relativistic effects. The current theoretical framework for cosmology is based on three main ingredients: a homogeneous and isotropic Friedmann-Lemaître-Robertson-Walker (FLRW) background, relativistic perturbation theory to describe fluctuations in the early universe and at very large scale, and Newtonian methods, notably N-body simulations, to study the evolution of fluctuations into the nonlinear regime of structure formation. Reconciling this framework with the observations requires the existence of dark components, cold dark matter (CDM) and a cosmological constant or some other form of dark energy. The resulting standard cosmological model, CDM, satisfies a vast class of observational constraints, in particular the high precision measurements of the cosmic microwave background anisotropies Ade et al. (2015). However, the existence and nature of these dark constituents are one of the most debated topics not only in modern cosmology, but also in theoretical physics. One aspect that has been the subject of intense debate is the question whether nonlinear relativistic “backreaction” effects due to formation of structures may play an important role in the average cosmic expansion Buchert (2000); Räsänen (2011); Buchert et al. (2015); Green and Wald (2015); Kolb et al. (2005).

Quantifying the systematic errors involved in the different modelling approximations, such as the use of Newtonian gravity for structure formation, is a crucial undertaking if one wishes to interpret correctly the data which will be produced by the upcoming precision surveys Amendola (2013); Maartens et al. (). Whilst some approaches have been introduced to estimate the role of relativistic corrections in -body simulations Bruni et al. (2014a); Thomas et al. (2015); Milillo et al. (2015); Adamek et al. (2013, 2014, 2016), the only viable avenue to an exact computation of the systematic errors resulting from the omission of these effects is the direct numerical integration of Einstein’s equation in the corresponding scenarios. Integrating the equations of general relativity, possibly coupled to stress-energy sources, is the field of numerical relativity, a framework strongly motivated by gravitational-wave–source modelling, but which has, over the years, developed in a number of parallel areas such as cosmology, mathematical relativity, and modified gravity Cardoso et al. (2015). Some of this work has already been aimed at studying inhomogeneous cosmologies Bentivegna and Korzynski (2012); Yoo et al. (2012, 2013); Bentivegna and Korzynski (2013); Yoo and Okawa (2014). While these numerical-relativity studies do not yet aspire to the level of realism achieved by -body simulations Bruni et al. (2014b); Fidler et al. (2015), they are useful testbeds to quantify the relativistic effects of nonlinear inhomogeneity on the cosmic expansion.

In this Letter, we integrate Einstein’s equation coupled to an inhomogeneous irrotational pressureless fluid (dust) with a three-dimensional density profile and no continuous symmetries. We choose initial data corresponding to a perturbed Einstein-de Sitter (EdS) model, i.e. a flat FLRW model with dust, with the aim of measuring, with no approximations, the departures of the fully nonlinear numerical solution from the idealised FLRW background and its perturbations. On the numerically-generated spacetimes, we measure a number of local and average properties of cosmological interest, such as the growth of overdensities and the formation of voids, the inhomogeneous and average expansion rate, and the backreaction term defined in the averaging framework Buchert (2000). The main results of this study are that (i) those perturbations that are large enough to collapse stop partaking in the cosmic expansion (i.e. reach the “turnaround” point) much earlier than expected from a spherical top-hat collapse model with the same initial density contrast; (ii) locally, the effects of nonlinear inhomogeneities can be substantial, leading to a departure from the average expansion rate of over at the underdensities; and (iii) the average expansion rate is hardly affected by the inhomogeneities, with a backreaction term which is never larger than .

Method. We integrate Einstein’s equation and the fluid conservation equation using a variant of the Baumgarte-Shapiro-Shibata-Nakamura formulation  Nakamura et al. (1987); Shibata and Nakamura (1995); Baumgarte and Shapiro (1998), along with the Wilson formulation for the hydrodynamical system Font (2008) and the conformal transverse traceless formulation for the Einstein constraints Lichnerowicz (1944); York (1971), an approach already used in cosmological settings Anninos and McKinney (1999); Giblin et al. (2016); Mertens et al. (2016). We choose to represent the spacetime in the synchronous-comoving gauge Landau and Lifshitz (1975), popular in cosmological perturbation theory Bruni et al. (2014b), which corresponds to the Lagrangian coordinates of the observers at rest with the matter.

To integrate this system, we use the Einstein Toolkit Loffler et al. (2012); ?, a free, open-source community infrastructure for numerical relativity. In particular, we use the McLachlan code mcl (); ? for the evolution of the gravitational variables, the Carpet Schnetter et al. (2004) package for handling adaptive mesh refinement, and the multigrid elliptic solver CT_MultiLevel Bentivegna (2014) to generate initial data; this is then coupled to a new module which evolves the hydrodynamical equations. All equations are discretized using fourth-order finite differencing.

The Einstein Toolkit is routinely used for simulations in relativistic astrophysics, and passes a variety of tests Loffler et al. (2012). Likewise, as will be presented elsewhere, the new module correctly reproduces several exact cosmological models with varying degrees of inhomogeneity. All results presented are convergent at the correct rate as the grid spacing is decreased, and we use this fact to extrapolate the continuum solution of the evolution system, and estimate the error bars resulting from its numerical integration at finite resolution. These are the quantities that appear in all plots.

Perturbations and averaging. We recall two approaches commonly used to solve the evolution system approximately, so that we can compare our solution to these schemes and check that we obtain the correct behavior in the appropriate regime.

For irrotational dust in the synchronous-comoving gauge, the line element can be written (with no loss of generality Landau and Lifshitz (1975)) as , where is the spatial metric. For spacetimes that are close enough to a FLRW model, one can use perturbation theory to follow the departures from the exact background solution. In the matter era, this is the spatially-flat EdS model, with metric , where the scale factor is a solution of Friedmann’s equations


the dot represents a time derivative, and we denote the EdS-background quantities with an overbar. The matter continuity equation gives for the background density. Starting from the inhomogeneous density , one can define the density contrast ; its growth in the synchronous-comoving gauge is governed, at first order, by:


The system of (1) and (2) is then solved by:


where and are the so-called growing and decaying modes. We will use these expressions below as a consistency check in the small-perturbation regime.

Another useful framework is that of cosmological averaging Buchert (2000), where Einstein’s equation is reduced from a set of partial differential equations for the fields to a set of ordinary differential equations in time for some of their averages over a given spatial region . Defining its volume as


where is the determinant of the spatial metric , one finds that the average scale factor satisfies a system similar to Friedmann’s (1), and in particular that




Here is the trace of the extrinsic curvature , , is the traceless part of , and denotes the average of a field over . Note that represents the local expansion rate, and in the FLRW background is the Hubble parameter. Whilst this setup is exact, the computation of itself requires tensorial quantities that do not satisfy ordinary differential equations, i.e. the system of ordinary differential equations for the averaged quantities is not closed. To circumvent this problem, one typically closes the system with a well-motivated ansatz for . One can, for instance, calculate its perturbative behavior: at first order, this term is identically zero, while at second order it scales as Kolb et al. (2005); Li and Schwarz (2008), but with a coefficient containing only surface terms of the averaging volume, which vanish for periodic domains. Beyond this order, the analytical approach becomes exceedingly difficult. A main goal of this Letter is to present an exact measurement of this quantity on an inhomogeneous spacetime.

Results. Our numerical investigation involves the evolution of a cubic domain of coordinate side , with periodic boundary conditions (since we set , will serve as the unit in which all other quantities, including mass and time, are measured). We discretize this domain with points (running two lower resolutions with and to quantify the error bars). We choose the initial density profile as that of the EdS model at the time when the Hubble horizon , plus a superimposed perturbation of initial amplitude (varying between and ) and comoving wavelength :


The ratio for is shown in Fig. 1. As decreases, we expect to recover a cubic domain of the EdS model. By increasing , we should then be able to observe the onset of nonperturbative effects.

Figure 1: Profile of the matter density ratio on the plane () for , when (left) and when (right).

We first need to solve the Einstein constraints; to simplify them, we choose a vanishing traceless part of the extrinsic curvature and a spatially constant . This corresponds to have, initially, a vanishing first-order perturbation of the expansion and a non-zero decaying mode in (4Bruni et al. (2014b). The momentum constraint is then identically satisfied, and the Hamiltonian constraint reduces to the nonlinear elliptic equation:


where . Using (9) and , we solve this equation with CT_MultiLevel Bentivegna (2014), obtaining the initial profile for (normalized to the EdS value) shown in Fig. 2.

Figure 2: Profile of on the plane () for , when (left) and when (right).

We then evolve the coupled gravitational and hydrodynamical equations, until the linear size of the domain has increased by roughly 100 times. We measure the departure of the volume expansion, represented by , from the EdS background model, for different initial amplitudes of the density contrast ; as clearly shown in Fig. 3, this difference is always small. We also monitor the density contrast at the overdensities and underdensities. As expected from linear perturbation theory, and shown in Fig. 4, for small values of the initial the density contrast grows linearly with , with a well-behaved evolution through . For , there is a clear departure from this behavior, with the overdensity becoming nonlinear already at , and eventually growing unbounded when .

Figure 3: Fractional difference of the scale factor of the simulation domain with respect to the EdS scale factor , as a function of the equal-time , for for (top to bottom). The numerical error bars, where visible, are included as shaded regions.
Figure 4: Growth of the density contrast at the overdensities (solid lines) and its negative at the underdensities (dashed lines), for (top to bottom). The linear-perturbation behavior is indicated by dotted lines.

In Fig. 5 we plot the fractional difference of (the local expansion rate) from the background value at the overdensities and underdensities. As expected, the expansion is larger at the underdensities and smaller at the overdensities. For , the departure from the expansion rate of the EdS background is substantial: again, the expansion is already visibly nonlinear at , and the overdensity reaches the turnaround point (signalled by ) at . At turnaround, the linearly extrapolated density contrast is only , much smaller than the standard value from spherical top-hat collapse,  Peacock (1999). For the same initial density contrast, the underdensity asymptotically approaches the expansion of the Milne model (a vacuum FLRW model with negative spatial curvature, represented by the solid gray line in Fig. 5), as predicted in Bruni et al. (1995a), with a fractional departure from EdS of over at . These are the first two important results of our calculations: even in this simple setup, with a perturbation wavelength initially four times larger than the EdS Hubble horizon, the onset of nonlinearity can occur very early, and inhomogeneities can affect the local expansion rate in a substantial, nonperturbative way.

In particular, the observed difference with respect to the spherical homogeneous top-hat collapse is due to the interplay of several factors, most notably the inhomogeneous character of the density, expansion rate, and 3-curvature, and the non-vanishing shear , absent in the top-hat case. Whilst, in the latter case, the perturbation is constrained to remain spatially constant, an inhomogeneous density and expansion accelerate the approach to turnaround at the peak, just like they do in the spherical Newtonian case Bruni et al. (1995b); Rubin and Loeb (2013). The shear also gives a small correction, which for is non negligible even in the initial perturbative regime Bouchet et al. (1992); Tellarini et al. (2015); Bruni et al. (2014b). These effects combine, leading to a negative contribution to the evolution of the local expansion rate , pushing it towards the turnaround (), and accelerating the collapse. The difference with the top-hat collapse is an important issue which we will investigate in detail in future work.

Figure 5: Fractional expansion rate at the overdensities (solid lines) and its negative at the underdensities (dashed lines), for (top to bottom). For the overdensity starts collapsing at . The underdensity with expands much faster than the background, asymptotically approaching the expansion of the Milne model (horizontal dark-gray line).

We then proceed to measure the backreaction quantity : the results are shown in Fig. 6. We extract a few relevant facts: first, given our initial conditions , it follows from the definition (8) that vanishes on the initial time slice. We also notice that, for smaller perturbations, remains zero within our error bars; for larger perturbations, it is clear from Fig. 6 that goes through a short transient phase before following the scaling for a period which is shorter for higher . Given that the only second-order contributions to are boundary terms that vanish on periodic domains like the one we used Li and Schwarz (2008), we conjecture that only higher-order terms are contributing to . Finally, enters the nonperturbative regime, where it is negative and its absolute value increases linearly with the scale factor. The effect is a very small deceleration of the expansion with respect to the EdS model. We conclude that the absolute value of remains generally quite small, but is not identically zero, as would follow from the assumptions of Green and Wald (2011).

Measuring the sign and scaling of the backreaction is a particularly relevant task, as many speculations on the effect of inhomogeneities on the average cosmic expansion rate are based on conjectures on these two properties. A back-of-the-envelope estimate involves the comparison of two competing effects, as quantities like the matter density at the overdensities quickly depart from the background value, but at the same time these regions take up a decreasing fractional volume and become proportionally less and less relevant to the average. Our results indicate that, at least for the specific configuration studied here, the former effect prevails, and the balance is towards an overall slowdown of the expansion rate.

In summary, within the limitations of our setup, in this work we found that, whilst local departures from the background density and expansion rate can be tangible, the average behavior of large volumes remains close to the FLRW background.

Figure 6: Absolute value of the backreaction as a function of the equal-time scale factor in Einstein-de Sitter space, for (top to bottom). The numerical error bars, where visible, are included as shaded regions. For comparison, we have superimposed dashed lines representing the scaling.
We are grateful to Alessio Notari and an anonymous referee for enlightening comments on this work. E.B. is supported by the project “Digitizing the universe: precision modelling for precision cosmology”, funded by the Italian Ministry of Education, University and Research (MIUR). M.B. is supported by the UK STFC Grant No. ST/K00090X/1 and ST/N000668/1. The simulations presented in this paper were carried out on the Sciama supercomputer at the Institute of Cosmology and Gravitation in Portsmouth.


  1. G. F. R. Ellis, General Relativity and Gravitation: A Centennial Perspective, eds. Ashtekar, Berger, Isenberg, MacCallum.  (Cambridge University Press, 2015).
  2. P. A. R. Ade et al. (Planck), “Planck 2015 results. XIII. Cosmological parameters,”  (2015), arXiv:1502.01589 .
  3. T. Buchert, “On average properties of inhomogeneous fluids in general relativity. 1. Dust cosmologies,” Gen. Rel. Grav. 32, 105–125 (2000).
  4. S. Räsänen, “Backreaction: directions of progress,” Class. Quant. Grav. 28, 164008 (2011).
  5. T. Buchert et al., “Is there proof that backreaction of inhomogeneities is irrelevant in cosmology?” Class. Quant. Grav. 32, 215021 (2015).
  6. S. R. Green and R. M. Wald, “Comments on Backreaction,”  (2015), arXiv:1506.06452 .
  7. E. W. Kolb, S. Matarrese, A. Notari,  and A. Riotto, “The Effect of inhomogeneities on the expansion rate of the universe,” Phys. Rev. D71, 023524 (2005).
  8. L. Amendola, “Cosmology and fundamental physics with the Euclid satellite,” Living Rev. Rel. 16, 6 (2013).
  9. R. Maartens, F. B. Abdalla, M. Jarvis,  and M. G. Santos (SKA Cosmology SWG Collaboration), “Cosmology with the SKA-overview,” arXiv:1501.04076  arXiv:1501.04076 .
  10. M. Bruni, D. B. Thomas,  and D. Wands, “Computing General Relativistic effects from Newtonian N-body simulations: Frame dragging in the post-Friedmann approach,” Phys. Rev. D89, 044010 (2014a).
  11. D. B. Thomas, M. Bruni,  and D. Wands, “The fully non-linear post-Friedmann frame-dragging vector potential: Magnitude and time evolution from N-body simulations,” Mon. Not. Roy. Astron. Soc. 452, 1727–1742 (2015).
  12. I. Milillo, D. Bertacca, M. Bruni,  and A. Maselli, “Missing link: A nonlinear post-Friedmann framework for small and large scales,” Phys. Rev. D92, 023519 (2015).
  13. J. Adamek, D. Daverio, R. Durrer,  and M. Kunz, “General Relativistic -body simulations in the weak field limit,” Phys. Rev. D88, 103527 (2013).
  14. J. Adamek, R. Durrer,  and M. Kunz, “N-body methods for relativistic cosmology,” Class. Quant. Grav. 31, 234006 (2014).
  15. J. Adamek, D. Daverio, R. Durrer,  and M. Kunz, “General relativity and cosmic structure formation,” Nature Phys. 12, 346–349 (2016).
  16. V. Cardoso, L. Gualtieri, C. Herdeiro,  and U. Sperhake, “Exploring New Physics Frontiers Through Numerical Relativity,” Living Rev. Relativity 18, 1 (2015).
  17. E. Bentivegna and M. Korzynski, “Evolution of a periodic eight-black-hole lattice in numerical relativity,” Class.Quant.Grav. 29, 165007 (2012).
  18. C.-M. Yoo, H. Abe, Y. Takamori,  and K.-i. Nakao, ‘‘Black Hole Universe: Construction and Analysis of Initial Data,” Phys.Rev. D86, 044027 (2012).
  19. C.-M. Yoo, H. Okawa,  and K.-i. Nakao, “Black Hole Universe: Time Evolution,” Phys.Rev.Lett. 111, 161102 (2013).
  20. E. Bentivegna and M. Korzynski, “Evolution of a family of expanding cubic black-hole lattices in numerical relativity,” Class.Quant.Grav. 30, 235008 (2013).
  21. C.-M. Yoo and H. Okawa, “Black Hole Universe with ,” Phys.Rev. D89, 123502 (2014).
  22. M. Bruni, J. C. Hidalgo, N. Meures,  and D. Wands, “Non-Gaussian Initial Conditions in ΛCDM: Newtonian, Relativistic, and Primordial Contributions,” Astrophys. J. 785, 2 (2014b).
  23. C. Fidler, C. Rampf, T. Tram, R. Crittenden, K. Koyama,  and D. Wands, “General relativistic corrections to -body simulations and the Zel’dovich approximation,” Phys. Rev. D92, 123517 (2015).
  24. T. Nakamura, K. Oohara,  and Y. Kojima, “General relativistic collapse to black holes and gravitational waves from black holes,”  Prog. Theor. Phys. Suppl. 90, 1–218 (1987) .
  25. M. Shibata and T. Nakamura, “Evolution of three-dimensional gravitational waves: Harmonic slicing case,”  Phys. Rev. D 52, 5428–5444 (1995) NoStop
  26. T. W. Baumgarte and S. L. Shapiro, “On the numerical integration of einstein’s field equations,”  Phys. Rev. D 59, 024007 (1998) .
  27. J. A. Font, “Numerical hydrodynamics and magnetohydrodynamics in general relativity,” Living Reviews in Relativity 11 (2008).
  28. A. Lichnerowicz, “L’integration des équations de la gravitation relativiste et le probleme des n corps,” J. Math. Pures Appl. 23, 37 (1944).
  29. J. York, James W., “Gravitational degrees of freedom and the initial-value problem,” Phys. Rev. Lett. 26, 1656–1658 (1971).
  30. P. Anninos and J. McKinney, “Relativistic hydrodynamics of cosmological sheets,” Phys. Rev. D 60, 064011 (1999).
  31. J. T. Giblin, J. B. Mertens,  and G. D. Starkman, “Departures from the Friedmann-Lemaitre-Robertston-Walker Cosmological Model in an Inhomogeneous Universe: A Numerical Examination,” Phys. Rev. Lett. 116, 251301 (2016).
  32. J. B. Mertens, J. T. Giblin,  and G. D. Starkman, “Integration of inhomogeneous cosmological spacetimes in the BSSN formalism,” Phys. Rev. D93, 124059 (2016).
  33. L. D. Landau and E. M. Lifshitz, The classical theory of fields (Oxford: Pergamon Press, 1975).
  34. F. Loffler, J. Faber, E. Bentivegna, T. Bode, P. Diener, et al., “The Einstein Toolkit: A Community Computational Infrastructure for Relativistic Astrophysics,” Class.Quant.Grav. 29, 115001 (2012).
  35. Einstein Toolkit, www.einsteintoolkit.org .
  36. McLachlan code, www.cct.lsu.edu/~eschnett/McLachlan .
  37. Kranc code, www.kranccode.org .
  38. E. Schnetter, S. H. Hawley,  and I. Hawke, “Evolutions in 3d numerical relativity using fixed mesh refinement,”  Class. Quant. Grav. 21, 1465–1488 (2004) .
  39. E. Bentivegna, “Solving the Einstein constraints in periodic spaces with a multigrid approach,” Class.Quant.Grav. 31, 035004 (2014).
  40. N. Li and D. J. Schwarz, “Scale dependence of cosmological backreaction,” Phys. Rev. D78, 083531 (2008).
  41. J. A. Peacock, Cosmological physics (Cambridge, UK: Univ. Pr., 1999).
  42. M. Bruni, S. Matarrese,  and O. Pantano, “Dynamics of silent universes,” Astrophys. J. 445, 958–977 (1995a).
  43. M. Bruni, S. Matarrese,  and O. Pantano, “A Local view of the observable universe,” Phys. Rev. Lett. 74, 1916–1919 (1995b)arXiv:astro-ph/9407054 [astro-ph] .
  44. D. Rubin and A. Loeb, “The virialization density of peaks with general density profiles under spherical collapse,” JCAP 1312, 019 (2013)arXiv:1311.5594 [astro-ph.CO] .
  45. F. R. Bouchet, R. Juszkiewicz, S. Colombi,  and R. Pellat, “Weakly nonlinear gravitational instability for arbitrary Omega,” Astrophys. J. 394, L5–L8 (1992).
  46. M. Tellarini, A. J. Ross, G. Tasinato,  and D. Wands, “Non-local bias in the halo bispectrum with primordial non-Gaussianity,” JCAP 1507, 004 (2015)arXiv:1504.00324 [astro-ph.CO] .
  47. S. R. Green and R. M. Wald, “A new framework for analyzing the effects of small scale inhomogeneities in cosmology,” Phys. Rev. D83, 084020 (2011).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description