# Replacing dark energy by silent virialisation

Boudewijn F. Roukema

► To cite this version:

Boudewijn F. Roukema. Replacing dark energy by silent virialisation. Astronomy & Astrophysics - A&A, 2018, 610, pp.A51. 10.1051/0004-6361/201731400. hal-01730086

HAL Id: hal-01730086

https://hal.science/hal-01730086v1

Submitted on 5 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Replacing dark energy by silent virialisation

Boudewijn F. Roukema<fnref n="1" />,<fnref n="2" />

Received 19 June 2017 / Accepted 22 October 2017

## ABSTRACT

**Context.** 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.

**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 $Q_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 (T³) 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_0^{EdS} = 37.7$ km s⁻<fnref n="1" />/Mpc and an effective Hubble constant $H_{eff,0} = 67.7$ km s⁻<fnref n="1" />/Mpc.

**Results.** From 2000 simulations at resolution 256³, 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.4}^{+0.1}$ Mpc/h$_{eff}$. Relativistically interpreted, this corresponds to strong average negative curvature evolution, with the mean (median) curvature functional $\Omega_R^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 ΛCDM expansion to first order. The software packages used here are free-licensed.

**Key words.** cosmology: theory – cosmological parameters – large-scale structure of Universe – dark energy

## 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 Λ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 et al. 2000, 2013; Buchert & Ostermann 2012). 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 Λ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 ($Q_D$) Zel'dovich approximation (QZA). This provides an analytical expression for kinematical backreaction $Q_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

<footnote n="1">¹ Toruń Centre for Astronomy, Faculty of Physics, Astronomy and Informatics, Grudziadzka 5, Nicolaus Copernicus University, ul. Gagarina 11, 87-100 Toruń, Poland
e-mail: boud@astro.uni.torun.pl</footnote>

