Does relativistic cosmology software handle emergent volume evolution?

Several software packages for relativistic cosmological simulations that do not fully implement the Einstein equation have recently been developed. Two of the free-licensed ones are inhomog and gevolution. A key question is whether globally emergent volume evolution that is faster than that of a Friedmannian reference model results from the averaged effects of structure formation. Checking that emergent volume evolution is correctly modelled by the packages is thus needed. We numerically replace the software's default random realisation of initial seed fluctuations by a fluctuation of spatially constant amplitude in a simulation's initial conditions. The average volume evolution of the perturbed model should follow that of a Friedmannian expansion history that corresponds to the original Friedmannian reference solution modified by the insertion of the spatially constant perturbation. We derive the equations that convert from the perturbed reference solution to the effective solution. We find that inhomog allows emergent volume evolution correctly at first order through to the current epoch. For initial conditions with a resolution of $N = 128^3$ particles and an initial non-zero extrinsic curvature invariant $I_i = 0.001$, inhomog matches an exact Friedmannian solution to -0.0058% (Einstein-de Sitter, EdS) or -0.0033% (LCDM). We find that gevolution models the decaying mode to fair accuracy, and excludes the growing mode by construction. For $N = 128^3$ and an initial scalar potential $\Phi$ = 0.001, gevolution is accurate for the decaying mode to 0.012% (EdS) or 0.013% (LCDM). We conclude that this special case of an exact non-linear solution for a perturbed Friedmannian model provides a robust calibration for relativistic cosmological simulations.

Abstract. Several software packages for relativistic cosmological simulations that do not fully implement the Einstein equation have recently been developed. Two of the free-licensed ones are inhomog and gevolution. A key question is whether globally emergent volume evolution that is faster than that of a Friedmannian reference model results from the averaged effects of structure formation. Checking that emergent volume evolution is correctly modelled by the packages is thus needed. We numerically replace the software's default random realisation of initial seed fluctuations by a fluctuation of spatially constant amplitude in a simulation's initial conditions. The average volume evolution of the perturbed model should follow that of a Friedmannian expansion history that corresponds to the original Friedmannian reference solution modified by the insertion of the spatially constant perturbation. We derive the equations that convert from the perturbed reference solution to the effective solution. We find that inhomog allows emergent volume evolution correctly at first order through to the current epoch. For initial conditions with a resolution of N = 128 3 particles and an initial non-zero extrinsic curvature invariant I i = 0.001, inhomog matches an exact Friedmannian solution to −0.0058% (Einstein-de Sitter, EdS) or −0.0033% (ΛCDM). We find that gevolution models the decaying mode to fair accuracy, and excludes the growing mode by construction. For N = 128 3 and an initial scalar potential Φ = 0.001, gevolution is accurate for the decaying mode to 0.012% (EdS) or 0.013% (ΛCDM). We conclude that this special case of an exact non-linear solution for a perturbed Friedmannian model provides a robust calibration for relativistic cosmological simulations.
We present a method of testing whether global volume evolution beyond that of an FLRW reference model, which could emerge from the inhomogeneous evolution of 3-Ricci curvature, is accurately modelled. We introduce a numerically uniform perturbation on a flat FLRW reference model. Depending on a particular code's structure and design, this should yield the global scale factor evolution of a curved or flat (depending on the perturbation) FLRW model, if the code accurately models the dynamical role of curvature and allows the emergence of global volume evolution beyond that of the reference model. This calibration test will not work on all codes: a simulation that checks for global mathematical consistency in every successive spatial hypersection of the cosmological pseudo-manifold would detect that the topology and curvature become inconsistent in this case. However, neither inhomog nor gevolution claim to check consistency between global spatial topology and curvature on the fly. They only establish this consistency in the initial conditions. Thus, the calibration test proposed here should test emergent volume evolution without causing the code to fail due to geometrytopology constraints.
We present our method in Sect. 2.4. We rewrite the standard FLRW relations in a convenient form in Sect. 2.4.1 and we present the generic way of relating the effective model to the reference model, given a perturbation, in Sect. 2.4.2. Time matching between a pair of original and expected FLRW models is discussed in Sect. 2.4.3. The relations for the FLRW scale factor solutions for the perturbed models for the specific codes are given for inhomog (Sect. 2.4.4) and gevolution (Sect. 2.4.5). The software details of inserting the perturbations in the codes are briefly described in Sect. 2.5.1 for inhomog and in Sect. 2.5.2 for gevolution. Results are presented in Sect. 3. Appendices include a derivation of the Hamiltonian constraint (Appendix A) and the Raychaudhuri equation (Appendix B) for the gevolution case, and the relations between the EdS reference and effective solutions illustrating decaying and growing modes (Appendix C). We discuss in Sect. 4 and conclude in Sect. 5.
This package aims to be fully reproducible (Akhlaghi et al., 2021), with a source package at zenodo.6794222♯ and live † † and archived † git repositories. This paper was produced with git commit cc5ca58 of the source package, which was configured, compiled and run on a Little Endian x86 64 architecture.

Method
To test the accuracy of the packages' modelling of emergent volume evolution, we perform the simulations with and without the addition of a spatially uniform perturbation to the initial conditions of the simulation. We call the model without (or prior to) the introduction of a perturbation the 'reference' model, while the 'effective' model (equivalently, the 'emergent' model) is the reference model with the addition of a weak, uniform perturbation. We additionally use the term 'expected' model to refer to the analytically expected effective model, as opposed to the numerically simulated effective model. The analytical calculations, presented in Sect. 2.4, represent the ideal, expected behaviour of the code. As stated above, the expected behaviour of a code that handles the perturbation correctly is to switch from the FLRW reference model to an effective model whose scale factor evolution is that of a new FLRW a(t) solution. This is, in principle, a very simple test. However, we are not aware of any cosmological simulation codes that, as cosmological time increases, replace their initial FLRW reference model a(t) solution by a new FLRW a(t) solution in a large (or global) spatial domain. The approximate methods of the codes tested in this paper are not guaranteed by construction to pass this calibration test. The degree to which these codes fail the test for a global perturbation should provide a guideline to the degree to which they may be inaccurate on more local domains.
The generic method of doing this testing is explained via a thought experiment in Sect. 2.1 and given algebraically in Sections 2.2, 2.4.1, 2.4.2 and 2.4.3. The specific elements of testing inhomog are given algebraically in Sect. 2.4.4 and in terms of software and calculations in Sect. 2.5.1. Similarly, the specific elements of testing gevolution are described in Sect. 2.4.5 (algebra; including the choice of perturbation mode -a solution of the Hamiltonian and Raychaudhuri equations) and in Sect. 2.5.2 (software and calculations). In each case, we start with an FLRW reference model defined by current-epoch values of the cosmological parameters. The expected result is that the average volume evolution of the model should follow the evolution of an 'emergent' Friedmannian model with current-epoch values of the cosmological parameters that are different to those of the reference model. We consider two FLRW models: an Einstein-de Sitter universe (EdS) model and a model with a non-zero cosmological constant Λ and cold dark matter particles (ΛCDM). Together, EdS and ΛCDM span a reasonable range of flat FLRW models. Both are of special interest to cosmology: EdS models are often used as a background model in works considering the possibility of dark energy being an effect of structure formation (e.g. Roukema et al. 2013; for a wider list of papers, see Sect. 1), while ΛCDM (with an undetermined topology) is the current standard model in cosmology.

