© ESO 2018

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 ≫ 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 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³ 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³ 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 thedensity 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³ 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¹ software). In Sect. 3.4 we explain how to numerically measure the invariants of what in Newtonian terms is the peculiar velocity gradient tensor. This is only done in the initial conditions for our main model. However, these invariants can also be measured at later time steps for fullN-body calculations.To the degree that the peculiar gradient tensor approximates the (general-relativistic) extrinsic curvature tensor, our approach should not deviate by much from a relativistic approach, even though the initial conditions are defined in Newtonian terms. The library form of the Delaunay Tessellation Field Estimator DTFE, slightly extended by the addition of a small number of functions and corrected for minor bugs, is used to estimate these invariants. The INHOMOG library is introduced in Sect. 3.5. This is used to evolve individual cubical domains using the QZA, each of comoving sidelength $L_D$. For convenience and in order to ease the use of this approach in full N-body simulations, this is done using RAMSES-SCALAV, an extension to RAMSES, as a front end to read in the initial conditions and call the DTFE and INHOMOG libraries. In Sect. 3.6 we define how to relate the domain-based average scale factors $a_D(t)$ to a global effective scale factor a_eff(t).

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

Since the estimation of $Q_D$ is crucial to this work, the accuracy 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 Λ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 relation between Newtonian and relativistic reasoning is discussed again in the light of these results, qualitatively presenting a perfect fluid effective model of virialisation, in Sect. 6. A summary and conclusions are given in Sect. 7. We set the spacetime unit conversion constant (Taylor & Wheeler 1992) c := 1 unless otherwise specified. Appendix A provides a computer algebra script for confirming the equivalence of expressions for invariants of the peculiar velocity gradient tensor (Eq. (11)). Appendix B provides the correct expression for the variance of the third invariant, which is wrong in Eqs. (C20)–(C22) of Buchert et al. (2000).

2 Virialisation versus spatial expansion

