paper.tex
\NeedsTeXFormat{LaTeX2e}
\documentclass[usenatbib]{mnras} % mnras.cls is available in e.g. Debian stable (9.10+)
%% (OPTIONAL) CONVENIENCE VARIABLE: Only relevant when you use Maneage's
%% '\includetikz' macro to build plots/figures within LaTeX using TikZ or
%% PGFPlots. If so, when the Figure files (PDFs) are already built, you can
%% avoid TikZ or PGFPlots completely by commenting/removing the definition
%% of '\makepdf' below. This is useful when you don't want to slow-down a
%% LaTeX-only build of the project (for example this happens when you run
%% './project make dist'). See the definition of '\includetikz' in
%% `tex/preamble-pgfplots.tex' for more.
%\newcommand{\makepdf}{}
%% Import the necessary preambles.
\input{tex/src/preamble-mnras-free.tex}
\input{tex/src/preamble-style.tex}
\input{tex/build/macros/project.tex}
\input{tex/src/preamble-pgfplots.tex}
\input{tex/src/preamble-journals-astro.tex}
%\input{tex/src/preamble-biblatex.tex}
\input{tex/src/preamble-shortcuts-astro.tex}
%% Details for people willing to re-LaTeX the source from ArXiv.
\newtoggle{showDetailedDerivations}
%\toggletrue{showDetailedDerivations}
\togglefalse{showDetailedDerivations} %% Typesetters/Proofreaders: please leave this set to \togglefalse.
\title[Relativistic cosmology and volume evolution]{\projecttitle}
\author[Borkowska \& Roukema]{Justyna Borkowska$^{1}$ \&
Boudewijn F. Roukema$^{1,2}$
\\
$^1$ Institute of Astronomy, Faculty of Physics,
Astronomy and Informatics, Nicolaus Copernicus
University, Grudziadzka 5, 87-100 Toru\'n, Poland
\\ %email: {\tt name@astro.uni.torun.pl} % avoid this for antispam reasons if we can
$^2$ Univ Lyon, Ens de Lyon, Univ Lyon1, CNRS, Centre de
Recherche Astrophysique de Lyon UMR5574, F--69007, Lyon,
France}
%Mohammad Akhlaghi$^{4,5}$
%$^4$ Instituto de Astrof\'isica de Canarias, Calle V\'ia L\'actea s/n, 38205 La Laguna, Tenerife, Spain.
%$^5$ Departamento de Astrof\'isica, Universidad de La Laguna, Avenida Astrof\'isico Francisco S\'anchez s/n, 38200 La Laguna, Tenerife, Spain.
\date{Accepted \ldots Received \ldots; in original form \ldots}
\pubyear{2021}
\begin{document}%\layout
\maketitle
%% mnras "...normally of not more than 250 words for Main Journal papers or 200 words for Letters"
%% Abstract
\begin{abstract}
%WCWC word count
%context
{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 {\inhomogname} and {\gevolutionname}.
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.}
%aims
{Checking that emergent volume evolution is correctly modelled by the packages is thus needed.}
%method
{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.}
%results
%solution of $\GevolutionAccuracyPercentEdSvalue$\% (EdS) or $\GevolutionAccuracyPercentLCDMvalue$\% ($\Lambda$CDM)
{We find that {\inhomogname} allows emergent volume evolution correctly at first order through to the current epoch.
For initial conditions with a resolution of $N=\Ncrootinhomogvalue^3$ particles and an initial non-zero extrinsic curvature invariant $\initial{\invI} = \InitPerturbInhomogPlusvalue$, {\inhomogname} matches an exact Friedmannian solution to $\InhomogAccuracyPercentEdSzdzzovalue$\% (Einstein--de~Sitter, EdS) or $\InhomogAccuracyPercentLCDMzdzzovalue$\% ($\Lambda$CDM).
%The Poisson-gauge formalism of {\gevolutionname} allows both growing and decaying modes for emergent volume evolution, and selects the decaying mode.
We find that {\gevolutionname} models the decaying mode to fair accuracy, and excludes the growing mode by construction.
For $N=\Ncrootgevolutionvalue^3$ and an initial scalar potential $\Phi = \InitPerturbGevnPlusvalue$, {\gevolutionname} is accurate for the decaying mode to $\GevolutionAccuracyPercentMaxEdSzdzzovalue$\% (EdS) or $\GevolutionAccuracyPercentMaxLCDMzdzzovalue$\% ($\Lambda$CDM).
}
%conclusions
{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.}
%WCWC word count
\end{abstract}
\begin{keywords}
methods: numerical, galaxies: evolution, cosmology: dark matter
\end{keywords}
\section{Introduction} \label{s-intro}
We introduce a method of investigating whether the recently-developed free-software--licensed packages for relativistic cosmological simulations are able to model global volume evolution that may emerge as the averaged effect of inhomogeneous non-linear spatial (3-Ricci) curvature evolution.
We use simulation packages that are both licensed as free software and distributed publicly, so the whole of the scientific community is able to test the software, modify it, introduce changes to improve scientific accuracy, add functionality or correct bugs, and freely redistribute these modified versions.
We apply our method to two fast codes that take shortcuts from full solutions of the Einstein equation.
These are {\gevolutionname} \citep{AdamekDDK15,AdamekDDK16code}\footnote{Full history (v1.0, v1.1, v1.2): \url{https://codeberg.org/boud/gevolution}}, which uses second-order expansions within the Poisson gauge and has become quite popular; and {\inhomogname} \citep{Roukema17silvir,RO19flatness}\footnote{\url{https://codeberg.org/boud/inhomog}}, which uses the relativistic Zel'dovich approximation to carry out very fast calculations that in terms of standard linear perturbation theory of structure formation are effectively non-perturbative calculations \citep{BKS00,BuchRZA1,BuchRZA2}.
The {\inhomogname} package is based on the relativistic Zel'dovich approximation (RZA) for the evolution of the kinematical backreaction $\CQ_\CD$ \citep{BuchRZA1,BuchRZA2}, i.e. the $\CQ_\CD$ Zel'dovich approximation (QZA) \citep{Roukema17silvir}.
This approach uses a flow-orthogonal foliation of spacetime and a Cartan co-frame formalism.
This allows for an approach that is effectively \enquote*{non-linear} in the context of linear perturbation theory.
For example, {\inhomogname} and the QZA can be used to derive the curvature and density \enquote*{$\Omega$} parameter dependence on the expansion rate during dark-matter halo collapse without needing restrictive simplifications such as Birkhoff's law or spherical symmetry \citep{RO19flatness,VignBuch19,Ostrowski19}.
While {\gevolutionname} is designed within a general-relativistic framework -- a Poisson gauge line element -- it linearises some expressions derived within the Poisson gauge and truncates others to second order.
The ability of {\inhomogname} and {\gevolutionname} to incorporate deviations from a flat Friedmann--Lema{\^i}tre--Robertson--Walker (FLRW) reference model will be investigated in this work.
Several packages that, in principle, solve the Einstein equation fully include the following.
The Einstein Toolkit (ET) approach of \citet{BentivegnaBruni15} and \citet{Macpherson17}'s {\sc FLRWSolver}\footnote{\url{https://github.com/hayleyjm/FLRWSolver_public}} uses the BSSN formalism \citep{ShibataNakamura95,BaumgarteShapiro99}.
In terms of the formulae justifying their calculations, the ET calculations aim to be fully relativistic, without truncating Maclaurin expansions.
Another free-licensed package for relativistic cosmological simulations using the BSSN formalism is {\sc cosmograph} \citep{GiblinMS17}\footnote{\url{https://github.com/cwru-pat/cosmograph}}.
The {\sc gramses} package \citep{Barrera2019} is not yet public, but should be distributable with a CECILL-compatible licence if it is derived from the {\sc ramses} code \citep{Teyssier02}, since the latter is licensed under CECILL.
Accurate modelling of non-linear 3-Ricci (spatial) curvature is a key challenge that all of these software packages need to meet if they are to be used to test the hypothesis that structure formation leads to recently emerging negative mean curvature that physically explains \enquote*{dark energy} \citep{Rasanen06negemergcurv,
NambuTanimoto05,KaiNambu07, Rasanen08peakmodel,
Larena09template,Chiesa14Larenatest, WiegBuch10,BuchRas12,
Wiltshire09timescape,DuleyWilt13,NazerW15CMB, ROB13,
Barbosa16viscos, BolejkoC10SwissSzek, LRasSzybka13,
Roukema17silvir, Sussman15LTBpertGIC,Chirinos16LTBav,
Bolejko17negcurv,Bolejko17styro,
Kras81kevolving,Kras82kevolving,Kras83kevolving,Stichel16,Stichel18,
Coley10scalars,KasparSvitek14Cartan,RaczDobos16}.
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 {\inhomogname} nor {\gevolutionname} 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 geometry--topology constraints.
We present our method in \SSS\ref{s-algebra-generic}.
We rewrite the standard FLRW relations in a convenient form in \SSS\ref{s-algebra-FLRW-rewritten} and we present the generic way of relating the effective model to the reference model, given a perturbation, in \SSS\ref{s-algebra-ref-to-eff}.
Time matching between a pair of original and expected FLRW models is discussed in \SSS{\ref{s-time-matching}.
The relations for the FLRW scale factor solutions for the perturbed models for the specific codes are given for {\inhomogname} (\SSS\ref{s-algebra-inhomog}) and {\gevolutionname} (\SSS\ref{s-algebra-gev}).
The software details of inserting the perturbations in the codes are briefly described in \SSS\ref{s-method-inhomog} for {\inhomogname} and in \SSS\ref{s-method-gev} for {\gevolutionname}.
Results are presented in \SSS\ref{s-results}.
We discuss in \SSS\ref{s-discuss} and conclude in \SSS\ref{s-conclu}.
This package aims to be fully reproducible \citep{Akhlaghi2020maneage}, with a source package at {\projectzenodoid}\footnote{{\projectzenodohrefShowURL}} and live\footnote{\projectgitrepository} and archived\footnote{\projectgitrepositoryarchived} {\sc git} repositories.
This paper was produced with git commit {\projectversion} of the source package, which was configured, compiled and run on a {\machinebyteorder} {\machinearchitecture} architecture.
\section{Method} \label{s-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.
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 \SSS\ref{s-thought-exper} and given algebraically in \SSSS\ref{s-norm-Hubble}, \ref{s-algebra-FLRW-rewritten}, \ref{s-algebra-ref-to-eff} and \ref{s-time-matching}.
The specific elements of testing {\inhomogname} are given algebraically in \SSS\ref{s-algebra-inhomog} and in terms of software and calculations in \SSS\ref{s-method-inhomog}.
Similarly, the specific elements of testing {\gevolutionname} are described in \SSS\ref{s-algebra-gev} (algebra; including the choice of perturbation mode -- a solution of the Hamiltonian and Raychaudhuri equations) and in \SSS\ref{s-method-gev} (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 \enquote*{emergent} Friedmannian model with current-epoch values of the cosmological parameters that are different to those of the reference model.
\begin{figure}
%% xfig -specialtext -latexfonts -startlatexFont default tex/xfig/gev_calibrate.fig % xfig is optional, and not currently shipped in Maneage
\input{tex/xfig/gev_calibrate.pstricks}
\caption{Schematic figure showing the relation between scale factor evolution in the biverse model (usually called the \enquote*{background} model plus perturbations), in which the two models are referred to here as a pair of \enquote*{reference} (superscript $\refmodel{{}}$) and \enquote*{effective} (superscript $\effmodel{{}}$) models, respectively.
The illustration assumes that the effective scale factor evolution, $\effmodel{a}(t)$, deviates significantly from that of the reference model, $\refmodel{a}(t)$.
The two models are matched {\em 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, $\effmodel{t}_0 < \refmodel{t}_0$).
If the effective model is the $\Lambda$CDM model fit to observations (i.e., $\effmodel{t}_0 \approx 13.8$~Gyr) and the reference model is an EdS model, and the two are matched {\em ab initio}, then there are two different Hubble constants: $\effmodel{H}_0 := \effmodel{H}(\effmodel{t}_0) > \refmodel{H}_0 := \refmodel{H}(\refmodel{t}_0)$.
Moreover, the reference model has $\refmodel{H}(\effmodel{t}_0) < \effmodel{H}_0$.
In this paper, by construction, the expected behaviour of the effective model is an FLRW solution.
\label{f-biverse-model}}
\end{figure}
\subsection{Thought experiment} \label{s-thought-exper}
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 the flow-orthogonal foliation of {\inhomogname} or the Poisson-gauge coordinate system of {\gevolutionname}, is big in volume, and homogeneous except for a compact domain $\CP$.
We assume that $\CP$ is small compared to the global volume $\CV$, but big compared to the volume scale $\CD$ of numerical gravitational averaging, i.e. $|\CD| \ll |\CP| \ll |\CV|$, where $|\cdot|$ gives the volume of a spatial domain.
We also assume that $\CP$ is itself numerically homogeneous and isotropic after averaging at the scale $\CD$, while the expansion rate or gravitational potential of $\CP$ is offset from that of the global part of the time slice.
If the simulation ignores the internal $\CD$-scale averaging and applies the Einstein equation at large scales, then the domain $\CP$ 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 $\CP$ and its complement that may gradually propagate to the interior of $\CP$.
To avoid these boundary effects and provide an unambiguous test method, we choose $\CP$ to be the global spatial section rather than just a small part of it, i.e. $\CP \equiv \CV$.
Moreover, we remove all numerical sources of spatial inhomogeneity.
Spatial inhomogeneity remains an implicit effect represented inside of the $\CD$ scale (if the code allows that), as we consider averaging to have already been performed inside the $\CD$ 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 pseudo-manifold and a perturbed pseudo-manifold mathematically co-exist as parallel universes, related to one another at early times by a diffeomorphism.
The reference pseudo-manifold is usually called a \enquote*{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 \enquote*{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 \enquote*{background}.
\subsection{Normalisation and Hubble constants} \label{s-norm-Hubble}
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 \enquote*{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 \citet{RaczDobos16} and \citet{RMBO16Hbg1}.
These papers used an EdS reference model, with the aim of approximately matching $\Lambda$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 \citet{RaczDobos16} and \citet{RMBO16Hbg1}, 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 $\Lambda$CDM model has the usual $H_0 \approx 70$~km/s/Mpc.
The reference and effective models reach unity scale factor after about 17.3~Gyr and 13.8~Gyr, respectively.
Figure~\ref{f-biverse-model} illustrates this schematically.
The reference $\refmodel{{}}$ and effective~$\effmodel{{}}$ models are matched {\em ab initio}, requiring $\refmodel{a}(t) \approx \effmodel{a}(t)$ at early times, where the time parameter $t$ is identical for the two models.
Using this notation, the matching \citep{RaczDobos16,RMBO16Hbg1} between a reference EdS model and an effective $\Lambda$CDM model yields an effective model that reaches $\effmodel{H}(\effmodel{t}_0) \approx \effmodel{H}(13.8~\mathrm{Gyr}) \approx 70$~km/s/Mpc for a reference EdS model with $\refmodel{H}_0 \approx \refmodel{H}(17.3~\mathrm{Gyr}) \approx 37.7$~km/s/Mpc.
This multiplicity of Hubble constants is contrary to the usual intuition of \enquote*{the} Hubble constant.
Confusion can arise if the different constants are not clearly distinguished.
Thus, in a similar spirit to \citet{RMBO16Hbg1}, 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.
\subsection{Emergent FLRW model parameters: expected behaviour} \label{s-algebra-generic}
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 ($\Omega_{X0}$ for the cosmological parameter of type $X$, and $H_0$, with the \enquote*{0} subscript indicating values when $a(t)=1$), 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 \citet{AdamekDDK16code} as a \enquote*{homogeneous mode}.
In \citet{AdamekDDK16code}'s terminology, the emergent (effective) model represents \enquote*{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 \SSS\ref{s-discuss-gevolution}).
The \enquote*{absorption} procedure is more likely to be necessary in the case of {\gevolutionname} rather than {\inhomogname}.
We first present the common steps shared by the two codes (\SSSS\ref{s-algebra-FLRW-rewritten}, \ref{s-algebra-ref-to-eff}, \ref{s-time-matching}), followed by the steps specific to the two codes in \SSSS\ref{s-algebra-inhomog} and \ref{s-algebra-gev}.
\subsubsection{Standard FLRW relations rewritten} \label{s-algebra-FLRW-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:
\begin{align}
\Omm &= \nrmlstn{\Omm} \,\nrmlstn{H}^2
{a}^{-1} {\dot{a}}^{-2} {\nrmlstn{a}}^{3}
%\left(\frac{a}{\nrmlstn{a}}\right)^{-1} \left(\frac{\dot{a}}{\nrmlstn{a}}\right)^{-2}
\nonumber\\
\Omk &= \nrmlstn{\Omk} \,\nrmlstn{H}^2 \left(\frac{\dot{a}}{\nrmlstn{a}}\right)^{-2}
\nonumber\\
\OmLam &= \nrmlstn{\OmLam} \,\nrmlstn{H}^2
\left(\frac{a}{\dot{a}}\right)^{2}
% \left(\frac{a}{\nrmlstn{a}}\right)^{2} \left(\frac{\dot{a}}{\nrmlstn{a}}\right)^{-2}
\nonumber\\
\dot{a} &= \nrmlstn{\dot{a}}
\sqrt{ \nrmlstn{\Omm} \left(\frac{a}{\nrmlstn{a}}\right)^{-1} +
\nrmlstn{\Omk} +
\nrmlstn{\OmLam} \left(\frac{a}{\nrmlstn{a}}\right)^{2} }
\,.
\label{e-FLRW-generalised}
\end{align}
The subscript $\nrmlstn{{}}$ indicates an epoch chosen for normalisation; the time-dependent cosmological parameters are $\Omm$, the non-relativistic matter density parameter, $\OmLam$, the dark energy parameter, and $\Omk$, the curvature parameter.
The derivative of the scale factor $a$ with respect to cosmological time $t$ is denoted $\dot{a}$.
The Hubble parameter is $H := \dot{a}/a$.
In the case of {\gevolutionname}, which uses conformal time in its internal calculations and output, conversions between conformal time $\tau$ and proper time $t$, via $\diffd t = a \diffd \tau$, are needed.
A convenient relation is ${\cal H} \equiv \dot{a}$, where ${\cal H}$ is \citet{AdamekDDK16code}'s conformal Hubble parameter.
We do not consider the $\dot{a} < 0$ case, which would require the negative square root for $\dot{a}$ in the expression for $\dot{a}$ in Eq.~\eqref{e-FLRW-generalised}.
As stated above, we retain the standard conventions of using the subscript $\currepoch{{}}$ for the epoch at which the scale factor attains unity, so that $\currepoch{a} :=1$.
The cosmological parameter values at an early epoch $\initial{t}$ are denoted by the subscript $\initial{{}}$.
Thus, we can rewrite the relations from Eq.~\eqref{e-FLRW-generalised} for the reference model (where the superscript $\refmodel{{}}$ is used to differentiate reference model variables from those of the effective model):
\begin{align}
\initial{\refmodel{\Omm}} &= \currepoch{\refmodel{\Omm}} \,{{\currepoch{\refmodel{H}}}^{2}} {\left({\initial{\refmodel{a}}}\right)^{-1}} {\left({\,\initial{\refmodel{\dot{a}}}}\right)^{-2}}
\nonumber\\
\initial{\refmodel{\Omk}} &= \currepoch{\refmodel{\Omk}} \,{{\currepoch{\refmodel{H}}}^{2}} {\left({\,\initial{{\refmodel{\dot{a}}}}}\right)^{-2}}
\nonumber\\
\initial{\refmodel{\OmLam}} &= \currepoch{\refmodel{\OmLam}} \,{{\currepoch{\refmodel{H}}}^{2}} {{\initial{\refmodel{a}}}^{2}} {\left({\,\initial{\refmodel{\dot{a}}}}\right)^{-2}}
\nonumber\\
\initial{\refmodel{\dot{a}}} &= \currepoch{\refmodel{H}}
\sqrt{ \currepoch{\refmodel{\Omm}} {{\initial{\refmodel{a}}}^{-1}} +
\currepoch{\refmodel{\Omk}} +
\currepoch{\refmodel{\OmLam}} {{\initial{\refmodel{a}}}^{2}} }
\,.
\label{e-FLRW-initial-ref}
\end{align}
\subsubsection{Reference FLRW to emergent FLRW conversion: common steps} \label{s-algebra-ref-to-eff}
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 $\mT^3$ (3-torus) spatial section, which is the case with both {\inhomogname} and {\gevolutionname}, is conserved, i.e.,
\begin{align}
M = \Omm H^2 a^3 \,,
\label{e-massconservation}
\end{align}
and allow a corrective factor depending on the assumptions of the modelling approach.
We assume the same cosmological constant $\Lambda$ for the reference and emergent models, since we do not consider dark energy models with more exotic equations of state.
This yields the conversion
\begin{align}
\initial{\effmodel{\OmLam}} &:= \refmodel{\OmLam}
\left( \frac{\initial{\effmodel{a}}}{\initial{\refmodel{a}}} \right)^2
\left( \frac{\initial{\refmodel{\dot{a}}}}{\initial{\effmodel{\dot{a}}}} \right)^2
\label{e-OmegaLami-eff} \,.
\end{align}
\begin{figure}
%% xfig -specialtext -latexfonts -startlatexFont default tex/xfig/gev_timeshift.fig % xfig is optional, and not currently shipped in Maneage
\input{tex/xfig/gev_timeshift.pstricks}
\caption{As for Fig.~\ref{f-biverse-model}, but taking into account that the \enquote*{initial} epoch at which the perturbation is added to the reference model is $\initial{\refmodel{t}} > 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 \SSS\protect\ref{s-time-matching}).
\label{f-timeshift}}
\end{figure}
We set the curvature parameter of the emergent model to satisfy the Hamiltonian constraint:
\begin{align}
\initial{\effmodel{\Omk}} &:= 1 - \initial{\effmodel{\Omm}} - \initial{\effmodel{\OmLam}}
\,.
\label{e-Omegaki-eff}
\end{align}
Provided that we have expressions for $\initial{\effmodel{a}}$, $\initial{\effmodel{\Omm}}$, and $\initial{\effmodel{\dot{a}}}$, we can now use the standard FLRW relations of Eq.~\eqref{e-FLRW-generalised} to write
\begin{align}
\currepoch{\effmodel{\Omm}} &= \initial{\effmodel{\Omm}}
\,\initial{\effmodel{H}}^2
\left({\currepoch{\effmodel{a}}}\right)^{-1}
\left({\currepoch{\effmodel{\dot{a}}}}\right)^{-2}
{\initial{\effmodel{a}}}^{3}
%\left(\frac{\currepoch{\effmodel{a}}}{\initial{\effmodel{a}}}\right)^{-1} \left(\frac{\currepoch{\effmodel{\dot{a}}}}{\initial{\effmodel{a}}}\right)^{-2}
\nonumber\\
\currepoch{\effmodel{\Omk}} &= \initial{\effmodel{\Omk}} \,\initial{\effmodel{H}}^2 \left(\frac{\currepoch{\effmodel{\dot{a}}}}{\initial{\effmodel{a}}}\right)^{-2}
\nonumber\\
\currepoch{\effmodel{\OmLam}} &= \initial{\effmodel{\OmLam}} \,\initial{\effmodel{H}}^2
%\left(\frac{\currepoch{\effmodel{a}}}{\initial{\effmodel{a}}}\right)^{2} \left(\frac{\currepoch{\effmodel{\dot{a}}}}{\initial{\effmodel{a}}}\right)^{-2}
\left(\frac{\currepoch{\effmodel{a}}}{\currepoch{\effmodel{\dot{a}}}}\right)^{2}
\nonumber\\
\currepoch{\effmodel{\dot{a}}} &= \initial{\effmodel{\dot{a}}}
\sqrt{ \initial{\effmodel{\Omm}} \left(\frac{\currepoch{\effmodel{a}}}{\initial{\effmodel{a}}}\right)^{-1} +
\initial{\effmodel{\Omk}} +
\initial{\effmodel{\OmLam}} \left(\frac{\currepoch{\effmodel{a}}}{\initial{\effmodel{a}}}\right)^{2} }
\,.
\label{e-FLRW-effective-currepoch}
\end{align}
As stated above, we adopt $\currepoch{\effmodel{a}} = 1$, so the cosmological time at the epoch at which this occurs is $\currepoch{\effmodel{t}}$, which, in general, is not the same as the epoch $\currepoch{\refmodel{t}} \approx 13.8$~Gyr in the reference model (when $\currepoch{\refmodel{a}} = 1$).
Three of the variables needed to successfully construct an emergent FLRW model using the above equations are $\initial{\effmodel{a}}$, $\initial{\effmodel{\dot{a}}}$ and $\initial{\effmodel{\Omm}}$.
Because of the differences in the analytical justifications and corresponding code structure, the expressions for these two variables differ between {\inhomogname} and {\gevolutionname}.
These expressions are given in the respective sections (\SSS\ref{s-algebra-inhomog} and \ref{s-algebra-gev}).
\subsubsection{FLRW time matching} \label{s-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 $\initial{\refmodel{t}}$ and $\initial{\effmodel{t}}$, 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.~\ref{f-biverse-model} is an analytical oversimplification; a numerical simulation does not start at the initial singularity.
This is illustrated in Fig.~\ref{f-timeshift}.
Functionally, we can write
\begin{align}
\initial{\refmodel{t}} &=
t_{\mathrm{FLRW}} \left(\initial{\refmodel{a}},
\refmodel{\currepoch{\Omm}},
\refmodel{\currepoch{\OmLam}},
\refmodel{\currepoch{H}} \right)
\\
\initial{\effmodel{t}} &=
t_{\mathrm{FLRW}} \left(\initial{\effmodel{a}},
\effmodel{\currepoch{\Omm}},
\effmodel{\currepoch{\OmLam}},
\effmodel{\currepoch{H}} \right) \,.
\end{align}
We match the two models at the switching epoch by setting
\begin{align}
\effmodel{t} - \initial{\effmodel{t}} &= \refmodel{t} - \initial{\refmodel{t}} \,.
\label{e-match-t}
\end{align}
At the switching epoch itself,
%\begin{align}
$\effmodel{t} - \initial{\effmodel{t}} = 0 =\refmodel{t} - \initial{\refmodel{t}}\,,$
%\end{align}
even though, in general, $\initial{\refmodel{t}} \ne \initial{\effmodel{t}}$.
The effective model, including Eq.~\eqref{e-match-t}, is only meaningful for $\effmodel{t} \ge \initial{\effmodel{t}}$.
The effective model universe age $\effmodel{t}$ 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 $\effmodel{\{ \currepoch{\Omm}, \currepoch{\OmLam}, \currepoch{H}\}}$, yields the expected values of $\effmodel{a}$.
For the purposes of FLRW conversions, which we calculate using {\cosmdistname}-{\cosmdistversion}\footnote{\url{https://codeberg.org/boud/cosmdist}}, we match reference and effective model \enquote*{redshifts} with scale factors in the usual way:
\begin{align}
\refmodel{a}(t) & =: \left(1+\refmodel{z}\right)^{-1} \\
\effmodel{a}(t) & =: \left(1+\effmodel{z}\right)^{-1} \,.
\end{align}
\subsubsection{Emergent FLRW model parameters: {\inhomogname}} \label{s-algebra-inhomog}
The intervention in the {\inhomogname} code is to insert a small non-zero global value of the extrinsic curvature invariant $\initavinvI$ (\SSS\ref{s-method-inhomog}).
As in \citet[][{\SSS}VI.A]{BuchRZA2}, \citet[][eq.~(20)]{Roukema17silvir}, this should not affect the scale factor at the epoch of intervention, but should affect the expansion rate, giving the relations for expected emergent model of
\begin{align}
\initial{\effmodel{a}} &= \initial{\refmodel{a}}
\label{e-ai-eff-inhomog}
\\
\initial{\effmodel{\dot{a}}} &=
\initial{\refmodel{\dot{a}}} \left(1 + \frac{\initavinvI}{3} \right)
\,.
\label{e-adoti-eff-inhomog}
\end{align}
Mass conservation (Eq.~\eqref{e-massconservation}), corrected to satisfy the relativistic Zel'dovich approximation initial conditions (\citealt{BuchRZA2}, eq.~(92); \citealt{Roukema17silvir}, eq.~(26)) gives
\begin{align}
\initial{\effmodel{\Omm}} &= \initial{\refmodel{\Omm}}
\left(1 + \frac{\initavinvI}{3} \right)^{-2}
\left(1 - \initavinvI\right) \,.
\label{e-Omegami-eff-inhomog}
\end{align}
Here, at the instant of inserting the perturbation in the {\inhomogname} simulation, we have implicitly made a choice of projection between curved and flat space.
In the case of the {\gevolutionname} simulation (see \SSS\ref{s-method-projection} below), this projection is implicit at every time step that evaluates quantities on a mesh.
The scalar averaging approach normally does not impose any particular projection to flat space.
Using Eqs~\eqref{e-FLRW-initial-ref} for the reference model, Eqs~\eqref{e-ai-eff-inhomog} and \eqref{e-adoti-eff-inhomog} to obtain $\initial{\effmodel{a}}$ and $\initial{\effmodel{\dot{a}}}$, Eqs~\eqref{e-Omegami-eff-inhomog}, \eqref{e-OmegaLami-eff}, and \eqref{e-Omegaki-eff} for the effective initial $\Omega$ parameters, and Eq.~\eqref{e-FLRW-effective-currepoch} for the effective cosmological parameters at $\currepoch{\effmodel{t}}$ , we obtain successive conversions from
$\refmodel{\{\currepoch{a}, \currepoch{\dot{a}}, \currepoch{\Omm}, \currepoch{\OmLam}, \currepoch{\Omk}, \currepoch{H}\}}$ to
$\refmodel{\{\initial{a}, \initial{\dot{a}}, \initial{\Omm}, \initial{\OmLam}, \initial{\Omk}, \initial{H}\}}$ to
$\effmodel{\{\initial{a}, \initial{\dot{a}}, \initial{\Omm}, \initial{\OmLam}, \initial{\Omk}, \initial{H}\}}$ to
$\effmodel{\{\currepoch{a}, \currepoch{\dot{a}}, \currepoch{\Omm}, \currepoch{\OmLam}, \currepoch{\Omk}, \currepoch{H}\}}$.
\subsubsection{Emergent FLRW model parameters: {\gevolutionname}} \label{s-algebra-gev}
The emergence of an effective scale factor different to that of the reference model in {\gevolutionname} is described in \citet[][\SSS5.3]{AdamekDDK16code} as the emergence of homogeneous modes.
\paragraph{Choice of projection} \label{s-method-projection}
The {\gevolutionname} package works with a three-dimensional grid of particles and data fields, scattered across a lattice created by a C++ library {\LATfieldtwo} \citep{Daverio16latfield2}.
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 the model of {\gevolutionname} is to evolve the perturbed model in a background-free way, i.e., independently of the reference model.
Whether or not this aim is achieved is a question which the current work partly investigates.
An obvious role of the reference-model dependence is that the scale factor evolution of the reference model, $\refmodel{a}(t)$ (or $\refmodel{a}(\tau)$), is an FLRW solution that in the code is evaluated independent of the evolution of structure in the simulation.
If deviations from $\refmodel{a}$ are small, i.e. if $\vert \effmodel{a} - \refmodel{a}\vert \ll 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 the structure-induced emergent dark energy hypothesis, then the use of low-order Maclaurin expansions may become inaccurate.
An implicit assumption that seems to be commonly made by numerical relativity codes -- including the case of {\gevolutionname} investigated here -- is that the line element contains all the information that is necessary to infer any relevant geometrical properties.
Analytically, this is correct.
However, if spatial curvature is to be accounted for geometrically, rather than just dynamically, then discretisation at a given mesh resolution requires (in the simplest approaches) a projection at each mesh cell from a coordinate patch of the perturbed 3-manifold onto a patch of a constant-curvature 3-manifold.
A choice of projection is used in the typical averaging method made at the scale of discretisation of a curved-space simulation.
In other words, at any time snapshot of the simulation, choosing a projection from each small spatial cell of the perturbed model -- which, in this case, is topologically a fragment of the $\mT^3$ hypersurface assumed by {\gevolutionname} -- onto a corresponding cubical cell of the Cartesian grid, is a simple numerical strategy, which does appear to be implicitly adopted in some parts of {\gevolutionname}.
Since the formal existence of the perturbed spatial section as a Riemannian 3-manifold requires the existence of a smooth invertible map -- a diffeomorphism -- from each curved-space \enquote*{cubical} cell of the perturbed manifold into the Euclidean 3-manifold $\mE^3$, one of the simplest options is to choose one of the many projections to a flat cubical cell rather than to a curved one.
The choice of projection necessarily implies a selection of which properties are conserved in the projection and which are distorted.
This can be illustrated using the example of the homogeneous two-dimensional manifolds $\mS^2$, the 2-sphere, and the infinite flat plane, $\mE^2$.
Any two-dimensional sky projection from $\mS^2$ into $\mE^2$ can conserve either areas (the Jacobian), angles, or geodesics, but not all simultaneously.
A similar situation holds in the three-dimensional case: the hypersphere $\mS^3$ (or a hyperbolic infinite space $\mH^3$) and the flat three-dimensional space $\mE^3$ are not isometric to one another.
A single projection from a positive-area patch of $\mS^3$ or $\mH^3$ into $\mE^3$ cannot conserve all geometrical properties simultaneously.
The distortion of geometrical properties affects (for a generic case rather than only our calibration case) the numerical treatment of continuous fields in curved space.
Here, we are interested in the particular case of {\gevolutionname} and {\LATfieldtwo}.
In the current form of these codes, when conceptually changing from flat-space mesh cell shapes to curved space mesh cell shapes, the geometrical change is not fully modelled numerically.
To illustrate the general issue in the two-dimensional case, $\mS^2$ cannot be tiled exactly by a Cartesian mesh of isometric quadrilaterals (it can be tiled exactly by six isometric squares, with straight equal-length edges and corner angles of 120$^\circ$, in a non-Cartesian mesh corresponding to the six faces of a single cube).
For example, attempting to define a mesh on $\mS^2$ using meridians and parallels of latitude would lead to the poles being surrounded by triangles, of which the edge on a parallel of constant latitude would be non-geodesic.
In other words, the triangles would each have one non-straight edge.
The assumption of {\LATfieldtwo} is that changes in cell shapes are not needed for the physical problem of interest; statistical combinations of pointwise estimates (\enquote*{finite differences}), with one point per mesh cell at a given resolution, are assumed to be sufficient to encode the spatially distributed geometry.
Thus, for {\gevolutionname} 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.
The projection will necessarily be an approximation that can preserve only one or a few, but not all, of the geometric characteristics of a cell.
The particle-to-mesh interpolation \citep[][app.~B]{AdamekDDK16code} 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 \citet[][app.~C, e.g. (C.1), (C.3)]{AdamekDDK16code}.
The inaccuracies introduced by the implicit choices of projections can be investigated numerically by carrying out simulations of successively higher resolution and checking for numerical convergence.
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 \enquote*{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 {\gevolutionname}'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 {\gevolutionname}, so we do not attempt to model them.
\paragraph{Reference FLRW to emergent FLRW conversion} \label{s-algebra-ref-to-eff-old}
The Poisson gauge line element used by default in v1.1 or later of {\gevolutionname} \citep[][eq.~(6)]{Adamek2017radiation}\footnote{For example, {\tt git diff 0906a1e..fe756e7 gevolution.hpp} in the \href{https://codeberg.org/boud/gevolution}{full history}.} is
\begin{align}
\diffd s^2 &= a^2 \left[-{\eeuler}^{2\Psi} \diffd \tau^2 -
2B_i \diffd x^i \diffd \tau +
\delta_{ij} \eeuler^{-2\Phi} \diffd x^i \diffd x^j + %\right. \nonumber \\ %&\left.\quad\quad
h_{ij} \diffd x^i \diffd x^j \right]
\,,
\label{e-line-element}
\end{align}
where $\tau$ is conformal time in the reference model ($\diffd t = a \diffd \tau$), Roman indices $i,j$ are spatial, giving $x^i$ as the spatial coordinates, $\Psi$ and $\Phi$ are scalar functions, $B_i$ and $h_{ij}$ are vector and tensor modes, respectively, $\delta_{ij}$ is the Kronecker $\delta$ ($\delta_{ij} = 1$ when $i=j$, $\delta_{ij} = 0$ otherwise; similarly for $\delta^{ij}$ and $\delta^i_j$ below).
For the homogeneous mode, $\chi := \Phi - \Psi$ is set to be equal by \citet{AdamekDDK16code}.
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 $\dot{\ } := \diffd/\diffd t$ and ${\ }' := \diffd/\diffd \tau$, so that
\begin{align}
\frac{\diffd\Phi}{\diffd t} =: \dotPhi = \frac{\Phi'}{\refmodel{a}} := \frac{\diffd \Phi}{\refmodel{a} \diffd\tau} \,.
\label{e-phi-t-tau-relation}
\end{align}
Continuing with biverse terminology for the standard model, the reference model has the line element of Eq.~\eqref{e-line-element} with $\Psi \equiv \Phi \equiv 0$, and the perturbed model only has the constraint $\Psi \equiv \Phi$.
At time $\initial{t}$, we now want to convert the set of parameters $\{\initial{a}, \initial{\dot{a}}, \initial{\Omm}, \initial{\OmLam}, \initial{\Omk}, \initial{H}\}^r$ of the reference model to the set $\{\initial{a}, \initial{\dot{a}}, \initial{\Omm}, \initial{\OmLam}, \initial{\Omk}, \initial{H}\}^e$ of an FLRW effective model that should correspond to the reference model perturbed by the spatially constant potential $\Phi$.
According to the perturbed FLRW metric (Eq.~\eqref{e-line-element}), distance scales proportionately to $\eeuler^{-\Phi}$, and the scale factor is linearly proportional to distance.
We compare the reference and perturbed model at the same epoch $\initial{t},$ giving
\begin{align}
\initial{\effmodel{a}} &:= \effmodel{a}(\initial{t}) = \initial{\refmodel{a}} \, \eeuler^{-\Phi} \,.
\label{e-ai-eff-gev}
\end{align}
Adopting the same mass between the reference and emergent models, i.e., $\effmodel{M} = \refmodel{M}$, using $a \propto \eeuler^{-\Phi}$ and $H := \dot{a} / a$, and adopting the cube-of-perturbed-side-length projection for calculating the perturbed volume as discussed in \SSS\ref{s-method-projection}, Eq.~\eqref{e-massconservation} yields
\begin{align}
\initial{\effmodel{\Omm}} &:= \initial{\refmodel{\Omm}}
\,\eeuler^{\Phi}
\left( \frac{\initial{\refmodel{\dot{a}}}}{\initial{\effmodel{\dot{a}}}} \right)^2
\,.
\label{e-Omegami-eff-gev}
\end{align}
\begin{figure}
\centering
\includegraphics[width=0.9\columnwidth]{inhomog_curvtest_EdS0_001}
\includegraphics[width=0.9\columnwidth]{inhomog_curvtest_LCDM0_001}
\caption{Emergent-to-reference model scale factor ratios
$\effmodel{a}/\refmodel{a}$ for an {\inhomogname} simulation without addition of a perturbation (dotted line) and
$\effmodel{a}/\refmodel{a}$ for an {\inhomogname} simulation with an initial invariant perturbation $\initial{\invI}=\InitPerturbInhomogPlusvalue$ (solid line).
The expected FLRW evolution $\expected{\effmodel{a}}/\refmodel{a}$ is also plotted as a dashed line, but is visually indistinguishable from the numerically measured value.
The residual $\delta$ (Eq.~\protect\eqref{e-resid-inhomog}) shows the numerical differences between these two curves.
{\em Top:} EdS reference model; {\em bottom:} $\Lambda$CDM.
Plain-text data for this figure and Fig.~\protect\ref{f-flat-FLRW-Iinv-test-minus} at \href{\projectzenodofilesbase/inhomog_scale_factors_EdS.dat}{\projectzenodoid/inhomog\_scale\_factors\_EdS.dat}, \href{\projectzenodofilesbase/inhomog_scale_factors_LCDM.dat}{\projectzenodoid/inhomog\_scale\_factors\_LCDM.dat}.
\label{f-flat-FLRW-Iinv-test-plus}}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=0.9\columnwidth]{inhomog_curvtest_EdSm0_001}
\includegraphics[width=0.9\columnwidth]{inhomog_curvtest_LCDMm0_001}
\caption{Emergent-to-reference model scale factor ratios, as for Fig.~\protect\ref{f-flat-FLRW-Iinv-test-plus}, for $\initial{\invI}=\InitPerturbInhomogMinusvalue$.
{\em Top:} EdS reference model; {\em bottom:} $\Lambda$CDM.
\label{f-flat-FLRW-Iinv-test-minus}}
\end{figure}
To evaluate the expansion rate in the emergent model, consider the thought experiment of a compact homogeneous domain $\CD$ within a global homogeneous spatial hypersurface, as presented above (\SSS\ref{s-thought-exper}).
We can make the conservative assumption that the line element (Eq.~\eqref{e-line-element}) correctly models both $\CD$ (where $\Phi \equiv \Psi \neq 0$) and the complement of $\CD$ (where $\Phi \equiv \Psi = 0$).
The {\gevolutionname} model has no information about the different expansion rates of $\CD$ and the complement of $\CD$ in this thought experiment, except via $\Phi$ in the spatial direction and $\Psi$ in the time direction.
The time part of the line element in the perturbed model differs from that in the reference model by a factor of $\eeuler^\Psi = \eeuler^\Phi$.
Returning to our reference and effective models, differentiating Eq.~\eqref{e-ai-eff-gev} and dividing by the change in effective time gives
\begin{align}
\initial{\effmodel{\dot{a}}}
&= \left( \initial{\refmodel{\dot{a}}} \, \eeuler^{-\Phi} + \initial{\refmodel{a}}\left(-\eeuler^{-\Phi}\right) \dotPhi \right) \eeuler^{-\Phi}
\nonumber \\
&= \eeuler^{-2\Phi} \initial{\refmodel{\dot{a}}} \left( 1 - \frac{\initial{\refmodel{a}}}{\initial{\refmodel{\dot{a}}}} \dotPhi \right)
\label{e-adoti-eff-gev}
\\
&= \eeuler^{-2\Phi} \initial{\refmodel{\dot{a}}} \left( 1 - \initial{\refmodel{a}} \initial{\left(\frac{\diffd \Phi}{\diffd \refmodel{a}}\right)} \right)
\,.
\label{e-adoti-eff-dphidaref-gev}
\end{align}
To first order in $\Phi$ and dropping terms in $\Phi \,{\diffd \Phi}/{\diffd \initial{\refmodel{a}}}$ or higher, Eq.~\eqref{e-adoti-eff-dphidaref-gev} reduces to
\begin{align}
\initial{\effmodel{\dot{a}}} &\approx
\initial{\refmodel{\dot{a}}}
\left(1 - 2\Phi - \initial{\refmodel{a}}
\initial{\left(\frac{\diffd \Phi}{\diffd \refmodel{a}}\right)} \right) \,,
\end{align}
as given in eq.~(5.11) of \citet{AdamekDDK16code}, where Eq.~\eqref{e-ai-eff-gev} corresponds to the second part of (5.11).
To evaluate Eq.~\eqref{e-adoti-eff-gev}, once a perturbation $\Phi$ is chosen, the Hamiltonian constraint in the {\gevolutionname} 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 App.~\ref{app-phiprime} (the Raychaudhuri equation is derived in App.~\ref{app-Raychaudhuri}).
Only one of these two solutions is allowed in {\gevolutionname}. This can be seen as follows.
Equations~(2.9) of \citet{AdamekDDK16code} for the older {\gevolutionname} form of the line element, and eq.~(9) of \citet{Adamek2017radiation} for the form used here, Eq.~\eqref{e-line-element}, are both linearised in $\Phi'$, while the non-linearised Einstein tensor time-time component $G^{00}$ is quadratic in $\Phi'$.
A solution to the non-linearised equation requires either the use of additional information about the system, or a selection of one of the two solutions.
The solution to the linearisation of $\alpha \xi^2 + \beta\xi + \gamma = 0$ corresponds to the limit of the negative square root in the solution of the quadratic, i.e. $\xi = \lim_{\alpha\rightarrow0} [(-\beta - \sqrt{\beta^2 - 4\alpha\gamma})/(2\alpha)] = -\gamma/\beta$.
The {\gevolutionname} code matches the linearised solution, Eq.~(9) of \citet{Adamek2017radiation}, i.e. the solution of Eq.~\eqref{e-Hamiltonian-constraint-quadratic} by linearisation is set to the negative square root branch of Eq.~\eqref{e-phiprime-solution-pair}.
This corresponds to the decaying mode of the Raychaudhuri equation, Eq.~\eqref{e-Raychaudhuri-Phi}, relating $\Phi''$ to $\Phi'$ and $\Phi$.
We return to the physical meaning of the growing mode in \SSS\ref{s-discuss-gevolution} below.
Thus, to check that {\gevolutionname} behaves as expected, using Eq.~\eqref{e-FLRW-initial-ref} for the reference model, the decaying branch in Eq.~\eqref{e-phiprime-solution-pair} (negative sign in $\pm$), Eqs~\eqref{e-ai-eff-gev} and \eqref{e-adoti-eff-gev} to obtain $\initial{\effmodel{a}}$ and $\initial{\effmodel{\dot{a}}}$, Eqs~\eqref{e-Omegami-eff-gev}, \eqref{e-OmegaLami-eff}, and \eqref{e-Omegaki-eff} for the effective initial $\Omega$ parameters, and Eq.~\eqref{e-FLRW-effective-currepoch} for the effective cosmological parameters at $\currepoch{\effmodel{t}}$, we obtain successive conversions from
$\refmodel{\{\currepoch{a}, \currepoch{\dot{a}}, \currepoch{\Omm}, \currepoch{\OmLam}, \currepoch{\Omk}, \currepoch{H}\}}$ to
$\refmodel{\{\initial{a}, \initial{\dot{a}}, \initial{\Omm}, \initial{\OmLam}, \initial{\Omk}, \initial{H}\}}$ to
$\effmodel{\{\initial{a}, \initial{\dot{a}}, \initial{\Omm}, \initial{\OmLam}, \initial{\Omk}, \initial{H}\}}$ to
$\effmodel{\{\currepoch{a}, \currepoch{\dot{a}}, \currepoch{\Omm}, \currepoch{\OmLam}, \currepoch{\Omk}, \currepoch{H}\}}$, i.e. we have the parameters of the expected FLRW scale factor solution.
\begin{figure}
\centering
\includegraphics[width=0.9\columnwidth]{gevn_curvtest_EdSm0_001}
\includegraphics[width=0.9\columnwidth]{gevn_curvtest_LCDMm0_001}
\caption{Scale factor ratios
$\effmodel{a}_{\initial{\Phi}=0}/\refmodel{a}$ for a standard {\gevolutionname} simulation (dotted line);
$\effmodel{a}_{\initial{\Phi}=\gevnparamZeromodeamplitudeEdSPlusvalue}/\refmodel{a}$ for a {\gevolutionname} simulation with an initial potential perturbation $\initial{\Phi}=\gevnparamZeromodeamplitudeEdSPlusvalue$ (solid line); and
the expected FLRW evolution of $\effmodel{a}_{\mathrm{FLRW}}/\refmodel{a}$ (as defined in \SSS\ref{s-algebra-gev}; dashed line), for the decaying mode (negative square root in Eq.~\protect\ref{e-phiprime-solution-pair}).
{\em Top:} EdS reference model; {\em bottom:} $\Lambda$CDM.
The sudden increase in value for the perturbed {\gevolutionname} simulation (solid curve) at an early time step is the effect of the $\initial{\Phi}$ 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.~\protect\ref{f-flat-FLRW-Phi-test-minus} at \href{\projectzenodofilesbase/gevolution_scale_factors_EdS.dat}{\projectzenodoid/gevolution\_scale\_factors\_EdS.dat}, \href{\projectzenodofilesbase/gevolution_scale_factors_LCDM.dat}{\projectzenodoid/gevolution\_scale\_factors\_LCDM.dat}.
\label{f-flat-FLRW-Phi-test-plus}}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=0.9\columnwidth]{gevn_curvtest_EdS0_001}
\includegraphics[width=0.9\columnwidth]{gevn_curvtest_LCDM0_001}
\caption{Scale factor ratios
$\effmodel{a}_{\initial{\Phi}=0}/\refmodel{a}$ for a standard {\gevolutionname} simulation, as for Fig.~\protect\ref{f-flat-FLRW-Phi-test-minus}, but with a negative initial potential perturbation;
$\effmodel{a}_{\initial{\Phi}=\gevnparamZeromodeamplitudeEdSMinusvalue}/\refmodel{a}$ for a {\gevolutionname} simulation with an initial potential perturbation $\initial{\Phi}=\gevnparamZeromodeamplitudeEdSMinusvalue$ (solid line); and
the expected FLRW decaying mode evolution of $\effmodel{a}_{\mathrm{FLRW}}/\refmodel{a}$ (as defined in \SSS\ref{s-algebra-gev}; dashed line).
{\em Top:} EdS reference model; {\em bottom:} $\Lambda$CDM.
The behaviour mirrors that for the simulation with a positive $\initial{\Phi}$.
Likewise, the initial point is removed from the residual plot for readability.
\label{f-flat-FLRW-Phi-test-minus}}
\end{figure}
\subsection{Software pipeline}
We would like the reader to be able to straightforwardly reproduce \citep{RougierHinsen2017Repdefn} 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 (\citealt{Akhlaghi2020maneage}; see also \citealt{Akhlaghi19IAU335,Infantesainz2020,PeperRouk2021}) which aims to provide a solution to the \enquote*{reproducibility crisis} in science \citep{Peng15,Baker16,Fanelli18}.
Using the free-licensed source package, available as a version of record at {\projectzenodohref} and in \href{\projectgitrepositoryPlainURL}{live} and \href{\projectgitrepositoryarchivedPlainURL}{archived} {\sc 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 {\projectversion} 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 $\Lambda$CDM observational proxy, and a $\Lambda$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.
Thus, the EdS reference model constants are $\currepoch{\refmodel{\Omm}} = \OmegaMEdSvalue$, $\currepoch{\refmodel{\OmLam}} = \OmegaLEdSvalue$, and $\currepoch{\refmodel{H}} = \HubbleEdSvalue$~km/s/Mpc, as explained above (\SSS\ref{s-method}; \citealt{RaczDobos16,RMBO16Hbg1}).
The $\Lambda$CDM reference model constants are $\currepoch{\refmodel{\Omm}} = \OmegaMLCDMvalue$, $\currepoch{\refmodel{\OmLam}} = \OmegaLLCDMvalue$, and $\currepoch{\refmodel{H}} = \HubbleLCDMvalue$~km/s/Mpc.
The aim in both codes is to minimally disrupt the existing code when inserting the spatially uniform perturbation.
\subsubsection{{\inhomogname}} \label{s-method-inhomog}
The {\inhomogname} package \citep{Roukema17silvir,RO19flatness} is a free-licensed C program and library that implements the relativistic Zel'dovich approximation (RZA) for the evolution of the kinematical backreaction $\CQ_\CD$ \citep{BuchRZA1,BuchRZA2}, i.e. the $\CQ_\CD$ Zel'dovich approximation (QZA) \citep{Roukema17silvir}.
The core functions of the library calculate the domain-averaged volume evolution represented by $a_{\CD}(t) := \scalaverageD{a}(t)$, using the RZA estimate for $\CQ_\CD(t)$, given an input of the initial values of the three invariants $\invI, \invII, \invIII$ of the extrinsic curvature in a comoving spatial domain in a flow-orthogonal time foliation.
This is used in the $N$-body simulation context by generating $N$-body simulation initial conditions from a power spectrum using {\mpgraficname} \citep{PrunetPichon08mpgrafic} and using the {\ramsesscalavname} extension of {\sc ramses} \citep{Teyssier02} as a front end.
The front end {\ramsesscalavname} reads the {\mpgraficname} realisation of initial conditions, calls {\dtfename} \citep{SchaapvdWeygaert00,vdWeygaertSchaap07,CW11,Kennel04kdtreeDTFE} to infer the Newtonian velocity gradients for non-overlapping subdomains that cover the full simulation volume, calculates approximations for $\initavinvI, \initavinvII, \initavinvIII$, and calls {\inhomogname} to calculate $a_{\CD}(t)$ for each of the subdomains.
The global mean volume evolution is calculated from the individual $a_{\CD}(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 the match with FLRW solutions.
To see if a global perturbation inserted into {\inhomogname}-{\inhomogversion} yields the expected FLRW expansion rate while minimally intervening in the use of {\inhomogname} via the {\ramsesscalavname} front end, the routine {\tt amr\_dtfe\_interface} is modified to read the value of $\initavinvI$ from a plain-text file, and $\initavinvII$ and $\initavinvIII$ are set to zero.
\begin{table*}
\centering
\caption{Conversion between reference and expected FLRW models. Initial first invariant $\initial{\invI}$ ({\inhomogname}) and initial potential perturbation $\initial{\Phi}$ ({\gevolutionname}); initial reference model epoch $\initial{\refmodel{a}}$; reference model $\refmodel{}$ and expected FLRW model $\effmodel{}$ current-epoch cosmological parameters; expected ($\expected{\effmodel{a}}(\refmodel{a} = 1)$) and numerically measured ($\effmodel{a}(\refmodel{a} = 1)$) final scale factor.
$\currepoch{\refmodel{H}}$ and $\currepoch{\effmodel{H}}$ are in units of km/s/Mpc; all other parameters are dimensionless.
Plain-text version available at \href{\projectzenodofilesbase/flrw_ref_eff_constants.dat}{\projectzenodoid/flrw\_ref\_eff\_constants.dat}
\label{t-ref-expected-measured}}
$\begin{array}{l r r c c c c c c c c c c}
\hline
\multicolumn{13}{c}{\mbox{\inhomogname}} \\
\mbox{model}\rule{0ex}{2.7ex}
& \multicolumn{1}{c}{\initial{\invI}}
& \initial{\refmodel{a}}
& \currepoch{\refmodel{\Omm}}
& \currepoch{\refmodel{\OmLam}}
& \currepoch{\refmodel{\Omk}}
& \currepoch{\refmodel{H}}
& \currepoch{\effmodel{\Omm}}
& \currepoch{\effmodel{\OmLam}}
& \currepoch{\effmodel{\Omk}}
& \currepoch{\effmodel{H}}
& \expected{\effmodel{a}}(\refmodel{a} = 1)
& \effmodel{a}(\refmodel{a} = 1)
\\
\hline
\mbox{EdS}\rule{0ex}{2.7ex}
& \InitPerturbInhomogPlusvalue
& \InhomogEdszZzoAInitRefvalue
& \InhomogEdszZzoOmegaMzRefvalue
& \InhomogEdszZzoOmegaLamzRefvalue
& \InhomogEdszZzoOmegaKzRefvalue
& \InhomogEdszZzoHzRefvalue
& \InhomogEdszZzoOmegaMzEffvalue
& \InhomogEdszZzoOmegaLamzEffvalue
& \InhomogEdszZzoOmegaKzEffvalue
& \InhomogEdszZzoADotZEffvalue
& \InhomogExpectedAratioEdSzdzzovalue & \InhomogPerturbAratioEdSzdzzovalue
\\
\mbox{EdS}
& \InitPerturbInhomogMinusvalue
& \InhomogEdsmzZzoAInitRefvalue
& \InhomogEdsmzZzoOmegaMzRefvalue
& \InhomogEdsmzZzoOmegaLamzRefvalue
& \InhomogEdsmzZzoOmegaKzRefvalue
& \InhomogEdsmzZzoHzRefvalue
& \InhomogEdsmzZzoOmegaMzEffvalue
& \InhomogEdsmzZzoOmegaLamzEffvalue
& \InhomogEdsmzZzoOmegaKzEffvalue
& \InhomogEdsmzZzoADotZEffvalue
& \InhomogExpectedAratioEdSmzdzzovalue & \InhomogPerturbAratioEdSmzdzzovalue
\\
\mbox{$\Lambda$CDM}
& \InitPerturbInhomogPlusvalue
& \InhomogLcdmzZzoAInitRefvalue
& \InhomogLcdmzZzoOmegaMzRefvalue
& \InhomogLcdmzZzoOmegaLamzRefvalue
& \InhomogLcdmzZzoOmegaKzRefvalue
& \InhomogLcdmzZzoHzRefvalue
& \InhomogLcdmzZzoOmegaMzEffvalue
& \InhomogLcdmzZzoOmegaLamzEffvalue
& \InhomogLcdmzZzoOmegaKzEffvalue
& \InhomogLcdmzZzoADotZEffvalue
& \InhomogExpectedAratioLCDMzdzzovalue & \InhomogPerturbAratioLCDMzdzzovalue
\\
\mbox{$\Lambda$CDM}
& \InitPerturbInhomogMinusvalue
& \InhomogLcdmmzZzoAInitRefvalue
& \InhomogLcdmmzZzoOmegaMzRefvalue
& \InhomogLcdmmzZzoOmegaLamzRefvalue
& \InhomogLcdmmzZzoOmegaKzRefvalue
& \InhomogLcdmmzZzoHzRefvalue
& \InhomogLcdmmzZzoOmegaMzEffvalue
& \InhomogLcdmmzZzoOmegaLamzEffvalue
& \InhomogLcdmmzZzoOmegaKzEffvalue
& \InhomogLcdmmzZzoADotZEffvalue
& \InhomogExpectedAratioLCDMmzdzzovalue & \InhomogPerturbAratioLCDMmzdzzovalue
\\
\hline
\multicolumn{12}{c}{\mbox{\gevolutionname}\rule{0ex}{2.7ex}} \\
\mbox{model}\rule{0ex}{2.7ex}
& \multicolumn{1}{c}{\initial{\Phi}}
& \initial{\refmodel{a}}
& \currepoch{\refmodel{\Omm}}
& \currepoch{\refmodel{\OmLam}}
& \currepoch{\refmodel{\Omk}}
& \currepoch{\refmodel{H}}
& \currepoch{\effmodel{\Omm}}
& \currepoch{\effmodel{\OmLam}}
& \currepoch{\effmodel{\Omk}}
& \currepoch{\effmodel{H}}
& \expected{\effmodel{a}}(\refmodel{a} = 1)
& \effmodel{a}(\refmodel{a} = 1)
\\
\hline
\mbox{EdS}\rule{0ex}{2.7ex}
& \GevolutionEdszZzoPhiInitvalue
& \GevolutionEdszZzoAInitRefvalue
& \GevolutionEdszZzoOmegaMzRefvalue
& \GevolutionEdszZzoOmegaLamzRefvalue
& \GevolutionEdszZzoOmegaKzRefvalue
& \GevolutionEdszZzoHzRefvalue
& \GevolutionEdszZzoOmegaMzEffvalue
& \GevolutionEdszZzoOmegaLamzEffvalue
& \GevolutionEdszZzoOmegaKzEffvalue
& \GevolutionEdszZzoADotZEffvalue
& \GevolutionExpectedAratioFinalEdSzdzzovalue & \GevolutionPhiAratioFinalEdSzdzzovalue
\\
\mbox{EdS}
& \GevolutionEdsmzZzoPhiInitvalue
& \GevolutionEdsmzZzoAInitRefvalue
& \GevolutionEdsmzZzoOmegaMzRefvalue
& \GevolutionEdsmzZzoOmegaLamzRefvalue
& \GevolutionEdsmzZzoOmegaKzRefvalue
& \GevolutionEdsmzZzoHzRefvalue
& \GevolutionEdsmzZzoOmegaMzEffvalue
& \GevolutionEdsmzZzoOmegaLamzEffvalue
& \GevolutionEdsmzZzoOmegaKzEffvalue
& \GevolutionEdsmzZzoADotZEffvalue
& \GevolutionExpectedAratioFinalEdSmzdzzovalue & \GevolutionPhiAratioFinalEdSmzdzzovalue
\\
\mbox{$\Lambda$CDM}
& \GevolutionLcdmzZzoPhiInitvalue
& \GevolutionLcdmzZzoAInitRefvalue
& \GevolutionLcdmzZzoOmegaMzRefvalue
& \GevolutionLcdmzZzoOmegaLamzRefvalue
& \GevolutionLcdmzZzoOmegaKzRefvalue
& \GevolutionLcdmzZzoHzRefvalue
& \GevolutionLcdmzZzoOmegaMzEffvalue
& \GevolutionLcdmzZzoOmegaLamzEffvalue
& \GevolutionLcdmzZzoOmegaKzEffvalue
& \GevolutionLcdmzZzoADotZEffvalue
& \GevolutionExpectedAratioFinalLCDMzdzzovalue & \GevolutionPhiAratioFinalLCDMzdzzovalue
\\
\mbox{$\Lambda$CDM}
& \GevolutionLcdmmzZzoPhiInitvalue
& \GevolutionLcdmmzZzoAInitRefvalue
& \GevolutionLcdmmzZzoOmegaMzRefvalue
& \GevolutionLcdmmzZzoOmegaLamzRefvalue
& \GevolutionLcdmmzZzoOmegaKzRefvalue
& \GevolutionLcdmmzZzoHzRefvalue
& \GevolutionLcdmmzZzoOmegaMzEffvalue
& \GevolutionLcdmmzZzoOmegaLamzEffvalue
& \GevolutionLcdmmzZzoOmegaKzEffvalue
& \GevolutionLcdmmzZzoADotZEffvalue
& \GevolutionExpectedAratioFinalLCDMmzdzzovalue & \GevolutionPhiAratioFinalLCDMmzdzzovalue
\\
\hline
\end{array}$
\end{table*}
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 {\cosmdistname}-{\cosmdistversion} where appropriate, independently from the {\inhomogname} code and its {\ramsesscalavname} front end.
Equations~\eqref{e-Omegami-eff-inhomog}, \eqref{e-OmegaLami-eff}, and \eqref{e-Omegaki-eff} relate the reference and emergent matter density parameters, $\refmodel{\Omm}$ and $\effmodel{\Omm}$, at the early epoch when the intervention occurs.
The equations in \SSSS\ref{s-algebra-FLRW-rewritten}, \ref{s-algebra-ref-to-eff} and \ref{s-time-matching} then show how the emergent FLRW model constants $\effmodel{\{ \currepoch{\Omm}, \currepoch{\OmLam}, \currepoch{\Omk}, \currepoch{H}\}}$ are calculated from those of the reference model, $\refmodel{\{ \currepoch{\Omm}, \currepoch{\OmLam}, \currepoch{\Omk}, \currepoch{H}\}}$.
\subsubsection{{\gevolutionname}} \label{s-method-gev}
The {\gevolutionname} software is a free-licensed C++ cosmological simulation code developed by \citet{AdamekDDK15,AdamekDDK16code} 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 \citep[e.g.][]{Hassani19,Reverberi19}.
The input parameter file for the simulation includes physical parameters such as $\refmodel{\Omm}$ and $\refmodel{H_0}$, 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.~\eqref{e-line-element}).
The code is highly modular.
The numerical core of {\gevolutionname} is {\LATfieldtwo}, 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 {\gevolutionname} is to determine whether the code is currently ({\gevolutionname}-{\gevolutionversion}) 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 {\tt generateIC\_basic()} in the file {\tt ic\_basic.hpp}, adding a step that allows the insertion of a uniform perturbation in the scalar gravitational potential $\Phi$.
The optional insertion of the perturbation and the perturbation's amplitude are controlled in the initial settings file.
The {\gevolutionname} 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\footnote{\projectgevpatchfileSWHhref}.
When we refer to {\gevolutionname}-{\gevolutionversion}, strictly speaking we mean {\gevolutionname}-{\gevolutionversion} 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 $\initial{\Phi}$ and $\initial{\dotPhi}$.
The reasoning is presented above in \SSS\ref{s-algebra-ref-to-eff-old}.
The conversion of initial scale factors to obtain $\initial{\effmodel{a}}$ is given in Eq.~\eqref{e-ai-eff-gev}.
The effective matter density $\initial{\effmodel{\Omm}}$ at the epoch of intervention is given in Eq.~\eqref{e-Omegami-eff-gev}.
The latter requires $\initial{\effmodel{\dot{a}}}$, which is provided by Eq.~\eqref{e-adoti-eff-gev}, provided that we know the value of $\dotPhi$ at this epoch.
The choice between growing and decaying modes for $\dotPhi$, corresponding to the positive and negative square roots in Eq.~\protect\ref{e-phiprime-solution-pair}, respectively (converted with Eq.~\ref{e-phi-t-tau-relation}), is set to the decaying mode, since {\gevolutionname} uses the linearised solution (see \SSS\ref{s-algebra-ref-to-eff-old}), which corresponds to the decaying perturbation.
As explained above (\SSS\ref{s-intro}, \SSS\ref{s-algebra-gev}), {\gevolutionname} and its backend {\LATfieldtwo} 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 $\Phi$ 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~\ref{app-phiprime} and mentioned above, this allows two possible flat perturbation modes satisfying the Hamiltonian constraint (choice of $\Phi'$, or equivalently, $\dotPhi$).
We compare the numerical evolution of {\gevolutionname} with the decaying solution, which is clearly the one chosen by the code.
The relations in \SSSS\ref{s-algebra-FLRW-rewritten}, \ref{s-algebra-ref-to-eff} and \ref{s-time-matching} are then again used to obtain the emergent FLRW model constants, and {\cosmdistname} is used to generate the expected emergent scale factor evolution.
\section{Results} \label{s-results}
\subsection{{\inhomogname}} \label{s-results-inhomog}
We ran the relativistic Zel'dovich approximation evolution package {\inhomogname}-{\inhomogversion} from $N$-body initial conditions, with a resolution of $N=\Ncrootinhomogvalue^3$ particles, fundamental domain comoving (reference model) size {\Lboxvalue}~Mpc/$h$, and the {\sc ramses} resolution parameter {\ramsesparamlevelmaxEdSname} $= \ramsesparamlevelmaxEdSvalue$, with and without a perturbation as described in \SSS\ref{s-method-inhomog}.
We compared the numerically measured scale factor evolution with the expected FLRW model by defining a residual $\delta$ and a fractional deviation $\epsilon$, defined
\begin{align}
\delta &:= \effmodel{a}/\refmodel{a} - \expected{\effmodel{a}}/\refmodel{a} \nonumber \\
\epsilon &:= \frac{\effmodel{a}/\refmodel{a} - \expected{\effmodel{a}}/\refmodel{a}}{
\expected{\effmodel{a}}/\refmodel{a}} \,.
\label{e-resid-inhomog}
\end{align}
The fractional deviation $\epsilon$ describes the numerical deviation as a fraction of the excess expansion that should be generated by the perturbation.
Figure~\ref{f-flat-FLRW-Iinv-test-plus} (for an initial density perturbation $\initial{\invI} = \InitPerturbInhomogPlusvalue$) and Fig.~\ref{f-flat-FLRW-Iinv-test-minus} (for $\initial{\invI} = \InitPerturbInhomogMinusvalue$) 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 $\epsilon = \InhomogAccuracyPercentEdSzdzzovalue$\% for $\initial{\invI} = \InitPerturbInhomogPlusvalue$, and $\epsilon = \InhomogAccuracyPercentEdSmzdzzovalue$\% for $\initial{\invI} = \InitPerturbInhomogMinusvalue$.
The corresponding $\Lambda$CDM reference model deviations are $\epsilon = \InhomogAccuracyPercentLCDMzdzzovalue$\% and $\epsilon = \InhomogAccuracyPercentLCDMmzdzzovalue$\%, respectively.
It is clear that {\sc 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~\ref{t-ref-expected-measured} 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 $\currepoch{\effmodel{\Omk}}$ appears in the expected FLRW model.
These are correctly handled by {\inhomogname}, while {\gevolutionname} instead returns to the reference model.
\subsection{{\gevolutionname}} \label{s-results-gevolution}
As in the {\inhomogname} case, we ran independent simulations for the EdS and $\Lambda$CDM models for {\gevolutionname}-{\gevolutionversion}, in each case with and without two cases of an initial perturbation: $\initial{\Phi}=\gevnparamZeromodeamplitudeEdSPlusvalue$ and $\initial{\Phi}=\gevnparamZeromodeamplitudeEdSMinusvalue$.
The simulations' particle resolution is $N={\Ncrootgevolutionvalue}^3$; the $\Lambda$CDM simulation was run with $\OmLamzeroref={\OmegaLLCDMvalue}$; and the fundamental domain comoving (reference model) size was {\Lboxvalue}~Mpc/$h$.
In the cases with $\initial{\Phi}=\gevnparamZeromodeamplitudeEdSPlusvalue$, 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.~\eqref{e-Omegami-eff-gev}.
Moreover, $\initial{\effmodel{\OmLam}}$ is also modified by the perturbation (Eq.~\eqref{e-OmegaLami-eff}), leaving $\initial{\effmodel{\Omk}}$ as the complement (Eq.~\eqref{e-Omegaki-eff}).
The evolution of the potential $\dotPhi$ enters into these relations via Eq.~\eqref{e-adoti-eff-gev}.
The linearisation in {\gevolutionname} implies the decaying mode of $\Phi$, which for $\Phi > 0$ yields $\dotPhi < 0$.
The formulae in \SSS\ref{s-algebra-generic} calculated using {\cosmdistname}-{\cosmdistversion} give a flat effective model ($\currepoch{\effmodel{\Omk}} = 0$) in these two cases, as shown in Table~\ref{t-ref-expected-measured}.
These formulae also give the expected behaviour of the $\effmodel{a}/\refmodel{a}$ ratio, which was calculated independently of {\gevolutionname}, using {\cosmdistname}.
%results
Figures~\ref{f-flat-FLRW-Phi-test-plus} and \ref{f-flat-FLRW-Phi-test-minus} show the behaviour of {\gevolutionname} with a perturbation in comparison to the expected FLRW behaviour equivalent to that of the decaying mode.
We find that {\gevolutionname} 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 {\gevolutionname}'s deviations $|\epsilon|$ from the expected behaviour for the EdS model are
$\epsilon = \GevolutionAccuracyPercentMaxEdSzdzzovalue$\% for $\initial{\Phi} = \gevnparamZeromodeamplitudeEdSPlusvalue$ and
$\epsilon = \GevolutionAccuracyPercentMaxEdSmzdzzovalue$\% for $\initial{\Phi} = \gevnparamZeromodeamplitudeEdSMinusvalue$; and for the $\Lambda$CDM model the deviations are
$\epsilon = \GevolutionAccuracyPercentMaxLCDMzdzzovalue$\% for $\initial{\Phi} = \gevnparamZeromodeamplitudeLCDMPlusvalue$ and
$\epsilon = \GevolutionAccuracyPercentMaxLCDMmzdzzovalue$\% for $\initial{\Phi} = \gevnparamZeromodeamplitudeLCDMMinusvalue$.
Those quantities exclude the initial spike seen in residual plots as an artefact of perturbation introduction.
In summary, we find that {\gevolutionname}-{\gevolutionversion} models the decaying mode well, but does not allow the growing mode of the quadratic equation in $\Phi'$ (we discuss this in \SSS\ref{s-discuss-gevolution}).
\section{Discussion} \label{s-discuss}
\subsection{{\inhomogname}} \label{s-discuss-inhomog}
As shown in Figs~\ref{f-flat-FLRW-Iinv-test-plus} and \ref{f-flat-FLRW-Iinv-test-minus}, {\inhomogname}-{\inhomogversion} is reasonably accurate in modelling the volume evolution implied by a density perturbation.
Nevertheless, the residual errors $\delta$ (Eq.~\eqref{e-resid-inhomog}) 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. $\delta \sim \left[\effmodel{a}(\refmodel{a}=1)\right]^2$.
This would appear to suggest a second-order 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 \citep[][\SSS{}V.A, V.B]{BuchRZA2}, 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\% \citep{RMBO16Hbg1} rather than 1\%.
\subsection{{\gevolutionname}} \label{s-discuss-gevolution}
In \citet[][\SSS5.3]{AdamekDDK16code}, 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.~\protect\eqref{e-phiprime-solution-pair}, i.e. the solution to the quadratic in $\Phi'$, the linearisation given in eqs~(2.9) of \citet{AdamekDDK16code} and eq.~(9) of \citet{Adamek2017radiation} restricts the evolution of a homogeneous mode to the decaying branch.
Equation~\eqref{e-Hamiltonian-constraint-quadratic} shows that this linearisation should be a good approximation when $\vert \Phi'\vert \ll \vert \refmodel{\dot{a}} \vert = \refmodel{\dot{a}}$, where the equality follows from only considering an expanding universe (not a contracting or stationary one).
For a decaying perturbation, if $\vert \Phi'\vert$ 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
\begin{align}
\vert \Phi' \vert > \refmodel{\dot{a}} &\Leftrightarrow \Phi' > \refmodel{\dot{a}}
\Leftrightarrow \dotPhi > \refmodel{\dot{a}}/\refmodel{a} \,.
\end{align}
Equation~\eqref{e-adoti-eff-gev}, with the subscript $\initial{{}}$ reinserted, shows that for small $\vert\Phi\vert$, this implies that $\initial{\effmodel{\dot{a}}} < 0 $.
As a global effect, this would represent a universe starting to collapse, which cannot correspond to any EdS solution nor to the $\Lambda$CDM solution.
Alternatively, as an effective model for a locally collapsing spatial region, average flatness in the scalar averaging sense does not strictly forbid gravitational turnaround from occurring in this way.
However, in practice, it is rare for a region that is on average flat, in an EdS or $\Lambda$CDM reference model evolving from a cosmologically typical initial power spectrum of perturbations, to slow in its expansion and collapse past its turnaround epoch.
This was shown using the relativistic Zel'dovich approximation \citep{RO19flatness}.
Thus, the relevance of the omitted growing mode is modest, since the mode is likely to be rare in practice.
Nevertheless, if the {\gevolutionname} convention of setting $\Phi = \Psi$ for the homogeneous mode were dropped, then more possibilites would be allowed.
An algebraic demonstration showing how both decaying and growing solutions are possible for uniform perturbations is given in App.~\ref{app-growing-decaying-modes} for the EdS case, for illustration.
Both the decaying and growing solutions are EdS solutions, i.e. they expand eternally.
Thus, allowing $\Phi \ne \Psi$ for the homogeneous mode could be expected to yield an extension of {\gevolutionname} that would enable the flat growing mode to be studied.
Given that {\gevolutionname} 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 $\Phi \ne \Psi$ 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 non-linearities 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 \SSS\ref{s-algebra-generic} 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 {\gevolutionname} 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.
\subsection{Numerical representation of geometrical information}
As explained above (\SSS\ref{s-method-projection}) in relation to {\gevolutionname}, 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.
\subsection{Free-licensed cosmology software ecosystem}
The results of this paper illustrate how the use of free-licensed, modular codes written in C ({\cosmdistname} and {\inhomogname}) or C++ ({\gevolutionname}) can be easily and efficiently checked against each other.
These should be easy to use in checks or extensions of the Einstein Toolkit \citep{BentivegnaBruni15,Macpherson17}, part of a much bigger free-licensed relativistic code project, or of {\sc cosmograph} \citep{GiblinMS17}.
\section{Conclusions} \label{s-conclu}
%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 {\inhomogname} 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 {\gevolutionname}, we found that out of the two solutions of a quadratic solution, the choice of solution is hardwired.
The choice made by {\gevolutionname} is the decaying mode that is implied by linearisation in the first conformal-time derivative of the potential, $\Phi_{,\tau}$.
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 {\gevolutionname}.
Work extending {\gevolutionname} 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 \projectversion.
\section*{Data Availability Statement} \label{s-data-availability}
The full reproducibility package for this paper, including the choice of simulation parameters (\enquote*{input data}), is available at {\projectzenodoid}\footnote{{\projectzenodohrefShowURL}} and in live\footnote{\projectgitrepository} and archived\footnote{\projectgitrepositoryarchived} {\sc git} repositories.
%The numerical output data for our main results is available as a DOI-identified record at \mbox{\href{\projectzenodofilesbase/voidgals_infall.dat}{\projectzenodoid/voidgals\_infall.dat}}.
The specific version of the source package used to produce this paper can be identified by its {\sc git} commit hash \projectversion.
%% Acknowledgements
\section*{Acknowledgments}
The authors wish to thank Julian Adamek for very detailed and rapid feedback, and Marius Peper, Krzysztof Bolejko, Jan Ostrowski, Pierre Mourier, Asta Heinesen and Quentin Vigneron
%...
for many useful suggestions.
Work on this paper has been supported by the ``A next-generation worldwide quantum sensor network with optical atomic clocks'' project of the TEAM IV programme of the Foundation for Polish Science co-financed by the European Union under the European Regional Development Fund.
Part of this work has been supported by the Polish MNiSW grant DIR/WK/2018/12.
Part of this work has been supported by the Pozna\'n Supercomputing and Networking Center (PSNC) computational grant 314.
Part of this work was supported by Universitas Copernicana Thoruniensis in Futuro under NCBR grant POWR.03.05.00-00-Z302/17.
\input{tex/build/macros/dependencies.tex}
%% Tell BibLaTeX to put the bibliography list here.
\bibliographystyle{\bibtexbstfile}
\bibliography{\bibtexbibfiletmp}
\begin{appendix}
\section{Evolution of the perturbation: Hamiltonian constraint} \label{app-phiprime}
In order to properly calculate the scale factor change in Eq.~\eqref{e-adoti-eff-gev}, a correct value of $\dot\Phi$ is needed -- whether as numerical data output by the code or calculated analytically.
As kindly pointed out by J.\/ Adamek, an analytical equation for the evolution of the homogeneous perturbation can be straightforwardly derived.
We start from the pointwise (rather than volume-averaged) Hamiltonian constraint in the $3+1$ decomposition of the Einstein equation \citep{Gourg07lecture}
%The Hamiltonian constraint in Poisson gauge for a homogeneous perturbation can be written as:
%
\begin{align}
{}^3{\CR} - K_{ij}K^{ij} + K^2 &= 16\pi G \mathbf{T}(\vec{n}, \vec{n}) = 16 \pi G T^{00}= 16 \pi G \rho\,,
\label{e-Hamiltonian-constraint}
\end{align}
%
where the (pointwise) quantities are the spatial curvature ${}^3{\CR}$, the extrinsic curvature at the same spacetime location $K_{ij}$, its trace $K$, the Newtonian gravitational constant $G$, the momentum--energy tensor $\mathbf{T}$ and its time--time component, which in this case is the matter density, $T^{00} = \rho$, and $\vec{n}$ is the normal to the spatial hypersurface at that point.
This equation needs to be expressed using the Poisson gauge formalism used in {\gevolutionname}.
The spatial part of the Ricci curvature can be calculated from the metric expressed as a line element, which in Poisson gauge (Eq.~\eqref{e-line-element}) and the case $B_{i} =0 = h_{ij}$ reduces to
%
\begin{align}
{}^3{\CR} = 2 \left(\refmodel{a}\right)^{-2} e^{2\Phi} \delta^{ij} (2 \Phi_{,ij} - \Phi_{,i}\Phi_{,j}) \,,
\label{e-3Ricci-Poisson}
\end{align}
%
where $\Phi_{,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
%
\begin{align}
N = \refmodel{a}{\mathrm e}^{\Psi}\,.
\label{e-lapse-function-Poisson}
\end{align}
Using $K_{ij} = -1/(2N)\,\partial/\partial \tau \ \gamma_{ij}$ [e.g., \citet[][(4.63); $\gamma_{ij}$ is the spatial part of the line element]{Gourg07lecture}], and Eq.~\eqref{e-phi-t-tau-relation}, gives
%
\begin{align}
K_{ij} = \refmodel{a} \left(\Phi' - \refmodel{\dot{a}} \right) \, \mathrm{e}^{-2\Phi-\Psi} \delta_{ij} \,,
\label{e-extrinsic-curvature}
\end{align}
%
making the switch between proper to coordinate time derivatives for convenience ($\refmodel{a}' = \refmodel{a}\refmodel{\dot{a}}$). The mixed, contravariant, and scalar extrinsic curvature are
%
\newcommand\detailedExCurvOne{%
\textbf{\em (A few intermediate details are given in red for extra checking only; these will be toggled off in the submitted version.)}
\begin{align*}%
\gamma^{ij} &= a^{-2} e^{2\Phi} \delta^{ij}
\end{align*}%
}%
\iftoggle{showDetailedDerivations}{{\color{red}\detailedExCurvOne}}{}%
\begin{align}
K^i_{\ j} &= \gamma^{ik} K_{kj} = \left(\refmodel{a}\right)^{-1} (\Phi' - \refmodel{\dot{a}}) \mathrm{e}^{-\Psi} \delta^i_{\ j} \,,\nonumber\\
K^{ij} &= \gamma^{ik} K^k_{\ j} = \left(\refmodel{a}\right)^{-3} (\Phi' - \refmodel{\dot{a}}) e^{2\Phi - \Psi} \delta^{ij} \,,\quad\mathrm{and}
\nonumber\\ %\label{e-extrinsic-curvature-inverse} \\
K &= K^i_{\ i} = 3 \left(\refmodel{a}\right)^{-1} (\Phi' - \refmodel{\dot{a}}) e^{-\Psi} \,,
\label{e-extrinsic-curvature-scalar}
\end{align}
respectively. %
\newcommand\detailedExCurvTwo{
\begin{align*}
K^{ij} K_{ij} = 3 a^{-2} (\Phi' - \refmodel{\dot{a}})^2 e^{-2\Psi}
\end{align*}
}%
\iftoggle{showDetailedDerivations}{{\color{red}\detailedExCurvTwo}}{}%
%\newcommand\detailedHamiltonianConstraint{
Thus, using Eqs~\eqref{e-extrinsic-curvature} and \eqref{e-extrinsic-curvature-scalar} in Eq.~\eqref{e-Hamiltonian-constraint}, we obtain
\begin{align}
&2 \left(\refmodel{a}\right)^{-2} e^{2\Phi} \delta^{ij} (2 \Phi_{,ij} - \Phi_{,i}\Phi_{,j}) \nonumber\\
&\phantom{2\ } - 3 a^{-2} (\Phi' - \refmodel{\dot{a}})^2 e^{-2\Psi} \nonumber\\
&\phantom{2\ }\, + 9 a^{-2} (\Phi' - \refmodel{\dot{a}})^2 e^{-2\Psi} %\nonumber\\
%&\phantom{2\ =}
= 16\pi G \rho \,.
\label{e-Hamiltonian-constraint-expanded}
\end{align}
%}
%\iftoggle{showDetailedDerivations}{{\color{red}\detailedHamiltonianConstraint}}{}
On the right-hand side of Eq.~\eqref{e-Hamiltonian-constraint-expanded} we have the density:
\begin{align}
G\,T^{00} = G\,\rho &= \frac{3G \currepoch{\refmodel{H}}^2 }{8\pi G\left(\refmodel{a}\right)^2}
\left( \currepoch{\refmodel{\Omm}}\, \left({\refmodel{a}}\right)^{-1} \mathrm{e}^{3\Phi} +
\currepoch{\refmodel{\OmLam}}\, \left({\refmodel{a}}\right)^2 \right)
\,,
\label{e-rho-in-Omegas}
\end{align}
from the reference model (FLRW) Hamiltonian constraint, which is $\left(\refmodel{\dot{a}}\right)^2 = \currepoch{\refmodel{H}}^2 \left( \currepoch{\refmodel{\Omm}} \left({\refmodel{a}}\right)^{-1} + \currepoch{\refmodel{\OmLam}} \left({\refmodel{a}}\right)^2 \right)$, together with a correction by a factor of $\mathrm{e}^{3\Phi}$ from the perturbed line element (Eq.~\eqref{e-line-element}).
Since we have a flat perturbation field with $\Phi_{,i} = 0 = \Phi_{,ij} \,\forall\, i,j$, Eq.~\eqref{e-Hamiltonian-constraint-expanded} simplifies to
\begin{align}
\mathrm{e}^{-2\Phi} \left(\refmodel{\dot{a}} - \Phi'\right)^2
&= \frac{8\pi G}{3} \left(\refmodel{a}\right)^2 \rho
\nonumber \\
&= \currepoch{\refmodel{H}}^2 \left( \currepoch{\refmodel{\Omm}}\, \left({\refmodel{a}}\right)^{-1} \mathrm{e}^{3\Phi} +
\currepoch{\refmodel{\OmLam}}\, \left({\refmodel{a}}\right)^2 \right)
\,,
\label{e-Hamiltonian-constraint-quadratic}
\end{align}
Subtracting the reference model Hamiltonian constraint gives
%
\begin{align}
\Phi' &= \refmodel{\dot{a}} \pm \mathrm{e}^{\Psi} \,
\sqrt{\left({\refmodel{\dot{a}}}\right)^2 +
\frac{\currepoch{\refmodel{H}}^2\, \currepoch{\refmodel{\Omm}}}
{\refmodel{a}}
\left(\mathrm{e}^{3\Phi} - 1\right)}
\,.
\label{e-phiprime-solution-pair}
\end{align}
The choice of the sign before the square root and the setting of $\Psi = \Phi$ may constrain the behaviour of the solution.
The negative square root together with $\Psi=\Phi$ give a decaying mode, in which over- and under-densities weaken with time.
The positive square root, possibly together with $\Psi \ne \Phi$, may give a growing mode; this is discussed in \SSS\ref{s-discuss-gevolution}.
A solution that linearises Eq.~\eqref{e-Hamiltonian-constraint-quadratic} in terms of $\Phi'$ (ignores the $\Phi'{}^2$ term) has only the decaying solution.
Examples of exact decaying and growing modes in the EdS case are given in \SSS\ref{app-growing-decaying-modes}.
\section{Evolution of the perturbation: Raychaudhuri equation} \label{app-Raychaudhuri}
As stated above, the shift function in our case is zero and the lapse has no spatial dependence, so
\begin{align}
S &= 0 \nonumber \\
S_{ij} &= 0 \nonumber \\
\mathrm{D}_i \mathrm{D}_j N &= 0 \,.
\end{align}
Thus, the Raychaudhuri equation \citep[][(4.64), cf (4.79 for the $N=1$ case)]{Gourg07lecture} becomes
\begin{align}
\frac{\partial}{\partial\tau} K_{ij}
&=
N \left( 0+ K K_{ij} -2 K_{ik} K^k_j -4\pi GT^{00} \gamma_{ij} \right)\,,
\end{align}
i.e., using Eqs~\eqref{e-lapse-function-Poisson}, \eqref{e-extrinsic-curvature}, \eqref{e-extrinsic-curvature-scalar}, and \eqref{e-rho-in-Omegas},
\newcommand\detailedRaychaudhuri{
\begin{align*}
&\left[ \refmodel{a}\refmodel{\dot{a}} \mathrm{e}^{-\Psi-2\Phi}
\left(\Phi' - \refmodel{\dot{a}} \right)
\right.
\nonumber\\
&\phantom{[} + \refmodel{a} \left(-\Psi'-2\Phi'\right) \mathrm{e}^{-\Psi-2\Phi}
\left(\Phi' - \refmodel{\dot{a}} \right)
\nonumber\\
&\phantom{[} + \left.
\refmodel{a} \left(-\Psi'-2\Phi'\right)
\left(\Phi'' - \refmodel{a}\refmodel{\dot{a}} \right)
\vphantom{e^{-\Psi}} \right] \delta_{ij}
\nonumber\\
&\phantom{[+} = \refmodel{a} \mathrm{e}^{\Psi}
\left[
\vphantom{\frac{3\currepoch{\refmodel{H}}^2}{8\pi \left(\refmodel{a}\right)^2}} % vphantom
3 \left(\refmodel{a}\right)^{-1} \mathrm{e}^{-\Psi}
\left(\Phi' - \refmodel{\dot{a}} \right) \,
\refmodel{a} \mathrm{e}^{-\Psi-2\Phi}
\left(\Phi' - \refmodel{\dot{a}} \right)
\right.
\nonumber\\
&\phantom{[+=} -2 \refmodel{a} \mathrm{e}^{-\Psi-2\Phi}
\left(\Phi' - \refmodel{\dot{a}} \right) \,
\left(\refmodel{a}\right)^{-1} \mathrm{e}^{-\Psi}
\left(\Phi' - \refmodel{\dot{a}} \right)
\nonumber\\
&\phantom{[+=}\left.
-4\pi \left(\refmodel{a}\right)^2 \mathrm{e}^{-2\Phi}
\frac{3\currepoch{\refmodel{H}}^2}{8\pi \left(\refmodel{a}\right)^2}
\left( \currepoch{\refmodel{\Omm}}\, \left({\refmodel{a}}\right)^{-1} \mathrm{e}^{3\Phi} +
\currepoch{\refmodel{\OmLam}}\, \left({\refmodel{a}}\right)^2 \right)
\right] \delta_{ij}
\nonumber\\
\\
& \mathrm{e}^{-\Psi-2\Phi}
\left[
\left(\Phi' - \refmodel{\dot{a}} \right)
\left(\refmodel{\dot{a}} -2\Psi' - \Phi' \right)
+ \Phi'' - \refmodel{a} \refmodel{\ddot{a}}\right]
\nonumber\\
&\phantom{e} =
\mathrm{e}^\Psi \left[ \vphantom{\frac{3}{2}}
\mathrm{e}^{-2\Psi-2\Phi}
\left(\Phi' - \refmodel{\dot{a}} \right)^2
\right.
\nonumber\\
&\phantom{e= \mathrm{e}^\Psi [}
\left.
- \frac{3}{2}\, \mathrm{e}^{-2\Phi}
\currepoch{\refmodel{H}}^2
\left( \currepoch{\refmodel{\Omm}}\, \left({\refmodel{a}}\right)^{-1} \mathrm{e}^{3\Phi} +
\currepoch{\refmodel{\OmLam}}\, \left({\refmodel{a}}\right)^2 \right)
\right]
\nonumber\\
\\
& (\Phi'-\refmodel{\dot{a}}) \left(\refmodel{\dot{a}} -\Psi' -2\Phi'\right)
+ \Phi'' - \refmodel{a}\, \refmodel{\ddot{a}}
\nonumber\\
&\phantom{xx} =
(\Phi'-\refmodel{\dot{a}})^2
- \frac{3}{2} \mathrm{e}^{2\Psi} H_0^2
\left(\Ommzero\, \refmodel{a}^{-1} \mathrm{e}^{3\Phi}+ \OmLamzero\, \left(\refmodel{a}\right)^{2} \right) \,.
\end{align*}}
%% The following line shows these intermediate details for internal drafts
%% only; set \toggletrue{showDetailedDerivations} at the top of the file
%% to enable it.
\iftoggle{showDetailedDerivations}{{\color{red}\detailedRaychaudhuri}}{}
\begin{align}
\Phi'' &= \refmodel{a}\, \refmodel{\ddot{a}}
- (\Phi'-\refmodel{\dot{a}}) \left(\refmodel{\dot{a}} -\Psi' -2\Phi'\right)
+ (\Phi'-\refmodel{\dot{a}})^2
\nonumber\\
&\phantom{=\ }- \frac{3}{2} \mathrm{e}^{2\Psi} H_0^2
\left(\Ommzero\, \refmodel{a}^{-1} \mathrm{e}^{3\Phi}+ \OmLamzero\, \left(\refmodel{a}\right)^{2} \right) \,.
\label{e-Raychaudhuri-Phi}
\end{align}
\section{Decaying and growing perturbations in effective solutions} \label{app-growing-decaying-modes}
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\^{\i}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.~\eqref{e-match-t}).
The reference and effective models are
\begin{align}
\refmodel{a} &= \left(\frac{\refmodel{t}}{\currepoch{\refmodel{t}}}\right)^{2/3} \,,\nonumber\\
\effmodel{a} &= \left(\frac{\effmodel{t}}{\currepoch{\effmodel{t}}}\right)^{2/3}
= \left(\frac{\refmodel{t} - \initial{\refmodel{t}} + \initial{\effmodel{t}}}{\currepoch{\effmodel{t}}}\right)^{2/3}
\,.
\label{e-EdS-illustration-2-free-params}
\end{align}
After setting $\currepoch{\refmodel{t}}$ for the reference model and $\initial{\refmodel{t}}$ for the insertion of the perturbation, the two parameters $\currepoch{\effmodel{t}}$ and $\initial{\effmodel{t}}$ both need to be chosen to provide the effective model.
A decaying perturbation occurs when $\currepoch{\refmodel{t}} = \currepoch{\effmodel{t}}$, since, as $t \rightarrow \infty$,
\begin{align}
\frac{\vert \effmodel{a} - \refmodel{a} \vert}{\refmodel{a}}
&\rightarrow
\left(\frac{\initial{\effmodel{t}} - \initial{\refmodel{t}}}{\refmodel{t}} \right)^{2/3}
\rightarrow 0 \,,
\end{align}
which will only be valid if $\initial{\effmodel{t}} - \initial{\refmodel{t}} > 0$.
When $\currepoch{\refmodel{t}} \neq \currepoch{\effmodel{t}}$, as $t \rightarrow \infty$,
\begin{align}
\frac{\vert \effmodel{a} - \refmodel{a} \vert}{\refmodel{a}}
&\rightarrow
\left| \left(\frac{\currepoch{\refmodel{t}}}{\currepoch{\effmodel{t}}}\right)^{2/3} -1 \right|
\rightarrow \infty \ \mathrm{or}\ 1 \,,
\end{align}
giving a perturbation that grows in amplitude.
%The situation is similar for the $\Lambda$CDM case.
\end{appendix}
%% Finish LaTeX
\end{document}
%% This file is part of a paper developed within the Maneage reproducible workflow
%% template:
%%
%% https://maneage.org
%
%% This file is free software: you can redistribute it and/or modify it
%% under the terms of the GNU General Public License as published by the
%% Free Software Foundation, either version 3 of the License, or (at your
%% option) any later version.
%
%% This file is distributed in the hope that it will be useful, but WITHOUT
%% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
%% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
%% for more details.
%
%% You should have received a copy of the GNU General Public License along
%% with this file. If not, see <https://www.gnu.org/licenses/>.