Thought experiment
This can be thought of in more detail via the following thought experiment, more complex than our simplified case. Let us suppose that a time slice, either in in the biverse model (usually called the 'background' model plus perturbations), in which the two models are referred to here as a pair of 'reference' (superscript r ) and 'effective' (superscript e ) models, respectively. The illustration assumes that the effective scale factor evolution, a e (t), deviates significantly from that of the reference model, a r (t). The two models are matched ab initio, in contrast to the typical observational matching of models, which requires models to match at the current epoch. Thus, the epochs at which the respective scale factors reach unity differ (in the illustration, t e 0 < t r 0 ). If the effective model is the ΛCDM model fit to observations (i.e., t e 0 ≈ 13.8 Gyr) and the reference model is an EdS model, and the two are matched ab initio, then there are two different Hubble constants: H e 0 := H e (t e 0 ) > H r 0 := H r (t r 0 ). Moreover, the reference model has H r (t e 0 ) < H e 0 . In this paper, by construction, the expected behaviour of the effective model is an FLRW solution.
the flow-orthogonal foliation of inhomog or the Poisson-gauge coordinate system of gevolution, is big in volume, and homogeneous except for a compact (finite volume) domain P. We assume that P is small compared to the global volume V, but big compared to the volume scale D of numerical gravitational averaging, i.e. |D| ≪ |P| ≪ |V|, where | · | gives the volume of a spatial domain. We also assume that P is itself numerically homogeneous and isotropic after averaging at the scale D, while the expansion rate or gravitational potential of P is offset from that of the global part of the time slice. If the simulation ignores the internal D-scale averaging and applies the Einstein equation at large scales, then the domain P should undergo scale factor evolution a(t) that can be approximated by an FLRW solution different to that of the reference model, apart from effects at the boundary between P and its complement that may gradually propagate to the interior of P.
To avoid these boundary effects and provide an unambiguous test method, we choose P to be the global spatial section rather than just a small part of it, i.e. P ≡ V. Moreover, we remove all numerical sources of spatial inhomogeneity. Spatial inhomogeneity remains an implicit effect represented inside of the D scale (if the code allows that), as we consider averaging to have already been performed inside the D scale. Thus, we expect FLRW-like effective a(t) behaviour after the perturbation has been applied, without any additional complications. This provides us with a robust calibration test.
The algebra and calculus for this test are, in principle, simple. However, the intuition is less simple, because of the habits of thinking in terms of the standard model. The standard model is mathematically a biverse model: a reference pseudomanifold and a perturbed pseudo-manifold mathematically co-exist as parallel universes, related to one another at early times by a diffeomorphism. The reference pseudomanifold is usually called a 'background' model. However, in order to have the property that positive and negative perturbations statistically cancel, it is usually assumed that averaging each spatial section (time slice) of the perturbed pseudo-manifold at successive cosmological times yields an alternative definition of the background that exactly matches the original reference model. Thus, the usual assumption is that the effective scale factor evolution is assumed to be identical to the reference model scale factor evolution. This makes the term 'background' ambiguous when we are considering relativistic cosmology, in which the expansion-structure decoupling hypothesis is dropped. Here, the primary question of interest is how the emergent average behaviour differs from the reference model. Since we wish to avoid ambiguity between the reference model and the average of the perturbed model, we avoid the ambiguous term 'background'.

Normalisation and Hubble constants
For both software packages, our goal is to study the emergence of super-or sub-Friedmannian volume evolution. This requires caution when normalising FLRW models. Normalisation usually refers to setting the scale factor to unity at the 'current epoch'. However, two different FLRW models that are normalised this way will, in general, have different current ages, so that at a given lookback time corresponding to an early epoch, the two models will disagree in scale factor evolution. In particular, a great enough lookback time in one model will correspond to a pre-big-bang epoch in the other. For a simulation that starts at an early time and evolves forwards in time with both a reference FLRW model and effective super-Friedmannian expansion, a different observational normalisation is needed. This was first shown independently by Rácz et al. (2017) and Roukema et al. (2017). These papers used an EdS reference model, with the aim of approximately matching ΛCDM behaviour (a proxy for observations) in the effective model. In this case, the initial scale factor growth in the reference and effective models is nearly identical as a function of a common cosmological time parameter, but later evolution of the scale factor in the effective model is super-Friedmannian. The effective model reaches a scale factor of unity after about 13.8 Gyr, while the reference model takes several more Gyr to reach a scale factor of unity. Thus, the constant values of the cosmological parameters used to parametrise the two models must differ. As found by Rácz et al. (2017) and Roukema et al. (2017), and with the usual subscript 0 for the epoch of the scale factor reaching unity, an observationally realistic EdS reference model matched this way has H 0 = 37.7 km/s/Mpc, while the effective ΛCDM model has the usual H 0 ≈ 70 km/s/Mpc. The reference and effective models reach unity scale factor after about 17.3 Gyr and 13.8 Gyr, respectively. Figure 1 illustrates this schematically. The reference r and effective e models are matched ab initio, requiring a r (t) ≈ a e (t) at early times, where the time parameter t is identical for the two models. Using this notation, the matching (Rácz et al., 2017;Roukema et al., 2017) Roukema et al. (2017), the following subsections present the physical reasoning and further notation. These make it possible to compare the scale factor evolution in the reference FLRW model to the emergent FLRW scale factor evolution that would be expected if the computer code were accurate. We ignore the contributions of radiation and neutrinos, but our test calibration could easily be extended to include them.