<footnote n="2">² Univ Lyon, ENS de Lyon, Univ Lyon1, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69007 Lyon, France</footnote>
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 Λ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, 2008; Buchert & Carfora 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 et al. 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, 2008; Nambu & Tanimoto 2005; Kai et al. 2007; 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. 2017; 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. (2017) using the free-licensed Einstein Toolkit, as well as those of Giblin et al. (2016a,b). All of these codes assume a $T^3$ 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 $T^3$ 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 et al. 2017). 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 $Q_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 $Q_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 Eqs. (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 $T^3$ spatial section is presented in Sect. 3.3. This is done using MPGRAFIC (this package and the others used
in this work are free-licensed<fnref n="1" /> 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 in-
variants can also be measured at later time steps for full *N*-body
calculations. To the degree that the peculiar gradient tensor ap-
proximates the (general-relativistic) extrinsic curvature tensor,
our approach should not deviate by much from a relativistic ap-
proach, even though the initial conditions are defined in Newto-
nian terms. The library form of the Delaunay Tessellation Field
Estimator *DTFE*, slightly extended by the addition of a small num-
ber 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*<sub>D</sub>. 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 *INHO-
MOG* libraries. In Sect. 3.6 we define how to relate the domain-
based average scale factors *a*<sub>D</sub>(*t*) to a global effective scale factor
*a*<sub>eff</sub>(*t*).

As a complement to the main model of this work, a method
of evolving initially cubical domains $\mathcal{D}$ in a full $N$-body simu-
lation, in which $Q_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 follow-
ing the Lagrangian volume defined by the particles of each initial
domain and again calculating an effective global expansion rate
$a_{\text{eff}}(t)$. This approach again uses RAMSES-SCALAV as a front end for
reading in initial conditions, and DTFE as the library for calculat-
ing the peculiar velocity gradient invariants. The INHOMOG library
is not used in this case.

Since the estimation of $Q_D$ is crucial to this work, the ac-
curacy of the DTFE numerical estimation of the peculiar velocity
gradient invariants (and thus $Q_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
$Q_D$ numerically in our $N$-body comparison method. The rela-
tion between Newtonian and relativistic reasoning is discussed
again in the light of these results, qualitatively presenting a per-
fect 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 oth-
erwise 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 pro-
vides the correct expression for the variance of the third in-
variant, 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 T³ 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 re-
minders in Kaiser 2017; Buchert 2018). In this case, if $Q_D$ evo-
lution in local domains is calculated within the constraints of
the same (expanding) T³ “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 cal-
culated in the local domains $\mathcal{D}$ using $Q_D$ and averaged. In terms
of kinematical backreaction, the global $Q_{T^3} = 0$ (see Sect. 3.4.2).
(Eq. (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 expan-
sion 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 $Q_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 $Q_D$ has
the same algebraic structure in the Newtonian (Buchert et al.
2000, Eq. (42); Buchert et al. 2013, Eq. (57); NZA) and relativis-
tic (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. (30), (31). This expression for $Q_D$
evolution is exact for some particular cases in the Newtonian
context (Buchert et al. 2000, III.D, III.E). In the relativistic con-
text, it is exact for plane symmetric collapse (Buchert et al. 2013,
Eq. (75)) and provides exact solutions for certain cases of spher-
ically symmetric collapse (Lemaître 1933; Tolman 1934; Bondi
1947; LTB) models (Buchert et al. 2013, V.B.1). The QZA re-
quires 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 T³ 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 T³ constraint; see
Sect. 3.4.2 for details.

How are virialisation, Newtonian and relativistic reasoning
related? The assumptions of the Newtonian collisionless, com-
pressible fluid approach for exact spherically symmetric col-
lapse, interpreted strictly, imply collapse to a singularity in a fi-
nite 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 col-
lapse. In other words, the validity of the approach is limited to
times before gravitational collapse events occur. More generi-
cally, 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 viri-
alisation of a set of particles is expected (and modelled) to oc-
cur instead of singularity formation. At a small enough scale,
baryon pressure also opposes gravitational collapse. Disregard-
ing baryon pressure, overdensities in a T³ *N*-body simulation
cannot be strictly spherically symmetric because they are com-
posed of particles; and they cannot evolve in an exact spherically

<footnote n="1">¹ https://www.gnu.org/licenses/licenses.html</footnote>
symmetric way because the global spatial section is T³. 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 T³ 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 T³ 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 viralised domains (presented in more detail for our main multi-domain partition modelling in Sect. 3.5.6). We thus define the Virialisation $Q_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_{EdS}$ and a reference-model expansion rate $H^{EdS}$ given by

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

where $H_1^{EdS} := H^{EdS}(a_{EdS} = 1)$. This EdS model is intended to approximately fit early epochs. In order to derive an observationally realistic value of $H_1^{EdS}$, an effective scale factor $a_{eff}$ which is nearly identical to the reference model scale factor at early times, that is, $a_{eff} \approx a_{EdS}$ for small $t$, is adopted. This is done by

normalising to $a_{eff} = 1$ at the present time $t_0 \equiv t_{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^{EdS} = 37.7 \text{ km s}^{-1}\text{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_{eff,0} = 67.74 \text{ km s}^{-1}\text{Mpc}$.

### 3.2. Scales

A standard T³ N-body simulation has three built-in length scales:

(i) $L_{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_{box}/N^{1/3}$, and

(v) $L_{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_{soft} := 2^{-\text{levelmax}}L_{box}$.

We insert two additional scales between $L_{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 back-reaction $Q_D$ and apply the scalar averaged evolution equations (Sect. 3.5). Since $Q_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_D$ – the scalar averaging initial domain size in comoving units (Sect. 3.5), and

(iii) $L_{DTFE}$ – the DTFE mesh size (Sect. 3.4).

We correspondingly define

$$n_D := L_{box}/L_D, \quad n_{DTFE} := L_{box}/L_{DTFE}, \quad n_{soft} := L_{box}/L_{soft}. \quad (2)$$

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

$$L_{box} \gg L_D \gg L_{DTFE} \gg L_N \gg L_{soft} \\ \Leftrightarrow 1 \ll n_D \ll n_{DTFE} \ll N^{1/3} \ll n_{soft} \quad (3)$$

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_{soft}$ is irrelevant. For the N-body comparison, we use $N = 128^3$ and set $n_{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)<fnref n="2" /> to generate cosmological initial conditions for an EdS model, with $H_1^{\text{EdS}} = 37.7$ 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_{\text{box}}$ (see Sect. 3.2) needs to be expressed in units of Mpc/$H_1^{\text{EdS}}/(100 \text{ km s}^{-1}/\text{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_{\text{box}}$ based on a given low-redshift length scale $L_{\text{DTFE}}$, so as above, we adopt the Ade et al. (2016) estimate of $H_{\text{eff},0} = 67.74$ km s⁻¹/Mpc. Thus, we have

$$ \frac{L_{\text{box}}}{\text{Mpc}/[H_1^{\text{EdS}}/(100 \text{ km s}^{-1}/\text{Mpc})]} = \frac{L_{\text{box}}}{\text{Mpc}/h_{\text{eff}}} \frac{H_1^{\text{EdS}}}{H_{\text{eff},0}} \\ = \frac{L_{\text{box}}}{L_{\text{DTFE}}} \frac{L_{\text{DTFE}}}{\text{Mpc}/h_{\text{eff}}} \frac{H_1^{\text{EdS}}}{H_{\text{eff},0}}, \quad (4) $$

where $h_{\text{eff}} := H_{\text{eff},0}/100 \text{ km s}^{-1}/\text{Mpc}$.

We match the Ade et al. (2016) power spectrum normalisation of $\sigma_8^{\Lambda\text{CDM}} = 0.8159$ to our EdS reference model. We evolve the effective Planck-parametrised $\Lambda\text{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_{\text{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 reference model expansion $a_{\text{EdS}}(t)$.

Delauayn 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)<fnref n="3" /> 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. 2010a,b; Sousbie 2011; Park et al. 2013; Pranav et al. 2017). DTFE uses the free-licensed library `cgal` (Computational Geometry Algorithms Library<fnref n="4" />) to make Delaunay tessellations.

A Delaunay tessellation of an N-body simulation snapshot is a division of the T<fnref n="3" /> 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_k^i, i = 1, 2, 3$, a linear interpolation of the peculiar velocity gradient $v_j^i$ relates velocity component changes from the 0th to the kth vertex for $k \neq 0$

$$ v_k^i - v_0^i = v_j^i (x_k^j - x_0^j), \quad (5) $$

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_k^j - x_0^j$ is invertible (which, in numerical practice, should almost always be the case), the velocity gradient can be estimated as

$$ v_{j,i}^i = (x_k^j - x_0^j)^{-1} (v_k^i - v_0^i). \quad (6) $$

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_{\text{DTFE}}^3$ cells within the simulation box, yielding an averaged interpolated estimate of $v_{j,i}^i$ in each DTFE grid cell.

#### 3.4.1. S¹-based estimates of velocity gradient uncertainties

We introduce a method of estimating the DTFE grid estimates of $v_{j,i}^i$ that does not seem to have been used previously. This new method is based on the numerical estimate of the integral of $v_{j,i}^i$ over any closed straight loop S¹ passing through the $n_{\text{DTFE}}$ points $x_k = \{x_k^j\}_j$ in an $n_{\text{DTFE}}^3$ cubical DTFE grid of velocity gradient estimates in a time slice of a simulation, where $k \in \{0, \dots, n_{\text{DTFE}} - 1\}$, and $x^{j1}, x^{j2}$ are two fixed values in the $j_1$th, $j_2$th directions, respectively, with $j_1 \neq j \neq j_2$. Analytically,

$$ \int_{S^1} v_{j,i}^i dx^j = 0 \quad (7) $$

should hold by the second fundamental theorem of calculus (Stokes' theorem) if $v_{j,i}^i$ 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_{j,i}^i}$ is Gaussian and identically and independently distributed at each of the DTFE grid points along a given S¹ loop. This is obviously an oversimplification. Collapsing or expanding structures on scales greater than $L_{\text{DTFE}}$ could reasonably lead to correlations in numerical errors in the estimates of $v_{j,i}^i$ across close pairs of DTFE grid points, and tensor invariants of $v_{j,i}^i$ 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

$$ \hat{\sigma}_{v_{j,i}^i} := \frac{\left| \sum_{k \in \{0, \dots, n_{\text{DTFE}} - 1\}} v_{j,i}^k \right|}{\sqrt{n_{\text{DTFE}} - 1}}, \quad (8) $$

<footnote n="2">² http://www2.iap.fr/users/pichon/mpgrafic.html, https://bitbucket.org/broukema/mpgrafic</footnote>

<footnote n="3">³ DTFE-1.1.1.Q: https://bitbucket.org/broukema/dtfE; DTFE-1.1.1: https://www.astro.rug.nl/%7Evoronoi/DTFE/dtfE.html</footnote>

<footnote n="4">⁴ https://www.cgal.org/</footnote>
where we use $n_{\text{DTFE}} - 1$ to increase the error estimate slightly at low $n_{\text{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_{\text{DTFE}}^2$ estimates for a fixed direction $j$. We define the rms (root mean square) of these $n_{\text{DTFE}}^2$ error estimates of the $i$th component summed in the $j$th direction:

$$ \hat{\sigma}_{v'_{j}}^{j} := \sqrt{\left\langle \left( \hat{\sigma}_{v'_{j}}^{j \neq j' \neq j_2} \right)^2 \right\rangle_{0 \le k_1 < n_{\text{DTFE}}, 0 \le k_2 < n_{\text{DTFE}}}}, \quad (9) $$

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 ESTIMATE SIGMA GRADIENT ONE DIREC. For propagating these uncertainties to estimates in uncertainties in the scalar invariants of $v'_{j}$ (Sect. 3.4.2), we estimate the rms over the diagonal and off-diagonal components of $v'_{j}$ separately:

$$ \hat{\sigma}_{v'_{j}}^{\text{diag}} := \sqrt{\left\langle \left( \hat{\sigma}_{v'_{j}}^{i} \right)^2 \right\rangle_i}, \quad \hat{\sigma}_{v'_{j}}^{\text{off}} := \sqrt{\left\langle \left( \hat{\sigma}_{v'_{j}}^{j} \right)^2 \right\rangle_{i,j,i \neq j}}. \quad (10) $$

These two error estimators summarise and simplify the uncertainty information from the $S^1$ 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. $S^1$ and $T^3$-based uncertainties in I and II

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

$$ I(v_{j,i}^i) := \text{tr}(v_{j,i}^i) = v_{i,i}^i = \nabla \cdot \mathbf{v} $$

$$ \begin{aligned} \text{II}(v_{j,i}^i) &:= \frac{1}{2} \left\{ \left[ \text{tr}(v_{j,i}^i) \right]^2 - \text{tr} \left[ \left( v_{j,i}^i \right)^2 \right] \right\} \\ &= \frac{1}{2} \left( (v_{i,i}^i)^2 - v_{j,i}^i v_{i,i}^j \right) \\ &= \frac{1}{2} \nabla \cdot \left( \mathbf{v}(\nabla \cdot \mathbf{v}) - (\mathbf{v} \cdot \nabla)\mathbf{v} \right) \end{aligned} $$

$$ \begin{aligned} \text{III}(v_{j,i}^i) &:= \det(v_{j,i}^i) \\ &= \frac{1}{3} v_{j,i}^i v_{i,k}^j v_{i,k}^k - \frac{1}{2} v_{i,i}^i (v_{j,i}^i v_{j,j}^j) + \frac{1}{6} (v_{i,i}^i)^3 \\ &= \frac{1}{3} \nabla \cdot \left\{ \frac{1}{2} \left[ \nabla \cdot (\mathbf{v}(\nabla \cdot \mathbf{v}) - (\mathbf{v} \cdot \nabla)\mathbf{v}) \right] \mathbf{v} \right. \\ & \qquad \left. - \left[ (\mathbf{v}(\nabla \cdot \mathbf{v}) - (\mathbf{v} \cdot \nabla)\mathbf{v}) \cdot \nabla \right] \mathbf{v} \right\}, \end{aligned} \quad (11) $$

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 $T^3$. Again, sums of I, II, or 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 $S^1$ 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 $S^1$ method of error estimation.

Again assuming statistically independent gaussian errors and using Eqs. (11) and (16), we obtain an $S^1$-based uncertainty for $\langle I \rangle_D, \langle II \rangle_D$, and $Q_D$ of

$$ \begin{aligned} \hat{\sigma}_{\langle I \rangle_D} &= \sqrt{3} \hat{\sigma}_{v'_{j}}^{\text{diag}} \\ \hat{\sigma}_{\langle II \rangle_D} &= \sqrt{\left\langle 2 \left( \hat{\sigma}_{v'_{j}}^{\text{diag}} v_{i,i}^i \right)^2 + \left[ \hat{\sigma}_{v'_{j}}^{\text{off}} v_{j,j}^i (1 - \delta_i^j) \right]^2 \right\rangle_x} \\ \hat{\sigma}_{Q_D} &= \sqrt{4 \hat{\sigma}_{\langle II \rangle_D}^2 + \frac{16}{9} \hat{\sigma}_{\langle I \rangle_D}^2 \langle v_{i,i}^i \rangle_x^2}, \end{aligned} \quad (12) $$

where $x$ represents all $n_{\text{DTFE}}^3$ DTFE grid points and $\delta_i^j \in \{0, 1\}$ is the usual Kronecker delta.

To test DTFE's numerical accuracy in estimating velocity field gradients, we define analytical velocity fields

$$ \begin{aligned} \mathbf{u}_1 &:= \left( \sin \frac{\pi n_x x^0}{L_{\text{box}}}, \sin \frac{\pi n_x x^1}{L_{\text{box}}}, \sin \frac{\pi n_x x^2}{L_{\text{box}}} \right), \\ \mathbf{u}_2 &:= \left( \sin \frac{\pi n_x x^1}{L_{\text{box}}}, \sin \frac{\pi n_x x^2}{L_{\text{box}}}, \sin \frac{\pi n_x x^0}{L_{\text{box}}} \right), \end{aligned} \quad (13) $$

for which I and 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 $I(\mathbf{u}_2) = 0 = \text{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_{\text{meas}}$ is the rms of the difference between the analytical expression for I or II and the numerical estimate?

(ii) What is the ratio of $\sigma_{\text{meas}}$ to the rms error predicted by the $S^1$-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. $Q_D$ Zel'dovich approximation (QZA)

Here we present our method of integrating the averaged Raychaudhuri evolution at the $L_D$ scale, used in our method (i), that we implement in the INHOMOG free-licensed (GPL-2+) library<fnref n="5" />. Averaging of $a_D$ within the full volume, that is, at the box scale $L_{\text{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, Eq. (9)) is:

$$ 3 \frac{\ddot{a}_D}{a_D} + 4\pi G \frac{M_{\text{D}i} a_{\text{D}i}^3}{V_{\text{D}i} a_D^3} - \Lambda = Q_D, \quad (14) $$

where $a_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 (i) values $a_D^3/a_{D_i}^3 = V_D(t)/V_{D_i}$,

<footnote n="5">⁵ https://bitbucket.org/broukema/inhomog, https://tracker.debian.org/inhomog</footnote>
$M_{\mathcal{D}_i}$ is the mass in the domain $\mathcal{D}$, and the averaged Hamiltonian constraint (Buchert et al. 2013, Eq. (10)) is:

$$ \left(\frac{\dot{a}_{\mathcal{D}}}{a_{\mathcal{D}}}\right)^2 - \frac{8\pi G}{3} \frac{M_{\mathcal{D}_i} a_{\mathcal{D}_i}^3}{V_{\mathcal{D}_i} a_{\mathcal{D}}^3} + \frac{\langle\mathcal{R}\rangle_{\mathcal{D}}}{6} - \frac{\Lambda}{3} = -\frac{Q_{\mathcal{D}}}{6} \quad (15) $$

(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 $\mathcal{D}$:

$$ Q_{\mathcal{D}} := 2 \langle \Pi \rangle_{\mathcal{D}} - \frac{2}{3} \langle \mathbf{I} \rangle_{\mathcal{D}}^2. \quad (16) $$

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

We use DTFE to numerically estimate the initial values of these invariants on a regular cubical mesh of $n_{\text{DTFE}}^3$ comoving positions in the EdS reference model. Thus, each domain $\mathcal{D}$ (again a cubical cell) contains $(n_{\text{DTFE}}/n_{\mathcal{D}})^3$ estimates of I and II. We average these to obtain $\langle \mathbf{I} \rangle_{\mathcal{D}}$ and $\langle \Pi \rangle_{\mathcal{D}}$. As stated above (Eq. (3)), it should be preferable to have $n_{\text{DTFE}}/n_{\mathcal{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)

$$ a_i := a_{\text{FL}}(t_i), \quad (17) $$

so that

$$ \varrho_{\text{FL}}(t) = \frac{\varrho_{\text{FL}}(t_i) a_i^3}{a(t)^3}, \quad (18) $$

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_{\text{FL}}(t_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 curvature tensor and the (normalised and zeroed) growth function $\xi$ (34) below) to have a time dimension ($[\Pi] = T^{-1}$, $[\Pi] = T^{-2}$, $[\Pi] = T^{-3}$, $[\xi] = T$, where $[x]$ denotes the dimension of $x$ and $T$ is the time dimension) or to be dimensionless. Buchert et al. (2013) uses the dimensional definitions throughout, except in Sect. VI.A. Here we adopt the dimensionless definitions. This requires the replacements

$$ \begin{aligned} \xi \rightarrow \frac{\dot{q}_i}{q_i} \xi, \qquad & \text{I} \rightarrow \left(\frac{q_i}{\dot{q}_i}\right) \text{I}, \\ \text{II} \rightarrow \left(\frac{q_i}{\dot{q}_i}\right)^2 \text{II}, \qquad & \text{III} \rightarrow \left(\frac{q_i}{\dot{q}_i}\right)^3 \text{III}, \end{aligned} \quad (19) $$

everywhere in the text of Buchert et al. (2013) except for Section VI.A, where the growing mode $q_{\text{FL}}(t)$ is given in Eqs. (35)

or (36). 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 RZA $a_{\mathcal{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 $\mathcal{D}$ in Lagrangian coordinates ($\mathcal{D}_i$),

$$ \begin{aligned} a_{\mathcal{D}_i} &= a_i \\ \dot{a}_{\mathcal{D}}(t_i) &= \dot{a}_{\text{FL}}(t_i) \left( 1 + \frac{1}{3} \langle \mathbf{I}_i \rangle_{\mathcal{D}_i} \right). \end{aligned} \quad (20) $$

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}_{\mathcal{D}}(t_i) > 0$, in order to select the positive square root for $\dot{a}_{\mathcal{D}}(t_i)$, yielding

$$ \begin{aligned} a_{\mathcal{D}_i} &= a_i \\ \dot{a}_{\mathcal{D}}(t_i) &= a_i \sqrt{\frac{8\pi G}{3} \frac{M_{\mathcal{D}_i}}{V_{\mathcal{D}_i} a_i^3} - \frac{\langle\mathcal{R}\rangle_{\mathcal{D}}(t_i)}{6} - \frac{Q_{\mathcal{D}}(t_i)}{6} + \frac{\Lambda}{3}}, \end{aligned} \quad (21) $$

where we postpone expressions for $\langle\mathcal{R}\rangle_{\mathcal{D}}(t_i)$ and $Q_{\mathcal{D}}(t_i)$ to (37) and (38) below.

### 3.5.3. Evolution equation

Writing the per-domain density $\rho_{\mathcal{D}} = (M_{\mathcal{D}_i}/V_{\mathcal{D}_i}) a_{\mathcal{D}_i}^3/a_{\mathcal{D}}^3$, Eq. (14) becomes

$$ 3 \frac{\ddot{a}_{\mathcal{D}}}{a_{\mathcal{D}}} + 4\pi G \langle \varrho \rangle_{\mathcal{D}} = Q_{\mathcal{D}} + \Lambda, \quad (23) $$

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

$$ 3 \frac{\ddot{a}_{\mathcal{D}}}{a_{\mathcal{D}}} + 4\pi G \frac{a_{\text{FL}}^3}{a_{\mathcal{D}}^3} \varrho_{\text{FL}} (1 - \langle \mathbf{I}_i \rangle_{\mathcal{D}_i}) = Q_{\mathcal{D}} + \Lambda, \quad (24) $$

and thus, using (18),

$$ 3 \frac{\ddot{a}_{\mathcal{D}}}{a_{\mathcal{D}}} + \frac{A}{a_{\mathcal{D}}^3} = Q_{\mathcal{D}} + \Lambda, \quad (25) $$

where the constant $A$ is defined

$$ \begin{aligned} A &:= 4\pi G \varrho_{\text{FL}}(t_i) a_{\mathcal{D}_i}^3 (1 - \langle \mathbf{I}_i \rangle_{\mathcal{D}_i}) \\ &= \frac{3}{2} H_{\text{FL}}^2(t_i) \Omega_m(t_i) a_{\mathcal{D}_i}^3 (1 - \langle \mathbf{I}_i \rangle_{\mathcal{D}_i}), \end{aligned} \quad (26) $$

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

$$ A = \frac{3}{2} H_{\text{FL}}^2(t_i) a_{\mathcal{D}_i}^3 (1 - \langle \mathbf{I}_i \rangle_{\mathcal{D}_i}). \quad (27) $$
This allows (22) to be rewritten

$$ \dot{a}_{\mathcal{D}}(t_i) = a_i \sqrt{\frac{2}{3} \frac{A}{a_i^3} - \frac{\langle \mathcal{R} \rangle_{\mathcal{D}}(t_i)}{6} - \frac{Q_{\mathcal{D}}(t_i)}{6} + \frac{\Lambda}{3}}, \quad (28) $$

and the Hamiltonian evolution equation becomes

$$ \dot{a}_{\mathcal{D}}(t) = \pm a_{\mathcal{D}} \sqrt{\frac{2}{3} \frac{A}{a^3} - \frac{\langle \mathcal{R} \rangle_{\mathcal{D}}(t)}{6} - \frac{Q_{\mathcal{D}}(t)}{6} + \frac{\Lambda}{3}}. \quad (29) $$

Typical collapsing solutions start with $\dot{a}_{\mathcal{D}} > 0$, the positive square root. Finding the first local minimum of $\dot{a}_{\mathcal{D}}^2$ in the positive square root case yields an estimate of the turnaround epoch $t_{\text{turn}}$ (when $\dot{a}_{\mathcal{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_{\text{turn}}$ and $a_{\mathcal{D}}(t)$, $\dot{a}_{\mathcal{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 $Q_{\mathcal{D}}$ and the mean curvature $\langle \mathcal{R} \rangle_{\mathcal{D}}$. The kinematical backreaction is given by Buchert et al. (2013, Eq. (50)) (or its Newtonian equivalent):

$$ Q_{\mathcal{D}} = \frac{\dot{\xi}^2 (\gamma_1 + \xi\gamma_2 + \xi^2\gamma_3)}{(1 + \xi \langle \mathbf{I}_i \rangle_I + \xi^2 \langle \mathbf{\Pi}_i \rangle_I + \xi^3 \langle \mathbf{\Pi}_i \rangle_I)^2}, \quad (30) $$

where

$$ \left\{ \begin{aligned} \gamma_1 &:= 2 \langle \mathbf{\Pi}_i \rangle_I - \frac{2}{3} \langle \mathbf{I}_i \rangle_I^2 \\ \gamma_2 &:= 6 \langle \mathbf{\Pi}_i \rangle_I - \frac{2}{3} \langle \mathbf{I}_i \rangle_I \langle \mathbf{I}_i \rangle_I \\ \gamma_3 &:= 2 \langle \mathbf{I}_i \rangle_I \langle \mathbf{\Pi}_i \rangle_I - \frac{2}{3} \langle \mathbf{\Pi}_i \rangle_I^2, \end{aligned} \right. \quad (31) $$

and the subscript $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, Eqs. (13), (54)) for the flat FLRW background case:

$$ \langle \mathcal{R} \rangle_{\mathcal{D}} = \frac{\dot{\xi}^2 (\tilde{\gamma}_1 + \xi \tilde{\gamma}_2 + \xi^2 \tilde{\gamma}_3)}{(1 + \xi \langle \mathbf{I}_i \rangle_I + \xi^2 \langle \mathbf{\Pi}_i \rangle_I + \xi^3 \langle \mathbf{\Pi}_i \rangle_I)^2}, \quad (32) $$

where

$$ \left\{ \begin{aligned} \tilde{\gamma}_1 &:= -2 \langle \mathbf{\Pi}_i \rangle_I - 12 \langle \mathbf{I}_i \rangle_I \frac{H}{\xi} - 4 \langle \mathbf{I}_i \rangle_I \frac{\xi}{\xi^2} \\ \tilde{\gamma}_2 &:= -6 \langle \mathbf{\Pi}_i \rangle_I - 24 \langle \mathbf{I}_i \rangle_I \frac{H}{\xi} - 8 \langle \mathbf{\Pi}_i \rangle_I \frac{\xi}{\xi^2} \\ \tilde{\gamma}_3 &:= -36 \langle \mathbf{\Pi}_i \rangle_I \frac{H}{\xi} - 12 \langle \mathbf{\Pi}_i \rangle_I \frac{\xi}{\xi^2}. \end{aligned} \right. \quad (33) $$

As detailed below in Sect. 3.5.5, we bypass the direct dependence on $\langle \mathcal{R} \rangle_{\mathcal{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, Eq. (32)),

$$ \xi(t) := \frac{[q_{\text{FL}}(t) - q_{\text{FL}}(t_i)]}{q_{\text{FL}}(t_i)}. \quad (34) $$

For the EdS reference model, the growing mode is

$$ q_{\text{FL}}(t) = a_{\text{FL}}(t) \quad (35) $$

(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_{m0} < 1$ and $\Omega_{\Lambda 0} := 1 - \Omega_{m0}$ is:

$$ q_{\text{FL}}(t) = \frac{5\Omega_{\text{m0}}}{2\Omega_{\Lambda 0}} \left\{ \frac{3\Omega_{\text{m0}} + \Omega_{\Lambda 0}a_{\text{FL}}}{\Omega_{\Lambda 0}a_{\text{FL}}} - \frac{3\Omega_{\text{m0}}\sqrt{\Omega_{\text{m0}} + \Omega_{\Lambda 0}a_{\text{FL}}}}{(\Omega_{\Lambda 0}a_{\text{FL}})^{3/2}} \operatorname{asinh}\left(\sqrt{\frac{\Omega_{\Lambda 0}a_{\text{FL}}}{\Omega_{\text{m0}}}}\right) \right\} \quad (36) $$

(Bildhauer et al. 1992, Eq. (21)).

Since $\xi_i = 0$ by definition (see (34)), the initial backreaction terms needed in (28) follow from (30)–(33):

$$ Q_{\mathcal{D}}(t_i) = \dot{\xi}_i^2 \left( 2 \langle \mathbf{\Pi}_i \rangle_I - \frac{2}{3} \langle \mathbf{I}_i \rangle_I^2 \right), \quad (37) $$

and

$$ \mathcal{R}_{\mathcal{D}}(t_i) = \dot{\xi}_i^2 \left( -2 \langle \mathbf{\Pi}_i \rangle_I - 12 \langle \mathbf{I}_i \rangle_I \frac{H}{\xi} - 4 \langle \mathbf{I}_i \rangle_I \frac{\xi}{\xi^2} \right). \quad (38) $$

For EdS, we have $a_{\text{FL}}(t) = a_i(t/t_i)^{2/3}$, so from (17) and (27) we have

$$ \begin{aligned} \frac{\dot{a}_{\text{FL}}(t_i)}{a_i} &= H_{\text{FL}}(t_i) = \frac{2}{3t_i} \\ A &= \frac{2a_i}{3t_i^2} (1 - \langle \mathbf{I}_i \rangle_{\mathcal{D}_i}). \end{aligned} \quad (39) $$

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

### 3.5.5. Numerical strategy

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

$$ \begin{aligned} \dot{a}_{\mathcal{D}} &= y \\ \dot{y} &= \frac{1}{3} \left[ (Q_{\mathcal{D}} + \Lambda) a_{\mathcal{D}} - \frac{A}{a_{\mathcal{D}}^2} \right], \end{aligned} \quad (40) $$

where $y \equiv \dot{a}_{\mathcal{D}}$ is defined to clarify the numerical strategy. This pair of first-order equations can be solved numerically by calculating $Q_{\mathcal{D}}(t)$ from (30) and (31) and the initial values of the invariants $\langle \mathbf{I}_i \rangle_{\mathcal{D}_i}$, $\langle \mathbf{\Pi}_i \rangle_{\mathcal{D}_i}$, and $\langle \mathbf{\Pi}_i \rangle_{\mathcal{D}_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. (in prep.), typical values of $\langle \mathbf{I}_i \rangle_{\mathcal{D}_i}$, $\langle \mathbf{\Pi}_i \rangle_{\mathcal{D}_i}$, and $\langle \mathbf{\Pi}_i \rangle_{\mathcal{D}_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. (B.2) in Appendix B, since it has not been published previously. However, in this work, we estimate $\langle \mathbf{I}_i \rangle_{\mathcal{D}_i}$, $\langle \mathbf{\Pi}_i \rangle_{\mathcal{D}_i}$, and $\langle \mathbf{\Pi}_i \rangle_{\mathcal{D}_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

$$
\begin{align}
\frac{\partial \dot{a}_{\mathcal{D}}}{\partial a_{\mathcal{D}}} &= 0 & \frac{\partial \dot{y}}{\partial a_{\mathcal{D}}} &= \frac{1}{3}(Q_{\mathcal{D}} + \Lambda) + \frac{2}{3} A a_{\mathcal{D}}^{-3} \nonumber \\
\frac{\partial \dot{a}_{\mathcal{D}}}{\partial y} &= 1 & \frac{\partial \dot{y}}{\partial y} &= 0 \nonumber \\
\frac{\partial \dot{a}_{\mathcal{D}}}{\partial t} &= 0 & \frac{\partial \dot{y}}{\partial t} &= \frac{1}{3} a_{\mathcal{D}} \frac{d(Q_{\mathcal{D}} + \Lambda)}{dt}. \tag{41}
\end{align}
$$

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),

$$
\begin{align}
\Omega_m^D &= \frac{8\pi G \langle\Omega\rangle_D}{3H_D^2} = \frac{H_{FL}^2(t_i) \Omega_m(t_i) a_{D_i}^3}{\dot{a}_D^2 a_D} (1 - \langle I_i \rangle_{D_i}) \nonumber \\
&= \frac{2 A \Omega_m(t_i)}{3 \dot{a}_D^2 a_D} \nonumber \\
\Omega_Q^D &= -\frac{Q_D}{6H_D^2} \nonumber \\
\Omega_\Lambda^D &= \frac{\Lambda}{3H_D^2} = \Omega_{\Lambda FL} \frac{H_{FL}^2}{H_D^2} \nonumber \\
\Omega_R^D &= 1 - \Omega_m^D - \Omega_\Lambda^D - \Omega_Q^D. \tag{42}
\end{align}
$$

3.5.6. Virialisation

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

$$
a_{\mathcal{D}}(t > t_{\text{coll}}) := a_{\text{EdS}}(t_{\text{coll}}) / (18\pi^2)^{1/3}, \quad (43)
$$

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, Eq. (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 impos-
ing 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,

$$
a_{\text{eff}}(t) := \left( \frac{\sum_{\mathcal{D}} a_{\mathcal{D}}^3(t)}{\sum_{\mathcal{D}} 1} \right)^{1/3} = \frac{\left( \sum_{\mathcal{D}} a_{\mathcal{D}}^3(t) \right)^{1/3}}{n_{\mathcal{D}}}, \quad (44)
$$

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 $Q_D$, so unless $\mathcal{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_{\text{eff}}(t) = a_{\text{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 Q<sub>D</sub>

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 code<fnref n="6" /> (Teyssier 2002;
Guillet & Teyssier 2011), available under CeCILL (GNU GPL
compatible) as Fortran 90 code, that we tentatively refer to as
ramses-scalav<fnref n="7" /> This serves as a front end to read in the ini-
tial conditions files (Sect. 3.3) and call the dTFE and inhomog li-
braries to estimate the per-domain initial invariants ⟨I⟩<sub>D<sub>i</sub></sub>, ⟨II⟩<sub>D<sub>i</sub></sub>,
and ⟨III⟩<sub>D<sub>i</sub></sub>, and to carry out VQZA evolution and volume aver-
aging as described in Sects. 3.5 and 3.6.

In an initial exploration of more numerical approaches, the
RAMSES-SCALAV front end also allows $Q_D$ to be calculated numer-
ically from the invariants at each major time step using Eq. (16)
instead of using the QZA analytical calculation of (30) 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 in-
tegrate using the classical (non-adaptive) Runge-Kutta (fourth-
order) algorithm. As in Sect. 3.6, we average the domain-level
scale factors $a_D$ to obtain $a_{\text{eff}}$ using Eq. (44). For any given do-
main $\mathcal{D}$, we store the list of particles initially present in $\mathcal{D}$, and
at time $t$ we calculate $\langle \text{I} \rangle_{\mathcal{D}}$ and $\langle \text{II} \rangle_{\mathcal{D}}$ by weighting I and II by the
volume element $d\mu_i$ of the $i$th particle. All parameters $I_i$, $\Pi_i$ and
$d\mu_i$ are calculated as (Euclidean) linear interpolations of the I, II
and $1/\rho$ values of the eight neighbouring dTFE cells surrounding
the particle. The global normalisation between $d\mu_i$ and $1/\rho_i$ is
arbitrary.

Since we do not change the list of particles in a domain be-
tween time steps, the domain shape, which we can conceptu-
ally imagine as the union of Voronoi tessellation tetrahedra sur-
rounding 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 T³ domain) in-
tegration, which is additive.

By default, a ramses cosmological simulation runs like a typ-
ical N-body simulation, with the reference model expansion in-
serted 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<sub>eff</sub>(t) (Eq. (44)) when cal-
culating 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.

<footnote n="6">⁶ https://bitbucket.org/rteyssie/ramses</footnote>

<footnote n="7">⁷ https://bitbucket.org/broukema/ramses-scalav</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Expanding-domain scale factor evolution $a_D^- / a_{EdS}$ for the biscale partition models illustrated in Fig. 1, comparing evolution of this domain according to the VQZA model (silent virialisation) and according to the standard N-body EdS “Newtonian” constraint (instantaneous feedback from virialisation). The VQZA scale factor evolution (+ symbols) correspond from bottom to top to those from bottom to top in Fig. 1. The standard-model scale-factor $a_D^{\text{Newt}}(t) / a_{EdS}(t)$ (solid curves; Eq. (48)) are indistinguishable following virialisation of the overdense domain, since by assumption, the expanding domain suddenly switches from hyperbolic (super-EdS) evolution to very-nearly flat (EdS) evolution ($a_D^{\text{Newt}}(t) / a_{EdS}(t) = 2^{1/3}$) when virialisation of the overdense domain occurs. A tiny difference from perfectly flat EdS evolution is present because the overdense domain has a small but non-zero fixed (stable clustering) volume.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Mean curvature functional $\Omega_R^{D^-}$ evolution of the expanding domain $D^-$ in the biscale partition models illustrated in Figs. 1 and 2. VQZA curvature evolution is shown with + symbols and the standard N-body EdS constraint (interpreted relativistically, not in Newtonian terms) is shown by curves. We hypothesise that the smooth curvature evolution (+) is relativistically more accurate than the sudden drop to flatness (solid curves) of the standard model.</figcaption></figure>

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 $T^3$ 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 Fig. 2 in the biscale case. This figure shows the evolution of the expanding domain $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

$$a_{D^{-}}^{\text{Newt}} := (2a_{\text{EdS}}^{3} - a_{D^{+}}^{3})^{1/3}, \qquad (48)$$

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: $\dot{a}_{D^{-}}^{\text{Newt}}$ has to suddenly become less positive or more negative when $\dot{a}_{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_R^{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,

$$\Omega_R^{D^-} := -\frac{a_{D^+}^3}{a_{D^-}^3} \Omega_R^{D^+} \frac{H_{D^+}^2}{H_{\text{FL}}^2}, \qquad (49)$$

where the $H_{D^-}^2/H_{\text{FL}}^2$ factor converts between globally normalised and per-domain normalised definitions of $\Omega$. Following collapse, the stable clustering assumption gives $H_{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_R^D$ growth to $\Omega_R^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.
<table><caption>**Table 1.** Software versions for numerical uncertainty tests on I, II (Sect. 5.1), and for VQZA (Sect. 5.2) and N-body-measured $Q_D$ (Sect. 5.4) simulations.</caption><thead><tr><th>Method</th><th>Package</th><th>Version and/or commit hash</th></tr></thead><tbody><tr><td>I, II errors</td><td>DTFE</td><td>1.1.1.Q, 1f2487e3</td></tr><tr><td>VQZA</td><td>MPGRAFIC</td><td>0.3.10</td></tr><tr><td>VQZA</td><td>RAMSES (ramses/trunk)</td><td>438f4659ae36</td></tr><tr><td>VQZA</td><td>RAMSES-SCALAV</td><td>53257c54</td></tr><tr><td>VQZA</td><td>DTFE</td><td>1.1.1.Q, 3d514657</td></tr><tr><td>VQZA</td><td>INHOMOG</td><td>0.0.67</td></tr><tr><td>N-body $Q_D$</td><td>MPGRAFIC</td><td>0.3.10</td></tr><tr><td>N-body $Q_D$</td><td>RAMSES (ramses/trunk)</td><td>438f4659ae36</td></tr><tr><td>N-body $Q_D$</td><td>RAMSES-SCALAV</td><td>1ba83730</td></tr><tr><td>N-body $Q_D$</td><td>DTFE</td><td>1.1.1.Q, 3d514657</td></tr><tr><td>N-body $Q_D$</td><td>INHOMOG</td><td>0.0.64</td></tr></tbody></table>

# 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.7}^{+0.3}$ 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_{\text{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 I and II

Measurement of $Q_D$ is crucial to calculating volume evolution, so we first consider numerical uncertainties. Application of the error estimation methods described in Sects. 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 $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_{\text{DTFE}} = 256$ points for both I and 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 $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.

<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Signal-to-noise ratios $S/N$ of I (above) and II (below) representing the ratio of rms signal to rms numerical noise $\sigma_{\text{meas}}$, for analytical model $u_1$ (which has non-zero invariants) in Eq. (13), sampled by realisations of $N = 128^3$ particles drawn from a uniform spatial random distribution for $8 \le n_{\text{DTFE}} \le 64$. The signal is negligibly small for $n_{\text{DTFE}} = 16, 128$ for $n_x = 8, 64$, respectively, since the test functions are sinusoidal and in phase with the box.</figcaption></figure>

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 $S^1$-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_{\text{DTFE}} = 2N$). These behaviours are reasonable.

The measured-to-predicted noise ratios for the second invariant $\Pi(u_2)$ (Fig. 5, lower panel, lower curves with + and × 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 $\Pi(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.
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** As in Fig. 4, ratio of measured noise $\sigma_{\text{meas}}$ to rms error predicted by S¹-based error estimators for I (above) and II (below), for both models $\mathbf{u}_1$ and $\mathbf{u}_2$. These ratios are within two orders of magnitude of unity in both cases.</figcaption></figure>

Global (T³) application of Stokes’ theorem gives consistent error propagation results to those for the S¹-based estimates, to within a few orders of magnitude. Figure 6 shows that $\langle I \rangle_{\text{box}}$ (global first invariant) errors are bounded above by the S¹-based errors with Gaussian error propagation. In contrast, $\langle \Pi \rangle_{\text{box}}$ errors are only bounded above by the S¹-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 $\langle \Pi \rangle_{\text{box}} = 0$ test at many $n_{\text{DTFE}}$ resolutions. The equivalent plot for $Q_{\text{box}}$ is not displayed, because the low amplitude of $|\langle I \rangle_{\text{box}}|$ estimates and Eq. (16) make the $\langle I \rangle_{\text{box}}^2$ contribution negligible, giving a plot visually identical to that for $\langle \Pi \rangle_{\text{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_{\text{box}} = 8L_{\mathcal{D}}$ and an MPGRAFIC grid of $N = 256^3$ particles, with $L_{\mathcal{D}}$ increasing within any of the three panels as labelled, and $L_{\mathcal{D}}/L_{\text{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_{\mathcal{D}} = 8$ Mpc/h$_{\text{eff}}$ and $L_{\mathcal{D}} = 16$ Mpc/h$_{\text{eff}}$. This is reasonable, because typical overdensities on these scales do not have time to collapse by the present.

On smaller averaging scales $L_{\mathcal{D}}$, silent virialisation (the VQZA, Definition 1, Sect. 4.2) starts to play a role in volume

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** As in Fig. 5, measured global $|\langle I \rangle_{\text{box}}|$ (above) and $|\langle \Pi \rangle_{\text{box}}|$ (below), which should be zero in the absence of numerical error (by Stokes’ theorem and Eq. (11)), shown as symbols, compared to S¹-based error estimators under independent Gaussian error propagation, $\hat{\sigma}_{(I)_{\text{box}}}$, $\hat{\sigma}_{(\Pi)_{\text{box}}}$, shown as curves, for both models $\mathbf{u}_1$ and $\mathbf{u}_2$.</figcaption></figure>

expansion. This can be seen in Fig. 7, where the scale factor growth is stronger as $L_{\mathcal{D}}$ decreases. The ΛCDM proxy scale factor evolution (× symbols) mostly lies between the $L_{\mathcal{D}} = 2$ and 4 Mpc/h$_{\text{eff}}$ curves. In other words, these curves provide the required super-EdS expansion required to reach a unity effective scale factor $a_{\text{eff}}(t)$ by $t = 13.8$ Gyr.

Figure 7 also indicates that the preferred $L_{\mathcal{D}}$ scale for providing sufficient super-EdS expansion is only weakly sensitive to the choice of $L_{\text{box}}$, $L_{\text{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_{\text{box}}$ sizes here are small. The box size $L_{\text{box}}$ can be increased further, but at the expense of leaving only small ratios $L_{\mathcal{D}}/L_{\text{DTFE}}$ and $L_{\text{DTFE}}/L_N$. Figure 8 shows that $2$ Mpc/h$_{\text{eff}} \lesssim L_{\mathcal{D}} \lesssim 4$ Mpc/h$_{\text{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_{\mathcal{D}}$ which gives $a_{\text{eff}}(13.8 \text{ Gyr}) = 1$. This is done here as follows. For a set of five simulations with different $L_{\mathcal{D}}$ values, as in one of the panels of Fig. 7, estimate $t_{a_{\text{eff}}=1}$ by a spline of $t$ against $a_{\text{eff}}(t)$ for each simulation. Least-squares fit a quadratic to $t_{a_{\text{eff}}=1}(L_{\mathcal{D}})$ and solve $t_{a_{\text{eff}}=1}(L_{\mathcal{D}}) = 13.8$ Gyr to obtain $L_{\mathcal{D}}$. We find that for $L_{\text{box}} = 8L_{\mathcal{D}} = 4L_{\text{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$_{\text{eff}}$. This is the scale at which $a_{\text{eff}}(t)$ reaches a unity scale factor at 13.8 Gyr, which is 16% above $a_{\text{EdS}}(13.8 \text{ Gyr})$. The error estimate is the standard
<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Full box cube-root-of-total-volume scale factor $a_{\text{eff}}(t)$ evolution for top: $L_{\text{box}} = 8L_{\text{D}} = 2L_{\text{DTFE}} = 16L_{\text{N}}$, middle: $L_{\text{box}} = 8L_{\text{D}} = 4L_{\text{DTFE}} = 8L_{\text{N}}$, bottom: $L_{\text{box}} = 8L_{\text{D}} = 8L_{\text{DTFE}} = 4L_{\text{N}}$, for $N = 256^3$ particles; where $h_{\text{eff}}$ is written as $h$ to reduce clutter. The $L_{\text{D}}$ scales are labelled, increasing in length *from top to bottom* within any panel. Cosmic variance is visible as non-uniformity in the $L_{\text{D}}$ dependence of the separation of the curves. Online: a fixed scale $L_{\text{D}}$ is shown with the same colour in all panels in this and later figures (for example, $L_{\text{D}} = 2.0 \text{ Mpc/h}_{\text{eff}}$ is purple).</figcaption></figure>

deviation of the 100 estimates. This error estimate represents the effect of cosmic variance for $4 \text{ Mpc/h}_{\text{eff}} \le L_{\text{box}} \le 64 \text{ Mpc/h}_{\text{eff}}$.

Corresponding values for smaller and bigger box sizes, $L_{\text{box}} = 8L_{\text{D}} = 2L_{\text{DTFE}} = 16L_{\text{N}}$ (Fig. 7, top panel) and $L_{\text{box}} = 8L_{\text{D}} = 8L_{\text{DTFE}} = 4L_{\text{N}}$ (Fig. 7, bottom panel), are

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** As in Fig. 7, full box cube-root-of-total-volume scale factor $a_{\text{eff}}(t)$, for $L_{\text{box}} = 64L_{\text{D}} = 2L_{\text{DTFE}} = 2L_{\text{N}}$, for $N = 256^3$ particles: the $L_{\text{D}} = 2 \text{ Mpc/h}_{\text{eff}}$ scale corresponds to a 128 Mpc/h_{eff} comoving box size.</figcaption></figure>

$L_{13.8} = 2.14 \pm 0.26 \text{ Mpc/h}_{\text{eff}}$ and $L_{13.8} = 2.59 \pm 0.25 \text{ Mpc/h}_{\text{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_{\text{box}} = 8L_{\text{D}} = 4L_{\text{DTFE}} = 8L_{\text{N}}$, as in Fig. 7, middle panel, for one quintuple of realisations): $L_{13.8} = 2.5_{-0.4}^{+0.1} \text{ Mpc/h}_{\text{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_{\text{D}}/L_{\text{DTFE}} = L_{\text{DTFE}}/L_{\text{N}} = 2$). This set of simulations gives a much more precise estimate: $L_{13.8} = 2.56 \pm 0.01 \text{ Mpc/h}_{\text{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.4}^{+0.1} \text{ Mpc/h}_{\text{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_{\text{D}}$. Using the FLRW critical density $\rho_{\text{crit}} = 1.52 \times 10^{11} M_{\odot}/\text{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_{\odot}$, corresponds to $1.10 \text{ Mpc/h}_{\text{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 $Q_{\text{D}}(t) = 0$ in all domains $\mathcal{D}$ at all times, then the super-EdS expansion rates should increase for any fixed $L_{\text{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}(Q_{\text{D}} := 0) \approx 8 \text{ Mpc/h}_{\text{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-$Q_{\text{D}}$ simulations.

Figure 10 helps to show the role of kinematical backreaction in our calculations, where the volume-weighted mean and
<figure><img src="image_9.png" /><figcaption>Fig. 9. As in Fig. 7, full box cube-root-of-total-volume scale factor $a_{\text{eff}}(t)$, but ignoring kinematical backreaction (setting $Q_D = 0 \forall D \forall t$), for (top) $L_{\text{box}} = 8L_D = 4L_{\text{DTFE}} = 8L_N$, (middle) $L_{\text{box}} = 8L_D = 8L_{\text{DTFE}} = 4L_N$, (bottom) $L_{\text{box}} = 8L_D = 16L_{\text{DTFE}} = 2L_N$.</figcaption></figure>

median values for the uncollapsed domains are defined

$$ \begin{align} \langle \Omega_Q^D \rangle_{\text{uncoll}} &:= -\frac{\langle Q_D \rangle_{\text{uncoll}}}{6H_{\text{eff}}^2} \nonumber \\ \langle \Omega_R^D \rangle_{\text{uncoll}} &:= -\frac{\langle R_D \rangle_{\text{uncoll}}}{6H_{\text{eff}}^2} \nonumber \\ \mu(\Omega_Q^D)_{\text{uncoll}} &:= -\frac{Q_{D_k}}{6H_{\text{eff}}^2} \nonumber \\ \mu(\Omega_R^D)_{\text{uncoll}} &:= -\frac{R_{D_k}}{6H_{\text{eff}}^2}, \tag{50} \end{align} $$

<figure><img src="image_10.png" /><figcaption>Fig. 10. Evolution of the volume-weighted mean (above) and median (below) globally normalised kinematical backreaction functional $\Omega_Q^D$ of uncollapsed domains, corresponding to the same scales as in the middle panel of Fig. 7.</figcaption></figure>

where $D_k$ is the $k$th domain among the sorted $i$th scale factors of uncollapsed domains satisfying

$$ \sum_{i \le k-1} a_{D_i}^3 \le \frac{1}{2} \sum_i a_{D_i}^3 \quad \text{and} \quad \sum_{i \ge k+1} a_{D_i}^3 \le \frac{1}{2} \sum_i a_{D_i}^3, \tag{51} $$

and Eqs. (30)–(33) are used to evaluate $Q$ and $R$ (ties in the median, if they occur, are upper weighted). The $\Omega_Q^D$ volume-weighted mean and median values for uncollapsed domains are positive, corresponding to negative $Q_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_{\text{eff}}(13.8 \text{ Gyr}) = 1$, as can be seen in Fig. 9.

The values of the mean and median $\Omega_Q^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 $T^3$ constraint on the global kinematical backreaction $Q_T^3$. This does appear to be the case in the upper panel of
<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Evolution of the volume-weighted mean (above) and median (below) globally normalised curvature functional $\Omega_R^D$ of uncollapsed domains, corresponding to the same scales as in the middle panel of Fig. 7.</figcaption></figure>

Fig. 10 at very early times, but the first virialisation events happen quickly, so $\langle \Omega_Q^D \rangle_{\text{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 $\langle \Omega_Q^D \rangle_{\text{uncoll}}$ as defined for this figure. Forcing a smooth $a_{\text{EdS}}(t)$ (or $a_{\Lambda\text{CDM}}(t)$) global volume evolution artificially hides this spikiness. Unsurprisingly, the volume-weighted medians $\mu(\Omega_Q^D)_{\text{uncoll}}$ evolve much more smoothly, since the collapse of a domain shifts the median within the central part of the $\Omega_Q^D$ distribution, which is denser than the tails. For $L_D = 2 \text{ Mpc}/h_{\text{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_{\text{eff}}(t)$ is accompanied by emerging negative average curvature, as shown in Fig. 11. At the $L_D = 2 \text{ Mpc}/h_{\text{eff}}$ scale (middle curve at early times), the volume-weighted mean $\Omega_R^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_R^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,

<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Super-EdS ratio $a_{\text{eff}}/a_{\text{EdS}}$ versus virialisation fraction $f_{\text{vir}}$ for 100 VQZA simulations with $L_D = 2 \text{ Mpc}/h_{\text{eff}}$ and $L_{\text{box}} = 8L_D = 4L_{\text{DTFE}} = 8L_N$, at an early (○), middle (×), and late (+) time, as indicated. At any fixed time, $a_{\text{eff}}/a_{\text{EdS}}$ correlates strongly with $f_{\text{vir}}$.</figcaption></figure>

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_{\text{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_{\text{coll}}$ of domains $\mathcal{D}$ that have gravitationally collapsed at any given time and normalising:

$$f_{\text{vir}} := \frac{N_{\text{coll}}}{n_{\mathcal{D}}^3} \qquad (52)$$

Figure 12 shows that at any fixed time, VQZA simulations run from independent realisations of initial conditions yield super-EdS expansion $a_{\text{eff}}/a_{\text{EdS}}$ that correlate strongly with $f_{\text{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_{\text{vir}}$ and $a_{\text{eff}}/a_{\text{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_{\text{vir}}(t)$ that results from a given realisation of initial conditions at a given time $t$ and the amount of super-EdS expansion $a_{\text{eff}}(t)/a_{\text{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_{\text{eff}}(t)/a_{\text{EdS}}(t) := 1$ by definition.
<figure><img src="image_13.png" /></figure>

**Fig. 13.** As in Fig. 7 (middle), for N-body evolution and Lagrangian tracing of initial domains, for $L_{\text{box}} = 8L_{\mathcal{D}} = 4L_{\text{DTFE}} = 4L_N$, with $N = 128^3$ particles, and $n_{\text{soft}}/(N^{1/3}) = 32$. The curves' behaviour is not monotonic with respect to $L_{\mathcal{D}}$ scales; the top curve is for $L_{\mathcal{D}} = 2.0$ Mpc/$h_{\text{eff}}$, the middle pair of curves are for $L_{\mathcal{D}} = 1.0$ Mpc/$h_{\text{eff}}$ (upper) and $L_{\mathcal{D}} = 4.0$ Mpc/$h_{\text{eff}}$ (lower), and the bottom pair are for $L_{\mathcal{D}} = 8.0$ Mpc/$h_{\text{eff}}$ (upper) and $L_{\mathcal{D}} = 0.5$ Mpc/$h_{\text{eff}}$ (lower).

<figure><img src="image_14.png" /></figure>

**Fig. 14.** As in Fig. 9 (top), for N-body evolution and Lagrangian tracing of initial domains, for $L_{\text{box}} = 8L_{\mathcal{D}} = 4L_{\text{DTFE}} = 4L_N$, with $N = 128^3$ particles, (artificially setting $Q_{\mathcal{D}} = 0 \ \forall \mathcal{D} \ \forall t$). The curves' behaviour is monotonic with respect to $L_{\mathcal{D}}$ scales, as labelled.

## 5.4. N-body measurement of $Q_{\mathcal{D}}$

What happens when $Q_{\mathcal{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? Figs. 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 $Q_{\mathcal{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_{\mathcal{D}} = 2, 4, 8$ Mpc/$h_{\text{eff}}$ has slightly stronger $a_{\text{eff}}(t)$ than for the VQZA case. In contrast, for the sequence of smaller scales ($L_{\mathcal{D}} = 2, 1, 0.5$ Mpc/$h_{\text{eff}}$) the $a_{\text{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 \lesssim L_{13.8}/\text{Mpc}/h_{\text{eff}} \lesssim 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 $Q_{\mathcal{D}}$ and their effects on $a_{\text{eff}}(t)$.

# 6. Discussion

## 6.1. Should global backreaction from $T^3$ 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 $T^3$, provided that the assumptions underlying the equations are satisfied. For flat space, given the definition of $Q_{\mathcal{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 $Q_{T^3}$ must be zero. Thus, Eq. (14) should yield $a_{\text{eff}}(t) \equiv a_{\text{EdS}}(t)$, no matter what scale of $\mathcal{D}$ is chosen, provided that the evolution of $Q_{\mathcal{D}}$ is calculated exactly. More realistically, in the present case $Q_{\mathcal{D}}$ is calculated using the QZA (Eqs. (30), (31)), so $a_{\text{eff}}(t) \approx a_{\text{EdS}}(t)$ should still be a good approximation, to the degree that Eqs. (30), (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 $Q_{\mathcal{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. 2017a,b; Ostrowski 2016; Ostrowski et al., in prep.). 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_{\mathcal{D}} \to 0^+$ in Eq. (14), $Q_{\mathcal{D}}$ via Eqs. (30) and (31) is not able to cancel self-gravity, so a literal interpretation of Eq. (14) would imply that $\dot{a}_{\mathcal{D}} \to -\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 $\mathcal{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 2018, 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, 2018; Bolejko & Lasky 2008; Bolejko & Ferreira 2012). 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 non-linear, 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

$$ \ddot{a}_{\mathcal{D}} \approx 0 \qquad (53) $$

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), (30), (31)) in one or more neighbouring domains, then this would lead to systematic error due to virialisation 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) Eqs. (14)–(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$_{\text{eff}}$, the calculations were performed (in particular see Fig. 7) down to $L_{\mathcal{D}} = 0.25$ Mpc/h$_{\text{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–$Q_{\mathcal{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 mod-
elled 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 $\mathcal{D}^{-}$ must undergo a sudden drop in acceleration $\dot{a}_{\mathcal{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}_{\mathcal{D}}$- for expanding domains $\mathcal{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 (Sects. 3.4.1, 3.4.2) a method of estimating numerical error in velocity gradient numerical estimates by applying Stokes' theorem over $S^1$. We found (Sect. 5.1) that for a simulation with $L_{\text{box}}/L_N = 128$, structure scales of $L_{\text{box}}/64$ or greater had $S/N$ ratios in the accuracy of DTFE determination of the velocity gradient invariants I and II of about ten or higher, except when attempting sub-mean-interparticle-separation sampling at $L_{\text{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_{\text{eff}}$ to be a reasonable representation of the averaging scale $\mathcal{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 ($T^3$ 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$ Mpc/$h_{\text{eff}}$, where the error is one standard deviation (not the standard error in the mean), but since these have low ratios $L_{\mathcal{D}}/L_{\text{DTFE}} = L_{\text{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_Q^\mathcal{D}$ and $\Omega_R^\mathcal{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 $\langle \Omega_Q^\mathcal{D} \rangle_{\text{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_R^\mathcal{D}$ mean and median values for un-collapsed 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_{\text{vir}}$, for a given choice of $L_{\mathcal{D}}$. The super-EdS scale factor ratio $a_{\text{eff}}/a_{\text{EdS}}$ is very strongly correlated with $f_{\text{vir}}$ (Sect. 5.3, Fig. 12).

(v) Artificially setting $Q_{\mathcal{D}} = 0$ gives $L_{13.8}(Q_{\mathcal{D}} := 0) \approx 8$ Mpc/$h_{\text{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 $Q_{\mathcal{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 $Q_{\mathcal{D}}$ scale of $L_{13.8} = 2.5^{+0.1}_{-0.4}$ Mpc/$h_{\text{eff}}$. This estimate is just a factor of about three lower than (more non-linear than) the 8 Mpc/$h_{\text{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 \text{ Mpc}/h_{\text{eff}} \lesssim r_0 \lesssim 10 \text{ Mpc}/h_{\text{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 8 Mpc/$h_{\text{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, Sects. 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.

*Acknowledgements.* Thank you to Jérémy Blaizot, Krzysztof Bolejko, Justyna Borkowska, Thomas Buchert, Mikolaj 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).

## References

Adamek, J., Clarkson, C., Durrer, R., & Kunz, M. 2015, Phys. Rev. Lett., 114, 051302

Adamek, J., Daverio, D., Durrer, R., & Kunz, M. 2016a, Nat. Phys., 12, 346

Adamek, J., Daverio, D., Durrer, R., & Kunz, M. 2016b, J. Cosmol. Astropart. Phys., 7, 053

Ade, P. A. R., Aghanim, N., Armitage-Caplan, C., et al. 2014, A&A, 571, A26

Ade, P. A. R., Aghanim, N., Arnaud, M., et al. 2016, A&A, 594, A13

Alles, A., Buchert, T., Al Roumi, F., & Wiegand, A. 2015, Phys. Rev. D, 92, 023512

Aragón-Calvo, M. A., Platen, E., van de Weygaert, R., & Szalay, A. S. 2010a, ApJ, 723, 364

Aragon-Calvo, M. A., Shandarin, S. F., & Szalay, A. 2010b, ArXiv e-prints [arXiv:1006.4178]

Aslanyan, G., & Manohar, A. V. 2012, J. Cosmol. Astropart. Phys., 6, 3

Aurich, R. 2015, MNRAS, 452, 1493

Bagla, J. S. 2005, Current Science, 88, 1088

Barbosa, R. M., Chirinos Isidro, E. G., Zimdahl, W., & Piattella, O. F. 2016, Gen. Relativity Gravitation, 48, 51

Bensby, T., Yee, J. C., Feltzing, S., et al. 2013, A&A, 549, A147

Bentivegna, E., & Bruni, M. 2016, Phys. Rev. Lett., 116, 251302

Bernardeau, F., & van de Weygaert, R. 1996, MNRAS, 279, 693

Bertschinger, E. 2001, ApJS, 137, 1

Bildhauer, S., Buchert, T., & Kasai, M. 1992, A&A, 263, 23

Blanchard, A., Douspis, M., Rowan-Robinson, M., & Sarkar, S. 2003, A&A, 412, 35

Bolejko, K. 2009, Gen. Relativity Gravitation, 41, 1585

Bolejko, K. 2017a, J. Cosmol. Astropart. Phys., 06, 025

Bolejko, K. 2017b, ArXiv e-prints [arXiv:1707.01800]

Bolejko, K. 2017c, Class. Quant. Grav., 35, 024003

Bolejko, K., & Célérier, M.-N. 2010, Phys. Rev. D, 82, 103510

Bolejko, K., & Ferreira, P. G. 2012, J. Cosmol. Astropart. Phys., 5, 003

Bolejko, K., & Lasky, P. D. 2008, MNRAS, 391, L59

Bondi, H. 1947, MNRAS, 107, 410

Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80

Buchert, T. 1994, MNRAS, 267, 811

Buchert, T. 2000a, in Proceedings of the ninth workshop on general relativity and gravitation, Hiroshima November 3-6, 1999, eds. Y. Eriguchi, T. Futamase, A. Hosoya, et al. (Hiroshima: Physics Department, Hiroshima University), 306

Buchert, T. 2000b, Gen. Rel. Grav., 32, 105

Buchert, T. 2001, Gen. Rel. Grav., 33, 1381

Buchert, T. 2005, Class. Quant. Grav., 22, L113

Buchert, T. 2008, Gen. Rel. Grav., 40, 467

Buchert, T. 2011, Class. Quant. Grav., 28, 164007

Buchert, T. 2018, MNRAS, 473, L46

Buchert, T., & Carfora, M. 2002, Class. Quant. Grav., 19, 6109

Buchert, T., & Carfora, M. 2003, Phys. Rev. Lett., 90, 031101

Buchert, T., & Carfora, M. 2008, Class. Quant. Grav., 25, 195001

Buchert, T., & Ehlers, J. 1997, A&A, 320, 1

Buchert, T., Kerscher, M., & Sicka, C. 2000, Phys. Rev. D, 62, 043525

Buchert, T., Nayet, C., & Wiegand, A. 2013, Phys. Rev. D, 87, 123503

Buchert, T., & Ostermann, M. 2012, Phys. Rev. D, 86, 023520

Buchert, T., & Räsänen, S. 2012, Ann. Rev. Nucl. Part. Sci., 62, 57

Cautun, M. C., & van de Weygaert, R. 2011, ArXiv e-prints [arXiv:1105.0370]

Chiesa, M., Maino, D., & Majerotto, E. 2014, J. Cosmol. Astropart. Phys., 12, 49

Chirinos Isidro, E. G., Barbosa, R. M., Piattella, O. F., & Zimdahl, W. 2017, Class. Quant. Grav., 34, 035001

Coley, A. A. 2010, Class. Quant. Grav., 27, 245017

Cox, D. P. G. 2007, Gen. Rel. Grav., 39, 87

Daverio, D., Dirian, Y., & Mitsou, E. 2017, Class. Quant. Grav., 34, 237001

Dehnen, W., & Read, J. I. 2011, Eur. Phys. J. Plus, 126, 55

Duley, J. A. G., Nazer, M. A., & Wiltshire, D. L. 2013, Class. Quant. Grav., 30, 175006

Eaton, J. W., Bateman, D., Hauberg, S., & Wehbring, R. 2014, GNU Octave version 3.8.1 manual: a high-level interactive language for numerical computations (CreateSpace Independent Publishing Platform)

Ehlers, J., & Buchert, T. 1997, Gen. Rel. Grav., 29, 733

Eisenstein, D. J., & Hu, W. 1998, ApJ, 496, 605

Eke, V. R., Cole, S., & Frenk, C. S. 1996, MNRAS, 282, 263

Ellis, G. F. R. 1984, in General Relativity and Gravitation, eds. B. Bertotti, F. De Felice, & A. Pascolini, 215

Ellis, G. F. R., & Gibbons, G. W. 2014, Class. Quant. Grav., 31, 025003

Ellis, G. F. R., & Stoeger, W. 1987, Class. Quant. Grav., 4, 1697

Ellis, G. F., & Tsagas, C. G. 2002, Phys. Rev. D, 66, 124015

Friedmann, A. 1923, Mir kak prostranstvo i vremya (The Universe as Space and Time) (Petrograd: Academia)

Friedmann, A. 1924, Zeitschr. Phys., 21, 326

Giblin, Jr, J. T., Mertens, J. B., & Starkman, G. D. 2016a, Phys. Rev. Lett., 116, 251301

Giblin, Jr, J. T., Mertens, J. B., & Starkman, G. D. 2016b, ApJ, 833, 247

Guillet, T., & Teyssier, R. 2011, J. Comput. Phys., 230, 4756

Hunt, P., & Sarkar, S. 2007, Phys. Rev. D, 76, 123504

Kai, T., Kozaki, H., Nakao, K., Nambu, Y., & Yoo, C. 2007, Prog. Theor. Phys., 117, 229

Kaiser, N. 2017, MNRAS, 469, 744

Kasai, M. 1995, Phys. Rev. D, 52, 5605

Kašpar, P., & Svíttek, O. 2014, Class. Quant. Grav., 31, 095012

Kennel, M. B. 2004, ArXiv e-prints [arXiv:physics/0408067]

Krasinski, A. 1981, Gen. Rel. Grav., 13, 1021

Krasinski, A. 1982, in The Birth of the Universe, eds. J. Audouze, & J. Tran Thanh van, 15

Krasinski, A. 1983, Gen. Rel. Grav., 15, 673

Lacey, C., & Cole, S. 1993, MNRAS, 262, 627

Larena, J., Alimi, J.-M., Buchert, T., Kunz, M., & Corasaniti, P.-S. 2009, Phys. Rev. D, 79, 083011

Lavinto, M., Räsänen, S., & Szybka, S. J. 2013, J. Cosmol. Astropart. Phys., 12, 51

Lemaître, G. 1927, Ann. de la Soc. Sc. de Brux., 47, 49

Lemaître, G. 1931, MNRAS, 91, 483

Lemaître, G. 1933, Ann. de la Soc. Sc. de Brux., 53, 51

Macpherson, H. J., Lasky, P. D., & Price, D. J. 2017, Phys. Rev. D, 95, 064028

Matarrese, S., Pantano, O., & Saez, D. 1994a, MNRAS, 271, 513

Matarrese, S., Pantano, O., & Saez, D. 1994b, Phys. Rev. Lett., 72, 320

Matarrese, S., & Terranova, D. 1996, MNRAS, 283, 400

Morita, M., Nakamura, K., & Kasai, M. 1998, Phys. Rev. D, 57, 6094

Nadathur, S., & Sarkar, S. 2011, Phys. Rev. D, 83, 063506

Nambu, Y., & Tanimoto, M. 2005, ArXiv e-prints, [arXiv:gr-qc/0507057]

Nazer, M. A., & Wiltshire, D. L. 2015, Phys. Rev. D, 91, 063519

Ostrowski, J. J. 2016, Ph.D. Thesis, Toruń Centre for Astronomy, Faculty of Physics, Astronomy and Informatics Nicolaus Copernicus University, Torun, Poland

Ostrowski, J. J., Buchert, T., & Roukema, B. F. 2017a, in Proceedings of the Fourteenth Marcel Grossmann Meeting on General Relativity, 12–18 July 2015, Rome, eds. M. Bianchi, R. T. Jantzen, & R. Ruffini (Singapore: World Scientific), [arXiv:1602.00302]

Ostrowski, J. J., Buchert, T., & Roukema, B. F. 2017b, Acta Physica Polonica B Proceedings Supplement, 10, 407

Park, C., Pranav, P., Chingangbam, P., et al. 2013, J. Korean Astron. Soc., 46, 125

Peebles, P. J. E. 1974, ApJ, 189, L51

Peebles, P. J. E. 1980, Large-Scale Structure of the Universe (Princeton University Press)

Peebles, P. J. E., & Hauser, M. G. 1974, ApJS, 28, 19

Pranav, P., Edelsbunner, H., van de Weygaert, R., et al. 2017, MNRAS, 465, 4281

Prunet, S., Pichon, C., Aubert, D., et al. 2008, ApJS, 178, 179

Rácz, G., Dobos, L., Beck, R., Szapudi, I., & Csabai, I. 2017, MNRAS, 469, L1

Räsänen, S. 2004, J. Cosmol. Astropart. Phys., 2, 003

Räsänen, S. 2006, J. Cosmol. Astropart. Phys., 11, 003

Räsänen, S. 2008, J. Cosmol. Astropart. Phys., 4, 026

Robertson, H. P. 1935, ApJ, 82, 284

Roukema, B. F., Ostrowski, J. J., & Buchert, T. 2013, J. Cosmol. Astropart. Phys., 10, 043

Roukema, B. F., France, M. J., Kazimierczak, T. A., & Buchert, T. 2014, MNRAS, 437, 1096

Roukema, B. F., Buchert, T., Ostrowski, J. J., & France, M. J. 2015, MNRAS, 448, 1660

Roukema, B. F., Buchert, T., Fujii, H., & Ostrowski, J. J. 2016, MNRAS, 456, L45

Roukema, B. F., Mourier, P., Buchert, T., & Ostrowski, J. J. 2017, A&A, 598, A111
B. F. Roukema: Dark energy as silent virialisation

Roy, X., Buchert, T., Carloni, S., & Obadia, N. 2011, Class. Quant. Grav., 28, 165004

Schaap, W. E., & van de Weygaert, R. 2000, A&A, 363, L29

Sikora, S., & Glód, K. 2017, Phys. Rev. D, 95, 063517

Sousbie, T. 2011, MNRAS, 414, 350

Steiner, F. 2016, ArXiv e-prints [arXiv:1608.03133]

Stichel, P. C. 2016, ArXiv e-prints [arXiv:1601.07030]

Sussman, R. A., Hidalgo, J. C., Dunsby, P. K. S., & German, G. 2015, Phys. Rev. D, 91, 063512

Szekeres, P. 1975, Communications in Mathematical Physics, 41, 55

Taylor, E. F., & Wheeler, J. A. 1992, Spacetime physics—Introduction to special relativity (2nd edn.) (W.H. Freeman and Co.)

Teyssier, R. 2002, A&A, 385, 337

Tolman, R. C. 1934, Proc. Nat. Acad. Sci., 20, 169

van de Weygaert, R., & Schaap, W. 2009, in Data Analysis in Cosmology, eds. V. J. Martínez, E. Saar, E. Martínez-González, & M.-J. Pons-Bordería (Berlin: Springer Verlag), Lect. Not. Phys., 665, 291

Villa, E., Matarrese, S., & Maino, D. 2011, J. Cosmol. Astropart. Phys., 8, 024

Weinberg, S. 1972, Gravitation and cosmology: Principles and applications of the general theory of relativity (New York: Wiley)

Wiegand, A., & Buchert, T. 2010, Phys. Rev. D, 82, 023523

Wiltshire, D. L. 2007a, New J. Phys., 9, 377

Wiltshire, D. L. 2007b, Phys. Rev. Lett., 99, 251101

Wiltshire, D. L. 2009, Phys. Rev. D, 80, 123512

Yodzis, P. 1974, Proceedings of the Royal Irish Academy Section A, 74, 61

Zel'dovich, Y. B. 1970, A&A, 5, 84
# Appendix A: Newtonian velocity gradient invariants

Here we provide a script (Appendix A Ostrowski 2016) that can be run using the free-licensed computer algebra system MAXIMA, version v.5.32.1⁸, 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 field */

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*(mtrace(vgrad)^2 -
  mtrace(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 */
angularAIR 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]) *
      angularAIR velocity[k];

vorticity_sq : 1/2 *
  mtrace(vorticity tensor .
    transpose(vorticity tensor) ), expand;

/* shear */
/* declare as a matrix */
shearaplus : zeromatrix(3,3);
for i : 1 thru 3 do
  for j : 1 thru 3 do
    shearaplus[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 */
shearaplus - transpose(shearaplus);

shearsq : 1/2 * mtrace(shearaplus .
  transpose(shearaplus) ),
  expand;

/* Vorticity-shear-expansion formula */

<footnote n="8">⁸ http://maxima.sourceforge.net</footnote>
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)
- 1/2 * vtr* mattrace(vgrad.vgrad)
+ 1/6 * (vtr)^3, expand;

/* Divergence formula */

IIItmp : IIb*v, expand;

IIb: 1/3 *
mydiv( IIItmp - my_v Dot_Del_w(IIItmp,v)),
expand;

/* Vorticity-shear-expansion formula */

IIIc: 1/3 * (1/9 * theta^3 +
theta *(vorticity_sq - shear_sq) +
mattrace(shear parensor.shear parensor.
 shear parensor)) +
angularVelocity.shear parensor.
transpose(angularVelocity),
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}[\langle \mathbf{III}_i \rangle_{\mathcal{B}_R}^2]
$$

As in Appendix B of BKS, for a Gaussian density field of power
spectrum Pᵢ(k) for a Fourier mode **k** of amplitude *k* := |**k**|, we
label independent modes **k**ᵢ, for 1 ≤ *i* ≤ 6, and estimate the en-
semble mean E[.] by using the spatial mean. The six-point cor-
relation function can be factorised by finding all distinct triples
of pairs, similarly to (C17) of BKS00,

$$
\begin{equation}
\begin{split}
\mathbb{E}[\tilde{\delta}_i(\mathbf{k}_1)\tilde{\delta}_i(\mathbf{k}_2)\tilde{\delta}_i(\mathbf{k}_3)\tilde{\delta}_i(\mathbf{k}_4)\tilde{\delta}_i(\mathbf{k}_5)\tilde{\delta}_i(\mathbf{k}_6)] =
\end{split}
\end{equation}
\begin{align*}
& P_i(\mathbf{k}_1)P_i(\mathbf{k}_3)P_i(\mathbf{k}_5)\delta^D(\mathbf{k}_1 + \mathbf{k}_2)\delta^D(\mathbf{k}_3 + \mathbf{k}_4)\delta^D(\mathbf{k}_5 + \mathbf{k}_6) \\
& + P_i(\mathbf{k}_1)P_i(\mathbf{k}_3)P_i(\mathbf{k}_4)\delta^D(\mathbf{k}_1 + \mathbf{k}_2)\delta^D(\mathbf{k}_3 + \mathbf{k}_5)\delta^D(\mathbf{k}_4 + \mathbf{k}_6) \\
& + P_i(\mathbf{k}_1)P_i(\mathbf{k}_3)P_i(\mathbf{k}_4)\delta^D(\mathbf{k}_1 + \mathbf{k}_2)\delta^D(\mathbf{k}_3 + \mathbf{k}_6)\delta^D(\mathbf{k}_4 + \mathbf{k}_5) \\
& + P_i(\mathbf{k}_1)P_i(\mathbf{k}_2)P_i(\mathbf{k}_5)\delta^D(\mathbf{k}_1 + \mathbf{k}_3)\delta^D(\mathbf{k}_2 + \mathbf{k}_4)\delta^D(\mathbf{k}_5 + \mathbf{k}_6) \\
& + P_i(\mathbf{k}_1)P_i(\mathbf{k}_2)P_i(\mathbf{k}_4)\delta^D(\mathbf{k}_1 + \mathbf{k}_3)\delta^D(\mathbf{k}_2 + \mathbf{k}_5)\delta^D(\mathbf{k}_4 + \mathbf{k}_6) \\
& + P_i(\mathbf{k}_1)P_i(\mathbf{k}_2)P_i(\mathbf{k}_4)\delta^D(\mathbf{k}_1 + \mathbf{k}_3)\delta^D(\mathbf{k}_2 + \mathbf{k}_6)\delta^D(\mathbf{k}_4 + \mathbf{k}_5) \\
& + P_i(\mathbf{k}_1)P_i(\mathbf{k}_2)P_i(\mathbf{k}_5)\delta^D(\mathbf{k}_1 + \mathbf{k}_4)\delta^D(\mathbf{k}_2 + \mathbf{k}_3)\delta^D(\mathbf{k}_5 + \mathbf{k}_6) \\
& + P_i(\mathbf{k}_1)P_i(\mathbf{k}_2)P_i(\mathbf{k}_3)\delta^D(\mathbf{k}_1 + \mathbf{k}_4)\delta^D(\mathbf{k}_2 + \mathbf{k}_5)\delta^D(\mathbf{k}_3 + \mathbf{k}_6) \\
& + P_i(\mathbf{k}_1)P_i(\mathbf{k}_2)P_i(\mathbf{k}_3)\delta^D(\mathbf{k}_1 + \mathbf{k}_4)\delta^D(\mathbf{k}_2 + \mathbf{k}_6)\delta^D(\mathbf{k}_3 + \mathbf{k}_5) \\
& + P_i(\mathbf{k}_1)P_i(\mathbf{k}_2)P_i(\mathbf{k}_4)\delta^D(\mathbf{k}_1 + \mathbf{k}_5)\delta^D(\mathbf{k}_2 + \mathbf{k}_3)\delta^D(\mathbf{k}_4 + \mathbf{k}_6) \\
& + P_i(\mathbf{k}_1)P_i(\mathbf{k}_2)P_i(\mathbf{k}_3)\delta^D(\mathbf{k}_1 + \mathbf{k}_5)\delta^D(\mathbf{k}_2 + \mathbf{k}_4)\delta^D(\mathbf{k}_3 + \mathbf{k}_6) \\
& + P_i(\mathbf{k}_1)P_i(\mathbf{k}_2)P_i(\mathbf{k}_3)\delta^D(\mathbf{k}_1 + \mathbf{k}_5)\delta^D(\mathbf{k}_2 + \mathbf{k}_6)\delta^D(\mathbf{k}_3 + \mathbf{k}_4) \\
& + P_i(\mathbf{k}_1)P_i(\mathbf{k}_2)P_i(\mathbf{k}_4)\delta^D(\mathbf{k}_1 + \mathbf{k}_6)\delta^D(\mathbf{k}_2 + \mathbf{k}_3)\delta^D(\mathbf{k}_4 + \mathbf{k}_5) \\
& + P_i(\mathbf{k}_1)P_i(\mathbf{k}_2)P_i(\mathbf{k}_3)\delta^D(\mathbf{k}_1 + \mathbf{k}_6)\delta^D(\mathbf{k}_2 + \mathbf{k}_4)\delta^D(\mathbf{k}_3 + \mathbf{k}_5) \\
& + P_i(\mathbf{k}_1)P_i(\mathbf{k}_2)P_i(\mathbf{k}_3)\delta^D(\mathbf{k}_1 + \mathbf{k}_6)\delta^D(\mathbf{k}_2 + \mathbf{k}_5)\delta^D(\mathbf{k}_3 + \mathbf{k}_4),
\end{split}
\tag{B.1}
\end{equation}
$$

where $\delta^D$ is the Dirac delta function. Use of (B.1) together with
(C9) and (C10) of BKS, where $\mathcal{B}_R$ is a ball of radius $R$ and $\tilde{W}_R$ is
the Fourier transform of the normalised top-hat window function
$W_R$, leads to

$$
\begin{equation}
\begin{split}
\mathbb{E}[\langle \mathrm{III}_i \rangle_{\mathcal{B}_R}^2] = & \\
& \frac{32 \pi^6 \int_{\mathbb{R}^3} d^3 k_1 \int_{\mathbb{R}^3} d^3 k_2 \int_{\mathbb{R}^3} d^3 k_3}{3 k_1^4 k_2^4 k_3^4} \\
& \times [k_1^2 k_2^2 k_3^2 - 3 (\mathbf{k}_1 \cdot \mathbf{k}_3)^2 k_1^2 + 2 (\mathbf{k}_1 \cdot \mathbf{k}_2) (\mathbf{k}_1 \cdot \mathbf{k}_3) (\mathbf{k}_2 \cdot \mathbf{k}_3)] \\
& \times [k_1^2 k_2^2 k_3^2 - (\mathbf{k}_1 \cdot \mathbf{k}_2)^2 k_3^2 - (\mathbf{k}_1 \cdot \mathbf{k}_3)^2 k_2^2 - (\mathbf{k}_2 \cdot \mathbf{k}_3)^2 k_1^2 \\
& \qquad + 2 (\mathbf{k}_1 \cdot \mathbf{k}_2) (\mathbf{k}_1 \cdot \mathbf{k}_3) (\mathbf{k}_2 \cdot \mathbf{k}_3)]
\end{split}
\tag{B.2}
\end{equation}
$$

rather than the original (C20)–(C22). Terms induced by (B.1)
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 integra-
tion with two sample points per dimension, performed with the
GNU Scientific Library.