Although the aim here, as in several related projects, is general-relativistic cosmology simulations, it is useful to consider Newtonian reasoning. However, defining a mathematically solid Newtonian cosmological model is not easy. Buchert & Ehlers (1997) showed that a 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 reminders in Kaiser 2017; Buchert 2018). In this case, if $Q_D$ evolution inlocal 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 calculated in the local domains $D$ using $Q_D$ and averaged. In terms of kinematical backreaction, the global $Q_{{\text{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 expansion is only valid until gravitational collapse and virialisation (or more generally, shell-crossing). Before discussing virialisation, let us first consider the QZA, that is, the algebraic expression for evaluating $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 relativistic (Buchert et al. 2013, Eq. (50)) cases, but with differing spatial definitions and interpretations (Buchert et al. 2013, IV.A); we cite this expression in Eqs. (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 context, it is exact for plane symmetric collapse (Buchert et al. 2013, Eq. (75)) and provides exact solutions for certain cases of spherically symmetric collapse (Lemaître 1933; Tolman 1934; Bondi 1947; LTB) models (Buchert et al. 2013, V.B.1). The QZA requires knowing the growth function for linear perturbations in a reference model (in our case, the EdS reference model), and of the initial invariants of the extrinsic curvature tensor. Here, we estimate the latter using the Newtonian peculiar velocity gradients, generated for a standard 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, compressible fluid approach for exact spherically symmetric collapse, interpreted strictly, imply collapse to a singularity in a finite time. A Lagrangian model is easier to relate to a relativistic model than an Eulerian model in the sense that the Lagrangian approach ties matter and space together directly. However, in this idealised case of spherical collapse, a positive Lagrangian volume drops suddenly to zero at the final moment of collapse. In other words, the validity of the approach is limited to times before gravitational collapse events occur. More generically, the assumptions fail when shell crossing first occurs. In N-body simulations or with physical particles, the usual way to avoid these Newtonian singularities is to adopt a hybrid fluid-particle model: the fluid model allows the model to be defined globally as shown by Buchert & Ehlers (1997, or approximated regionally in an adaptive mesh refinement method), while virialisation of a set of particles is expected (and modelled) to occur instead of singularity formation. At a small enough scale, baryon pressure also opposes gravitational collapse. Disregarding baryon pressure, overdensities in a T³ N-body simulation cannot be strictly spherically symmetric because they are composed of particles; and they cannot evolve in an exact spherically symmetric way because the global spatial section is 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 (Ω_Λ). 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 virialised 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_{\text{EdS}}:={\left(3H_1^{\text{EdS}}t/2\right)}^{2/3},H^{\text{EdS}}:={\dot{a}}_{\text{EdS}}/a_{\text{EdS}}=2/(3t),$$(1)

where $H_1^{\text{EdS}}:=H^{\text{EdS}}(a_{\text{EdS}}=1)$. This EdS model is intended to approximately fit early epochs. In order to derive an observationally realistic value of $H_1^{\text{EdS}}$, an effective scale factor a_eff which is nearly identical to the reference model scale factor at early times, that is, a_eff ≈ a_EdS for small t, is adopted. This is done by normalising to a_eff = 1 at the present time $t_0\equiv t_{a_{\text{eff}}=1}$, where the Planck Surveyor (Table 4, sixth data column, Ade et al. 2016) parametrisation of the Λ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^{\text{EdS}}=37.7$km s^-1 ∕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.74km s^-1 ∕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^−levelmaxL_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 backreaction $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_{\text{box}}/L_D,n_{\text{DTFE}}:=L_{\text{box}}/L_{\text{DTFE}},n_{\text{soft}}:=L_{\text{box}}/L_{\text{soft}}.$$(2)

Using this terminology, numerical scalar-averaging modelling of cosmological expansion requires, in general, $$\begin{array}{ll}{L}_{\text{box}}\gg \hfill & L_D\gg L_{\text{DTFE}}\gg L_N\gg L_{\text{soft}}\hfill \\ \iff \hfill & 1\ll n_D\ll n_{\text{DTFE}}\ll N^{1/3}\ll n_{\text{soft}}\hfill \end{array}$$(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³ or 4096³, respectively. Given these heavy requirements in computer resources, we limit this initial study to numerical exploration of the modest value N = 256³. 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³ 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)² to generate cosmological initial conditions for an EdS model, with $H_1^{\text{EdS}}=37.7$km s^-1 ∕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_box (see Sect. 3.2) needs to be expressed in units of $\text{Mpc}/[H_1^{\text{EdS}}/(100 \text{km}{\text{s}}^{-1}\text{}/\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_box based on a given low-redshift length scale L_DTFE, so as above, we adopt the Ade et al. (2016) estimate of H_eff,0 = 67.74km s^-1 ∕Mpc. Thus, we have $$\begin{array}{ll}\frac{L_{\text{box}}}{\text{Mpc}/[H_1^{\text{EdS}}/(100 \text{km}{\text{s}}^{-1}\text{}/\text{Mpc})]}\hfill & =\frac{L_{\text{box}}}{\text{Mpc/}{h}_{\text{eff}}}\frac{H_1^{\text{EdS}}}{H_{\text{eff},0}}\hfill \\ \hfill & =\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}},\hfill \end{array}$$(4)

where h_eff := H_eff,0∕100 km s^-1 ∕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 ΛCDM model backwards to z = 1000, that is, to t = 547 Myr, and then forwards using our EdS model, obtaining σ₈ = 1.0395 when a_EdS = 1, which we use for MPGRAFIC initial conditions.

The value of L_N is set in MPGRAFIC, and later read in automatically by RAMSES.

3.4 Estimating “velocity gradient” invariants

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

Delaunay and Voronoi tessellation are two methods that Bernardeau & van de Weygaert (1996) argued provide close to optimal tracing of the (Newtonian) peculiar velocity field. Here, we use the former, as implemented in the Delaunay Tessellation Field Estimator (DTFE)³ 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⁴) to make Delaunay tessellations.

A Delaunay tessellation of an N-body simulation snapshot is a division of the T³ 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≠0 $$\begin{array}{ll}{v}_k^i-v_0^i\hfill & =v_{,j}^i(x_k^j-x_0^j),\hfill \end{array}$$(5)

for the usualEinstein 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 $$\begin{array}{ll}{v}_{,j}^i\hfill & ={(x_k^j-x_0^j)}^{-1}(v_k^i-v_0^i).\hfill \end{array}$$(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$ 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$ 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$ over any closed straight loop S¹ passing through the n_DTFE points $x_k={\left\{x_k^j\right\}}_j$ in an $n_{\text{DTFE}}{}^3$ cubical DTFE grid of velocity gradient estimates in a time slice of a simulation, where k ∈ {0, …, n_DTFE − 1}, and x^j₁, x^j₂ are two fixed values in thej₁th, j₂ th directions, respectively,with j₁ ≠ j ≠ j₂. Analytically, $$\begin{array}{ll}{\displaystyle {\int }_{S^1}{v}_{,j}^i}\text{d}{x}^j\hfill & =0\hfill \end{array}$$(7)

should hold by the second fundamental theorem of calculus (Stokes’ theorem) if $v_{,j}^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}$ 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_DTFE could reasonably lead to correlations in numerical errors in the estimates of $v_{,j}^i$ across close pairs of DTFE grid points, and tensor invariants of $v_{,j}^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 $$\begin{array}{ll}{\widehat{\sigma }}_{v_{,j}^i}^{j_1\ne j\ne j_2}\hfill & :=\frac{\left|{\displaystyle \sum _{k\in \{0,\dots ,n_{\text{DTFE}}-1\}}{v}_{,j}^i}\right|}{\sqrt{n_{\text{DTFE}}-1}},\hfill \end{array}$$(8)

where we use n_DTFE − 1 to increase the error estimate slightly at low n_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 ith component summed in the jth direction: $$\begin{array}{ll}{\widehat{\sigma }}_{v_{,j}^i}^j\hfill & :=\sqrt{{\langle {\left({\widehat{\sigma }}_{v_{,j}^i}^{j_1\ne j\ne j_2}\right)}^2\rangle }_{0\le k_1<n_{\text{DTFE}},0\le k_2<n_{\text{DTFE}}}},\hfill \end{array}$$(9)

where k₁, k₂ label the grid positions in the j₁th, j₂ th directions, respectively. We introduce this into DTFE with the function ESTIMATESIGMAGRADIENTONEDIREC. For propagating these uncertainties to estimates in uncertainties in the scalar invariants of $v_{,j}^i$ (Sect. 3.4.2), we estimate the rms over the diagonal and off-diagonal components of $v_{,j}^i$ separately: $$\begin{array}{ll}{\widehat{\sigma }}_{v_{,j}^i}^{\text{diag}}\hfill & :={\sqrt{\langle {\left({\widehat{\sigma }}_{v_{,i}^i}^i\right)}^2\rangle }}_i,{\widehat{\sigma }}_{v_{,j}^i}^{\text{off}}:=\sqrt{{\langle {\left({\widehat{\sigma }}_{v_{,j}^i}^j\right)}^2\rangle }_{i,j,i\ne j}}.\hfill \end{array}$$(10)

These two error estimators summarise and simplify the uncertainty information from the S¹ 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¹ and T³-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) $$\begin{array}{ll}\text{I}(v_{,j}^i):=\hfill & \text{tr}(v_{,j}^i)=v_{,i}^i=\nabla \cdot v\hfill \\ \text{II}(v_{,j}^i):=\hfill & \frac{1}{2}\left\{{\left[\text{tr}\left(v_{,j}^i\right)\right]}^2-\text{tr}\left[{\left(v_{,j}^i\right)}^2\right]\right\}\hfill \\ =\hfill & \frac{1}{2}\left({(v_{,i}^i)}^2-v_{,j}^i{v}_{,i}^j\right)\hfill \\ =\hfill & \frac{1}{2}\nabla \cdot (v(\nabla \cdot v)-(v\cdot \nabla )v)\hfill \\ \text{III}(v_{,j}^i):=\hfill & \text{det}(v_{,j}^i)\hfill \\ =\hfill & \frac{1}{3}{v}_{,j}^i{v}_{,k}^j{v}_{,i}^k-\frac{1}{2}{v}_{,i}^i(v_{,j}^i{v}_{,i}^j)+\frac{1}{6}{(v_{,i}^i)}^3\hfill \\ =\hfill & \frac{1}{3}\nabla \cdot \{\frac{1}{2}\left[\nabla \cdot (v(\nabla \cdot v)-(v\cdot \nabla )v)\right]v\hfill \\ \hfill & -\left[(v(\nabla \cdot v)-(v\cdot \nabla )v)\cdot \nabla \right]v\},\hfill \end{array}$$(11)

where v := (v⁰, v¹, v²). Thus, all three invariants can be expressed as divergences, so that Stokes’ theorem again applies, this time over the full T³. 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¹ 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¹ method of error estimation.

Again assuming statistically independent gaussian errors and using Eqs. (11) and (16), we obtain an S¹ -based uncertainty for ${\langle \text{I}\rangle }_D,{\langle \text{II}\rangle }_D$, and $Q_D$ of $$\begin{array}{ll}{\widehat{\sigma }}_{{\langle \text{I}\rangle }_D}\hfill & =\sqrt{3}{\widehat{\sigma }}_{v_{,j}^i}^{\text{diag}}\hfill \\ {\widehat{\sigma }}_{{\langle \text{II}\rangle }_D}\hfill & =\sqrt{{\langle 2{\left({\widehat{\sigma }}_{v_{,j}^i}^{\text{diag}}{v}_{,i}^i\right)}^2+{\left[{\widehat{\sigma }}_{v_{,j}^i}^{\text{off}}{v}_{,j}^i(1-{\delta }_i^j)\right]}^2\rangle }_x}\hfill \\ {\widehat{\sigma }}_{Q_D}\hfill & =\sqrt{4{\widehat{\sigma }}_{{\langle \text{II}\rangle }_D}{}^2+\frac{16}{9}{\widehat{\sigma }}_{{\langle \text{I}\rangle }_D}{}^2{\langle v_{,i}^i\rangle }_x{}^2},\hfill \end{array}$$(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{array}{ll}{u}_1\hfill & :=\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),\hfill \\ u_2\hfill & :=\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),\hfill \end{array}$$(13)

for which I and II are straightforward to calculate analytically; the only non-zero gradient components are aligned with the sinusoidal directions in u ₁, and orthogonal to them in u₂ (giving I(u₂) = 0 = II(u₂)). We study the two following questions:

• (i)

  What is the signal-to-noise ratio (S∕N) of the rms of u₁ or u ₂, where the signal is defined by Eq. (13) and the noise σ_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 σ_meas to the rms error predicted by the S¹-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⁵. Averaging of $a_D$ within the full volume, that is, at the box scale L_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 Λ.

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_{D_i}{a}_{D_i}^3}{V_{D_i}{a}_D^3}-\Lambda =Q_D,$$(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}$, $M_{D_i}$ is the mass in the domain $D$, and the averaged Hamiltonian constraint (Buchert et al. 2013, Eq. (10)) is: $${\left(\frac{{\dot{a}}_D}{a_D}\right)}^2-\frac{8\pi G}{3}\frac{M_{D_i}{a}_{D_i}^3}{V_{D_i}{a}_D^3}+\frac{{\langle \text{R}\rangle }_D}{6}-\frac{\Lambda }{3}=-\frac{Q_D}{6}$$(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 $D$: $$\begin{array}{ll}{Q}_D\hfill & :=2{\langle \text{II}\rangle }_D-\frac{2}{3}{\langle \text{I}\rangle }_D{}^2.\hfill \end{array}$$(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 $D$ (again a cubical cell) contains ${(n_{\text{DTFE}}/n_D)}^3$ estimates of I and II. We average these to obtain ${\langle \text{I}\rangle }_D$ and ${\langle \text{II}\rangle }_D$. As stated above (Eq. (3)), it should be preferable to have $n_{\text{DTFE}}/n_D\gg 1$.

We only consider the Λ = 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),$$(17)

so that $${\varrho }_{\text{FL}}(t)=\frac{{\varrho }_{\text{FL}}(t_i)a_i{}^3}{a{(t)}^3},$$(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_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 cuvature tensor and the (normalised and zeroed) growth function ξ ((34) below) to have a time dimension ([I] = T⁻¹, [II] = T⁻², [III] = T⁻³, [ξ] = 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{array}{lll}\hfill & \xi \to \frac{{\dot{q}}_i}{q_i}\xi ,\hfill & \text{I}\to \left(\frac{q_i}{{\dot{q}}_i}\right)\text{I},\hfill \\ \hfill & \text{II}\to {\left(\frac{q_i}{{\dot{q}}_i}\right)}^2\text{II},\hfill & \text{III}\to {\left(\frac{q_i}{{\dot{q}}_i}\right)}^3\text{III},\hfill \end{array}$$(19)

everywhere in the text of Buchert et al. (2013) except for Section VI.A, where the growing mode q_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 ${ℝ}^3$ or T³ (some other flat FLRW spatial sections also allow Fourier analysis). Here we adopt a standard ${ℝ}^3$ power spectrum (Eisenstein & Hu 1998).

Buchert et al. (2013) VI.A gives the initial conditions for ${}^{\text{RZA}}{a}_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 ξ and the invariants, and a bold i subscript toindicate that the invariants are calculated at the initial time on the domain $D$ in Lagrangian coordinates ($D_i$), $$\begin{array}{ll}{a}_{D_i}\hfill & =a_i\hfill \\ {\dot{a}}_D\left(t_i\right)\hfill & ={\dot{a}}_{\text{FL}}\left(t_i\right)\left(1+\frac{1}{3}{\langle {\text{I}}_i\rangle }_{{\text{D}}_i}\right).\hfill \end{array}$$(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}}_D(t_i)>0$, in order to select the positive square root for ${\dot{a}}_D(t_i)$, yielding $$\begin{array}{ll}{a}_{D_i}\hfill & =a_i\hfill \\ {\dot{a}}_D\left(t_i\right)\hfill & =a_i\sqrt{\frac{8\pi G}{3}\frac{M_{D_i}}{V_{D_i}{a}_i{}^3}-\frac{{\langle \text{R}\rangle }_D\left(t_i\right)}{6}-\frac{Q_D\left(t_i\right)}{6}+\frac{\Lambda }{3}},\hfill \end{array}$$

where we postpone expresssions for ${\langle \text{R}\rangle }_D\left(t_i\right)$ and $Q_D\left(t_i\right)$ to (37) and (38) below.

3.5.3 Evolution equation

Writing the per-domain density ${\rho }_D=\left(M_{D_i}/V_{D_i}\right)a_{D_i}^3/a_D^3$, Eq. (14) becomes $$3\frac{{\ddot{a}}_D}{a_D}+4\pi G{\langle \varrho \rangle }_D=Q_D+\Lambda ,$$(23)

and with Buchert et al. (2013, Eq. (92)) (with a dimensionless definition of ξ and the invariants) can be expressed as $$3\frac{{\ddot{a}}_D}{a_D}+4\pi G\frac{a_{\text{FL}}^3}{a_D^3}{\varrho }_{\text{FL}}\left(1-{\langle {\text{I}}_i\rangle }_{D_i}\right)=Q_D+\Lambda ,$$(24)

and thus, using (18), $$3\frac{{\ddot{a}}_D}{a_D}+\frac{A}{a_D^3}=Q_D+\Lambda ,$$(25)

where the constant A is defined $$\begin{array}{lll}A\hfill & :=\hfill & 4\pi G{\varrho }_{\text{FL}}(t_i)a_{D_i}^3\left(1-{\langle {\text{I}}_i\rangle }_{{\text{D}}_i}\right)\hfill \\ \hfill & =\hfill & \frac{3}{2}{H}_{\text{FL}}^2(t_i){\text{Ω}}_{\text{m}}(t_i)a_{D_i}^3\left(1-{\langle {\text{I}}_i\rangle }_{D_i}\right),\hfill \end{array}$$(26)

with H_FL := ȧ_FL∕a_FL and ${\text{Ω}}_{\text{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_{D_i}^3\left(1-{\langle {\text{I}}_i\rangle }_{D_i}\right).$$(27)

This allows (22) to be rewritten $$\begin{array}{ll}{\dot{a}}_D\left(t_i\right)\hfill & =a_i\sqrt{\frac{2}{3}\frac{A}{a_i{}^3}-\frac{{\langle \text{R}\rangle }_D\left(t_i\right)}{6}-\frac{Q_D\left(t_i\right)}{6}+\frac{\Lambda }{3}},\hfill \end{array}$$(28)

and the Hamiltonian evolution equation becomes $$\begin{array}{ll}{\dot{a}}_D\left(t\right)\hfill & =\pm a_D\sqrt{\frac{2}{3}\frac{A}{a^3}-\frac{{\langle \text{R}\rangle }_D\left(t\right)}{6}-\frac{Q_D\left(t\right)}{6}+\frac{\Lambda }{3}}\cdot \hfill \end{array}$$(29)

Typical collapsing solutions start with ${\dot{a}}_D>0$, the positivesquare root. Finding the first local minimum of ${\dot{a}}_D^2$ in the positive square root case yields an estimate of the turnaround epoch t_turn (when ${\dot{a}}_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_turn and $a_D(t),{\dot{a}}_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_D$ and the mean curvature ${\langle \text{R}\rangle }_D$. The kinematical backreaction is given by Buchert et al. (2013, Eq. (50)) (or its Newtonian equivalent): $$\begin{array}{ll}{Q}_D=\hfill & \frac{{\dot{\xi }}^2\left({\gamma }_1+\xi {\gamma }_2+{\xi }^2{\gamma }_3\right)}{{\left(1+\xi {\langle {\text{I}}_i\rangle }_I+{\xi }^2{\langle {\text{II}}_i\rangle }_I+{\xi }^3{\langle {\text{III}}_i\rangle }_I\right)}^2},\hfill \end{array}$$(30)

where $$\{\begin{array}{l}{\gamma }_1:=2{\langle {\text{II}}_i\rangle }_I-\frac{2}{3}{\langle {\text{I}}_i\rangle }_I{}^2\hfill \\ {\gamma }_2:=6{\langle {\text{III}}_i\rangle }_I-\frac{2}{3}{\langle {\text{II}}_i\rangle }_I{\langle {\text{I}}_i\rangle }_I\hfill \\ {\gamma }_3:=2{\langle {\text{I}}_i\rangle }_I{\langle {\text{III}}_i\rangle }_I-\frac{2}{3}{\langle {\text{II}}_i\rangle }_I{}^2,\hfill \end{array}$$(31)

and the subscript $I$ indicates that a fully relativistic calculation would integrate the curved, Riemannian volume even at thisearly time when perturbations are weak, and the dimensionless growth rate ξ is given in Eq. (34). The mean curvature is given by Buchert et al. (2013, Eqs. (13), (54)) for the flat FLRW background case: $$\begin{array}{ll}{\langle \text{R}\rangle }_D=\hfill & \frac{{\dot{\xi }}^2\left({\tilde{\gamma }}_1+\xi {\tilde{\gamma }}_2+{\xi }^2{\tilde{\gamma }}_3\right)}{{\left(1+\xi {\langle {\text{I}}_i\rangle }_I+{\xi }^2{\langle {\text{II}}_i\rangle }_I+{\xi }^3{\langle {\text{III}}_i\rangle }_I\right)}^2},\hfill \end{array}$$(32)

where $$\{\begin{array}{l}{\tilde{\gamma }}_1:=-2{\langle {\text{II}}_i\rangle }_I-12{\langle {\text{I}}_i\rangle }_I\frac{H}{\dot{\xi }}-4{\langle {\text{I}}_i\rangle }_I\frac{\ddot{\xi }}{{\dot{\xi }}^2}\hfill \\ {\tilde{\gamma }}_2:=-6{\langle {\text{III}}_i\rangle }_I-24{\langle {\text{II}}_i\rangle }_I\frac{H}{\dot{\xi }}-8{\langle {\text{II}}_i\rangle }_I\frac{\ddot{\xi }}{{\dot{\xi }}^2}\hfill \\ {\tilde{\gamma }}_3:=-36{\langle {\text{III}}_i\rangle }_I\frac{H}{\dot{\xi }}-12{\langle {\text{III}}_i\rangle }_I\frac{\ddot{\xi }}{{\dot{\xi }}^2}.\hfill \end{array}$$(33)

As detailed below in Sect. 3.5.5, we bypass the direct dependence on ${\langle \text{R}\rangle }_D$ by integrating Eq. (25) instead. Evaluation of these expressions requires the growth rate ξ(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)}\cdot$$(34)

For the EdS reference model, the growing mode is $$q_{\text{FL}}(t)=a_{\text{FL}}(t)$$(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 ΛCDM) with present density parameter Ω_m0 < 1 and Ω_Λ0 := 1 − Ω_m0 is: $$\begin{array}{ll}{q}_{\text{FL}}(t)=\hfill & \frac{5{\text{Ω}}_{\text{m0}}}{2{\text{Ω}}_{\Lambda 0}}\{\frac{3{\text{Ω}}_{\text{m0}}+{\text{Ω}}_{\Lambda 0}{a}_{\text{FL}}}{{\text{Ω}}_{\Lambda 0}{a}_{\text{FL}}}\hfill \\ \hfill & -\frac{3{\text{Ω}}_{\text{m0}}\sqrt{{\text{Ω}}_{\text{m0}}+{\text{Ω}}_{\Lambda 0}{a}_{\text{FL}}}}{{({\text{Ω}}_{\Lambda 0}{a}_{\text{FL}})}^{3/2}}\text{asinh}\left(\sqrt{\frac{{\text{Ω}}_{\Lambda 0}{a}_{\text{FL}}}{{\text{Ω}}_{\text{m0}}}}\right)\}\hfill \end{array}$$(36)

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

Since ξ_i = 0 by definition (see (34)), the initial backreaction terms needed in (28) follow from (30)–(33): $$\begin{array}{ll}{Q}_D(t_i)\hfill & ={\stackrel{.}{{\xi }_i}}^2\left(2{\langle {\text{II}}_i\rangle }_I-\frac{2}{3}{\langle {\text{I}}_i\rangle }_I{}^2\right),\hfill \end{array}$$(37)

and $$\begin{array}{ll}{R}_D(t_i)\hfill & ={\stackrel{.}{{\xi }_i}}^2\left(-2{\langle {\text{II}}_i\rangle }_I-12{\langle {\text{I}}_i\rangle }_I\frac{H}{\dot{\xi }}-4{\langle {\text{I}}_i\rangle }_I\frac{\ddot{\xi }}{{\dot{\xi }}^2}\right)\cdot \hfill \end{array}$$(38)

For EdS, we have $a_{\text{FL}}(t)=a_i{(t/t_i)}^{2/3}$, so from (17) and (27) we have $$\begin{array}{lll}\frac{{\dot{a}}_{\text{FL}}(t_i)}{a_i}\hfill & =\hfill & H_{\text{FL}}(t_i)=\frac{2}{3t_i}\hfill \\ A\hfill & =\hfill & \frac{2a_i}{3t_i{}^2}\left(1-{\langle {\text{I}}_i\rangle }_{D_i}\right).\hfill \end{array}$$(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{array}{lll}{\dot{a}}_D\hfill & =\hfill & y\hfill \\ \dot{y}\hfill & =\hfill & \frac{1}{3}\left[(Q_D+\Lambda )a_D-\frac{A}{a_D^2}\right],\hfill \end{array}$$(40)

where $y\equiv {\dot{a}}_D$ is defined toclarify the numerical strategy. This pair of first-order equations can be solved numerically by calculating $Q_D(t)$ from (30) and (31) and the initial values of the invariants ${\langle {\text{I}}_i\rangle }_{D_i}$, ${\langle {\text{II}}_i\rangle }_{D_i}$, and ${\langle {\text{III}}_i\rangle }_{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 {\text{I}}_i\rangle }_{D_i}$, ${\langle {\text{II}}_i\rangle }_{D_i}$, and ${\langle {\text{III}}_i\rangle }_{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 {\text{I}}_i\rangle }_{D_i}$, ${\langle {\text{II}}_i\rangle }_{D_i}$, and ${\langle {\text{III}}_i\rangle }_{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{array}{lll}\frac{\partial {\dot{a}}_D}{\partial a_D}=0\hfill & \hfill & \frac{\partial \dot{y}}{\partial a_D}=\frac{1}{3}(Q_D+\Lambda )+\frac{2}{3}A{a}_D^{-3}\hfill \\ \frac{\partial {\dot{a}}_D}{\partial y}=1\hfill & \hfill & \frac{\partial \dot{y}}{\partial y}=0\hfill \\ \frac{\partial {\dot{a}}_D}{\partial t}=0\hfill & \hfill & \frac{\partial \dot{y}}{\partial t}=\frac{1}{3}{a}_D\frac{D(Q_D+\Lambda )}{Dt}\cdot \hfill \end{array}$$(41)

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{array}{ll}{\text{Ω}}_{\text{m}}^D=\frac{8\pi G{\langle \varrho \rangle }_D}{3H_D{}^2}\hfill & =\frac{H_{\text{FL}}^2(t_i){\text{Ω}}_{\text{m}}(t_i)a_{D_i}^3}{{\dot{a}}_D^2{a}_D}\left(1-{\langle {\text{I}}_i\rangle }_{D_i}\right)\hfill \\ \hfill & =\frac{2}{3}\frac{A{\text{Ω}}_{\text{m}}(t_i)}{{\dot{a}}_D^2{a}_D}\hfill \\ {\text{Ω}}_Q^D\hfill & =-\frac{Q_D}{6H_D^2}\hfill \\ {\text{Ω}}_{\Lambda }^D\hfill & =\frac{\Lambda }{3H_D^2}={\text{Ω}}_{\Lambda }{}_{\text{FL}}\frac{H_{\text{FL}}^2}{H_D^2}\hfill \\ {\text{Ω}}_R^D\hfill & =1-{\text{Ω}}_{\text{m}}^D-{\text{Ω}}_{\Lambda }^D-{\text{Ω}}_Q^D.\hfill \end{array}$$(42)

3.5.6 Virialisation

After turnaround of a given domain $D$, that is, when ${\dot{a}}_D$ becomes negative, we continue integration towards numerical collapse ($a_D\ll 1$) at foliation time t_coll. The initial solution $\{a_D(t),{\dot{a}}_D(t):t<t_{\text{coll}}\}$ is used as an approximation that is improved iteratively. We assume virialisation – we set $$a_D(t>t_{\text{coll}}):=a_{\text{EdS}}(t_{\text{coll}})/{(18{\pi }^2)}^{1/3},$$(43)

where 18π² 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-Λ 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 imposing the reference model expansion rate and abandoning the aim of estimating structure-generated effective expansion.

3.6 Regional versus global averaging

As in Räsänen (2006, 3.1) and Rácz et al. (2017), we cubically average the domain-wise scale factors, $$a_{\text{eff}}(t):={\left(\frac{{\displaystyle \sum _D{a}_D^3}(t)}{{\displaystyle \sum _{D}1}}\right)}^{1/3}=\frac{{\left({\displaystyle \sum _D{a}_D^3}(t)\right)}^{1/3}}{n_D},$$(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 $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_eff (t) = a_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_D$

The two main classes of faster simulation techniques are tree codes and particle-mesh (PM) based methods (e.g., Bagla 2005; Dehnen & Read 2011). In this work, we provide a patch to the adaptive mesh refinement (AMR) RAMSES code⁶ (Teyssier 2002; Guillet & Teyssier 2011), available under CeCILL (GNU GPL compatible) as Fortran 90 code, that we tentatively refer to as RAMSES-SCALAV⁷ This serves as a front end to read in the initial conditions files (Sect. 3.3) and call the DTFE and INHOMOG libraries to estimate the per-domain initial invariants ${\langle \text{I}\rangle }_{D_i}$, ${\langle \text{II}\rangle }_{D_i}$, and ${\langle \text{III}\rangle }_{D_i}$, and to carry outVQZA evolution and volume averaging 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 numerically 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 integrate 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_eff using Eq. (44). For any given domain $D$, we store the list of particles initially present in $D$, and at time t we calculate ${\langle \text{I}\rangle }_D$ and ${\langle \text{II}\rangle }_D$ by weighting I and II by the volume element dμ_i of the ith particle. All parameters I_i, II _i and d μ_i are calculated as (Euclidean) linear interpolations of the I, II and 1∕ρ values of the eight neighbouring DTFE cells surrounding the particle. The global normalisation between d μ_i and 1∕ρ_i is arbitrary.

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

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

4 VQZA super-EdS expansion: biscale partition

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

4.1 Newtonian conservation of EdS expansion prior to virialisation

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

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

We consider several subcases – a class that includes planar collapse: $${\langle \text{II}\rangle }_{D_i^{-}}=0,{\langle \text{III}\rangle }_{D_i^{-}}=0;$$(46)

a class that includes spherically symmetric expansion or collapse: $${\langle \text{II}\rangle }_{D_i^{-}}={\langle \text{I}\rangle }_{D_i^{-}}{}^2/3,{\langle \text{III}\rangle }_{D_i^{-}}={\langle \text{I}\rangle }_{D_i^{-}}{}^3/27;$$(47)

and cases where either ${\langle \text{II}\rangle }_{D_i^{-}}$ or ${\langle \text{III}\rangle }_{D_i^{-}}$ is set to zero and the non-zero invariants are set to values that give terms in Eqs. (30) and (31) of roughly similar orders of magnitude. Specific values are indicated in the caption of Fig. 1.

The planar collapse class of solutions (Buchert et al. 2000, III.D) constitute an exact solution family, which is clearly supported numerically in Fig. 1 prior to the collapse of the collapsing domain at about 6.2 Gyr. Prior to 6.2 Gyr, the global super-EdS scale factor ratio is unity to within a precision of |1 − a_eff(t)∕a_EdS| < 2.3 × 10⁻⁵, where a_eff is defined in Eq. (44). In other words, prior to virialisation, the QZA in the planar class of solutions is fully consistent with the Newtonian T³ cosmology that is normally used to calculate structure formation in the flat FLRW models (EdS and ΛCDM).

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

Fig. 1 Pre- and post-virialisation EdS-normalised scale factor evolution aeff∕aEdS for a biscale partition evolved using the VQZA as motivated in Sect. 2 and defined in Sect. 4. From bottom to top at t ≳ 7 Gyr, the arbitrarily chosen dimensionless values of the average initial invariants (at initial scale factor$a_{\text{EdS}}{}_i=0.005$) in the expanding domain are $({\langle \text{I}\rangle }_{D_i^{-}}$, ${\langle \text{II}\rangle }_{D_i^{-}}$, ${\langle \text{III}\rangle }_{D_i^{-}})=(0.01,0,0)$ (“planar” case); (0.01, 10−4∕3, 10−6∕27) (“spherical” $D^{-}$ case); (0.01, 10−4, 0); and (0.01, 0, 10−6); respectively. The sharp transitions from pre-virialisation to post-virialisation epochs are clearly visible and occur at roughly 7 Gyr to 2 Gyr, respectively. Prior to virialisation, the QZA is fullycompatible with EdS global evolution (which can be thought of as Newtonian T3 cosmology), especially in the planar case, in which the QZA is exact in the Newtonian case.

4.2 Aftervirialisation: the VQZA

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

• (i)

  divide a large volume to be studied into contiguous domains $D$;

• (ii)

  evolve the scale factor $a_D(t)$ in each domain according to the Raychaudhuri equation (Eqs. (25), (26)), where the kinematical backreaction $Q_D$ is approximated by the $Q_D$ Zel’dovich approximation (QZA) (Eqs. (11), (30), (31)), unless gravitational collapse and virialisation (iii) occur at some time t_coll;

• (iii)

  if virialisation (Sect. 3.5.6) occurs in a domain $D$, QZA evolution ceases at t = t_coll and the stable clustering approximation (Eq. (43) in the EdS case) for $a_D(t)$ is adopted for t ≥ t_coll for the domain $D$;

• (iv)

  the global scale factor evolution is estimated by cubical averaging of the per-domain scale factors $a_D(t)$ (Eq. (44)).

Is this model Newtonian or relativistic? Our aim is to model spatial expansion in a close approximation to general-relativistic behaviour. As discussed above, a fully Newtonian cosmology is difficult to define. The Buchert & Ehlers (1997) fluid model justification of T³ 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⁻⁵ 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 Λ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}}:={\left(2a_{\text{EdS}}{}^3-a_{D^{+}}^3\right)}^{1/3},$$(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: ${\ddot{a}}_{D^{-}}^{\text{Newt}}$ has to suddenly become less positive or more negative when ${\dot{a}}_{D^{+}}$ very rapidlyswitches 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 isgiven in Eq. (42) in the VQZA case. In the EdS constraint case, we again interpret ${\text{Ω}}_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},$$(49)

where the $H_{D^{+}}^2/H_{\text{FL}}^2$ factor converts between globally normalised and per-domain normalised definitions of Ω. 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 ${\text{Ω}}_R^D$ growth to ${\text{Ω}}_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.

Fig. 2 Expanding-domain scale factor evolution $a_{D^{-}}/a_{\text{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_{\text{EdS}}(t)$ (solid curves; Eq. (48)) are indistinguishable following virialisation of the overdense domain, since by assumption, the expanding domain suddenly switchesfrom hyperbolic (super-EdS) evolution to very-nearly flat (EdS) evolution ($a_{D^{-}}^{\text{Newt}}(t)/a_{\text{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.

Fig. 3 Mean curvature functional ${\text{Ω}}_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.

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 Λ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_eff that reaches unity just a little earlier than t = 15 Gyr, rather than the ΛCDM current age of the Universe estimate of t₀ = 13.80 Gyr. Nevertheless, it is premature to make detailed observational tests using the VQZA approach, so in this work we consider the ΛCDM proxy estimate of the age of the Universe as a reasonable first-order approximation.