Software pipeline
We would like the reader to be able to straightforwardly reproduce (Rougier et al., 2017) our results, i.e., check the method and obtain identical results using precisely identified versions of high-level software, lower level software libraries, compilers and input values. We provide the full calculation pipeline for this paper using the Maneage template (Akhlaghi et al. 2021; see also Akhlaghi 2019;Infante-Sainz et al. 2020;Peper and Roukema 2021) which aims to provide a solution to the 'reproducibility crisis' in science (Peng, 2015;Baker, 2016;Fanelli, 2018). Using the free-licensed source package, available as a version of record at zenodo.6794222 and in live and archived git repositories, the aim is that the reader should be able to produce the results of this paper on any Unix-like operating system. Git commit hashes are used to identify contributing software sources and the full package itself. In particular, git commit cc5ca58 of the source package was used to produce this paper.
For each code, we test both an EdS reference model, in which emergent volume growth may match that of the ΛCDM observational proxy, and a ΛCDM reference model itself. While the amount of extra volume growth is not intended to be matched to observations in this paper, since the aim is to test the accuracy of the codes against the known FLRW solution, it is still preferable to use realistic values of cosmological density parameters Ω and Hubble constants H 0 (Ω X0 for the cosmological parameter of type X, and H 0 , with the '0' subscript indicating values when a(t) = 1). Thus, the EdS reference model constants are Ω m

Emergent FLRW model parameters: expected behaviour
The generic method of parametrisation of a perturbed simulation is common to the two codes considered here (and, we expect, to other codes). We define the following terminology. We need to calculate the cosmological parameter constants Ω X0 for the effective FLRW a(t) model. The effective model and its constants are determined by the reference model modified by the addition of a uniform perturbation. This kind of perturbation is referred to by Adamek et al. (2016b) as a 'homogeneous mode'. In Adamek et al. (2016b)'s terminology, the emergent (effective) model represents 'absorbing the homogeneous mode', and, in principle, could be used to replace the original reference model. The replacement of the reference model by the effective model would enable the simulation to continue with the homogeneous mode removed (we briefly return to this in Sect. 4.2). The 'absorption' procedure is more likely to be necessary in the case of gevolution rather than inhomog.
We first present the common steps shared by the two codes (Sections 2.4.1, 2.4.2, 2.4.3), followed by the steps specific to the two codes in Sections 2.4.4 and 2.4.5.

Standard FLRW relations rewritten
We first write the standard FLRW relations in a general form, allowing the cosmological parameters (constants) to be evaluated at an arbitrary epoch (not necessarily the current epoch). These relations represent, respectively, comoving conservation of rest-frame mass, normalisation by the expansion rate squared, constancy of dark energy modelled as a cosmological constant, and the Raychaudhuri equation: (1) The subscript n indicates an epoch chosen for normalisation; the time-dependent cosmological parameters are Ω m , the non-relativistic matter density parameter, Ω Λ , the dark energy parameter, and Ω k , the curvature parameter. The derivative of the scale factor a with respect to cosmological time t is denotedȧ. The Hubble parameter is H :=ȧ/a. In the case of gevolution, which uses conformal time in its internal calculations and output, conversions between conformal time τ and proper time t, via dt = adτ , are needed. A convenient relation is H ≡ȧ, where H is Adamek et al. (2016b)'s conformal Hubble parameter. We do not consider theȧ < 0 case, which would require the negative square root forȧ in the expression forȧ in Eq. (1).
As stated above, we retain the standard conventions of using the subscript 0 for the epoch at which the scale factor attains unity, so that a 0 := 1. The cosmological parameter values at an early epoch t i are denoted by the subscript i . Thus, we can rewrite the relations from Eq. (1) for the reference model (where the superscript r is used to differentiate reference model variables from those of the effective model): (2)

Reference FLRW to emergent FLRW conversion: common steps
Steps to convert from the reference FLRW a(t) solution to the effective FLRW solution that are common to the two codes include the following. We assume by default that the mass of the spatial section in an FLRW model with a T 3 (3-torus) spatial section, which is the case with both inhomog and gevolution, is conserved, i.e., and allow a corrective factor depending on the assumptions of the modelling approach. We assume the same cosmological constant Λ for the reference and emergent models, since we do not consider dark energy models with more exotic equations of state. This yields the conversion (4) We set the curvature parameter of the emergent model to satisfy the Hamiltonian constraint:  Figure 2: As for Fig. 1, but taking into account that the 'initial' epoch at which the perturbation is added to the reference model is t r i > 0 in the reference model, not t = 0. The FLRW parametrisation of the expected effective model needs to take into account the difference between the two FLRW time parametrisations (see Sect. 2.4.3).
As stated above, we adopt a e 0 = 1, so the cosmological time at the epoch at which this occurs is t e 0 , which, in general, is not the same as the epoch t r 0 ≈ 13.8 Gyr in the reference model (when a r 0 = 1). Three of the variables needed to successfully construct an emergent FLRW model using the above equations are a e i ,ȧ e i and Ω m e i . Because of the differences in the analytical justifications and corresponding code structure, the expressions for these two variables differ between inhomog and gevolution. These expressions are given in the respective sections (Sect. 2.4.4 and 2.4.5).

FLRW time matching
In the method presented here, our intervention in the reference model, which is expected to switch it to an effective model, is carried out at an epoch that is early, but is later than the initial singularity. We expect the expansion behaviour of the effective model to correspond to an exact FLRW a(t) solution. However, when we switch between the reference FLRW a(t) model and the effective one, the universe ages t r i and t e i , as calculated from their corresponding FLRW parameter values at the corresponding reference and effective scale factors, will, in general, differ at the switching epoch. In other words, the scheme in Fig. 1 is an analytical oversimplification; a numerical simulation does not start at the initial singularity. This is illustrated in Fig. 2. Functionally, we can write We match the two models at the switching epoch by setting At the switching epoch itself, t e −t e i = 0 = t r −t r i , even though, in general, t r i = t e i . The effective model, including Eq. (9), is only meaningful for t e ≥ t e i . The effective model universe age t e does not represent the integral of proper time since the initial singularity; it only represents the time parameter of the FLRW a(t) solution that, together with the parameters {Ω m0 , Ω Λ0 , H 0 } e , yields the expected values of a e .
For the purposes of FLRW conversions, which we calculate using cosmdist-0.3.12 ‡, we match reference and effective model 'redshifts' with scale factors in the usual way:

Emergent FLRW model parameters:
inhomog In this work, as described below in Sect. 2.5.1, inhomog uses a power spectrum of Gaussian initial conditions to generate N-body simulation initial conditions, infer initial invariants of the extrinsic curvature in spatial subdomains D, and use the relativistic Zel'dovich approximation (RZA) to calculate the evolution of the kinematical backreaction Q D (Buchert and Ostermann, 2012;Buchert et al., 2013), i.e. the Q D Zel'dovich approximation (QZA) (Roukema, 2018). The evolution of Q D yields the evolution of effective scale factors a D (t) and density parameters in each domain D. The intervention in inhomog is to replace the value of the extrinsic curvature invariant I i in each domain D by a single small nonzero global value (Sect. 2.5.1). As in Buchert et al. (2013, Sect. VI.A), Roukema (2018, eq. (20)), this should not affect the scale factor at the epoch of intervention, but should affect the expansion rate, giving the relations between the reference model and the emergent model of Mass conservation (Eq. (3)), corrected to satisfy the relativistic Zel'dovich approximation initial conditions (Buchert et al. 2013, eq. (92); Roukema 2018, eq. (26)) gives Here, at the instant of inserting the perturbation in the inhomog simulation, we have implicitly made a choice of projection between curved and flat space. In the case of the gevolution simulation (see Sect. 2.4.5 below), this projection is implicit at every ‡ https://codeberg.org/boud/cosmdist time step that evaluates quantities on a mesh. The scalar averaging approach normally does not impose any particular projection to flat space.
Using Eqs (2) for the reference model, Eqs (12) and (13) to obtain a e i anḋ a e i , Eqs (14), (4), and (5) for the effective initial Ω parameters, and Eq. (6) for the effective cosmological parameters at t e 0 , we obtain successive conversions from

Emergent FLRW model parameters: gevolution
The emergence of an effective scale factor different to that of the reference model in gevolution is described in Adamek et al. (2016b, Sect. 5.3) as the emergence of homogeneous modes.
The gevolution package works with a three-dimensional grid of particles and data fields, scattered across a lattice created by a C++ library latfield2 (Daverio et al., 2016). The biverse nature of the model is built as in standard perturbation theory, as stated above, in the Poisson gauge. In the Poisson gauge, there is an unperturbed background, which we refer to here as the reference model to avoid ambiguity, and in parallel there is a nearly identical, but perturbed, universe. However, the aim of gevolution is to evolve the perturbed model in a background-free way, i.e., independently of the reference model, using the reference model only to set the initial conditions of the simulation. Whether or not this aim is achieved is a question that the current work partly investigates.
The scale factor evolution of the reference model, a r (t) (or a r (τ )), is a solution from the FLRW family. In the code, it is evaluated independently of the evolution of structure in the simulation. If deviations from a r are small, i.e. if |a e − a r | ≪ 1 at all times, then the Maclaurin expansions used to justify approximations in the code should remain accurate. On the contrary, if there is significant non-Friedmannian effective volume evolution -as expected, for example, for the structure-induced emergent dark energy hypothesis -then the use of low-order Maclaurin expansions may become inaccurate.
Choice of projection An implicit assumption commonly made by numerical relativity codes -including the case of gevolution -is that the line element contains all the information that is necessary to infer any relevant geometrical properties. Analytically, this is correct. However, numerically, some of the geometrical information is lost when performing a projection from a patch of a curved space to a corresponding patch of a flat space, and discretising to a mesh. The choice of projection necessarily implies a selection of which properties are conserved in the projection and which are distorted; this selection affects the mesh discretisation. For example, any two-dimensional sky projection from the 2-sphere S 2 into the flat 2-plane E 2 can conserve either areas (via the Jacobian), angles, or geodesics, but not all simultaneously. Two examples of non-conservation of properties are the following. The 2-sphere S 2 cannot be tiled exactly by a Cartesian mesh of isometric quadrilaterals; though it can be tiled exactly in a non-Cartesian mesh by six isometric squares, with straight equal-length edges and corner angles of 120 • (not 90 • ; the corner angles are 'distorted' if S 2 is projected to the ordinary cube). A projection of S 2 to E 2 followed by a mesh using parallels of latitude and meridians would lead to the poles being surrounded by small triangles on S 2 , appearing as rectangles in E 2 , with edges along parallels of latitude that are geodesic on E 2 but non-geodesic on S 2 ; geodesics on S 2 would in general be non-geodesic on E 2 .
For more generic cases, such as a perturbed 3-manifold modelling a spatial slice of the Universe, the projection and discretisation will affect the numerical treatment of continuous fields in curved space. Thus, most numerical relativity simulations necessarily make a projection and associated discretisation, where the discretisation is itself associated with the numerical resolution of the simulation. In principle, the inaccuracies introduced by the implicit choices of projection and discretisation can be investigated numerically by carrying out simulations of successively higher resolution and checking for numerical convergence.
Here, we are interested in the particular case of gevolution and latfield2. The assumption of latfield2 is that changes in cell shapes are not needed for the physical problem of interest; statistical combinations of pointwise estimates ('finite differences'), are assumed to be sufficient to encode the spatially distributed geometry. Thus, for gevolution to handle curvature on the scale of discretisation into a mesh, prior to using a Fourier transform, a choice of projection has to be made. As stated above, this will result in an approximation that cannot conserve all the geometric characteristics of a cell. The particle-to-mesh interpolation (Adamek et al., 2016b, app. B) uses a combination of regular-cubical-pyramidal cloud-in-cell and nearest-grid-point averaging in coordinate space, i.e. the computationally simplest projection is adopted. In the (deprojected) curved 3-manifold itself, the cubical pyramids are, in general, irregular. For the Hamiltonian constraint, discrete averages based on sampling that appears to be uniform in flat space, rather than being intrinsically uniform on the curved 3-manifold, are used Adamek et al. (2016b, app. C, e.g. (C.1), (C.3)).
In our work, we wish to minimise changes to the internal structure of each code. In this context, we assume that the volume of a 'cubical' cell in the perturbed model is given by the cube of the side length of the cell, where the side length is given by the (perturbed) line element in the Poisson gauge. A more detailed model for future calibration tests could perform sub-cell discretisation uniformly in the curved space rather than in the flat space. The choice of cubing the Poisson-gauge line-element side length translates to an assumption that gevolution's approximation of curvature only affects the dynamics in the domain, not the geometry on sub-cell scales that is implicitly represented by the particle distribution. Alternative choices could be based on constant curvature models of the cells surrounding a point. However, none of these alternative choices of geometrical projection appear to be included in gevolution, so we do not attempt to model them.
Reference FLRW to emergent FLRW conversion The Poisson-gauge line elementused by default in the following form in v1.1 or later of gevolution (Adamek et al., 2017, eq. (6)) § -is where τ is conformal time in the reference model (defined by dτ = a −1 dt), Roman indices i, j are spatial, giving x i as the spatial coordinates. The parameters that define the model are Ψ, Φ, B i and h ij , where Ψ and Φ are scalar functions, and B i and h ij are vector and tensor modes, respectively. We use the Kronecker δ ij (δ ij = 1 when i = j, δ ij = 0 otherwise; similarly for δ ij and δ i j below). For this work, as we are modelling a perturbation effectively corresponding to the homogeneous mode, we set χ := Φ − Ψ to zero, as in Adamek et al. (2016b, Sect. 5.3) for the homogeneous mode. Although the vector and tensor modes appear in the line element, for the purposes of this test their homogeneous modes have been set to zero, so they do not contribute to the emergent FLRW a(t) model that is expected. We write˙:= d/dt and ′ := d/dτ , so that Continuing with biverse terminology, the reference model has the line element of Eq. (15) with Ψ ≡ Φ ≡ 0, while the perturbed model only has the constraint Ψ ≡ Φ (as stated above). At time t i , we now want to convert the set of model-defining parameters {a i ,ȧ i , Ω mi , Ω Λi , Ω ki , H i } r of the reference model to the set {a i ,ȧ i , Ω mi , Ω Λi , Ω ki , H i } e of an FLRW effective model. The latter should correspond to the reference model perturbed by the spatially constant potential Φ.
According to the perturbed metric (Eq. (15)), distance scales proportionately to e −Φ , and the scale factor is linearly proportional to distance. We compare the reference and perturbed model at the same epoch t i , giving To relate the density parameters, we assume that the total mass content of a given volume of the Universe is the same for the reference and emergent models, i.e. M e = M r , and use the relations a ∝ e −Φ and H :=ȧ/a. We calculate the perturbed volume by assuming that it can be described as the cube of the side length described by the perturbed metric (as described in Sect. 2.4.5, strictly valid only in flat space). Equation To evaluate the expansion rateȧ e in the emergent model, consider the thought experiment of a compact homogeneous domain D within a global homogeneous spatial hypersurface, as presented above (Sect. 2.1). We can make the conservative assumption that the line element (Eq. (15)) correctly models both D (where Φ ≡ Ψ = 0) and the § For example, git diff 0906a1e..fe756e7 gevolution.hpp in the full history. differentiating Eq. (17) and dividing by the change in effective time giveṡ To first order in Φ and dropping terms in Φ dΦ/da r i or higher, Eq. (20) reduces tȯ as given in eq. (5.11) of Adamek et al. (2016b), where Eq. (17) corresponds to the second part of (5.11). To evaluate Eq. (19), once a perturbation Φ is chosen, the Hamiltonian constraint in the gevolution context, i.e. without allowing a mean per-mesh-cell curvature term, allows a flat perturbation that is either decaying or growing, as derived in Appendix A (the Raychaudhuri equation is derived in Appendix B).
Only one of these two solutions is allowed in gevolution. This can be seen as follows. Both equation ( 15)) -are linearised in Φ ′ (they are approximations that ignore some terms). The non-linearised Einstein tensor time-time component G 00 is quadratic in Φ ′ . Thus, solving the Hamiltonian constraint for Φ ′ -with the assumption of Φ being spatially homogeneous -yields a pair of solutions. Algebraically, these correspond to the choice of either a plus or a minus sign in solving a quadratic equation (see Appendix A for details). Approximation by linearisation ignores one of the two solutions of the quadratic. For example, if we write a general quadratic equation as αξ 2 +βξ+γ = 0, then approximating |ξ| 2 ≪ |(β/α)ξ| gives the solution ξ ≈ −γ/β. In reality, the quadratic equation has (for a positive discriminant) two real solutions. In the case of gevolution, the calculations in the code match the linearised solution (Eq. (9) of Adamek et al. (2017)) that omits terms in (Φ ′ ) 2 or of higher order. Thus, the linearisation effectively chooses the sign of the solution to the non-linearised equation (i.e. the solution of Eq. (A.8) by linearisation is set to the negative branch of the square root in Eq. (A.9)). This choice corresponds to the decaying mode of the Raychaudhuri equation, which relates Φ ′′ to Φ ′ and Φ (Eq. (B.3)). Without linearisation, the growing mode is not excluded. We return to the physical meaning of the growing mode, and whether the approximation adopted by gevolution could imply that the code fails to correctly model structure evolution, in Sect. 4.2 below.
Thus, to check that gevolution behaves as expected, using Eq. (2) for the reference model, the decaying branch in Eq. (A.9) (negative sign in ±), Eqs (17) and (19) to obtain a e i andȧ e i , Eqs (18), (4), and (5) for the effective initial Ω parameters, and Eq. (6) for the effective cosmological parameters at t e 0 , we obtain successive conversions from , Ω k0 , H 0 } e , i.e. we have the parameters of the expected FLRW scale factor solution.

