Replacing dark energy by silent virialisation

TL;DR

This paper proposes a 'silent' virialisation model in cosmology that explains the observed accelerated expansion without dark energy by accounting for relativistic effects during structure formation.

Contribution

It introduces a novel approach using the QZA to model volume evolution in non-collapsed regions, matching ΛCDM expansion without dark energy.

Findings

Average expansion reaches 16% above EdS at 13.8 Gyr

Strong correlation between virialisation fraction and super-EdS expansion

Model reproduces ΛCDM expansion to first order

Abstract

Standard cosmological $N$-body simulations have background scale factor evolution that is decoupled from non-linear structure formation. Prior to gravitational collapse, kinematical backreaction ($Q_D$) justifies this approach in a Newtonian context. However, the final stages of a gravitational collapse event are sudden; a globally imposed expansion rate thus forces at least one expanding region to suddenly decelerate. This is relativistically unrealistic. Instead, we allow non-collapsed domains to evolve in volume according to the $Q_D$ Zel'dovich Approximation (QZA). We study the inferred average expansion under this "silent" virialisation hypothesis. We set standard (mpgrafic) EdS cosmological $N$-body initial conditions. Using RAMSES, we call DTFE to estimate the initial values of the three invariants of the extrinsic curvature tensor in Lagrangian domains $D$. We integrate the Raychaudhuri equation in each domain using inhomog, adopt the stable clustering hypothesis (VQZA), and average spatially. We adopt an early-epoch--normalised EdS reference-model Hubble constant $H_1^{bg} = 37.7$ km/s/Mpc and an effective Hubble constant $H_0^{eff} = 67.7$ km/s/Mpc. From 2000 simulations at resolution $256^3$, a unity effective scale factor is reached at 13.8~Gyr (16% above EdS) for an averaging scale of $L_{13.8}=2.5^{+0.1}_{-0.4}$ Mpc/$h^{eff}$. Relativistically interpreted, this corresponds to strong average negative curvature evolution. The virialisation fraction and super-EdS expansion correlate strongly at fixed cosmological time. Thus, starting from EdS initial conditions and averaging on a typical non-linear structure formation scale, the VQZA dark-energy--free average expansion matches $\Lambda$CDM expansion to first order. The software packages used here are free-licensed.

Tables

Table 1. Table 1: Software versions for numerical uncertainty tests on I I {\mathrm{I}} , II II {\mathrm{II}} (Sect.  5.1 ), and for VQZA (Sect.  5.2 ) and N 𝑁 N -body–measured 𝒬 𝒟 subscript 𝒬 𝒟 {\cal Q}_{\cal D} (Sect.  5.4 ) simulations.

method	package	version and/or
		commit hash
$I$, $\mathrm{II}$ errors	dtfe	1.1.1.Q, 1f2487e3
VQZA	mpgrafic	0.3.10
VQZA	ramses (ramses/trunk)	438f4659ae36
VQZA	ramses-scalav	53257c54
VQZA	dtfe	1.1.1.Q, 3d514657
VQZA	inhomog	0.0.67
$N$-body ${𝒬}_{𝒟}$	mpgrafic	0.3.10
$N$-body ${𝒬}_{𝒟}$	ramses (ramses/trunk)	438f4659ae36
$N$-body ${𝒬}_{𝒟}$	ramses-scalav	1ba83730
$N$-body ${𝒬}_{𝒟}$	dtfe	1.1.1.Q, 3d514657
$N$-body ${𝒬}_{𝒟}$	inhomog	0.0.64

Equations

$$a_{\mathrm{EdS}}:=(3H_1^{\mathrm{EdS}}t/2{)}^{2/3},H^{\mathrm{EdS}}:={\dot{a}}_{\mathrm{EdS}}/a_{\mathrm{EdS}}=2/(3t),$$

$$n_{\mathcal{D}}:=L_{\mathrm{box}}/L_{\mathcal{D}},n_{\mathrm{DTFE}}:=L_{\mathrm{box}}/L_{\mathrm{DTFE}},n_{\mathrm{soft}}:=L_{\mathrm{box}}/L_{\mathrm{soft}}.$$

$$L_{\mathrm{box}}\gg$$

$$\iff$$

$$\displaystyle\frac{L_{\mathrm{box}}}{\mbox{Mpc/[$H_{1}^{\mathrm{EdS}}$/(100\leavevmode\nobreak\ km/s/Mpc)]}}$$

$$\displaystyle=\frac{L_{\mathrm{box}}}{L_{\mathrm{DTFE}}}\frac{L_{\mathrm{DTFE}}}{\mbox{Mpc/$h_{\mathrm{eff}}$}}\frac{H_{1}^{\mathrm{EdS}}}{H_{\mathrm{eff},0}},$$

$$v_k^i-v_0^i$$

$$v_{,j}^i$$

$${\int }_{S^1}{v}_{,j}^i\mathrm{d}{x}^j$$

$${\hat{\sigma }}_{v_{,j}^i}^{j_1\ne j\ne j_2}$$

$${\hat{\sigma }}_{v_{,j}^i}^j$$

$${\hat{\sigma }}_{v_{,j}^i}^{\mathrm{diag}}$$

$$\mathrm{I}(v_{,j}^i):=$$

$$\mathrm{II}(v_{,j}^i):=$$

$$=$$

$$\mathrm{III}(v_{,j}^i):=$$

$$=$$

$$\displaystyle\;-\left[\;\Big{(}\mathbf{v}(\nabla\cdot\mathbf{v})-(\mathbf{v}\cdot\nabla)\mathbf{v}\Big{)}\cdot\nabla\right]\;\mathbf{v}\Bigg{\}},$$

$${\hat{\sigma }}_{{⟨\mathrm{I}⟩}_{\mathcal{D}}}$$

$${\hat{\sigma }}_{{⟨\mathrm{II}⟩}_{\mathcal{D}}}$$

$${\hat{\sigma }}_{{\mathcal{Q}}_{\mathcal{D}}}$$

$${\mathbf{u}}_1$$

$${\mathbf{u}}_2$$

$$3\frac{{\ddot{a}}_{\mathcal{D}}}{a_{\mathcal{D}}}+4\pi G\frac{M_{{\mathcal{D}}_{\mathbf{i}}}{a}_{{\mathcal{D}}_{\mathbf{i}}}^3}{V_{{\mathcal{D}}_{\mathbf{i}}}{a}_{\mathcal{D}}^3}-\Lambda ={\mathcal{Q}}_{\mathcal{D}},$$

$${\left(\frac{{\dot{a}}_{\mathcal{D}}}{a_{\mathcal{D}}}\right)}^2-\frac{8\pi G}{3}\frac{M_{{\mathcal{D}}_{\mathbf{i}}}{a}_{{\mathcal{D}}_{\mathbf{i}}}^3}{V_{{\mathcal{D}}_{\mathbf{i}}}{a}_{\mathcal{D}}^3}+\frac{{⟨\mathcal{R}⟩}_{\mathcal{D}}}{6}-\frac{\Lambda }{3}=-\frac{{\mathcal{Q}}_{\mathcal{D}}}{6}$$

$${\mathcal{Q}}_{\mathcal{D}}$$

$$a_{\mathbf{i}}:=a_{\mathrm{FL}}(t_{\mathbf{i}})$$

$${\varrho }_{\mathrm{FL}}(t)=\frac{{\varrho }_{\mathrm{FL}}(t_{\mathbf{i}}){a_{\mathbf{i}}}^3}{a(t{)}^3},$$

$$\xi \to \frac{{\dot{q}}_{\mathbf{i}}}{q_{\mathbf{i}}}\xi ,$$

$$\mathrm{II}\to {\left(\frac{q_{\mathbf{i}}}{{\dot{q}}_{\mathbf{i}}}\right)}^2\mathrm{II},$$

$$a_{{\mathcal{D}}_{\mathbf{i}}}$$

$${\dot{a}}_{\mathcal{D}}\left(t_{\mathbf{i}}\right)$$

$$a_{{\mathcal{D}}_{\mathbf{i}}}$$

$${\dot{a}}_{\mathcal{D}}\left(t_{\mathbf{i}}\right)$$

$$3\frac{{\ddot{a}}_{\mathcal{D}}}{a_{\mathcal{D}}}+4\pi G{⟨\varrho ⟩}_{\mathcal{D}}={\mathcal{Q}}_{\mathcal{D}}+\Lambda ,$$

$$3\frac{{\ddot{a}}_{\mathcal{D}}}{a_{\mathcal{D}}}+4\pi G\frac{a_{\mathrm{FL}}^3}{a_{\mathcal{D}}^3}{\varrho }_{\mathrm{FL}}\left(1-{⟨{\mathrm{I}}_{\mathbf{i}}⟩}_{{\mathcal{D}}_{\mathbf{i}}}\right)={\mathcal{Q}}_{\mathcal{D}}+\Lambda$$

Full text

*Aims. *However, the final stages of a gravitational collapse event are sudden; a globally imposed smooth expansion rate forces at least one expanding region to suddenly and instantaneously decelerate in compensation for the virialisation event. This is relativistically unrealistic. A more conservative hypothesis is to allow non-collapsed domains to continue their volume evolution according to the ${\cal Q}_{\cal D}$ Zel’dovich approximation (QZA). We aim to study the inferred average expansion under this ‘silent’ virialisation hypothesis.

*Methods. *We set standard (mpgrafic) EdS 3-torus (T3) cosmological $N$-body initial conditions. Using ramses, we partitioned the volume into domains and called the dtfe library to estimate the per-domain initial values of the three invariants of the extrinsic curvature tensor that determine the QZA. We integrated the Raychaudhuri equation in each domain using the inhomog library, and adopted the stable clustering hypothesis to represent virialisation (VQZA). We spatially averaged to obtain the effective global scale factor. We adopted an early-epoch–normalised EdS reference-model Hubble constant $H_{1}^{\mathrm{EdS}}=37.7$ $\mathrm{km/s/Mpc}$ and an effective Hubble constant $H_{\mathrm{eff},0}=67.7$ $\mathrm{km/s/Mpc}$.

*Results. *From 2000 simulations at resolution $256^{3}$, we find that reaching a unity effective scale factor at 13.8 Gyr (16% above EdS), occurs for an averaging scale of $L_{13.8}=2.5^{+0.1}_{-0.4}$ Mpc/$h_{\mathrm{eff}}$. Relativistically interpreted, this corresponds to strong average negative curvature evolution, with the mean (median) curvature functional $\Omega_{{\cal R}}^{{\cal D}}$ growing from zero to about 1.5–2 by the present. Over 100 realisations, the virialisation fraction and super-EdS expansion correlate strongly at fixed cosmological time.

*Conclusions. *Thus, starting from EdS initial conditions and averaging on a typical non-linear structure formation scale, the VQZA dark-energy–free average expansion matches $\Lambda$CDM expansion to first order. The software packages used here are free-licensed.

1 Introduction

Cosmological $N$-body simulations normally assume that cosmological expansion is gravitationally decoupled from structure formation: scale factor evolution is inserted into the simulation independently of the non-linear structure growth calculated in the simulation. This is a simplification primarily due to the limitations of both analytical and numerical methods of calculation. However, in practice, this decoupling assumption is usually considered to be a corollary of the homogeneity and isotropy assumption underlying the Friedmann–Lemaître–Robertson–Walker (FLRW; Friedmann 1923, 1924; Lemaître 1927, 1931; Robertson 1935) family of cosmological solutions of the Einstein equation, which includes the present standard model of cosmology, the $\Lambda$CDM (cosmological-constant–dominated cold dark matter) model. The validity of the decoupling assumption has long been debated (e.g. Ellis & Stoeger 1987; Buchert 2011, and references therein).

Here, our aim is to take into account the role of virialisation and calculate the effective scale factor evolution in an $N$-body simulation context, guided by general-relativistic constraints defined by scalar averaging (Buchert 2000b, 2001, 2011) that are closely analogous to Newtonian equivalents but are justified relativistically (Buchert, Kerscher, & Sicka 2000; Buchert & Ostermann 2012; Buchert, Nayet, & Wiegand 2013). By providing tools for calculating effective background expansion from structure formation itself, we also encourage the community to explore alternatives to the conservative approach adopted in this paper. In Sect. 2, we explain why the standard treatment of virialisation appears to imply, via an externally imposed background expansion rate, a relativistically unrealistic constraint on spatial expansion. Specific calculations to illustrate the argument are given in Sect. 4.

Given that at early redshifts $z\gg 10$, the observational distinction between the $\Lambda$CDM model and the Einstein-de Sitter (EdS) model (a flat FLRW model with a zero cosmological constant) is very weak, we followed Occam’s razor and adopted an EdS model at early times, setting up standard (mpgrafic) cosmological $N$-body initial conditions. We evolved the initial conditions by a primarily analytical approach – which we refer to here as the kinematical-backreaction (${\cal Q}_{\cal D}$) Zel’dovich approximation (QZA). This provides an analytical expression for kinematical backreaction ${\cal Q}_{\cal D}$ evolution in the Newtonian (Buchert et al. 2000) and general-relativistic contexts (Kasai 1995; Morita et al. 1998 in the form Buchert & Ostermann 2012; Buchert et al. 2013; Alles et al. 2015; see also Matarrese & Terranova 1996; Villa et al. 2011), with closely related but differing definitions and meaning.

It has long been accepted observationally that non-linear structure, that is, statistical patterns of overdensities and underdensities with respect to an ‘average’ (mean) density, can be characterised by a length scale of several megaparsecs (Peebles 1974; Peebles & Hauser 1974). If the average expansion history in an initially EdS simulation is distinct from EdS evolution, then it can reasonably be expected that the difference will depend on the choice of length scale corresponding to non-linear structure. Moreover, if the treatment of volume evolution proposed here is a fair approximation of the evolution of the real Universe, then the length scale that yields the appropriate amount of super-EdS scale factor evolution to match a $\Lambda$CDM approximation of observational data should be reasonably close to the several-megaparsec scale.