Introduction of perturbations in the simulation codes
In this section, we provide the detailed method of how to apply the analytical methods described in section Sect. 2.4. For both codes, this was done with the aim of minimally disrupting the existing code when inserting the spatially uniform perturbation.  Figure 5: Scale factor ratios a e Φ i =0 /a r for a standard gevolution simulation (dotted line); a e Φ i =0.001 /a r for a gevolution simulation with an initial potential perturbation Φ i = 0.001 (solid line); and the expected FLRW evolution of a e FLRW /a r (as defined in Sect. 2.4.5; dashed line), for the decaying mode (negative square root in Eq. A.9). Top: EdS reference model; bottom: ΛCDM. The sudden increase in value for the perturbed gevolution simulation (solid curve) at an early time step is the effect of the Φ i offset being inserted into the simulation, and this datapoint is removed from residuals to improve readability. The expected and actual curves match prior to this initial peak in the curve, by construction. Plain-text data for this figure and Fig. 6

inhomog
The inhomog package (Roukema, 2018;Roukema and Ostrowski, 2019) is a free-licensed C program and library that implements the relativistic Zel'dovich approximation (RZA) for the evolution of the kinematical backreaction Q D (Buchert and Ostermann, 2012;Buchert et al., 2013), i.e. the Q D Zel'dovich approximation (QZA) (Roukema, 2018). The core functions of the library calculate the domain-averaged volume evolution represented by a D (t) := a D (t), using the RZA estimate for Q D (t), given an input of the initial values of the three invariants I, II, III of the extrinsic curvature in a comoving spatial domain in a flow-orthogonal time foliation.  Figure 6: Scale factor ratios a e Φ i =0 /a r for a standard gevolution simulation, as for Fig. 6, but with a negative initial potential perturbation; a e Φ i =−0.001 /a r for a gevolution simulation with an initial potential perturbation Φ i = −0.001 (solid line); and the expected FLRW decaying mode evolution of a e FLRW /a r (as defined in Sect. 2.4.5; dashed line). Top: EdS reference model; bottom: ΛCDM. The behaviour mirrors that for the simulation with a positive Φ i . Likewise, the initial point is removed from the residual plot for readability. This is used in the N-body simulation context by generating N-body simulation initial conditions from a power spectrum using mpgrafic (Prunet et al., 2008) and using the ramses-scalav extension of ramses (Teyssier, 2002) as a front end. The front end ramses-scalav reads the mpgrafic realisation of initial conditions, calls dtfe (Schaap and van de Weygaert, 2000;van de Weygaert and Schaap, 2009;Cautun and van de Weygaert, 2011;Kennel, 2004) to infer the Newtonian velocity gradients for non-overlapping subdomains that cover the full simulation volume, calculates approximations for I i , II i , III i , and calls inhomog to calculate a D (t) for each of the subdomains. The global mean volume evolution is calculated from the individual a D (t) evolution of the individual subdomains. For the purposes of this paper, we avoid the issue of subdomain collapse and virialisation, since these are unnecessary for testing To see if a global perturbation inserted into inhomog-0.1.10-1da3bed yields the expected FLRW expansion rate while minimally intervening in the use of inhomog via the ramses-scalav front end, the routine amr dtfe interface is modified to read the value of I i from a plain-text file, and II i and III i are set to zero.
Our expectation is that this should yield an emergent FLRW solution for the expansion history different to that of the reference model. The properties of the expected emergent solution are calculated using a shell script, using cosmdist-0.3.12 where appropriate, independently from the inhomog code and its ramses-scalav front end. Equations (14), (4), and (5)  , Ω Λ0 , Ω k0 , H 0 } r .

gevolution
The gevolution software is a free-licensed C++ cosmological simulation code developed by Adamek et al. (2016a,b) aiming to be relativistically more accurate than traditional codes. It has gained popularity as a tool for making FLRW observational predictions and for testing beyond-FLRW cosmological scenarios (e.g. Hassani et al., 2019;Reverberi and Daverio, 2019). The input parameter file for the simulation includes physical parameters such as Ω m r and H 0 r , the comoving spatial side length of the 3-torus domain being simulated, the initial distribution of the particles prior to their perturbation by a power spectrum and starting and output redshifts. The user selects desired outputs such as particle distributions at specified redshifts or power spectra of different perturbations in the Poisson gauge line element (Eq. (15)). The code is highly modular. The numerical core of gevolution is latfield2, a library containing tools for solving equations of classical lattice field theory. The code runs efficiently either on a personal computer of limited computational power or on a large computer cluster, depending primarily on the particle resolution.
The goal of testing gevolution is to determine whether the code is currently (gevolution-1.2-0404c0b) able to handle a global perturbation while intervening only minimally in the internal code structure. In order to test the effect of a spatially uniform global perturbation alone, we first provide an input power spectrum of perturbations with zero amplitude. In the initial cycle of generating the particle distribution, we modify the function generateIC basic() in the file ic basic.hpp, adding a step that allows the insertion of a uniform perturbation in the scalar gravitational potential Φ. The optional insertion of the perturbation and the perturbation's amplitude are controlled in the initial settings file. The gevolution code outputs the reference model scale factors and other parameters. A patch file with our exact set of modifications is provided in the reproducibility package of this paper and is archived at Software Heritage ¶. When we refer to gevolution-1.2-0404c0b, strictly speaking we mean gevolution-1.2-0404c0b as modified using this patch file.
We again calculate the expected FLRW scale factor evolution using a shell script independent from the simulation code being tested, apart from the shared use of Φ i andΦ i . The reasoning is presented above in Sect. 2.4.5. The conversion of initial scale factors to obtain a e i is given in Eq. (17). The effective matter density Ω m e i at the epoch of intervention is given in Eq. (18). The latter requiresȧ e i , which is provided by Eq. (19), provided that we know the value ofΦ at this epoch. The choice between ¶ https://archive.softwareheritage.org/ swh:1:cnt:043550246bb1c9c26664a0dbebedb334693033b1 growing and decaying modes forΦ, corresponding to the positive and negative square roots in Eq. A.9, respectively (converted with Eq. 16), is set to the decaying mode, since gevolution uses the linearised solution (see Sect. 2.4.5), which corresponds to the decaying perturbation.
As explained above (Sect. 1, Sect. 2.4.5), gevolution and its backend latfield2 do not check the consistency between global spatial topology and curvature except in the sense that this is assumed to be implemented in the initial conditions. Moreover, parameters representing the geometry per mesh cell are not stored. Thus, the insertion of a spatially constant perturbation in Φ cannot correspond to offsetting the mean curvature in each mesh cell; instead it can only lead to a modified spatial section that is again of zero spatial curvature. As shown in Appendix A and mentioned above, this allows two possible flat perturbation modes satisfying the Hamiltonian constraint (choice of Φ ′ , or equivalently,Φ). We compare the numerical evolution of gevolution with the decaying solution, which is clearly the one chosen by the code. The relations in Sections 2.4.1, 2.4.2 and 2.4.3 are then again used to obtain the emergent FLRW model constants, and cosmdist is used to generate the expected emergent scale factor evolution.

inhomog
We ran the relativistic Zel'dovich approximation evolution package inhomog-0.1.10-1da3bed from N-body initial conditions, with a resolution of N = 128 3 particles, fundamental domain comoving (reference model) size 40 Mpc/h, and the ramses resolution parameter levelmax EdS HTwEight = 10, with and without a perturbation as described in Sect. 2.5.1. We compared the numerically measured scale factor evolution with the expected FLRW model by defining a residual δ and a fractional deviation ǫ, defined δ := a e /a r − a e FLRW /a r ǫ := a e /a r − a e FLRW /a r a e FLRW /a r .
The fractional deviation ǫ describes the numerical deviation as a fraction of the excess expansion that should be generated by the perturbation. Figure 3 (for an initial density perturbation I i = 0.001) and Fig. 4 (for I i = −0.001) show that the code behaves as expected to fair numerical accuracy, correctly taking into account the effect of spatial curvature and matching the emergent Friedmannian model. For the EdS reference model, the deviation from the expected model increases in amplitude with time, with the fractional deviation at the final output time being ǫ = −0.0058% for I i = 0.001, and ǫ = −0.0059% for I i = −0.001. The corresponding ΛCDM reference model deviations are ǫ = −0.0033% and ǫ = −0.0034%, respectively.
It is clear that inhomog allows the emergence of scale factor evolution beyond that of the reference model, i.e. it allows the emergence of the homogeneous mode. Table 1 lists the scale factors at the final reference model epoch, and the reference and effective FLRW model current epoch constants. In all eight cases, a dynamical non-null curvature parameter Ω k e 0 appears in the expected FLRW model. These are correctly handled by inhomog, while gevolution instead returns to the reference model.