The approach presented here adds detail to previous work in modelling inhomogeneous curvature and expansion (Räsänen 2004; Buchert 2008; Wiegand & Buchert 2010; Buchert & Räsänen 2012) that is constrained by the Raychaudhuri equation and the Hamiltonian constraint (whose spatial averages generalise the Friedmannian expansion and acceleration laws). The recent emergence of average negative scalar curvature in cosmic history (Buchert 2005; Buchert & Carfora 2008; Buchert 2008; see also Wiltshire 2007a, b), in which curvature inhomogeneity evolves together with kinematical backreaction (according to a combined conservation law; Buchert 2000b), is found by several of these groups to push the effective scale factor to evolve in a way that potentially mimics ‘dark energy’ on large scales (Roy et al. 2011). Observationally, the non-negligible peculiar expansion rate of voids (Roukema, Ostrowski, & Buchert 2013) and the 10%-level environmental dependence of the baryon acoustic oscillation (BAO) peak location (Roukema et al. 2015, 2016) are consistent with a structure formation alternative to dark energy (Roukema et al. 2017). However, precise observational analysis in a more accurate universe model than the FLRW model requires precise modelling. Specific implementations of emerging negative curvature models have been developed using many different approximations (Räsänen 2006; Nambu & Tanimoto 2005; Kai et al. 2007; Räsänen 2008; Larena et al. 2009; Chiesa et al. 2014; Wiegand & Buchert 2010; Buchert & Räsänen 2012; Wiltshire 2009; Duley et al. 2013; Nazer & Wiltshire 2015; Roukema et al. 2013; Barbosa et al. 2016; Bolejko & Célérier 2010; Lavinto et al. 2013; see also Sussman et al. 2015; Chirinos Isidro et al. 2016; Bolejko 2017a, b; Krasinski 1981, 1982, 1983; Stichel 2016; Coley 2010; Kašpar & Svítek 2014; Rácz et al. 2017).

There have been very recent attempts to run “fully” relativistic cosmological simulations, which nevertheless require numerical shortcuts and strategies in order to obtain results in reasonable amounts of computing time, and, in some methods, adopt standard flat-space tools such as Fourier transforms. However, these have so far failed to show compatibility with the simpler derivations of emerging average negative curvature models cited above. Two possible reasons for this may be that these codes are too tightly coupled to a background whose average curvature is not allowed to evolve, and they do not appear to take virialisation into account. As explained below in Sect. 2, we adopt the stable clustering hypothesis (Peebles 1980) for virialisation, and the QZA for non-virialised domains, without tightly constraining the average curvature. This approach thus allows the emergence of average negative curvature associated with super-EdS expansion rates.

In the terminology of Adamek et al. (2016a, b), the approach we present here aims to model the growth of what these authors call “homogeneous modes” in their hybrid particle–mesh simulations that aim towards consistency with the Einstein equation. Pure mesh-based cosmological simulation codes aiming at consistency with the Einstein equation include those of Bentivegna & Bruni (2016) and Macpherson et al. (2016) using the free-licensed Einstein Toolkit, as well as those of Giblin, Mertens, & Starkman (2016a, b). All of these codes assume a T3 spatial section of their universe model, with the aim being numerical convenience, even though this choice of topology is physically reasonable in terms of the Einstein equation (Aslanyan & Manohar 2012; Ade et al. 2014; Roukema et al. 2014; Aurich 2015; Steiner 2016, and references therein). Our approach only imposes the dynamical constraint of a T3 uniformly flat spatial section in the initial conditions. The mesh-based codes do not yet claim to resolve gravitational collapse (virialisation). The Adamek et al. (2016b) code aims to do this, but is currently tied closely to an FLRW background, and is complementary to the approach presented here (see also Daverio, Dirian, & Mitsou 2016). For approximations related to virialisation and/or gradients of effective pressure in the context of inhomogeneous cosmology, see Buchert (2001), Bolejko & Lasky (2008), Bolejko & Ferreira (2012), and Adamek et al. (2015).

The method adopted here to some degree resembles the Räsänen (2006) approach, in the sense that collapsing and expanding regions are evolved individually, and the Rácz et al. (2017) approach of partitioning $N$-body realisations into small domains and evolving these in order to generate an effective scale factor evolution. However, both of those approaches assume spherical symmetry to calculate the dynamical evolution of their domains, and neglect spatial continuity by assuming that the kinematical backreaction ${\cal Q}_{\cal D}$ is zero, which is only strictly correct in special cases, such as that of spherically symmetric collapse against a spatially flat background [see also Bolejko (2009) regarding Szekeres (1975) models]. We do not assume spherical symmetry of domains, we have spatial continuity, and we include the role of ${\cal Q}_{\cal D}$. Räsänen (2006) does not consider domains whose disjoint union constitutes the global spatial section. Rácz et al. (2017) do consider the density of domains whose union gives the full spatial section, but the domains are in reality non-spherical; the authors’ assumption of spherical symmetry is an approximation that assumes kinematical backreaction to be zero. Their domains appear to be non-Lagrangian: shifting of matter across domains does not seem to be taken into account in their Equations (3) and (4).

Thus, the method of this paper differs from Räsänen (2006) and Rácz et al. (2017) in that we use generic initial condition realisations over a contiguous spatial section, we calculate kinematical backreaction explicitly and use it in the Raychaudhuri equation for evolving each domain (thus bypassing any need to assume or impose spherical symmetry), and we adopt the stable clustering hypothesis for virialisation. After presentation of definitions and method, the method proposed in this paper is summarised in Definition 1 (Sect. 4.2). Future development of our approach will be to make it free from an assumed reference model, as in the Wiegand & Buchert (2010) volume-partitioning approach or the Wiltshire (2007a, b) biscale partition. The latter divides the spatial section contiguously, though not smoothly, into negatively curved void domains and spatially flat collapsing and virialised structures (‘walls’), with a statistical justification for matching walls to the expansion history.