gevolution
As in the inhomog case, we ran independent simulations for the EdS and ΛCDM models for gevolution-1.2-0404c0b, in each case with and without two cases of an initial perturbation: Φ i = 0.001 and Φ i = −0.001. The simulations' particle resolution is N = 128 3 ; the ΛCDM simulation was run with Ω ref Λ0 = 0.6911; and the fundamental domain comoving (reference model) size was 40 Mpc/h.
In the cases with Φ i = 0.001, the perturbed spatial section is smaller than in the reference model, corresponding to a spatial section overdense in comparison to the background model. However, this is insufficient to determine the expected spatial curvature, since density has to be normalised by the squared expansion rate in order to obtain the emergent density parameter, as shown in Eq. (18). Moreover, Ω Λ e i is also modified by the perturbation (Eq. (4)), leaving Ω k e i as the complement (Eq. (5)). The evolution of the potentialΦ enters into these relations via Eq. (19). The linearisation in gevolution implies the decaying mode of Φ, which for Φ > 0 yieldsΦ < 0. The formulae in Sect. 2.4 calculated using cosmdist-0.3.12 give a flat effective model (Ω k e 0 = 0) in these two cases, as shown in Table 1. These formulae also give the expected behaviour of the a e /a r ratio, which was calculated independently of gevolution, using cosmdist. Figures 5 and 6 show the behaviour of gevolution with a perturbation in comparison to the expected FLRW behaviour equivalent to that of the decaying mode. We find that gevolution provides a fair approximation to the expected FLRW behaviour, yielding decaying modes that return to the original FLRW reference solutions, though with a small difference between the numerical and analytical solutions.
Numerically, maximal gevolution's deviations |ǫ| from the expected behaviour for the EdS model are ǫ = 0.012% for Φ i = 0.001 and ǫ = 0.012% for Φ i = −0.001; and for the ΛCDM model the deviations are ǫ = 0.013% for Φ i = 0.001 and ǫ = 0.013% for Φ i = −0.001. Those quantities exclude the initial spike seen in residual plots as an artefact of perturbation introduction.
In summary, we find that gevolution-1.2-0404c0b models the decaying mode well, but does not allow the growing mode of the quadratic equation in Φ ′ (we discuss this in Sect. 4.2).

inhomog
As shown in Figs 3 and 4, inhomog-0.1.10-1da3bed is reasonably accurate in modelling the volume evolution implied by a density perturbation. Nevertheless, the residual errors δ (Eq. (22)) appear to grow systematically, especially at late times. The residual errors at the current epoch of the reference model are of the order of the square of the change in scale factor, i.e. δ ∼ [a e (a r = 1)] 2 . This would appear to suggest a secondorder error in the QZA model. Tightening of the integration accuracy parameters did not significantly affect this, although the visually identical nature of the patterns of the residuals for positive versus negative fluctuations for a fixed reference model was strengthened with improved numerical precision. The pair of residual curves for a fixed reference model differ numerically, despite the visual resemblance.
Thus, although the QZA model is in several ways non-linear in the sense of traditional perturbation theory, and is exact in some special cases (Buchert et al., 2013, Sect. V.A, V.B), it appears that for the evolution of average volume in a perturbed FLRW model, further improvements are likely to be needed, especially if the additional scale factor component is of the order of 16% (Roukema et al., 2017) rather than 1%.

gevolution
In Adamek et al. (2016b, Sect. 5.3), the possibility that a homogeneous mode could emerge is discussed, with suggestions of how it could be handled. However, here we find that out of the two possible branches given in Eq. (A.9), i.e. the solution to the quadratic in Φ ′ , the linearisation given in eqs (2.9) of Adamek et al. (2016b) and eq. (9) of Adamek et al. (2017) restricts the evolution of a homogeneous mode to the decaying branch. Equation (A.8) shows that this linearisation should be a good approximation when |Φ ′ | ≪ |ȧ r | =ȧ r , where the equality follows from only considering an expanding universe (not a contracting or stationary one). For a decaying perturbation, if |Φ ′ | also decreases with increasing time, then an initial perturbation satisfying the condition is likely to lead to a solution in which the linearisation remains accurate.
The growing mode is the situation where an emergent model leads to Equation (19), with the subscript i reinserted, shows that for small |Φ|, this implies thatȧ e i < 0. As a global effect, this would represent a universe starting to collapse, which cannot correspond to any EdS solution nor to the ΛCDM solution. The possible relevance of the growing mode could be that it is an effective model for a locally collapsing spatial region, since average flatness in the scalar averaging sense does not strictly forbid gravitational turnaround from occurring in this way. In practice, in an EdS or ΛCDM reference model evolving from a cosmologically typical initial power spectrum of perturbations, it is rare for a region that is on average flat to slow in its expansion and collapse past its turnaround epoch. This was shown using the relativistic Zel'dovich approximation (Roukema and Ostrowski, 2019). Thus, the relevance of the omitted growing mode in the context of our test for homogeneous modes is modest, since the mode is likely to be rare in practice. It would be more likely to be relevant in the more general case of an inhomogeneous perturbation, in which case local gravitational collapse is expected. In that case, the Hamiltonian constraint will still include squares of Φ ′ , so the selection of the growing mode branch of the solutions would appear to be necessary for full modelling of structure formation within the gevolution approach. An appropriate condition would have to be derived analytically and evaluated numerically to decide when to switch between the two modes. This would be particularly interesting for the formation of cosmic voids, whose role in effective volume expansion remains a key question of study.
Nevertheless, returning to homogeneous modes, if the gevolution convention of setting Φ = Ψ for the homogeneous mode were dropped, then more possibilities would be allowed. An algebraic demonstration showing how both decaying and growing solutions are possible for uniform perturbations is given in App. Appendix C for the EdS case, for illustration. Both the decaying and growing solutions are EdS solutions, i.e. they expand eternally. Thus, allowing Φ = Ψ for the homogeneous mode could be expected to yield an extension of gevolution that would enable the flat growing mode to be studied.
Given that gevolution is free-licensed, it is likely that we, the original authors, or other cosmologists will extend the code to allow for the emergence of the growing homogeneous mode. This would quite likely involve the inclusion of more second-order terms in the algebraic justification and in the coding itself, together with allowing Φ = Ψ for the homogeneous mode. This future extension would not be completely trivial, since the code does not currently include unit tests or regression tests, and the built-in dependence on an implicitly fixed FLRW reference model of a(t) could yield nonlinearities that are numerically hard to control. Nevertheless, this would be a promising followup of the current analysis.
Another option, but possibly unstable numerically, would be to use the formulae in Sect. 2.4 to calculate a new effective model at each time step and use this instead of the reference model in all appropriate formulae. The source package for this paper includes an optional patch (not applied by default) that illustrates how this approach might be attempted.
One of the issues not considered here is the treatment of radiation and other particle species (eg. non-cold dark matter). The gevolution package has data structures for modelling these, but connecting these to the issue of global non-zero curvature significantly increases the complexity of the problem, so this is left for future work. Dependence of inhomog and gevolution accuracy parameter ǫ (defined in Eq. (22)) as a function of simulation resolution N 3 .

Role of simulation resolution
In typical cosmological simulations, increasing the spatial resolution generally leads to numerical convergence for the calculation of quantities of interest. However, in the current work, our calibration test is that of a homogeneous perturbation. Increasing the resolution in the context of a homogeneous perturbation, for a fixed comoving reference volume and identical starting epoch, should lead to statistically equivalent numerical effects in each of a greater number of smaller regions, and is unlikely to show numerical convergence. We ran our full set of simulations for the resolutions listed in Table 2, which shows the dependence of ǫ on the particle resolution N. The disagreement between the expected and numerically calculated scale factor evolution is worse for higher N in the case of inhomog. This is consistent with the automatic choice of starting epoch set by mpgrafic. For a fixed comoving reference domain size, 40 Mpc/h in the current case, higher spatial resolution typically implies that a fixed amplitude of perturbations can be started in the linear regime at an earlier epoch. For example, the N = 32 3 EdS simulations are started at a r i = 0.060, while the N = 64 3 and N = 128 3 EdS simulations are started earlier, at a r i = 0.044 and a r i = 0.033, respectively. The earlier starting epochs give longer evolution times, but start from the same initial homogeneous perturbation size. Since the expected and numerically calculated effective scale factors diverge as time increases (see Figs 3, 4), the longer evolution time yields greater errors.
In the case of gevolution, which does not couple the starting epoch to the resolution, increasing N has a negligible effect on ǫ.

Numerical representation of geometrical information
As explained above (Sect. 2.4.5) in relation to gevolution, the full representation of geometrical information in a numerical simulation is not trivial. Most cosmological codes conceptually use three-dimensional meshes composed of polyhedral cells (typically cubes) as their numerical backbone, or particles that discretely represent a smooth mass distribution. These do not completely represent the geometrical information present in the line element of the spacetime's metric, and instead use projections. Given that not only spacetime curvature in general, but also spatial curvature in particular, is a crucial element of a general-relativistic spacetime, the issue of geometric projections cannot be avoided if the treatment of curvature in a code is to be correctly understood and analysed. Codes claiming to be fully relativistic will need to correctly handle these geometrical issues in order to study their influence on the evolution of the Universe.

Free-licensed cosmology software ecosystem
The results of this paper illustrate how the use of free-licensed, modular codes written in C (cosmdist and inhomog) or C++ (gevolution) can be easily and efficiently checked against each other. These should be easy to use in checks or extensions of the Einstein Toolkit (Bentivegna and Bruni, 2016;Macpherson et al., 2017), part of a much bigger free-licensed relativistic code project, or of cosmograph (Giblin et al., 2017).

Conclusions
The question of finding exact inhomogeneous cosmological solutions that are appropriate for the 3-torus spatial topology common to N-body and other cosmological simulations, in order to calibrate them, has long been a challenge.
The curved FLRW models do not provide a full solution to this challenge, but if we consider only the dynamical role of the curved models in volume evolution rather than their directly geometrical and topological role, as has been the tradition in the N-body simulation community, then it is clear that they can provide a robust calibration for relativistic cosmological simulations. The QZA code inhomog passes this calibration test to first order in the change in the scale factor, and the higher order deviations hint at possibilities for improvements.
The flat FLRW models, appropriately recalibrated to match the numerical code, have shown their utility in illustrating the risk of linearisation. In the case of gevolution, we found that out of the two solutions of a quadratic solution, the choice of solution is hardwired. The choice made by gevolution is the decaying mode that is implied by linearisation in the first conformal-time derivative of the potential, Φ ,τ . Thus, inhomogeneous initial conditions in a cosmological simulation that lead to an emergent homogeneous mode that is growing rather than decaying are likely to be excluded by construction in gevolution. Work extending gevolution to allow the growing mode might provide a useful extension of the code.
The checks presented here have been done aiming at long-term reproducibility; the git commit hash of the source package of this paper, for those wishing to confirm that attempted reproductions are using the same version of the source package, is cc5ca58.

Data Availability Statement
The full reproducibility package for this paper, including the choice of simulation parameters ('input data'), is available at zenodo.6794222 + and in live * and archived♯ git repositories. The specific version of the source package used to produce this paper can be identified by its git commit hash cc5ca58.
where Φ ,i denotes a spatial derivative in the direction i. The extrinsic curvature K ij can be calculated from the spatial line element and the lapse and shift functions. Since in our example we assume that vector and tensor perturbations are zero, our shift covector B i is zero; the lapse function is N = a r e Ψ . (A.3) Using K ij = −1/(2N) ∂/∂τ γ ij [e.g., Gourgoulhon (2007, (4.63); γ ij is the spatial part of the line element)], and Eq. (16), gives making the switch between proper to coordinate time derivatives for convenience (a r′ = a rȧr ). The mixed, contravariant, and scalar extrinsic curvature are K i j = γ ik K kj = (a r ) −1 (Φ ′ −ȧ r )e −Ψ δ i j , K ij = γ ik K k j = (a r ) −3 (Φ ′ −ȧ r )e 2Φ−Ψ δ ij , and The choice of the sign before the square root and the setting of Ψ = Φ may constrain the behaviour of the solution. The negative square root together with Ψ = Φ give a decaying mode, in which over-and under-densities weaken with time. The positive square root, possibly together with Ψ = Φ, may give a growing mode. The possible physical relevance of this case is discussed in Sect. 4.2. A solution that linearises Eq. (A.8) in terms of Φ ′ (ignores the Φ ′2 term) has only the decaying solution. Examples of exact decaying and growing modes in the EdS case are given in Sect. Appendix C.
Appendix B. Evolution of the perturbation: Raychaudhuri equation As stated above, the shift function in our case is zero and the lapse has no spatial dependence, so Thus, the Raychaudhuri equation (Gourgoulhon, 2007, (4.64), cf (4.79 for the N = 1 case)) becomes When finding the FLRW effective scale factor history that corresponds to a reference model combined with a uniform perturbation, the EdS solutions illustrate how both decaying and growing modes are possible. Although the EdS model for a fixed initial singularity has only one free parameter (either the Hubble-Lemaître constant or, equivalently, the current universe age), the effective model is not required to have its pre-switching evolution match that of the reference model, since the effective model is only meaningful after the switching point. In other words, the implied big bang moment of the effective model is a virtual value rather than a physical constraint. Thus, a second free parameter is provided by the time offset (see Eq. (9)). The reference and effective models are After setting t r 0 for the reference model and t r i for the insertion of the perturbation, the two parameters t e 0 and t e i both need to be chosen to provide the effective model. A decaying perturbation occurs when t r 0 = t e 0 , since, as t → ∞, |a e − a r | a r → t e i − t r