The method is presented in Sect. 3, starting with EdS reference model parameters in Sect. 3.1 and with terminology for the different scales present in $N$-body simulations in general, and for the method presented here in particular, in Sect. 3.2. The generation of initial conditions on a mesh covering the fundamental domain of the T3 spatial section is presented in Sect. 3.3. This is done using mpgrafic (this package and the others used in this work are free-licensed111https://www.gnu.org/licenses/licenses.html software). In Sect. 3.4 we explain how to numerically measure the invariants of what in Newtonian terms is the peculiar velocity gradient tensor. This is only done in the initial conditions for our main model. However, these invariants can also be measured at later time steps for full $N$-body calculations. To the degree that the peculiar gradient tensor approximates the (general-relativistic) extrinsic curvature tensor, our approach should not deviate by much from a relativistic approach, even though the initial conditions are defined in Newtonian terms. The library form of the Delaunay Tessellation Field Estimator dtfe, slightly extended by the addition of a small number of functions and corrected for minor bugs, is used to estimate these invariants. The inhomog library is introduced in Sect. 3.5. This is used to evolve individual cubical domains using the QZA, each of comoving sidelength $L_{{\cal D}}$. For convenience and in order to ease the use of this approach in full $N$-body simulations, this is done using ramses-scalav, an extension to ramses, as a front end to read in the initial conditions and call the dtfe and inhomog libraries. In Sect. 3.6 we define how to relate the domain-based average scale factors $a_{{\cal D}}(t)$ to a global effective scale factor $a_{\mathrm{eff}}(t)$.

As a complement to the main model of this work, a method of evolving initially cubical domains ${\cal D}$ in a full $N$-body simulation, in which ${\cal Q}_{\cal D}$ is estimated numerically at each major time step instead of evolved using the QZA, is presented in Sect. 3.7. In this case, the Raychaudhuri equation is integrated by following the Lagrangian volume defined by the particles of each initial domain and again calculating an effective global expansion rate $a_{\mathrm{eff}}(t)$. This approach again uses ramses-scalav as a front end for reading in initial conditions, and dtfe as the library for calculating the peculiar velocity gradient invariants. The inhomog library is not used in this case.

Since the estimation of ${\cal Q}_{\cal D}$ is crucial to this work, the accuracy of the dtfe numerical estimation of the peculiar velocity gradient invariants (and thus ${\cal Q}_{\cal D}$) is presented in Sect. 5.1. The main results of this work, the comparison of effective scale factor evolution to EdS evolution and $\Lambda$CDM evolution, are presented in Sect. 5.2. In Sect. 5.3, the close relation between virialisation and super-EdS global effective expansion is shown. In Sect. 5.4, we briefly compare our main results with those of calculating ${\cal Q}_{\cal D}$ numerically in our $N$-body comparison method. The relation between Newtonian and relativistic reasoning is discussed again in the light of these results, qualitatively presenting a perfect fluid effective model of virialisation, in Sect. 6. A summary and conclusions are given in Sect. 7. We set the spacetime unit conversion constant (Taylor & Wheeler 1992) $c:=1$ unless otherwise specified. Appendix A provides a computer algebra script for confirming the equivalence of expressions for invariants of the peculiar velocity gradient tensor (Eq. (11)). Appendix B provides the correct expression for the variance of the third invariant, which is wrong in Eqs (C20)–(C22) of Buchert et al. (2000).

2 Virialisation versus spatial expansion

Although the aim here, as in several related projects, is general-relativistic cosmology simulations, it is useful to consider Newtonian reasoning. However, defining a mathematically solid Newtonian cosmological model is not easy. Buchert & Ehlers (1997) showed that a T3 Newtonian cosmological model defined in terms of a dust fluid and gravitational potential, which avoids the divergent infinite sum problem in the pointwise two-point attraction approach, can provide a self-consistent solution (see also Ehlers & Buchert 1997; Ellis & Gibbons 2014 for a recent review and an interesting family of solutions; and recent reminders in Kaiser 2017; Buchert 2017). In this case, if ${\cal Q}_{\cal D}$ evolution in local domains is calculated within the constraints of the same (expanding) T3 ‘absolute’ space, then the global scale factor evolution of an initially EdS model retains EdS behaviour (proportionality to $t^{2/3}$), even when scale factor evolution is calculated in the local domains ${\cal D}$ using ${\cal Q}_{\cal D}$ and averaged. In terms of kinematical backreaction, the global ${\cal Q}_{\mathrm{T}^{3}}=0$ (see Sect. 3.4.2). [Equation (61) of Buchert et al. (2000) is incorrect; the middle element of the double equation has to be omitted to make the equation correct.] However, this justification for exact EdS global expansion is only valid until gravitational collapse and virialisation (or more generally, shell-crossing). Before discussing virialisation, let us first consider the QZA, that is, the algebraic expression for evaluating ${\cal Q}_{\cal D}$ time evolution as presented in Buchert et al. (2000) and Buchert et al. (2013) in the spirit of the Zel’dovich approximation (Zel’dovich 1970).

The QZA evolution of the kinematical backreaction ${\cal Q}_{\cal D}$ has the same algebraic structure in the Newtonian [Buchert et al. 2000, eq. (42); Buchert et al. 2013, eq. (57); NZA] and relativistic [Buchert et al. 2013, eq. (50)] cases, but with differing spatial definitions and interpretations (Buchert et al. 2013, IV.A); we cite this expression in Eqs (LABEL:e-resultQ2), (31). This expression for ${\cal Q}_{\cal D}$ evolution is exact for some particular cases in the Newtonian context (Buchert et al. 2000, III.D, III.E). In the relativistic context, it is exact for plane symmetric collapse [Buchert et al. 2013, eq. (75)] and provides exact solutions for certain cases of spherically symmetric collapse (Lemaître 1933; Tolman 1934; Bondi 1947, LTB) models (Buchert et al. 2013, V.B.1). The QZA requires knowing the growth function for linear perturbations in a reference model (in our case, the EdS reference model), and of the initial invariants of the extrinsic curvature tensor. Here, we estimate the latter using the Newtonian peculiar velocity gradients, generated for a standard T3 cosmological $N$-body simulation. The initial invariants in non-global domains should collectively satisfy a partitioning formula [Buchert & Carfora 2008, Sect. 3.2; Wiegand & Buchert 2010, eqs (16), (17)] that yields globally zero mean invariants due to the T3 constraint; see Sect. 3.4.2 for details.

How are virialisation, Newtonian and relativistic reasoning related? The assumptions of the Newtonian collisionless, compressible fluid approach for exact spherically symmetric collapse, interpreted strictly, imply collapse to a singularity in a finite time. A Lagrangian model is easier to relate to a relativistic model than an Eulerian model in the sense that the Lagrangian approach ties matter and space together directly. However, in this idealised case of spherical collapse, a positive Lagrangian volume drops suddenly to zero at the final moment of collapse. In other words, the validity of the approach is limited to times before gravitational collapse events occur. More generically, the assumptions fail when shell crossing first occurs. In $N$-body simulations or with physical particles, the usual way to avoid these Newtonian singularities is to adopt a hybrid fluid–particle model: the fluid model allows the model to be defined globally as shown by Buchert & Ehlers (1997) (or approximated regionally in an adaptive mesh refinement method), while virialisation of a set of particles is expected (and modelled) to occur instead of singularity formation. At a small enough scale, baryon pressure also opposes gravitational collapse. Disregarding baryon pressure, overdensities in a T3 $N$-body simulation cannot be strictly spherically symmetric because they are composed of particles; and they cannot evolve in an exact spherically symmetric way because the global spatial section is T3. Moreover, Gaussian initial conditions will generically lead to at least moderately anisotropic collapse. Finally, excessively high two-particle close-encounter accelerations are avoided by a small-length-scale softening parameter. Thus, the limitations of the Newtonian fluid approach are dealt with in $N$-body simulations by considering singularity formation to be unrealistic (unless a simulation is detailed enough to include active galactic nucleus formation and star formation through to supermassive and stellar black hole formation).

This avoidance of singularity formation is generally modelled without any overt compensation in volume evolution at the global T3 volume evolution level. Below we argue that in the standard approach (whether analytical or $N$-body), there is, in fact, hidden compensation in volume evolution at the global level.

In this work, we are not interested in modelling the details of virialisation. What is relevant for a relativistic interpretation is whether or not compensation between volume loss and gain in collapsing and expanding (respectively) Lagrangian domains is included in the model. It is clear that, when interpreting observations in terms of FLRW models, the virialisation of high-mass collapsed objects correlates closely with the cosmologically recent appearance of a non-negligible value of the dark energy parameter ($\Omega_{\Lambda}$). This observational clue is quantified in Roukema et al. (2013). Relativistically, the assumption of no sudden compensation between virialising and expanding domains can be described as a ‘silent virialisation’ approximation (following Matarrese et al. 1994b, a for the somewhat different “silent universe” approximation; see also Ellis & Tsagas 2002; Bolejko 2017b).

In Sect. 4, we present a two-domain partition of the T3 volume that illustrates our claim that the standard approach of passing from pre- to post-virialisation is relativistically unrealistic. In Sect. 4.1, we define the biscale partitions and study volume evolution prior to the collapse and virialisation event. In Sect. 4.2, we present the dilemma between choosing the standard approach versus using the QZA for the expanding domain and the stable clustering approximation for virialised domains (presented in more detail for our main multi-domain partition modelling in Sect. 3.5.6). We thus define the Virialisation ${\cal Q}_{\cal D}$ Zel’dovich Approximation.

3 Method

3.1 EdS reference model extrapolated from early epochs

The use of an Einstein–de Sitter (EdS) reference model at early times (that we previously called a “background” EdS model; Roukema et al. 2013, 2017) or on large spatial scales in inhomogeneous cosmological models that aim to be more accurate than the FLRW model risks leading to confusion when referring to ‘the’ Hubble constant, since several different values compete for this name. Here, as in Roukema et al. (2017), we set up the EdS reference model in terms of a reference-model scale factor $a_{\mathrm{EdS}}$ and a reference-model expansion rate $H^{\mathrm{EdS}}$ given by

[TABLE]

where $H_{1}^{\mathrm{EdS}}:=H^{\mathrm{EdS}}(a_{\mathrm{EdS}}=1)$. This EdS model is intended to approximately fit early epochs. In order to derive an observationally realistic value of $H_{1}^{\mathrm{EdS}}$, an effective scale factor $a_{\mathrm{eff}}$ which is nearly identical to the reference model scale factor at early times, that is, $a_{\mathrm{eff}}\approx a_{\mathrm{EdS}}$ for small $t$, is adopted. This is done by normalising to $a_{\mathrm{eff}}=1$ at the present time ${t_{0}}\equiv{t_{a_{\mathrm{eff}}=1}}$, where the Planck Surveyor (Table 4, sixth data column, Ade et al. 2016) parametrisation of the $\Lambda$CDM model is used as a proxy for cosmological observations. As explained in Roukema et al. (2017), this requires that we adopt an early-epoch–normalised EdS–reference-model Hubble constant $H_{1}^{\mathrm{EdS}}=37.7$ $\mathrm{km/s/Mpc}$ (Rácz et al. 2017 adopt this value too) and an effective Hubble constant – the limiting value that should be observed in the local few hundred Mpc according to an FLRW fit of the data – $H_{\mathrm{eff},0}=67.74$ $\mathrm{km/s/Mpc}$.

3.2 Scales

A standard T3 $N$-body simulation has three built-in length scales:

• (i)

$L_{\mathrm{box}}$ – the side length of the fundamental domain, often called the box size,

• (iv)

$L_{N}$ – the mean interparticle separation along one of the three fundamental directions between ‘adjacent’ particles among a set of $N^{1/3}$ particles sorted along that axis, where the simulation contains $N$ particles; that is, $L_{N}:=L_{\mathrm{box}}/N^{1/3}$, and

• (v)

$L_{\mathrm{soft}}$ – the softening length of Newtonian two-point instantaneous attraction; in the case of ramses (Sect. 3.7), this can be considered to be the maximum level of resolution for calculating the Newtonian gravitational potential in an adaptively defined cell, $L_{\mathrm{soft}}:=2^{-{\mathrm{levelmax}}}L_{\mathrm{box}}$.

We insert two additional scales between $L_{\mathrm{box}}$ and $L_{N}$ (we order the definitions from (i) to (v) to monotonically match length scales). Non-linear structure formation has characteristic scales (Peebles 1974; Peebles & Hauser 1974), so as in earlier work, we need to set a scale at which we estimate the kinematical backreaction ${\cal Q}_{\cal D}$ and apply the scalar averaged evolution equations (Sect. 3.5). Since ${\cal Q}_{\cal D}$ is an average of the peculiar velocity gradient tensor invariants, the latter should be estimated at a smaller scale. However, to reduce the Poisson noise effects of using a finite number of particles, the scale at which these invariants are estimated (Sect. 3.4) should be greater than $L_{N}$. Thus, our intermediate scales are

• (ii)

$L_{{\cal D}}$ – the scalar averaging initial domain size in comoving units (Sect. 3.5), and

• (iii)

$L_{\mathrm{DTFE}}$ – the dtfe mesh size (Sect. 3.4).

We correspondingly define

[TABLE]

Using this terminology, numerical scalar-averaging modelling of cosmological expansion requires, in general,

[TABLE]

for an $N$-particle simulation, where the validity of the minimal ratios in this multiple inequality needs to be studied both numerically and analytically. Setting the first three ratios of these values to $8$ or $16$ would require $N=512^{3}$ or $4096^{3}$, respectively. Given these heavy requirements in computer resources, we limit this initial study to numerical exploration of the modest value $N=256^{3}$. For our main calculations, in which the QZA uses only the initial values of the invariants of the peculiar velocity gradient (in a Newtonian interpretation), $n_{\mathrm{soft}}$ is irrelevant. For the $N$-body comparison, we use $N=128^{3}$ and set $n_{\mathrm{soft}}/(N^{1/3})=32$.

3.3 Initial conditions

We use the free-licensed (GNU General Public License, version 2 or later; GPL-2$+$) package mpgrafic-0.3.10 (Bertschinger 2001; Prunet et al. 2008)222http://www2.iap.fr/users/pichon/mpgrafic.html, https://bitbucket.org/broukema/mpgrafic to generate cosmological initial conditions for an EdS model, with $H_{1}^{\mathrm{EdS}}=37.7$ $\mathrm{km/s/Mpc}$ – which, in the absence of structure formation, would give a unity scale factor (for the reference model) at the unrealistically late foliation time of 17.3 Gyr [see Eq. (1)].

The comoving fundamental domain size $L_{\mathrm{box}}$ (see Sect. 3.2) needs to be expressed in units of Mpc/[$H_{1}^{\mathrm{EdS}}$/(100 km/s/Mpc)], since the simulation starts with the EdS reference model. However, in order for length scales to be interpretable in terms of standard descriptions of low-redshift observations, it is more useful to choose $L_{\mathrm{box}}$ based on a given low-redshift length scale $L_{\mathrm{DTFE}}$, so as above, we adopt the Ade et al. (2016) estimate of $H_{\mathrm{eff},0}=67.74$ $\mathrm{km/s/Mpc}$. Thus, we have

[TABLE]

where $h_{\mathrm{eff}}:=H_{\mathrm{eff},0}/$100 $\mathrm{km/s/Mpc}$.

We match the Ade et al. (2016) power spectrum normalisation of $\sigma_{8}^{\Lambda\mathrm{CDM}}=0.8159$ to our EdS reference model. We evolve the effective Planck-parametrised $\Lambda$CDM model backwards to $z=1000$, that is, to $t=547$ Myr, and then forwards using our EdS model, obtaining $\sigma_{8}=1.0395$ when $a_{\mathrm{EdS}}=1$, which we use for mpgrafic initial conditions.

The value of $L_{N}$ is set in mpgrafic, and later read in automatically by ramses.

3.4 Estimating ‘velocity gradient’ invariants

In order to model what from a general-relativistic point of view is the extrinsic curvature tensor, we numerically estimate what in Newtonian terms is the peculiar velocity gradient tensor. That is, we make the standard assumption that world lines of particles can have their time derivatives separated into ‘peculiar velocities’ and a refence model expansion $a_{\mathrm{EdS}}(t)$.

Delaunay and Voronoi tessellation are two methods that Bernardeau & van de Weygaert (1996) argued provide close to optimal tracing of the (Newtonian) peculiar velocity field. Here, we use the former, as implemented in the Delaunay Tessellation Field Estimator (dtfe)333DTFE-1.1.1.Q: https://bitbucket.org/broukema/dtfe;

DTFE-1.1.1:

https://www.astro.rug.nl/%7Evoronoi/DTFE/dtfe.html which is a free-licensed code that has been found to be highly competitive in comparison with other codes that aim to extract peculiar velocity fields from cosmological $N$-body simulations (Schaap & van de Weygaert 2000; van de Weygaert & Schaap 2009; Cautun & van de Weygaert 2011; Kennel 2004), with good tracing of the details of the cosmic web of virialised structure (Aragón-Calvo et al. 2010; Sousbie 2011; Aragon-Calvo et al. 2010; Park et al. 2013; Pranav et al. 2017). Dtfe uses the free-licensed library cgal (Computational Geometry Algorithms Library444https://www.cgal.org/) to make Delaunay tessellations.

A Delaunay tessellation of an $N$-body simulation snapshot is a division of the T3 volume into tetrahedra using the particles’ positions, in such a way that none of the particles lies inside the circum-2-sphere of any of the tetrahedra. Given one such tetrahedron with vertices $k=0,1,2,3$ and simulation peculiar velocity vectors at these vertices $v^{i}_{k}$, $i=1,2,3$, a linear interpolation of the peculiar velocity gradient $v^{i}_{,j}$ relates velocity component changes from the 0-th to the $k$-th vertex for $k\not=0$

[TABLE]

for the usual Einstein summation notation over corresponding upper/lower indices and a comma in the subscript to indicate a partial derivative. If the matrix $x^{j}_{k}-x^{j}_{0}$ is invertible (which, in numerical practice, should almost always be the case), the velocity gradient can be estimated as

[TABLE]

We use dtfe to carry out this linear interpolation using dtfe “method 1”. This consists of Monte Carlo sampling with a quasi-random spatial distribution within a Delaunay tetrahedron, giving interpolated values that accumulate in cubical cells of a regular grid of ${n_{\mathrm{DTFE}}}^{3}$ cells within the simulation box, yielding an averaged interpolated estimate of $v^{i}_{,j}$ in each DTFE grid cell.

3.4.1 S1-based estimates of velocity gradient uncertainties

We introduce a method of estimating the dtfe grid estimates of $v^{i}_{,j}$ that does not seem to have been used previously. This new method is based on the numerical estimate of the integral of $v^{i}_{,j}$ over any closed straight loop S1 passing through the $n_{\mathrm{DTFE}}$ points ${\mathbf{x}}_{k}=\{x^{j}_{k}\}_{j}$ in an ${n_{\mathrm{DTFE}}}^{3}$ cubical DTFE grid of velocity gradient estimates in a time slice of a simulation, where $k\in\{0,\ldots,n_{\mathrm{DTFE}}-1\}$, and $x^{j_{1}}$, $x^{j_{2}}$ are two fixed values in the $j_{1}$-th, $j_{2}$-th directions, respectively, with $j_{1}\neq j\neq j_{2}$. Analytically,

[TABLE]

should hold by the second fundamental theorem of calculus (Stokes’ theorem) if $v^{i}_{,j}$ is continuous and real-valued. Provided that we allow the usual S-shaped peculiar velocity profiles across filaments and virialised objects, this should be satisfied by a standard interpretation of the information numerically represented by the simulation particles.

Numerically, we make the assumption that the numerical uncertainty $\sigma_{v^{i}_{,j}}$ is Gaussian and identically and independently distributed at each of the dtfe grid points along a given S1 loop. This is obviously an oversimplification. Collapsing or expanding structures on scales greater than $L_{\mathrm{DTFE}}$ could reasonably lead to correlations in numerical errors in the estimates of $v^{i}_{,j}$ across close pairs of dtfe grid points, and tensor invariants of $v^{i}_{,j}$ have represent symmetries that may exacerbate these correlations. It is also quite realistic for the per-grid-point errors to be non-gaussian. Nevertheless, this order of magnitude estimate of the error depends on very few assumptions, and provides a consistency check that can be interpreted in terms of expected failures of these minimal assumptions.

Thus, we introduce

[TABLE]

where we use $n_{\mathrm{DTFE}}-1$ to increase the error estimate slightly at low $n_{\mathrm{DTFE}}$, and the absolute value becomes irrelevant when this formula is used in Eq. (9) below. There is no constraint $i=j$: the velocity components’ derivatives need to correspond to the direction of integration, but the velocity components themselves do not need to. This provides $n_{\mathrm{DTFE}}^{2}$ estimates for a fixed direction $j$. We define the rms (root mean square) of these $n_{\mathrm{DTFE}}^{2}$ error estimates of the $i$-th component summed in the $j$-th direction:

[TABLE]

where $k_{1},k_{2}$ label the grid positions in the $j_{1}$-th, $j_{2}$-th directions, respectively. We introduce this into dtfe with the function estimateSigmaGradientOneDirec. For propagating these uncertainties to estimates in uncertainties in the scalar invariants of $v^{i}_{,j}$ (Sect. 3.4.2), we estimate the rms over the diagonal and off-diagonal components of $v^{i}_{,j}$ separately:

[TABLE]

These two error estimators summarise and simplify the uncertainty information from the S1 constraint from all velocity gradient components in all three orthogonal directions (of the simulation interpreted as a Newtonian simulation) over the whole fundamental domain. The motivation for separating diagonal from off-diagonal components is that in a physically realistic simulation, this might help to separate different types of error.

3.4.2 S1 and T3-based uncertainties in

${\mathrm{I}}$ and ${\mathrm{II}}$

In flat space, the tensor invariants of the peculiar velocity gradient can be written (Buchert 1994; Ehlers & Buchert 1997; see App. A for a derivation)

[TABLE]

where $\mathbf{v}:=(v^{0},v^{1},v^{2})$. Thus, all three invariants can be expressed as divergences, so that Stokes’ theorem again applies, this time over the full T3. Again, sums of ${\mathrm{I}}$, ${\mathrm{II}}$, or ${\mathrm{III}}$ over the full simulation volume should analytically be zero, but numerically should indicate the level of numerical noise. However, assuming that the dtfe grid over the full box gives statistically independent and identically distributed errors is likely to be a worse approximation than assuming this within any S1 straight loop. Nevertheless, global means of the invariants and of kinematical backreaction [see Eq. (16) below] can be checked for consistency with zero based on the S1 method of error estimation.

Again assuming statistically independent gaussian errors and using Eqs (11) and (16), we obtain an S1-based uncertainty for $\left\langle{\mathrm{I}}\right\rangle_{\cal D},\left\langle{\mathrm{II}}\right\rangle_{\cal D},$ and ${{\cal Q}}_{{\cal D}}$ of

[TABLE]

where $\mathbf{x}$ represents all $n_{\mathrm{DTFE}}^{3}$ dtfe grid points and $\delta^{j}_{i}\in\{0,1\}$ is the usual Kronecker delta.

To test dtfe’s numerical accuracy in estimating velocity field gradients, we define analytical velocity fields

[TABLE]

for which ${\mathrm{I}}$ and ${\mathrm{II}}$ are straightforward to calculate analytically; the only non-zero gradient components are aligned with the sinusoidal directions in $\mathbf{u}_{1}$, and orthogonal to them in $\mathbf{u}_{2}$ (giving ${\mathrm{I}}(\mathbf{u}_{2})=0={\mathrm{II}}(\mathbf{u}_{2})$). We study the two following questions:

• (i)

What is the signal-to-noise ratio ($S/N$) of the rms of $\mathbf{u}_{1}$ or $\mathbf{u}_{2}$, where the signal is defined by Eq. (13) and the noise $\sigma_{\mathrm{meas}}$ is the rms of the difference between the analytical expression for ${\mathrm{I}}$ or ${\mathrm{II}}$ and the numerical estimate?

• (ii)

What is the ratio of $\sigma_{\mathrm{meas}}$ to the rms error predicted by the S1-based error estimators indicated above [Eqs (8)–(10), (12)]?

The results of these calculations using test_vgrad_DTFE.cpp in dtfe are presented in Sect. 5.1.

3.5 ${\cal Q}_{\cal D}$ Zel’dovich approximation (QZA)

Here we present our method of integrating the averaged Raychaudhuri evolution at the $L_{{\cal D}}$ scale, used in our method (i), that we implement in the inhomog free-licensed (GPL-2$+$) library555https://bitbucket.org/broukema/inhomog, https://tracker.debian.org/inhomog. Averaging of $a_{{\cal D}}$ within the full volume, that is, at the box scale $L_{\mathrm{box}}$, is discussed below in Sect. 3.6. For the reader’s convenience, we make this subsection slightly more general than is needed for the specific method adopted here: we include the case of a non-zero cosmological constant $\Lambda$.

3.5.1 Scalar averaging equations

The averaged Raychaudhuri equation in the relativistic case [Buchert et al. 2013, (9)] is:

[TABLE]

where $a_{{\cal D}}$ is defined in the usual way [Buchert et al. 2013, eqs (2), (3)], we write the per-domain scale factor and volume evolution in terms of their initial ($\mathbf{i}$) values $a_{{\cal D}}^{3}/a_{{{{\cal D}}_{\mathbf{i}}}}^{3}=V_{{\cal D}}(t)/V_{{{{\cal D}}_{\mathbf{i}}}}$, $M_{{{{\cal D}}_{\mathbf{i}}}}$ is the mass in the domain ${\cal D}$, and the averaged Hamiltonian constraint [Buchert et al. 2013, eq. (10)] is:

[TABLE]

(Buchert 2000a, b, 2001), where the kinematical backreaction is defined in terms of the invariants of the extrinsic curvature tensor, approximated in terms of the invariants of the velocity gradient tensor on a domain ${\cal D}$:

[TABLE]

For a more physically motivated definition in the Newtonian case, see Buchert et al. (2000, II.B., Eq. (5), App. A).

We use dtfe to numerically estimate the initial values of these invariants on a regular cubical mesh of ${n_{\mathrm{DTFE}}}^{3}$ comoving positions in the EdS reference model. Thus, each domain ${\cal D}$ (again a cubical cell) contains $(n_{\mathrm{DTFE}}/n_{{\cal D}})^{3}$ estimates of ${{\mathrm{I}}}$ and ${{\mathrm{II}}}$. We average these to obtain $\left\langle{\mathrm{I}}\right\rangle_{\cal D}$ and $\left\langle{\mathrm{II}}\right\rangle_{\cal D}$. As stated above [Eq. (3)], it should be preferable to have $n_{\mathrm{DTFE}}/n_{{\cal D}}\gg 1$.

We only consider the $\Lambda=0$ (dark-energy–free) case here, since Occam’s razor favours first calculating a DE-free model with expansion generated consistently with structure formation, which is the main theme of this paper. Non-flat FLRW backgrounds are more difficult to model correctly, since Fourier analysis is invalid in these cases; the power spectrum needs to be expressed in terms of an orthonormal basis for a constant non-zero curvature 3-manifold of the appropriate topology.

3.5.2 Initial invariants and other initial conditions

The inhomog software allows the reference model scale factor $a(t)$ to be normalised at unity either at the initial time or at the present, setting (cf Buchert et al. 2000 C.1 last paragraph)

[TABLE]

so that

[TABLE]

where the subscript ‘FL’ indicates the FLRW reference model. In this paper, the reference model is the EdS model described in Sect. 3.1, and $a_{\mathrm{FL}}({{t}_{\mathbf{i}}})$ is set as described in Sect. 3.3. In Buchert et al. (2013), the reference model is termed a background.

Buchert et al. (2000) internally alternates between defining the invariants of the extrinsic cuvature tensor and the (normalised and zeroed) growth function $\xi$ [(34) below] to have a time dimension ($[{\mathrm{I}}]=\mathrm{T}^{-1}$, $[{\mathrm{II}}]=\mathrm{T}^{-2}$, $[{\mathrm{III}}]=\mathrm{T}^{-3}$, $[\xi]=\mathrm{T}$, where $[x]$ denotes the dimension of $x$ and $\mathrm{T}$ is the time dimension) or to be dimensionless. Buchert et al. (2013) uses the dimensional definitions throughout, except in Section VI.A. Here we adopt the dimensionless definitions. This requires the replacements

[TABLE]

everywhere in the text of Buchert et al. (2013) except for Section VI.A, where the growing mode $q_{\mathrm{FL}}(t)$ is given in Eqs (35) or (LABEL:e-growth-flat). In many formulae, the replacements cancel so that the published formulae are unaffected by this change of convention.

Numerical evaluation of these invariants for the RZA version of the QZA formalism developed in Buchert et al. (2000); Buchert & Ostermann (2012); Buchert et al. (2013) requires using a power spectrum in the reference model, so the spatial section of the reference model must be $\mathbb{R}^{3}$ or $T^{3}$ (some other flat FLRW spatial sections also allow Fourier analysis). Here we adopt a standard $\mathbb{R}^{3}$ power spectrum (Eisenstein & Hu 1998).

Buchert et al. (2013) VI.A gives the initial conditions for ${}^{{\rm RZA}}a_{{\cal D}}$ (for brevity we drop the ‘RZA’ pre-superscript), following from (A3) and (49) in Buchert et al. (2000), where here we use the dimensionless definition of the growth rate $\xi$ and the invariants, and a bold $\mathbf{i}$ subscript to indicate that the invariants are calculated at the initial time on the domain ${\cal D}$ in Lagrangian coordinates (${{{\cal D}}_{\mathbf{i}}}$),

[TABLE]

An alternative choice would be to use the Hamiltonian constraint (15) and assume that the initial epoch is early enough that turnaround has not yet been reached, $\dot{a}_{{\cal D}}({{t}_{\mathbf{i}}})>0$, in order to select the positive square root for $\dot{a}_{{\cal D}}({{t}_{\mathbf{i}}})$, yielding

[TABLE]

where we postpone expresssions for $\left\langle{\cal R}\right\rangle_{\cal D}\left({{t}_{\mathbf{i}}}\right)$ and ${\cal Q}_{{\cal D}}\left({{t}_{\mathbf{i}}}\right)$ to (37) and (38) below.

3.5.3 Evolution equation

Writing the per-domain density $\rho_{{\cal D}}=\left(M_{{{{\cal D}}_{\mathbf{i}}}}/V_{{{{\cal D}}_{\mathbf{i}}}}\right)\,a_{{{{\cal D}}_{\mathbf{i}}}}^{3}/a_{{\cal D}}^{3}$, equation (14) becomes

[TABLE]

and with Buchert et al. (2013, eq. (92)) (with a dimensionless definition of $\xi$ and the invariants) can be expressed as

[TABLE]

and thus, using (18),

[TABLE]

where the constant $A$ is defined

[TABLE]

with $H_{\mathrm{FL}}:=\dot{a}_{\mathrm{FL}}/a_{\mathrm{FL}}$ and $\Omega_{\mathrm{m}}:=8\pi G\rho_{\mathrm{FL}}/(3H_{\mathrm{FL}}^{2})$ the usual FLRW matter density parameter. For the EdS reference model we have

[TABLE]

This allows (22) to be rewritten

[TABLE]

and the Hamiltonian evolution equation becomes

[TABLE]

Typical collapsing solutions start with $\dot{a}_{{\cal D}}>0$, the positive square root. Finding the first local minimum of $\dot{a}_{{\cal D}}^{2}$ in the positive square root case yields an estimate of the turnaround epoch $t_{\mathrm{turn}}$ (when $\dot{a}_{{\cal D}}$ drops to zero), enabling a switch to the negative square root for an initial approximate solution that continues through to collapse. Iterative improvement of the estimates of $t_{\mathrm{turn}}$ and $a_{{\cal D}}(t),\dot{a}_{{\cal D}}(t)$ yields improved accuracy.

3.5.4 Auxiliary equations

To numerically solve Eq. (15) would require estimating the evolution of both the kinematical backreaction ${\cal Q}_{{\cal D}}$ and the mean curvature $\left\langle{\cal R}\right\rangle_{\cal D}$. The kinematical backreaction is given by Buchert et al. (2013) (50) (or its Newtonian equivalent):

[TABLE]

where

[TABLE]

and the subscript ${\cal I}$ indicates that a fully relativistic calculation would integrate the curved, Riemannian volume even at this early time when perturbations are weak, and the dimensionless growth rate $\xi$ is given in Eq. (34). The mean curvature is given by Buchert et al. (2013) (13), (54) for the flat FLRW background case:

[TABLE]

where

[TABLE]

As detailed below in Sect. 3.5.5, we bypass the direct dependence on $\left\langle{\cal R}\right\rangle_{\cal D}$ by integrating Eq. (25) instead. Evaluation of these expressions requires the growth rate $\xi(t)$, for which we use the dimensionless form of Buchert et al. (2013) (32),

[TABLE]

For the EdS reference model, the growing mode is

[TABLE]

(for example, Chapter 15, Weinberg 1972; eq. (21) Bildhauer et al. 1992). For completeness, the growing mode for low-density flat FLRW backgrounds (such as $\Lambda$CDM) with present density parameter $\Omega_{\mathrm{m0}}<1$ and $\Omega_{\Lambda 0}:=1-\Omega_{\mathrm{m0}}$ is:

[TABLE]

[Bildhauer et al. 1992, (21)].

Since ${{\xi}_{\mathbf{i}}}=0$ by definition (see (34)), the initial backreaction terms needed in (28) follow from (LABEL:e-resultQ2)–(33):

[TABLE]

and

[TABLE]

For EdS, we have $a_{\mathrm{FL}}(t)={{a}_{\mathbf{i}}}(t/{{t}_{\mathbf{i}}})^{2/3}$, so from (17) and (27) we have

[TABLE]

The standard FLRW expressions can be used for the low-density flat background case.

3.5.5 Numerical strategy

The Raychaudhuri equation (25) is a second-order ordinary differential equation (ODE), which can be reduced to two first-order equations

[TABLE]

where $y\equiv\dot{a}_{{\cal D}}$ is defined to clarify the numerical strategy. This pair of first-order equations can be solved numerically by calculating ${\cal Q}_{\cal D}(t)$ from (LABEL:e-resultQ2) and (31) and the initial values of the invariants $\left\langle{{{\mathrm{I}}}_{\mathbf{i}}}\right\rangle_{{\cal D}_{\mathbf{i}}},\left\langle{{{\mathrm{II}}}_{\mathbf{i}}}\right\rangle_{{\cal D}_{\mathbf{i}}},$ and $\left\langle{{{\mathrm{III}}}_{\mathbf{i}}}\right\rangle_{{\cal D}_{\mathbf{i}}}$, and setting the scale factor initial conditions from Eq. (20) and the constant $A$ from Eqs (27) or (39), and the growth function from Eq. (35). In principle, as in Ostrowski et al. (2018), typical values of $\left\langle{{{\mathrm{I}}}_{\mathbf{i}}}\right\rangle_{{\cal D}_{\mathbf{i}}},\left\langle{{{\mathrm{II}}}_{\mathbf{i}}}\right\rangle_{{\cal D}_{\mathbf{i}}},$ and $\left\langle{{{\mathrm{III}}}_{\mathbf{i}}}\right\rangle_{{\cal D}_{\mathbf{i}}}$, could be evaluated under the assumption of Gaussianity of the initial invariants, using eqs (C5), (C14) and (C18) of Buchert et al. (2000) and the corrected equivalent of eqs (C20)–(C22) of Buchert et al. (2000), which we provide here in Eq. (LABEL:e-BKS00-C20-corrected) in App. B, since it has not been published previously. However, in this work, we estimate $\left\langle{{{\mathrm{I}}}_{\mathbf{i}}}\right\rangle_{{\cal D}_{\mathbf{i}}},\left\langle{{{\mathrm{II}}}_{\mathbf{i}}}\right\rangle_{{\cal D}_{\mathbf{i}}},$ and $\left\langle{{{\mathrm{III}}}_{\mathbf{i}}}\right\rangle_{{\cal D}_{\mathbf{i}}}$ from the initial conditions generated by mpgrafic, which only assumes Gaussianity of the density fluctuations. Thus, the distributions of the second and third invariants of the peculiar velocity gradients are not constrained to be Gaussian.

Numerical solvers for ODEs often require Jacobians. In the case of Eq. (40), these are

[TABLE]

We solve these equations using the embedded Runge–Kutta–Fehlberg $(4,5)$ method of the GNU Scientific Library (GSL) routine gsl_odeiv_step_rkf45.

The averaged density parameters then follow from (18) and (19) in Buchert et al. (2013),

[TABLE]

3.5.6 Virialisation

After turnaround of a given domain ${\cal D}$, that is, when $\dot{a}_{{\cal D}}$ becomes negative, we continue integration towards numerical collapse ($a_{{\cal D}}\ll 1$) at foliation time $t_{\mathrm{coll}}$. The initial solution $\{a_{{\cal D}}(t),\dot{a}_{\cal D}(t):t<t_{\mathrm{coll}}\}$ is used as an approximation that is improved iteratively. We assume virialisation – we set

[TABLE]

where $18\pi^{2}$ is the usual EdS stable clustering (Peebles 1980) virialisation overdensity (Lacey & Cole 1993, eq. (A16)). For completeness, the inhomog package provides the corresponding value for the flat-$\Lambda$ case if that is chosen, using the fitting formula of Bryan & Norman (1998, (6)) to Eke et al. (1996)’s calculation.

As motivated below in Sect. 4, we do not introduce any behaviour in other domains that compensates for virialisation events. Indeed, it is not obvious how such compensating volume evolution could reasonably be approximated, apart from imposing the reference model expansion rate and abandoning the aim of estimating structure-generated effective expansion.

3.6 Regional versus global averaging

As in Räsänen (2006, 3.1) and Rácz et al. (2017), we cubically average the domain-wise scale factors,

[TABLE]

where the volumes of collapsed domains are bound below by their virialisation volumes at the time of collapse, as defined in Eq. (43). Although this represents volume averaging at the global level, virialisation is not represented in ${\cal Q}_{\cal D}$, so unless ${\cal D}$ is chosen large enough for no virialisation events to occur, Eq. (14) is expected to fail at the global level, since the single fluid stream assumption fails at virialisation. Thus, $a_{\mathrm{eff}}(t)=a_{\mathrm{EdS}}(t)$ is not expected in the presence of virialisation. See Sect. 6.1 for more discussion of this point.

3.7 ramses front end and $N$-body measurement of ${\cal Q}_{\cal D}$

The two main classes of faster simulation techniques are tree codes and particle-mesh (PM) based methods (e.g., Bagla 2005; Dehnen & Read 2011). In this work, we provide a patch to the adaptive mesh refinement (AMR) ramses code666https://bitbucket.org/rteyssie/ramses (Teyssier 2002; Guillet & Teyssier 2011), available under CeCILL (GNU GPL compatible) as Fortran 90 code, that we tentatively refer to as ramses-scalav777https://bitbucket.org/broukema/ramses-scalav This serves as a front end to read in the initial conditions files (Sect. 3.3) and call the dtfe and inhomog libraries to estimate the per-domain initial invariants $\left\langle{\mathrm{I}}\right\rangle_{{\cal D}_{\mathbf{i}}}$, $\left\langle{\mathrm{II}}\right\rangle_{{\cal D}_{\mathbf{i}}}$, and $\left\langle{\mathrm{III}}\right\rangle_{{\cal D}_{\mathbf{i}}}$, and to carry out VQZA evolution and volume averaging as described in Sections 3.5 and 3.6.

In an initial exploration of more numerical approaches, the ramses-scalav front end also allows ${\cal Q}_{\cal D}$ to be calculated numerically from the invariants at each major time step using Eq. (16) instead of using the QZA analytical calculation of (LABEL:e-resultQ2) and (31), which is based on the initial values of the invariants and the growth rate of the reference model. We then integrate Eq. (25), again rewriting it in the form of Eq. (40). In this case, we integrate using the classical (non-adaptive) Runge–Kutta (fourth-order) algorithm. As in Sect. 3.6, we average the domain-level scale factors $a_{{\cal D}}$ to obtain $a_{\mathrm{eff}}$ using Eq. (44). For any given domain ${\cal D}$, we store the list of particles initially present in ${\cal D}$, and at time $t$ we calculate $\left\langle{\mathrm{I}}\right\rangle_{\cal D}$ and $\left\langle{\mathrm{II}}\right\rangle_{\cal D}$ by weighting ${\mathrm{I}}$ and ${\mathrm{II}}$ by the volume element $\mathrm{d}\mu_{i}$ of the $i$-th particle. All parameters ${{\mathrm{I}}}_{i}$, ${{\mathrm{II}}}_{i}$ and $\mathrm{d}\mu_{i}$ are calculated as (Euclidean) linear interpolations of the ${{\mathrm{I}}}$, ${{\mathrm{II}}}$ and $1/\rho$ values of the eight neighbouring dtfe cells surrounding the particle. The global normalisation between $\mathrm{d}\mu_{i}$ and $1/\rho_{i}$ is arbitrary.

Since we do not change the list of particles in a domain between time steps, the domain shape, which we can conceptually imagine as the union of Voronoi tessellation tetrahedra surrounding the particles, is unlikely to remain cubical. It is likely to become irregular and disconnected. The domain at late times will contain holes that exclude particles that have entered from neighbouring domains, and will include small islands of space around particles that have escaped the main body of the domain. In principle, this is not a problem, since Eq. (25) is based on Riemannian (or Lebesgue in the case of a global T3 domain) integration, which is additive.

By default, a ramses cosmological simulation runs like a typical $N$-body simulation, with the reference model expansion inserted at each time step from a pre-calculated model, decoupled from structure formation. For future work, our ramses-scalav front end makes it easy to adopt a more consistent approach, that is, to use the cubically averaged $a_{\mathrm{eff}}(t)$ [Eq. (44)] when calculating the Newtonian gravitational potential that is used for inferring accelerations of particles. In the present work, we limit calculations to those using the reference model EdS scale factor evolution, in order to compare with the VQZA approach.

4 VQZA super-EdS expansion: biscale partition

Räsänen (2006) proposed that the volume average of expanding and contracting domains defined against an FLRW reference model may evolve differently from, and in particular, at later times, evolve faster than the background. As stated in Sect. 2, prior to virialisation, this claim fails in a strictly Newtonian T3 case for contiguous space: the kinematical backreaction ${\cal Q}_{\cal D}$ in the expanding and contracting domains is a non-linear effect that compensates for what otherwise appears to be dynamical asymmetry in non-linear evolution of initially small density perturbations. This was established in (Buchert & Ehlers 1997).

4.1 Newtonian conservation of EdS expansion prior to virialisation

We illustrate this here by dividing a T3 volume of an EdS model into exactly two domains ${\cal D}^{+}$ and ${\cal D}^{-}$ of initially equal volume: T${}^{3}={\cal D}^{+}\cup{\cal D}^{-}$; ${\cal D}^{+}\cap{\cal D}^{-}=\emptyset$; ${{V}_{{\cal D}^{-}_{\mathbf{i}}}}={{V}_{{\cal D}^{+}_{\mathbf{i}}}}$. We set ${\cal D}^{-}$ to be the expanding domain and ${\cal D}^{+}$ the contracting domain. Thanks to Stokes’ theorem, the global (Newtonian) averages of the three invariants of the peculiar velocity gradient (approximating what in the relativistic case is the peculiar expansion tensor; see III.B.1 in Buchert et al. 2013), where the invariants are given in Sect. 3.4.2, are each zero on this compact spatial 3-manifold. With volume partitioning [see Wiegand & Buchert 2010, eqs (16), (17)] and our choice to give the two domains initially equal volumes, we have per-domain average initial invariants related by

[TABLE]

Using Eqs (LABEL:e-resultQ2), (31), (34), and (35), the QZA on this biscale volume split should provide approximate confirmation that the EdS scale factor evolution is conserved during a Newtonian pre-virialisation period in cases in which the QZA is approximate, and exact (within numerical accuracy) confirmation in cases in which it is exact.

We consider several subcases – a class that includes planar collapse:

[TABLE]

a class that includes spherically symmetric expansion or collapse:

[TABLE]

and cases where either $\left\langle{\mathrm{II}}\right\rangle_{{\cal D}^{-}_{\mathbf{i}}}$ or $\left\langle{\mathrm{III}}\right\rangle_{{\cal D}^{-}_{\mathbf{i}}}$ is set to zero and the non-zero invariants are set to values that give terms in Eqs (LABEL:e-resultQ2) and (31) of roughly similar orders of magnitude. Specific values are indicated in the caption of Fig. 1.

The planar collapse class of solutions (Buchert et al. 2000, III.D) constitute an exact solution family, which is clearly supported numerically in Fig. 1 prior to the collapse of the collapsing domain at about 6.2 Gyr. Prior to 6.2 Gyr, the global super-EdS scale factor ratio is unity to within a precision of $|1-a_{\mathrm{eff}}(t)/a_{\mathrm{EdS}}|<2.3\times 10^{-5}$, where $a_{\mathrm{eff}}$ is defined in Eq. (44). In other words, prior to virialisation, the QZA in the planar class of solutions is fully consistent with the Newtonian $T^{3}$ cosmology that is normally used to calculate structure formation in the flat FLRW models (EdS and $\Lambda$CDM).

The spherical symmetric class of solutions is also exact (Buchert et al. 2000, III.E), with $\gamma_{1}=\gamma_{2}=\gamma_{3}=0$. However, if Eq. (47) is satisfied so that the expanding domain ${\cal D}^{-}$ is a member of this class, then Eq. (45) prevents the contracting domain from satisfying Eq. (47) (with ${\cal D}^{-}$ replaced by ${\cal D}^{+}$), since the squared second invariant is necessarily positive in both cases. Thus, a small deviation from global EdS evolution is visible in Fig. 1: the global super-EdS scale factor ratio is close to unity but drops a few percent below unity just prior to the collapse at 4.3 Gyr. The other two cases show similar levels of deviation from an exact global EdS expansion just prior to collapse of the collapsing domain.

4.2 After virialisation: the VQZA

In Fig. 1, the global scale factor evolution follows from assuming the stable clustering hypothesis [Eq. (43)] for the collapsed domain and by assuming that the expanding domain continues its expansion as given in the QZA in Eqs (LABEL:e-resultQ2), (31), (34), and (35). This conservative assumption motivates the following definition.

**Definition **\thetheorem

Virialisation ${\cal Q}_{\cal D}$ Zel’dovich approximation (VQZA):

• (i)

Divide a large volume to be studied into contiguous domains ${\cal D}$;

• (ii)

evolve the scale factor $a_{{\cal D}}(t)$ in each domain according to the Raychaudhuri equation [Eqs (25), (26)], where the kinematical backreaction ${\cal Q}_{\cal D}$ is approximated by the ${\cal Q}_{\cal D}$ Zel’dovich approximation (QZA) [Eqs (11), (LABEL:e-resultQ2), (31)], unless gravitational collapse and virialisation (iii) occur at some time $t_{\mathrm{coll}}$;

• (iii)

if virialisation [Sect. 3.5.6] occurs in a domain ${\cal D}$, QZA evolution ceases at $t=t_{\mathrm{coll}}$ and the stable clustering approximation [Eq. (43) in the EdS case] for $a_{{\cal D}}(t)$ is adopted for $t\geq t_{\mathrm{coll}}$ for the domain ${\cal D}$;

• (iv)

the global scale factor evolution is estimated by cubical averaging of the per-domain scale factors $a_{{\cal D}}(t)$ [Eq. (44)].

Is this model Newtonian or relativistic? Our aim is to model spatial expansion in a close approximation to general-relativistic behaviour. As discussed above, a fully Newtonian cosmology is difficult to define. The Buchert & Ehlers (1997) fluid model justification of T3 Newtonian cosmology does not apply past shell-crossing. When one domain collapses and the others continue expanding, should there be a sudden compensation in the behaviour of the expanding domain(s) in response to the virialisation of the collapsed domain? We hypothesise that this particular feedback effect would not occur in a relativistic model. In other words, we hypothesise silent virialisation: that the details of collapse have a negligible effect at large distances from the collapsing domain (Matarrese et al. 1994b, a; Ellis & Tsagas 2002; Bolejko 2017b). In this sense, the VQZA is a relativistic approximation. It is not exactly relativistic. For example, the sum of rest masses is conserved via the constants $A$ used in per domain Raychaudhuri integration [Eqs. (25), (26)]; peculiar velocities of the order of $10^{-2.5}$ would imply that the system mass (in the Minkowski spacetime sense) should be about $10^{-5}$ higher than the sum of the rest masses. In work closely analogous to the present paper, Bolejko (2017b, c) uses a silent universe family of cosmological solutions of the Einstein equation based on four scalar quantities (density, expansion rate, shear and Weyl curvature). That model numerically evolves an ensemble of worldlines starting with $\Lambda$CDM initial conditions, making similar assumptions to those of the VQZA: stable clustering and silent virialisation. The emergence of strong negative mean curvature is inferred in those works.

The dynamical difference between the VQZA and the standard approach is also shown in Figure 2 in the biscale case. This figure shows the evolution of the expanding domain ${\cal D}^{-}$ in the VQZA ($+$ symbols) compared with evolution that follows the standard (implicit) dynamical assumption made about expanding domains in cosmological $N$-body simulations and typical analytical calculations (solid curves). The latter is defined by

[TABLE]

where the factor of two represents division of the spatial section into two domains. The expanding domain has smooth VQZA behaviour independently of the collapse event. In contrast to the smooth volume evolution under the VQZA, the standard model implicitly requires the expanding domain(s) to be suddenly pulled more strongly by the virialised object once collapse and virialisation have finished, in order to preserve a globally EdS expansion rate: $\ddot{a}_{{\cal D}^{-}}^{\mathrm{Newt}}$ has to suddenly become less positive or more negative when $\dot{a}_{{\cal D}^{+}}$ very rapidly switches from a strong negative value to zero.

Interpreted relativistically, the expanding domain prior to collapse corresponds to negatively curved space that is below the critical density. This is shown in Fig. 3 for the same biscale examples. The curvature functional is given in Eq. (42) in the VQZA case. In the EdS constraint case, we again interpret $\Omega_{{\cal R}}^{{\cal D}}$ relativistically (an EdS model is intended as a relativistic, not Newtonian, model), using the partitioning formula, eq. (10) of Wiegand & Buchert (2010) applied to a globally flat spatial section,

[TABLE]

where the ${H_{{\cal D}^{+}}^{2}}/{H{{}_{\mathrm{FL}}}^{2}}$ factor converts between globally normalised and per-domain normalised definitions of $\Omega$. Following collapse, the stable clustering assumption gives $H_{{\cal D}^{+}}:=0$, and both domains are interpreted as flat. Thus, in the EdS constraint case, the expanding domain in the biscale case needs to suddenly switch from smooth $\Omega_{{\cal R}}^{{\cal D}}$ growth to $\Omega_{{\cal R}}^{{\cal D}}\approx 0$. Interpreted from a Newtonian point of view, this sudden switch is allowed; here, we hypothesise that it is relativistically inaccurate.

Ironically, it is the standard model that assumes a form of feedback (backreaction) in this situation; by not adopting this virialisation-induced feedback, the VQZA is more conservative than the standard approach.

5 Results

In studying the question of whether dark energy can be replaced by a non-linear structure formation scale at a first level of approximation, we treat the $\Lambda$CDM model as an observational proxy, as in Roukema et al. (2017). In principle, given the microlensed Galactic Bulge star likely upper age limit of $T\approx 14.7^{+0.3}_{-0.7}$ Gyr (Bensby et al. 2013; Roukema et al. 2017) and EdS CMB fits indicating an age of the Universe slightly higher than this (Blanchard et al. 2003; Hunt & Sarkar 2007; Nadathur & Sarkar 2011), it could reasonably argued that we should consider the structure formation scale that yields an effective scale factor $a_{\mathrm{eff}}$ that reaches unity just a little earlier than $t=15$ Gyr, rather than the $\Lambda$CDM current age of the Universe estimate of $t_{0}=13.80$ Gyr. Nevertheless, it is premature to make detailed observational tests using the VQZA approach, so in this work we consider the $\Lambda$CDM proxy estimate of the age of the Universe as a reasonable first-order approximation.

5.1 Uncertainties in peculiar velocity gradient invariants ${\mathrm{I}}$ and ${\mathrm{II}}$

Measurement of ${\cal Q}_{\cal D}$ is crucial to calculating volume evolution, so we first consider numerical uncertainties. Application of the error estimation methods described in Sections 3.4.1 and 3.4.2 to the test functions [Eq. (13)] are illustrated in Figs 4–6, using the dtfe version indicated in Table 1.

Figure 4 shows values mostly greater than unity. In other words, given that we know the true test-function velocity field $\mathbf{u}_{1}$ [Eq. (13)], numerical dtfe estimates of the first and second tensor invariants of the peculiar velocity gradient tensor typically have errors of smaller amplitude than the signal. This is consistent with the claims in the literature regarding dtfe’s accuracy (Schaap & van de Weygaert 2000; van de Weygaert & Schaap 2009; Cautun & van de Weygaert 2011).

The $n_{x}=64,n_{\mathrm{DTFE}}=256$ points for both ${\mathrm{I}}$ and ${\mathrm{II}}$ (upper and lower plots, respectively) are for dtfe sampling at a higher resolution than the $N=128^{3}$ particle resolution. That is, one full sinusoidal wavelength of $\mathbf{u}_{1}$ in one direction is sampled by about two particles, and dtfe estimates the velocity gradient at four grid points along this direction. To obtain any reasonable velocity gradient in this situation would obviously push the extremes of the available information content, so obtaining $S/N\sim 1$, as shown in both panels of the figure, is clearly reasonable.

However, the error estimates in Fig. 4 are based on prior knowledge of the true velocity field, which, in general, is not the case. The S1-based error estimation method introduced in Sect. 3.4.1 bypasses this requirement. Figure 5 shows that this new method appears to be viable, in the sense that the ratios of rms measured noise to rms predicted noise are within at most two orders of magnitude of unity. The assumptions of statistical independence and Gaussianity used in the derivation of Eq. (8) are unlikely to be fully satisfied, so better agreement between modelled and measured noise would be surprising. Structures sized just twice that of the particle resolution ($n_{x}=64=N/2$, asterisks) tend to give the worst (highest) ratios. Most of the ratios worsen when the dtfe resolution is higher than the particle resolution ($n_{\mathrm{DTFE}}=2N$). These behaviours are reasonable.

The measured-to-predicted noise ratios for the second invariant ${\mathrm{II}}(\mathbf{u}_{2})$ (Fig. 5, lower panel, lower curves with $+$ and $\times$ symbols) mostly are lower than unity. That is, the measured errors are lower than that given by Eq. (8). This is presumably related to the fact that ${\mathrm{II}}(\mathbf{u}_{2})=0$, and possibly also to symmetries in sampling this function which cancel each other very effectively, so that the assumption of statistically independent distributions overestimates the noise level.

Global (T3) application of Stokes’ theorem gives consistent error propagation results to those for the S1-based estimates, to within a few orders of magnitude. Figure 6 shows that $\left\langle{\mathrm{I}}\right\rangle_{\mathrm{box}}$ (global first invariant) errors are bounded above by the S1-based errors with Gaussian error propagation. In contrast, $\left\langle{\mathrm{II}}\right\rangle_{\mathrm{box}}$ errors are only bounded above by the S1-based errors for large structures ($n_{x}=1$ or $8$). Structures whose wavelength is twice that of the particle resolution ($n_{x}=64=N/2$, asterisks) numerically fail the analytically required $\left\langle{\mathrm{II}}\right\rangle_{\mathrm{box}}=0$ test at many $n_{\mathrm{DTFE}}$ resolutions. The equivalent plot for ${\cal Q}_{\mathrm{box}}$ is not displayed, because the low amplitude of $|\left\langle{\mathrm{I}}\right\rangle_{\mathrm{box}}|$ estimates and Eq. (16) make the $\left\langle{\mathrm{I}}\right\rangle_{\mathrm{box}}^{2}$ contribution negligible, giving a plot visually identical to that for $\left\langle{\mathrm{II}}\right\rangle_{\mathrm{box}}$, but raised by a factor of two.

5.2 VQZA super-EdS expansion

Figure 7 shows effective scale factors for VQZA simulations with $L_{\mathrm{box}}=8L_{{\cal D}}$ and an mpgrafic grid of $N=256^{3}$ particles, with $L_{{\cal D}}$ increasing within any of the three panels as labelled, and $L_{{\cal D}}/L_{\mathrm{DTFE}}$ increasing from two in the top panel to eight in the bottom panel. (See Table 1 for software versions.) The scale factor growth is very close to EdS for L_{{\cal D}}=8\,\mbox{Mpc/h_{\mathrm{eff}}} and L_{{\cal D}}=16\,\mbox{Mpc/h_{\mathrm{eff}}}. This is reasonable, because typical overdensities on these scales do not have time to collapse by the present.

On smaller averaging scales $L_{{\cal D}}$, silent virialisation (the VQZA, Definition 1, Sect. 4.2) starts to play a role in volume expansion. This can be seen in Figure 7, where the scale factor growth is stronger as $L_{{\cal D}}$ decreases. The $\Lambda$CDM proxy scale factor evolution ($\times$ symbols) mostly lies between the $L_{{\cal D}}=2$ and $4\,$Mpc/$h_{\mathrm{eff}}$ curves. In other words, these curves provide the required super-EdS expansion required to reach a unity effective scale factor $a_{\mathrm{eff}}(t)$ by $t=13.8$ Gyr.

Figure 7 also indicates that the preferred $L_{{\cal D}}$ scale for providing sufficient super-EdS expansion is only weakly sensitive to the choice of $L_{\mathrm{box}}$, $L_{\mathrm{DTFE}}$, and $L_{N}$. The spacing between the curves in any individual panel is not uniform: this is very likely a sign of cosmic variance: the $L_{\mathrm{box}}$ sizes here are small. The box size $L_{\mathrm{box}}$ can be increased further, but at the expense of leaving only small ratios $L_{{\cal D}}/L_{\mathrm{DTFE}}$ and $L_{\mathrm{DTFE}}/L_{N}$. Figure 8 shows that 2\,{\mbox{Mpc/h_{\mathrm{eff}}}}\,\lower 2.58334pt\hbox{\buildrel<\over{\sim}}\,L_{{\cal D}}\,\lower 2.58334pt\hbox{\buildrel<\over{\sim}}\,4\,{\mbox{Mpc/h_{\mathrm{eff}}}} is still favoured in this case.

By running a large number of these simulations, we can estimate $L_{13.8}$, the value of $L_{{\cal D}}$ which gives $a_{\mathrm{eff}}(13.8\leavevmode\nobreak\ \mathrm{Gyr})=1$. This is done here as follows. For a set of five simulations with different $L_{{\cal D}}$ values, as in one of the panels of Fig. 7, estimate $t_{a_{\mathrm{eff}}=1}$ by a spline of $t$ against $a_{\mathrm{eff}}(t)$ for each simulation. Least-squares fit a quadratic to $t_{a_{\mathrm{eff}}=1}(L_{{\cal D}})$ and solve $t_{a_{\mathrm{eff}}=1}(L_{{\cal D}})=13.8$ Gyr to obtain $L_{{\cal D}}$. We find that for $L_{\mathrm{box}}=8L_{{\cal D}}=4L_{\mathrm{DTFE}}=8L_{N}$ (as in Fig. 7, middle panel), 100 sets of 5 simulations yield an averaging scale of $L_{13.8}=2.51\pm 0.22$ Mpc/$h_{\mathrm{eff}}$. This is the scale at which $a_{\mathrm{eff}}(t)$ reaches a unity scale factor at 13.8 Gyr, which is 16% above $a_{\mathrm{EdS}}(13.8\leavevmode\nobreak\ \mathrm{Gyr})$. The error estimate is the standard deviation of the 100 estimates. This error estimate represents the effect of cosmic variance for 4\,\mbox{Mpc/h_{\mathrm{eff}}}\leq L_{\mathrm{box}}\leq 64\,\mbox{Mpc/h_{\mathrm{eff}}}.

Corresponding values for smaller and bigger box sizes, $L_{\mathrm{box}}=8L_{{\cal D}}=2L_{\mathrm{DTFE}}=16L_{N}$ (Fig. 7, top panel) and $L_{\mathrm{box}}=8L_{{\cal D}}=8L_{\mathrm{DTFE}}=4L_{N}$ (Fig. 7, bottom panel), are $L_{13.8}=2.14\pm 0.26$ Mpc/$h_{\mathrm{eff}}$ and $L_{13.8}=2.59\pm 0.25$ Mpc/$h_{\mathrm{eff}}$, respectively, for 100 sets of 5 simulations of $N=256^{3}$ particles in each case. Thus, the uncertainty in changing the box size and resolution parameters is of about the same magnitude as the uncertainty associated with cosmic variance, and about ten times the standard error in the mean of $L_{13.8}$ for any fixed set of size and resolution parameters. We summarise these by writing the (dominating) systematic error and taking the median of the three sets of calculations ($L_{\mathrm{box}}=8L_{{\cal D}}=4L_{\mathrm{DTFE}}=8L_{N}$, as in Fig. 7, middle panel, for one quintuple of realisations): $L_{13.8}=2.5^{+0.1}_{-0.4}$ Mpc/$h_{\mathrm{eff}}$.

We also carried out 100 sets of five simulations on still bigger box scales (as in Fig. 8), at the risk of greater numerical error due to the small factors between the shorter length scales ($L_{{\cal D}}/L_{\mathrm{DTFE}}=L_{\mathrm{DTFE}}/L_{N}=2$). This set of simulations gives a much more precise estimate: $L_{13.8}=2.56\pm 0.01$ Mpc/$h_{\mathrm{eff}}$, most likely because of the reduced cosmic variance. However, since there is a greater risk of numerical error, this estimate is not likely to be more accurate (the resolution-related systematic error is likely to be higher). Thus we retain $L_{13.8}=2.5^{+0.1}_{-0.4}$ Mpc/$h_{\mathrm{eff}}$ as a conservative summary of our main quantitative result.

Räsänen (2006)’s early model was only a toy model, without a proposal for a specific characteristic length scale $L_{{\cal D}}$. Using the FLRW critical density $\rho_{\mathrm{crit}}=1.52\times 10^{11}M_{\odot}/\mathrm{Mpc}^{3}$, Rácz et al. (2017)’s preferred mass scale for matching the observational distance-modulus–redshift relation, $2.03\times 10^{11}M_{\sun}$, corresponds to $1.10$ Mpc/$h_{\mathrm{eff}}$, which is a little below half our scale. We do not expect our main result to closely match that of Rácz et al. (2017), since our methods differ. If we artificially suppress kinematical backreaction in order to compare with Rácz et al. (2017), that is, if we set ${\cal Q}_{{\cal D}}(t)=0$ in all domains ${\cal D}$ at all times, then the super-EdS expansion rates should increase for any fixed $L_{{\cal D}}$. Figure 9 shows that this is indeed the case. Since kinematical backreaction is ignored, stronger super-EdS global volume evolution results, giving L_{13.8}({\cal Q}_{{\cal D}}:=0)\approx 8\leavevmode\nobreak\ \mbox{Mpc/h_{\mathrm{eff}}}, that is, an overestimate by a factor of about three. This is a factor of about seven greater than Rácz et al. (2017)’s estimate: the latter’s model does not quite match our zero-${\cal Q}_{\cal D}$ simulations.

Figure 10 helps to show the role of kinematical backreaction in our calculations, where the volume-weighted mean and median values for the uncollapsed domains are defined

[TABLE]

where ${{\cal D}}_{k}$ is the $k$-th domain among the sorted $i$-th scale factors of uncollapsed domains satisfying

[TABLE]

and Eqs (LABEL:e-resultQ2), (31), (LABEL:e-resultR2), (33) are used to evaluate ${\cal Q}$ and ${\cal R}$ (ties in the median, if they occur, are upper weighted). The $\Omega_{{\cal Q}}^{{\cal D}}$ volume-weighted mean and median values for uncollapsed domains are positive, corresponding to negative ${\cal Q}_{\cal D}$, which tends to cause earlier collapse of overdensities and slower expansion of voids. Thus, if kinematical backreaction were artificially set to zero in all domains, as appears to be the case in Rácz et al. (2017) and is shown here in Fig. 9, then it would be reasonable for the faster expansion of voids to contribute to much stronger scale factor growth. In that case, larger scale (lower amplitude) initial overdensities would be sufficient to obtain $a_{\mathrm{eff}}(13.8\leavevmode\nobreak\ {\mathrm{Gyr}})=1$, as can be seen in Fig. 9.

The values of the mean and median $\Omega_{{\cal Q}}^{{\cal D}}$ values shown in Fig. 10 only sample uncollapsed domains, and are normalised globally [Wiegand & Buchert 2010, eq. (10)], rather than per-domain [Eq. (42)]. This implies that prior to virialisation, the mean should be constrained to be close to zero by the Stokes’ theorem T3 constraint on the global kinematical backreaction ${\cal Q}_{T^{3}}$. This does appear to be the case in the upper panel of Fig. 10 at very early times, but the first virialisation events happen quickly, so $\left\langle\Omega_{{\cal Q}}^{{\cal D}}\right\rangle_{\mathrm{uncoll}}$ quickly grows away from zero.

The spikiness of Fig. 10 (the vertical axis is logarithmic) is physically realistic: it emphasises the discontinuous nature of virialisation. Each time a domain collapses, it no longer contributes to $\left\langle\Omega_{{\cal Q}}^{{\cal D}}\right\rangle_{\mathrm{uncoll}}$ as defined for this figure. Forcing a smooth $a_{\mathrm{EdS}}(t)$ (or $a_{\Lambda\mathrm{CDM}}(t)$) global volume evolution artificially hides this spikiness. Unsurprisingly, the volume-weighted medians $\mu(\Omega_{{\cal Q}}^{{\cal D}})_{\mathrm{uncoll}}$ evolve much more smoothly, since the collapse of a domain shifts the median within the central part of the $\Omega_{{\cal Q}}^{{\cal D}}$ distribution, which is denser than the tails. For $L_{{\cal D}}=2\,$Mpc/$h_{\mathrm{eff}}$, these medians peak rapidly at about $0.03$ and then gradually decay to the present.

As expected by Buchert (2005; see Sect. 1 above for more literature), the super-EdS $a_{\mathrm{eff}}(t)$ is accompanied by emerging negative average curvature, as shown in Fig. 11. At the $L_{{\cal D}}=2\,$Mpc/$h_{\mathrm{eff}}$ scale (middle curve at early times), the volume-weighted mean $\Omega_{{\cal R}}^{{\cal D}}$ for uncollapsed domains grows from zero to about $1$ by about half the age of the Universe and then grows more slowly to about $1.5$ by the present. The volume-weighted median $\Omega_{{\cal R}}^{{\cal D}}$ also grows quickly to about $1.5$ by $t\sim 4$ Gyr, the same epoch, and then grows slowly to about $2$ by the present. New observational tests will be needed to constrain recent emerging curvature such as that modelled here, since constraints that assume comoving non-evolving uniform curvature do not constrain evolving non-uniform mean curvature. Moreover, in this work, no distinction is made between volume-average (“bare”) $\Omega$ parameters and those likely to be observed from our Galaxy (“dressed”; Buchert & Carfora 2002, 2003), so more work will be needed to calculate parameters such as the detailed estimations for the Timescape model given in Nazer & Wiltshire (2015, Table II).

5.3 Dependence of super-EdS expansion on virialisation fraction $f_{\mathrm{vir}}$

Roukema et al. (2013) argued that virialisation and non-rigidity of voids constitute a key element to an emerging negative average curvature alternative to dark energy. The VQZA provides an ab initio model of this argument, so it should be possible to verify the relationship between super-EdS expansion and virialisation. The virialisation fraction (fraction of mass in virialised objects) can be estimated by counting the number $N_{\mathrm{coll}}$ of domains ${\cal D}$ that have gravitationally collapsed at any given time and normalising:

[TABLE]

Figure 12 shows that at any fixed time, VQZA simulations run from independent realisations of initial conditions yield super-EdS expansion $a_{\mathrm{eff}}/a_{\mathrm{EdS}}$ that correlate strongly with $f_{\mathrm{vir}}$. Spearman’s rank correlation $\rho$ tests overwhelmingly confirm this, giving a probability $p=3\times 10^{-9},10^{-9},$ and $10^{-6}$, of no correlation between $f_{\mathrm{vir}}$ and $a_{\mathrm{eff}}/a_{\mathrm{EdS}}$ at the early, middle, and late times (as listed in Fig. 12), respectively.

Thus, unsurprisingly, silent virialisation leads to a strong correlation between the degree of virialisation $f_{\mathrm{vir}}(t)$ that results from a given realisation of initial conditions at a given time $t$ and the amount of super-EdS expansion $a_{\mathrm{eff}}(t)/a_{\mathrm{EdS}}(t)$ at that same epoch. Enforcement of Newtonian global expansion, as in standard cosmological $N$-body simulations, would destroy this correlation, since in that case $a_{\mathrm{eff}}(t)/a_{\mathrm{EdS}}(t):=1$ by definition.

5.4 $N$-body measurement of ${\cal Q}_{\cal D}$

What happens when ${\cal Q}_{\cal D}$ is calculated from the Lagrangian domains in the evolving $N$-body simulation, as detailed in Sect. 3.7 above, where numerical noise plays a stronger role? Figures 13 and 14 show the equivalent of Figs 7 (middle) and 9 (top), respectively. (See Table 1 for software versions.)

Figures 9 (top) and 14, in which ${\cal Q}_{\cal D}$ is artificially set to zero, appear, by inspection, to be almost identical. This close correspondence provides a check on the numerical integration and coordination between software modules, since the former integration is carried out by calls from the inhomog library to the GSL routine gsl_odeiv_step_rkf45, while the latter is carried out by the ramses-scalav routine update_av_scalars_global_module with a hardwired version of the classical (non-adaptive) Runge–Kutta algorithm. The initial conditions were not set to be identical between these two panels, so the small differences indicate both cosmic variance and numerical error.

For more realistic calculations, when kinematical backreaction is not set to zero, Fig. 13 shows more complex behaviour than that in Fig. 7 (middle). The sequence of bigger scales $L_{{\cal D}}=2,4,8$ Mpc$/h_{\mathrm{eff}}$ has slightly stronger $a_{\mathrm{eff}}(t)$ than for the VQZA case. In contrast, for the sequence of smaller scales ($L_{{\cal D}}=2,1,0.5$ Mpc$/h_{\mathrm{eff}}$) the $a_{\mathrm{eff}}(t)$ growth is successively weaker. Cosmic variance is strong on these scales, so seeking an explanation from a small number of small simulations is not strongly justified: many competing numerical effects such as close-encounter scattering and softening may be involved. More importantly for the present work, a value of 1\,\lower 2.58334pt\hbox{\buildrel<\over{\sim}}\,L_{13.8}/\mbox{Mpc/h_{\mathrm{eff}}}\,\lower 2.58334pt\hbox{\buildrel<\over{\sim}}\,4 indicates consistency with the VQZA calculations. Future work on higher numbers of higher $N$ simulations will enable a statistically more significant comparison between the analytical and numerical models of ${\cal Q}_{\cal D}$ and their effects on $a_{\mathrm{eff}}(t)$.

6 Discussion

6.1 Should global backreaction from T3 initial conditions be zero?

The scalar averaging equations, given in Sect. 3.5.1 together with Eq. (44), can be applied on any domain size, including the global volume, which in this work is T3, provided that the assumptions underlying the equations are satisfied. For flat space, given the definition of ${\cal Q}_{\cal D}$ in Eq. (16) and the flat space expression of the peculiar velocity gradient invariants as exact divergences in Eq. (11), it would appear that ${\cal Q}_{\mathrm{T}^{3}}$ must be zero. Thus, Eq. (14) should yield $a_{\mathrm{eff}}(t)\equiv a_{\mathrm{EdS}}(t)$, no matter what scale of ${\cal D}$ is chosen, provided that the evolution of ${\cal Q}_{\cal D}$ is calculated exactly. More realistically, in the present case ${\cal Q}_{\cal D}$ is calculated using the QZA [Eqs (LABEL:e-resultQ2), (31)], so $a_{\mathrm{eff}}(t)\approx a_{\mathrm{EdS}}(t)$ should still be a good approximation, to the degree that Eqs (LABEL:e-resultQ2), (31) remain accurate at $t$, modulo Poisson noise contributed by combining many domains. Indeed, Fig. 1 shows that prior to virialisation, this approximation holds reasonably well in illustrative cases.

However, Eq. (14) and the definition of ${\cal Q}_{\cal D}$ normally used in the literature are derived for a pressureless perfect fluid, that is, ‘dust’, with no allowance for shell crossings or avoidance of collapse to a singularity. This is not a problem for studying non-linear gravitational collapse times of overdense objects, in which case these equations yield realistic estimates in a more elegant way than the usual alternatives (Ostrowski et al. 2016; Ostrowski 2016; Ostrowski et al. 2017, 2018). On the other hand, Eq. (14) (in its form stated above) is no longer valid once the final stages of gravitational collapse and virialisation have taken place. As $a_{{\cal D}}\rightarrow 0^{+}$ in Eq. (14), ${\cal Q}_{\cal D}$ via Eqs (LABEL:e-resultQ2) and (31) is not able to cancel self-gravity, so a literal interpretation of Eq. (14) would imply that $\dot{a}_{{\cal D}}\rightarrow-\infty$.

In this sense, the global volume-averaging expressed in Eq. (44) is used in this work as an extension beyond the usual set of scalar averaging equations; the metric volumes are summed, but the choice of the ${\cal D}$ scale affects the result via virialisation.

Obviously, if the aim of calculations of this sort were to study EdS cosmology and to substitute a perfect fluid by a set of collisionless particles in which shell-crossings were allowed, then it would have to be possible to modify Eq. (14) by adding a virialisation term to the kinematical backreaction in order to prevent collapse to a singularity. For example, this could take the form of the dynamical backreaction proposed in eqs (13a), (13d) of Buchert (2001) (see also Yodzis 1974; Buchert 2017, footnote 7). To retain a global EdS expansion history, the new term would not only have to prevent collapsing domains from forming singularities, but would also have to correspondingly slow down the expansion rate of expanding domains, in order to retain a net EdS expansion rate during the post-virialisation phase.

The main aim of this work is a step towards relativistic cosmology rather than testing the consistency of Newtonian cosmology. The latter requires instant reaction at a distance, but the former is constrained by causality, making instantaneous reaction to virialisation unphysical. Thus, the approach defined in Definition 1 (Sect. 4.2) is retained for the present work. In other words, here, the QZA formulae for expanding regions are not modified nor supplemented in reaction to virialisation, and the resulting global volume evolution is not constrained to be Newtonian.

Moreover, since stable clustering is adopted for virialised domains, if we consider their volume-weighted curvature contributions to be negligible (either because the volumes are tiny, the curvature is weak, or both), then the negative curvature of the expanding domains will tend to yield an overall negative mean curvature, invalidating the derivation of the divergence expressions in Eq. (11) at late times, since a flat manifold is no longer a good approximation.

6.2 How can virialisation provide an order-unity dark-energy expansion?

How is it possible that stable clustering, which is typically thought of as a Newtonian statistical gravitational effect, can provide a reasonable relativistic approximation in this context? The reasoning above can be reworded as follows. Prior to virialisation, the QZA formulae provide a good cancellation between expansion and collapse compared to EdS evolution of the scale factor. At a virialisation event, the spatially contracting contribution of a virialised domain to the general cosmological volume evolution is cancelled, either by the stable clustering assumption or (in future work) using more detailed analytical models (Buchert 2001; Bolejko & Lasky 2008; Bolejko & Ferreira 2012; Buchert 2017). To better than an order of magnitude, about half the mass of any large volume today is located in virialised dark matter haloes. Thus, among the ensemble of domains on the nonlinear, few-megaparsec scale, the domains that contribute about half the total mass have their gravitational roles cancelled, in the sense that the averaged Raychaudhuri equation, Eq. (14), is replaced by

[TABLE]

for those domains. Under an EdS constraint, expanding regions are forced to slow their expansion down in response to this virialisation. Relativistically interpreted, this would mean that the expanding, negatively curved regions are forced to slow down and flatten. Without the EdS constraint, that is, without modifying the QZA formula for expanding domains (or in other words, with the hypothesis that silent virialisation is a fair approximation), about half the mass that contributes to deceleration of the global, effective scale factor is effectively cancelled, in terms of its contribution to average deceleration. Thus, the cancelling of about half the deceleration results in volume acceleration. This effect (related to that of “finite infinity” regions: Ellis 1984; Wiltshire 2007a; Cox 2007) is overlooked in the standard approach because Eq. (44) (cubical averaging of domain-wise scale factors) is not normally used.

In this sense, the ‘new physics’ that contributes significantly to dark energy is old physics – Newtonian virialisation – that is assumed to not have instantaneous effects at a distance. This could be described as non-compensation for virialisation. The global acceleration effect emerges indirectly from a statistical Newtonian particle effect via the combination of Eq. (14) for pre-virialisation average volume evolution (with $\Lambda\equiv 0$), Eq. (53) for virialised domains, and Eq. (44) to silently prevent the latter from affecting the former.

6.3 Caveats of the VQZA, the fluid model and vorticity

The accuracy of the VQZA depends on the accuracy of the assumptions stated in Definition 1, Sect. 4.2. For example, if propagation of relativistic constraints across domain boundaries implied that virialisation in one domain induced inaccuracy in the QZA formulae [Eqs (11), (LABEL:e-resultQ2), (31)] in one or more neighbouring domains, then this would lead to systematic error due to violation of the silent virialisation assumption of (ii) and (iv) of Definition 1. Other simplifying assumptions either stated above or implicit to the VQZA approach presented here include:

• (i)

a spacetime foliation that matches the EdS foliation at early times is assumed (proper time and coordinate time are not distinguished);

• (ii)

the difference between the sum of rest masses and the system mass is assumed to be negligible (see Adamek et al. 2016a, b, for an example of modelling Lorentzian effects);

• (iii)

equations (14), (15), and (16) are justified in the relativistic case under the assumption of zero vorticity;

• (iv)

there is assumed to be no shell-crossing prior to virialisation;

• (v)

although the scale found to be of most interest is found to be $L_{13.8}=2.5^{+0.1}_{-0.4}$ Mpc/$h_{\mathrm{eff}}$, the calculations were performed (in particular see Fig. 7) down to $L_{{\cal D}}=0.25$ Mpc/$h_{\mathrm{eff}}$, where it is clear that the zero vorticity and lack of shell-crossing assumptions could potentially lead to non-negligible correction terms.

Thus, as with each of the other approaches towards general-relativistic cosmological simulations published so far (see Sect. 1), the approximations and simplifying assumptions of this approach have to be studied further in order to judge how much they affect the results.

7 Conclusion

This work presents the first in a series of cosmological expansion modelling improvements. These aim to incorporate successively more elements towards obtaining general-relativistically valid $N$-body simulations in which the effective expansion rate is calculated from local structure formation rather than inserted arbitrarily. A key element of the present work in obtaining super-EdS expansion is described in Sect. 2 and illustrated quantitatively in Sect. 4. This is the virialisation–${\cal Q}_{\cal D}$ Zel’dovich approximation (Definition 1, Sect. 4.2), in which locally Newtonian expansion modelled by the QZA in combination with virialisation modelled under the stable clustering hypothesis together imply super-EdS evolution of the effective expansion rate.

The standard FLRW structure formation alternative to the VQZA can be summarised as follows:

• (i)

prior to a collapse/virialisation event, the Lagrangian per-domain average scale factor evolution is implicitly modelled to be smooth for both the collapsing domain and the expanding domain;

• (ii)

the final stages of a gravitational collapse and virialisation event are very sudden in comparison with cosmological time scales;

• (iii)

the global scale factor evolution is assumed to be smooth (FLRW) passing from times prior to the collapse/virialisation event through to times following the event;

• (iv)

assumptions (i) + (ii) + (iii) together imply that at least one expanding domain ${\cal D}^{-}$ must undergo a sudden drop in acceleration $\ddot{a}_{{\cal D}^{-}}$ when the collapse/virialisation event occurs, in order to compensate for the latter; smooth expansion would violate either (ii) or (iii); in the case of a biscale model, the mean spatial 3-Ricci curvature (in a relativistic interpretation) of the expanding domain must suddenly switch from negative curvature to zero curvature.

The hypothesis adopted here is that (iv) is relativistically unrealistic (as shown by the solid curves in Fig. 2 for biscale examples). This motivates the VQZA, in which (i) and (ii) hold, but (iii) is violated to the degree that any collapse/virialisation process is accepted to be a sudden event, and to be silent — without a significant sudden effect on expanding domains. Thus, the VQZA avoids (iv) the non-differentiability of $\dot{a}_{{\cal D}^{-}}$ for expanding domains ${\cal D}^{-}$. In other words, we postulate that smoothness in volume growth of the expanding domain, as approximated in the QZA, is relativistically more accurate than smoothness in global mean volume growth in situations in which collapse terminates suddenly in a small, already greatly contracted domain. In Sect. 4, this argument is illustrated analytically and numerically with a biscale model that shows the sudden change required to happen in the expanding domain in the standard model.

The consequences of the VQZA were examined here for a more generic case, by partitioning general $N$-body simulation initial conditions into initially cubical domains.

• (i)

We introduced (Sections 3.4.1, 3.4.2) a method of estimating numerical error in velocity gradient numerical estimates by applying Stokes’ theorem over S1. We found (Sect. 5.1) that for a simulation with $L_{\mathrm{box}}/L_{N}=128$, structure scales of $L_{\mathrm{box}}/64$ or greater had $S/N$ ratios in the accuracy of dtfe determination of the velocity gradient invariants ${\mathrm{I}}$ and ${\mathrm{II}}$ of about ten or higher, except when attempting sub–mean-interparticle-separation sampling at $L_{\mathrm{DTFE}}=L_{N}/2$.

• (ii)

Based on our VQZA simulations (Figs 7, 8), we consider $L_{13.8}=2.5^{+0.1}_{-0.4}$ Mpc/$h_{\mathrm{eff}}$ to be a reasonable representation of the averaging scale ${\cal D}$ which provides the optimal amount of super-EdS scale factor evolution to match observations as represented by an $\Lambda$CDM proxy. The error represents our estimate of systematic error associated with changing the box (T3 fundamental domain), dtfe and particle resolution scales, and running 100 VQZA simulations of size $N=256^{3}$ at each relevant set of scales. Our largest box size simulations yield L_{13.8}=2.56\pm 0.01\leavevmode\nobreak\ {\mbox{Mpc/h_{\mathrm{eff}}}}, where the error is one standard deviation (not the standard error in the mean), but since these have low ratios $L_{{\cal D}}/L_{\mathrm{DTFE}}=L_{\mathrm{DTFE}}/L_{N}=2$, it is likely that systematic, resolution-related error is greater than this random error.

• (iii)

This super-EdS scale factor evolution is associated with kinematical backreaction and curvature normalised functionals $\Omega_{{\cal Q}}^{{\cal D}}$ and $\Omega_{{\cal R}}^{{\cal D}}$, respectively, whose volume-weighted means and medians for uncollapsed domains are mostly positive throughout cosmic history (Figs 10, 11). Kinematical backreaction is strongest at early times, with $\left\langle\Omega_{{\cal Q}}^{{\cal D}}\right\rangle_{\mathrm{uncoll}}$ reaching a median value of about 0.03 at $t\approx 2$ Gyr and then dropping gradually towards zero. The recent emergence of average negative curvature is confirmed in these calculations, with the $\Omega_{{\cal R}}^{{\cal D}}$ mean and median values for uncollapsed domains growing to about 1.5 to 2, respectively, by the present.

• (iv)

Cosmic variance, as represented by independent realisations of initial conditions, gives different numbers of collapsed objects at a given cosmological time, and thus virialisation fractions $f_{\mathrm{vir}}$, for a given choice of $L_{{\cal D}}$. The super-EdS scale factor ratio $a_{\mathrm{eff}}/a_{\mathrm{EdS}}$ is very strongly correlated with $f_{\mathrm{vir}}$ (Sect. 5.3, Fig. 12).

• (v)

Artificially setting ${\cal Q}_{\cal D}=0$ gives $L_{13.8}({\cal Q}_{\cal D}:=0)\approx 8$ Mpc/$h_{\mathrm{eff}}$ (Figs 9, 14). This is the result that we would expect from Rácz et al. (2017)’s approach.

• (vi)

Small numbers of small exploratory simulations with numerical estimation of ${\cal Q}_{\cal D}$ at every major time step give rough agreement with the VQZA $L_{13.8}$ scale (Fig. 13).

Thus, from EdS initial conditions, dropping the decoupling hypothesis and estimating spatial expansion using the VQZA provides sufficient dark-energy–free average expansion by the standard age of the Universe estimate, for an averaging ${\cal Q}_{\cal D}$ scale of $L_{13.8}=2.5^{+0.1}_{-0.4}$ Mpc/$h_{\mathrm{eff}}$. This estimate is just a factor of about three lower than (more non-linear than) the 8\mbox{Mpc/h_{\mathrm{eff}}} scale at which the $\sigma_{8}$ normalisation of the initial power spectrum of density perturbations in flat space $P_{k}$ is normally defined. The historically important non-linearity scale, the $r_{0}$ scale in the two-point spatial auto-correlation function of galaxies, at which the excess probability of finding a pair of galaxies is unity, is close in value. This was estimated over four decades ago (Peebles 1974) as 5\mbox{Mpc/h_{\mathrm{eff}}}\,\lower 2.58334pt\hbox{\buildrel<\over{\sim}}\,r_{0}\,\lower 2.58334pt\hbox{\buildrel<\over{\sim}}\,10\mbox{Mpc/h_{\mathrm{eff}}} (Peebles & Hauser 1974, Eq. (48)). In this presentation of the VQZA, we set $\sigma_{8}$ appropriate for an early-epoch Planck-parametrised $\Lambda$CDM model (Sect. 3.3). Since silent virialisation is the key driver of the super-EdS expansion that we find here, it is unsurprising that (i) $L_{13.8}$, (ii) the 8Mpc/$h_{\mathrm{eff}}$ scale at which $\sigma_{8}$ is not far from unity, and (iii) the galaxy two-point auto-correlation function zero point $r_{0}$ are close in value to one another.

This work provides a step towards correcting standard $N$-body simulations for the general-relativistic constraints imposed by scalar averaging. This work is too preliminary to provide observational predictions at comparable precision to that of the $\Lambda$CDM model. Future work will need to include: inhomogeneous spatial curvature for calculating averages; curvature effects of particle movement; feedback of the global expansion rate on the effective background model growth rate $\xi(t)$; improvement beyond the QZA, since the QZA is only a first-order approximation (Buchert et al. 2013, Sections IV, V; see Alles et al. 2015), in particular by taking into account an effective background with evolving average curvature; and ray-tracing for calculating luminosity distances and angular diameter distances through inhomogeneously curved spatial sections (see, for example, Sikora & Głód 2017, in which the effective scale factor of the scalar-averaging averaging approach is found to be a good match to detailed ray-tracing calculations in a specific cosmological inhomogeneous metric).

Nevertheless, it seems reasonable to consider this work to be more accurate than the standard model: sudden, acausal switches in curvature and expansion rates in expanding regions that contain no singular spacetime events would be relativistically surprising.

Acknowledgments

Thank you to Jérémy Blaizot, Krzysztof Bolejko, Justyna Borkowska, Thomas Buchert, Mikołaj Korzyński, Bartosz Lew, Jan Ostrowski, Pierre Mourier, Yann Rasera, Joakim Rosdahl, Quentin Vigneron, David Wiltshire, and an anonymous referee for useful comments. A part of this project was funded by the National Science Centre, Poland, under grant 2014/13/B/ST9/00845. Part of this work consists of research conducted within the scope of the HECOLS International Associated Laboratory, supported in part by the Polish NCN grant DEC-2013/08/M/ST9/00664. A part of this project has made use of computations made under grant 314 of the Poznań Supercomputing and Networking Center (PSNC). This work has used the free-licensed GNU Octave package (Eaton et al. 2014).

Appendix A Newtonian velocity gradient invariants

Here we provide a script (App. A Ostrowski 2016) that can be run using the free-licensed computer algebra system maxima, version v.5.32.1,888http://maxima.sourceforge.net in order to confirm the equivalence of the expressions in Eq. (11). The three printf statements at the end should yield “0” to show algebraic equivalences.

/* Maxima script for checking Euclidean velocity gradient invariant relations

(C) 2016,2017 J.J. Ostrowski, B.F. Roukema, GPL-3+

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/gpl-3.0.html */

v: [v_x,v_y,v_z]; /* declare velocity fld */

depends(v,[x,y,z]); /* each component may depend on any coordinate */

derivabbrev : true; /* select subscript style of showing derivatives */

/* 1st order partial derivs of velocity */ dvdx: diff(v,x,1); dvdy: diff(v,y,1); dvdz: diff(v,z,1);

/* gradient of the velocity field */ vgrad: transpose(matrix(dvdx,dvdy,dvdz));

/* packages for matrix manipulations */ load("nchrpl"); load("eigen");

/* divergence and curl of vector field */ mydiv(v) := diff(v[1],x,1) + diff(v[2],y,1) + diff(v[3],z,1); mycurl(v) := [diff(v[3],y,1)-diff(v[2],z,1), diff(v[1],z,1)-diff(v[3],x,1), diff(v[2],x,1)-diff(v[1],y,1)];

/* gradient of scalar field */ mygrad(phi) := [diff(phi,x,1), diff(phi,y,1), diff(phi,z,1)];

/* directional derivative of vector w in direction v */ my_v_dot_del_w(v,w) := [mygrad(w[1]).transpose(v), mygrad(w[2]).transpose(v), mygrad(w[3]).transpose(v)];

// / divergence formula for second invariant / //

/* Invariant II by definition: / IIdefn : 1/2(mattrace(vgrad)^2 - mattrace(vgrad.vgrad)),expand;

/* Divergence formula */ IItmp : v * mydiv(v) - my_v_dot_del_w(v,v), expand; IIb : 1/2 * mydiv(IItmp), expand;

/* expansion rate */ theta : mydiv(v);

/* for Levi-Civita symbol */ load(itensor);

/* vorticity / angular_velocity : 1/2 * mycurl(v); / declare as a matrix */ vorticity_tensor : zeromatrix(3,3); for i : 1 thru 3 do for j : 1 thru 3 do for k : 1 thru 3 do vorticity_tensor[i,j] : vorticity_tensor[i,j] - levi_civita([i,j,k]) * angular_velocity[k];

vorticity_sq : 1/2 * mattrace( vorticity_tensor . transpose(vorticity_tensor) ), expand;

/* shear / / declare as a matrix */ shear_tensor : zeromatrix(3,3); for i : 1 thru 3 do for j : 1 thru 3 do shear_tensor[i,j] : vgrad[i,j] - 1/3 * kdelta([i],[j]) * theta - vorticity_tensor[i,j], expand;

/* verify that shear is symmetric, i.e. zero antisymmetric part */ shear_tensor - transpose(shear_tensor);

shear_sq : 1/2 * mattrace( shear_tensor . transpose(shear_tensor) ), expand;

/* Vorticity-shear-expansion formula */ IIc : vorticity_sq - shear_sq + 1/3 * theta^2, expand;

// / divergence formula for third invariant / //

/* Invariant III by definition: */ IIIdefn: determinant(vgrad);

/* Invariant III as sum of sums over indices: / vtr : mattrace(vgrad); IIIa : 1/3 * mattrace( vgrad.vgrad.vgrad ) - 1/2 * vtr mattrace(vgrad.vgrad) + 1/6 * (vtr)^3, expand;

/* Divergence formula / IIItmp : IIbv, expand; IIIb : 1/3 * mydiv( IIItmp - my_v_dot_del_w(IItmp,v)), expand;

/* Vorticity-shear-expansion formula */ IIIc : 1/3 * (1/9 * theta^3 + theta *(vorticity_sq - shear_sq) + mattrace(shear_tensor.shear_tensor. shear_tensor)) + angular_velocity.shear_tensor. transpose(angular_velocity), expand;

// / Check that definitions match alternative expressions / //

/* Defn of invariant vs divergence formula */ printf(false,"IIdefn - IIb: ~a", expand(IIdefn-IIb));

/* Defn of invariant vs vorticity-shear-expansion */ printf(false,"IIdefn - IIc: ~a", expand(IIdefn-IIc));

/* Defn of invariant vs sum of sums */ printf(false,"IIIdefn - IIIa: ~a", expand(IIIdefn-IIIa));

/* Defn of invariant vs divergence formula */ printf(false,"IIIdefn - IIIb: ~a", expand(IIIdefn-IIIb));

/* Defn of invariant vs vorticity-shear-expansion */ printf(false,"IIIdefn - IIIc: ~a", expand(IIIdefn-IIIc));

Appendix B Variance of the third invariant

$\mathbb{E}[\left<\mathbf{III}_{\mathbf{i}}\right>^{2}_{{\cal B}_{R}}]$

As in App. B of BKS, for a Gaussian density field of power spectrum $P_{\mathbf{i}}(k)$ for a Fourier mode $\mathbf{k}$ of amplitude $k:=|\mathbf{k}|$, we label independent modes $\mathbf{k}_{i}$, for $1\leq i\leq 6$, and estimate the ensemble mean $\mathbb{E}[.]$ by using the spatial mean. The six-point correlation function can be factorised by finding all distinct triples of pairs, similarly to (C17) of BKS00,

[TABLE]

where $\delta^{D}$ is the Dirac delta function. Use of (54) together with (C9) and (C10) of BKS, where ${\cal B}_{R}$ is a ball of radius $R$ and $\widetilde{W}_{R}$ is the Fourier transform of the normalised top-hat window function $W_{R}$, leads to

[TABLE]

rather than the original (C20)–(C22). Terms induced by (54) that involve $\tilde{W}_{R}({{k}}_{1})^{2}$ cancel exactly.

These expressions have been derived with the help of GNU Octave and the free-licensed computer algebra system maxima and tested numerically using numerical 18-dimensional integration with two sample points per dimension, performed with the GNU Scientific Library.