A numerical relativity scheme for cosmological simulations

Let us consider the metric in ADM form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\al}{\alpha} \newcommand{\ed}{{\rm d}} \newcommand{\bo}{\square} \displaystyle \ed s^2 = - \al^2 \ed t^2 + \ga_{ij} \left(\ed x^i + \be^i \ed t \right) \left(\ed x^{\,j} + \be^{\,j} \ed t \right) \,, \nonumber \end{align} \tag{ 2.1 }$$

where $\newcommand{\al}{\alpha} \al$ is the lapse function, $\newcommand{\bo}{\square} \newcommand{\et}{\eta} \newcommand{\e}{{\rm e}} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\be}{\beta} \be^i$ is the shift vector and $\newcommand{\ga}{\gamma} \ga_{ij}$ is the induced 3-metric, so that $\newcommand{\pa}{\partial} \newcommand{\bo}{\square} \newcommand{\al}{\alpha} \newcommand{\et}{\eta} \newcommand{\e}{{\rm e}} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\be}{\beta} n := \al^{-1} \left(\pa_t - \be^i \pa_i \right)$ is the unit vector normal to the spatial hypersurfaces. The gravitational equations of motion of GR in standard ADM form then read ($8\pi G = c = 1$)

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Lie}{\mathcal{L}} \newcommand{\la}{\lambda} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\al}{\alpha} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \left(\pa_t - \Lie_{\be} \right) \ga_{ij} = -2 \al K_{ij} \,, \label{eq:EOMq} \nonumber \end{align} \tag{ 2.2 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Lie}{\mathcal{L}} \newcommand{\La}{\Lambda} \newcommand{\la}{\lambda} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\al}{\alpha} \newcommand{\na}{\nabla} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \left(\pa_t - \Lie_{\be} \right) K_{ij} = - \al \left[ 2 K_{ik} K_j^k - K_{ij} K - R_{ij} + \La \ga_{ij} + S_{ij} - \frac{1}{2}\, \ga_{ij} \left(S - E \right) \right] - \na_i \na_j \al \,, \label{eq:EOMK} \nonumber \end{align} \tag{ 2.3 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\La}{\Lambda} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle 0 = E + \frac{1}{2} \left(K_{ij} K^{ij} - K^2 - R \right) + \La \,, \label{eq:Hconst} \nonumber \end{align} \tag{ 2.4 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\na}{\nabla} \newcommand{\bo}{\square} \displaystyle 0 = P_i - \na_j K_i^{\,j} + \na_i K \,, \label{eq:Mconst} \nonumber \end{align} \tag{ 2.5 }$$

where we use $\newcommand{\ga}{\gamma} \ga_{ij}$ to displace the indices, $K_{ij}$ is the extrinsic curvature, $\newcommand{\na}{\nabla} \na$ and $R_{ij}$ are the connection and Ricci tensor of $\newcommand{\ga}{\gamma} \ga_{ij}$, $\newcommand{\La}{\Lambda} \La$ is the cosmological constant and E, P_i, $S_{ij}$ are the canonical matter energy, momentum and stress densities, respectively. These are the components of the energy-momentum tensor $T_{\mu\nu}$ in the n-frame

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\si}{\sigma} \newcommand{\ro}{\rho} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \label{eq:EPSdef} E := n^{\mu} n^{\nu} T_{\mu\nu},\quad P^{\mu} := - h^{\mu\nu} n^{\ro} T_{\nu\ro}, \quad S^{\mu\nu} := h^{\mu\ro} h^{\nu\si} T_{\ro\si}, \nonumber \end{align} \tag{ 2.6 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle h_{\mu\nu} := g_{\mu\nu} + n_{\mu} n_{\nu} \,, \nonumber \end{align} \tag{ 2.7 }$$

is the associated projector. Incidentally, $\newcommand{\e}{{\rm e}} P^0, S^{0\mu} \equiv 0$, while $\newcommand{\ga}{\gamma} P_i := \ga_{ij} P^{\,j}$ and $\newcommand{\ga}{\gamma} S_{ij} := \ga_{ik} \ga_{jl} S^{kl}$. Assuming the presence of an independently conserved PF, we can split $\newcommand{\e}{{\rm e}} T_{\mu\nu} \equiv \mathcal{T}_{\mu\nu} + T^{\rm extra}_{\mu\nu}$, where the curly letters ($\mathcal{E}$, $\mathcal{P}_i$, $\mathcal{S}_{ij}$) denote the PF part, while the 'extra' superscript labels the extra matter content. We thus have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ro}{\rho} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \label{eq:TPF} \mathcal{T}^{\mu\nu} := \left(\ro + p \right) U^{\mu} U^{\nu} + p\, g^{\mu\nu} \,, \nonumber \end{align} \tag{ 2.8 }$$

where $U^{\mu}$ is the fluid's 4-velocity and $\newcommand{\ro}{\rho} \ro$ and p are the energy density and pressure in the fluid's rest-frame U. Using (2.6) for the PF sector we then get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ro}{\rho} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \label{eq2:9} \mathcal{E} = W^2 \left(\ro + p \right) - p \,, \nonumber \end{align} \tag{ 2.9 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ro}{\rho} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \label{qe2:10} \mathcal{P}^i = W^2 \left(\ro + p \right) v^i \,, \nonumber \end{align} \tag{ 2.10 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ro}{\rho} \newcommand{\la}{\lambda} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \mathcal{S}^{ij} = W^2 \left(\ro + p \right) v^i v^{\,j} + p \ga^{ij} \equiv \frac{\mathcal{P}^i \mathcal{P}^{\,j}}{\mathcal{E} + p} + p \ga^{ij} \,, \label{eq:Hij} \nonumber \end{align} \tag{ 2.11 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\al}{\alpha} \newcommand{\bo}{\square} \displaystyle v^i := \al^{-1} \left(\frac{U^i}{U^0} - \frac{n^i}{n^0} \right),\quad W := -n_{\mu} U^{\mu} \equiv \left[ 1 - v_i v^i \right]^{-1/2} \,, \nonumber \end{align} \tag{ 2.12 }$$

are the proper time 3-velocity measured by the canonical observer $n^{\mu}$ and the generalization of the corresponding Lorentz factor, respectively. In particular, we can invert

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ro}{\rho} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \label{eq:roE} \ro = \mathcal{E} - \frac{\mathcal{P}_i \mathcal{P}^i}{\mathcal{E} + p},\,\, v^i = \frac{\mathcal{P}^i}{\mathcal{E} + p} \, . \nonumber \end{align} \tag{ 2.13 }$$

Using the above $\newcommand{\ro}{\rho} \ro\, (\mathcal{E}, \mathcal{P}, p)$ relation, along with the equation of state $\newcommand{\ro}{\rho} p = p(\ro)$, we can express p in terms of $\mathcal{E}$ and $\mathcal{P}_i$. For instance, in the case of constant equation of state parameter $\newcommand{\ro}{\rho} w := p/\ro$, we find

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle p = \frac{1}{2} \left[ \sqrt{\left(1 + w \right)^2 \mathcal{E}^2 - 4w \mathcal{P}_i \mathcal{P}^i} - \left(1 - w \right) \mathcal{E} \right] \, . \nonumber \end{align} \tag{ 2.14 }$$

Equation (2.11) then implies that all PF quantities can be expressed algebraically in terms of $\mathcal{E}$ and $\mathcal{P}_i$, i.e. the degrees of freedom of a PF, a fact that is important to keep in mind in what follows.

In the most common NR schemes, the so-called 'free evolution' schemes, and in particular in the cosmological simulations performed so far, one solves numerically the evolution equations, i.e. the ones that are differential equations in time. Here this corresponds to solving equations (2.2) and (2.3) for the gravitational sector. For a generic matter sector, one needs to specify the matter action and derive the corresponding equations of motion. However, in the case of PF matter, the fact that there are only four degrees of freedom $\mathcal{E}$ and $\mathcal{P}_i$ means that their evolution is entirely determined by the energy-momentum conservation equation

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\na}{\nabla} \newcommand{\bo}{\square} \displaystyle \label{eq:EMTcons} \na_{\mu} \mathcal{T}^{\mu\nu} = 0 \,, \nonumber \end{align} \tag{ 2.15 }$$

since we have assumed that the PF is minimally coupled and thus conserved independently of $T_{\mu\nu}^{\rm extra}$. Indeed, using the above translations between n-frame and U-frame components, we can express (2.15) in terms of $\mathcal{E}$ and $\mathcal{P}^i$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Lie}{\mathcal{L}} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\al}{\alpha} \newcommand{\na}{\nabla} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \left(\pa_t - \Lie_{\be} \right) \mathcal{E} = - \na_i(\al \mathcal{P}^i) - \mathcal{P}^i \na_i \al + \al \left(K_{ij} \mathcal{S}^{ij} + K \mathcal{E} \right) \,, \label{eq:evolE} \nonumber \end{align} \tag{ 2.16 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Lie}{\mathcal{L}} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\al}{\alpha} \newcommand{\na}{\nabla} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \left(\pa_t - \Lie_{\be} \right) \mathcal{P}^i = - \na_j (\al \mathcal{S}^{ij}) - \mathcal{E} \na^i \al + \al \left(2 K^i_j \mathcal{P}^{\,j} + K \mathcal{P}^i \right) \,, \label{eq:evolP} \nonumber \end{align} \tag{ 2.17 }$$

which form a closed set thanks to (2.11) and the equation of state. With the evolution equations being the ones that are solved numerically, the constraint equations (2.4) and (2.5) need only to be solved at the initial data surface and are then monitored at each time-step⁶, thus providing a first and valuable stability check of the simulation.

Here we propose the opposite choice for dealing with the PF field. Instead of using the four energy-momentum conservation equations (2.16) and (2.17) to evolve $\mathcal{E}$ and $\mathcal{P}_i$ in time, we choose to use the four constraint equations (2.4) and (2.5) to determine $\mathcal{E}$ and $\mathcal{P}_i$ algebraically out of the rest of the involved fields at each time-step

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\La}{\Lambda} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \mathcal{E} = - \frac{1}{2} \left(K_{ij} K^{ij} - K^2 - R \right) - \La - E^{\rm extra} \,, \label{eq:H} \nonumber \end{align} \tag{ 2.18 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\na}{\nabla} \newcommand{\bo}{\square} \displaystyle \mathcal{P}_i = \na_j K_i^{\,j} - \na_i K - P^{\rm extra}_i \, . \label{eq:Hi} \nonumber \end{align} \tag{ 2.19 }$$

With this algebraic substitution, the PF fields are totally eliminated and the resulting system is unconstrained⁷. In the standard ADM approach in which we work here, our scheme would therefore correspond to solving (2.2), (2.3) and the evolution equations of the 'extra' matter fields, where (2.18) and (2.19) have been used to eliminate the PF fields. These last two equations can then be used to obtain the PF information out of the fields that are actually been evolved at any time. As a consequence, it is the PF evolution equations (2.16) and (2.17) that are now redundant, i.e. they are no longer needed for closing the system. In analogy with the free evolution schemes, one could instead monitor these equations to check numerical stability. As already mentioned in the introduction, however, the study of the accuracy of the PF dynamics using this indirect method is postponed to future work. The only check we will perform here is related to the FLRW solution.

From now on let us focus on the case $T^{\rm extra}_{\mu\nu} = 0$ for simplicity. We see that $\newcommand{\ga}{\gamma} \ga_{ij}$, $K_{ij}$ parametrize the constraint hypersurface in phase space of the full GR  +  PF system, i.e. we only have physical states by construction. The numerical constraint violation of free evolution schemes is now translated into an error on the PF quantities $\mathcal{E}$ and $\mathcal{P}_i$ and therefore to a violation of energy-momentum conservation $\newcommand{\na}{\nabla} \na_{\mu} \mathcal{T}^{\mu\nu} \neq 0$, i.e. the relation between two subsequent states is not necessarily physical. It therefore seems that we have only traded one kind of error for another. However, what we did gain is that now the PF fields no longer participate in the robustness test, since they are not evolved in the integration to begin with. They are only physically relevant combinations of the fields that we are actually solving numerically. This way we avoid the growing modes discussed in the introduction (e.g. the matter density contrast), that are present already at the analytical linear level. Since the gravitational perturbations are bounded in time in the appropriate gauge, we are now able to set up a robustness test around a FLRW solution.

Let us finally get a bit more insight into how the gravitational and PF degrees of freedom (DOF) coexist inside $\newcommand{\ga}{\gamma} \ga_{ij}$ and $K_{ij}$. To that end, we linearize (2.2) and (2.3) around a flat FLRW background, using (2.11), (2.18) and (2.19). Being only interested in the spectrum, we look at small enough scales/periods and can thus effectively set the scale factor to 1. Moreover, we can also integrate out the $K_{ij}$ variables, thus expressing the equations in second-order form. We perform a scalar-vector-tensor (SVT) decomposition

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\al}{\alpha} \newcommand{\bo}{\square} \displaystyle \label{e220} \al = 1 + \psi \,, \nonumber \end{align} \tag{ 2.20 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ch}{\chi} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \label{e221} \be^i = \pa^i \chi + \chi^i \,, \nonumber \end{align} \tag{ 2.21 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\si}{\sigma} \newcommand{\la}{\lambda} \newcommand{\de}{\delta} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \label{e222} \ga_{ij} = \left(1 - 2 \si \right) \de_{ij} + 2 \pa_i \pa_j h + 2 \pa_{(i} h_{j)} + 2 h_{ij} \,, \nonumber \end{align} \tag{ 2.22 }$$

i.e. all vectors/tensors have zero divergence and trace, and we displace the indices using $\newcommand{\de}{\delta} \de_{ij}$. We get two wave equations

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\si}{\sigma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \ddot{\si} = \bar{c}_s^2 \pa^2 \si,\quad \ddot{h}_{ij} = \pa^2 h_{ij} \,, \nonumber \end{align} \tag{ 2.23 }$$

where $\newcommand{\ed}{{\rm d}} \newcommand{\ro}{\rho} \newcommand{\e}{{\rm e}} c_s^2 := \ed p/ \ed \ro$ is the speed of sound and appears through the perturbation of p in $\mathcal{S}_{ij}$, using the fact that $\newcommand{\ro}{\rho} \ro = \mathcal{E}$ at the linear level (2.13). We see that the gravitational potential $\newcommand{\si}{\sigma} \si$ is now dynamical and carries the sound DOF. In vacuum GR we have that $\bar{c}_s^2 \to 1$, but $\newcommand{\si}{\sigma} \si$ has no DOF because of the Hamiltonian and (the divergence of) the momentum constraints. As for the h and h_i components, it is convenient to integrate in two auxiliary fields π and $\pi_i$ (with $\newcommand{\pa}{\partial} \newcommand{\e}{{\rm e}} \pa_i \pi_i \equiv 0$) in order to express the equations in first-order form

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ch}{\chi} \newcommand{\ph}{\phi} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \dot{h} = \pi + \chi \,,\quad \dot{\pi} = - \phi + \psi \,, \nonumber \end{align} \tag{ 2.24 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ch}{\chi} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \dot{h}_i = \pi_i + \chi_i \,,\quad \dot{\pi}_i = 0 \, . \nonumber \end{align} \tag{ 2.25 }$$

The π and $\pi_i$ fields are linear combinations of $\newcommand{\ga}{\gamma} \ga_{ij}$ and $K_{ij}$ perturbations, whose precise form is irrelevant for our purposes here. Indeed, the relevant feature is the structure of the equations. In particular, we see that there are no DOF in the $(h, \pi)$ couple, because we can trivialize the evolution of both fields by choosing appropriately χ and ψ. As for $(h_i, \pi_i)$, the h_i components can be trivialized with $\newcommand{\ch}{\chi} \chi_i$, so we only have two DOF in $\pi_i$, corresponding to the rotational velocity of the fluid. In vacuum GR the latter are neutralized by the rotational part of the momentum constraint. Thus, the four DOF of the fluid $\mathcal{E}$ and $\mathcal{P}_i$ lie in the lower helicity modes of the gravitational field that are now unconstrained. We are therefore basically evolving the PF through the gravitational fields.

One should also observe that this 'de-constraining' procedure is only possible in a cosmological context, because we need $\mathcal{E}+p > 0$ in (2.11). This, however, does not mean that we cannot handle local voids $\mathcal{E} = p = 0$. Indeed, assuming the weak energy condition $\newcommand{\ro}{\rho} p \geqslant -\ro$, we get from (2.13) that $(\mathcal{E} + p){\hspace{0pt}}^2 \geqslant \mathcal{P}_i \mathcal{P}^i$, so the numerator in $\mathcal{P}_i \mathcal{P}_j/(\mathcal{E} + p)$ tends faster to zero than the denominator. One can thus safely set $\newcommand{\ga}{\gamma} \mathcal{S}_{ij} \to \ga_{ij} p$ by hand as soon as $\mathcal{E}$ goes below some threshold value.

Finally, note that any generally-covariant theory of gravity can be deconstrained in some analogous manner, simply because $\mathcal{E}$ and $\mathcal{P}_i$ will always appear linearly in the diffeomorphism constraints. However, depending on the precise structure of the theory and extra matter content, this might not be enough to get rid of all growing linear modes in the fields that need to be evolved.

We start by rewriting the unconstrained GR  +  PF system, i.e. (2.2) and (2.3) with (2.11), (2.18) and (2.19), in the form which provides a better numerical behavior. We first decompose $\newcommand{\ga}{\gamma} \ga_{ij}$ and $K_{ij}$ in the conformal/traceless fashion

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \displaystyle \ga_{ij} \equiv {\rm e}^{4\ph} \ti{\ga}_{ij},\quad K_{ij} \equiv {\rm e}^{4\ph} \left(\ti{A}_{ij} + \frac{1}{3}\, \ti{\ga}_{ij} K \right) \,,\, \nonumber \end{align} \tag{ 3.1 }$$

where $\newcommand{\ti}{\tilde} \newcommand{\ga}{\gamma} \newcommand{\de}{\delta} \newcommand{\e}{{\rm e}} \det \ti{\ga}_{ij} \equiv 1$, $\newcommand{\ti}{\tilde} \newcommand{\ga}{\gamma} \newcommand{\e}{{\rm e}} \ti{\ga}^{ij} \ti{A}_{ij} \equiv 0$ and all indices are displaced using $\newcommand{\ti}{\tilde} \newcommand{\ga}{\gamma} \ti{\ga}_{ij}$ from now on. The above separation, also used in the BSSN scheme [54, 55], is especially important in cosmology because it distinguishes the fields with a non-trivial background $\newcommand{\ph}{\phi} \ph$, K from the rest $\newcommand{\ti}{\tilde} \newcommand{\ga}{\gamma} \ti{\ga}_{ij}$, $\newcommand{\ti}{\tilde} \ti{A}_{ij}$. We then consider the following slicing

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\al}{\alpha} \newcommand{\ed}{{\rm d}} \newcommand{\bo}{\square} \displaystyle \al = a^2 {\rm e}^{-2\ph},\quad a(t) := \left(\frac{\int \ed^3 x \, {\rm e}^{6 \ph}}{\int \ed^3 x} \right)^{1/3} \,, \nonumber \end{align} \tag{ 3.2 }$$

where the integrals are carried over the full spatial slice. On a FLRW space-time, both a and $\newcommand{\ph}{\phi} {\rm e}^{2\ph}$ (and thus $\newcommand{\al}{\alpha} \al$) reduce to the scale factor, so this is a non-linear generalization of conformal time. This choice is mostly out of convenience, because the linear propagating fluctuations have a constant frequency with respect to that time.

Next, we consider the gauge $\newcommand{\pa}{\partial} \newcommand{\ti}{\tilde} \newcommand{\ga}{\gamma} \newcommand{\Ga}{\Gamma} \ti{\Ga}^i := - \pa_j \ti{\ga}^{ij} = 0$, whose conservation $\newcommand{\ti}{\tilde} \newcommand{\Ga}{\Gamma} \dot{\ti{\Ga}}^i = 0$ leads to an elliptic equation for $\newcommand{\bo}{\square} \newcommand{\et}{\eta} \newcommand{\e}{{\rm e}} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\be}{\beta} \be^i$. In practice, however, the modified condition $\newcommand{\ti}{\tilde} \newcommand{\al}{\alpha} \newcommand{\Ga}{\Gamma} \newcommand{\la}{\lambda} \dot{\ti{\Ga}}^i - \al \la K \ti{\Ga}^i = 0$ with $\newcommand{\la}{\lambda} \la > 0$ helps the convergence of numerical deviations from $\newcommand{\ti}{\tilde} \newcommand{\Ga}{\Gamma} \ti{\Ga}^i = 0$, and the detailed equation reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\la}{\lambda} \newcommand{\Ga}{\Gamma} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\al}{\alpha} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \label{eq:beeq} \ti{\pa}^2 \be^i + \frac{1}{3}\, \ti{\ga}^{ij} \pa_j \pa_k \be^k - 2 \al \left(\pa_j \ti{A}^{ij} - 2 \ti{A}^{ij} \pa_j \ph \right) - \al \la K \ti{\Ga}^i = 0 \,, \nonumber \end{align} \tag{ 3.3 }$$

where $\newcommand{\pa}{\partial} \newcommand{\ti}{\tilde} \newcommand{\ga}{\gamma} \ti{\pa}^2 := \ti{\ga}^{ij} \pa_i \pa_j$ and we have used $\newcommand{\ti}{\tilde} \newcommand{\Ga}{\Gamma} \ti{\Ga}^i = 0$ to simplify the $\newcommand{\ti}{\tilde} \newcommand{\Ga}{\Gamma} \newcommand{\si}{\sigma} \sim \dot{\ti{\Ga}}^i$ part. The evolution equations (2.2) and (2.3) now read

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\al}{\alpha} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \dot{\ph} = - \frac{1}{6} \, \al K + \be^i \pa_i \ph + \frac{1}{6}\, \pa_i \be^i \,, \label{eq:phidot} \nonumber \end{align} \tag{ 3.4 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\lambda} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\al}{\alpha} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \dot{\ti{\ga}}_{ij} = - 2 \al \ti{A}_{ij} + \be^k \pa_k \ti{\ga}_{ij} + \ti{\ga}_{kj} \pa_i \be^k + \ti{\ga}_{ik} \pa_j \be^k - \frac{2}{3}\, \ti{\ga}_{ij} \pa_k \be^k \,, \label{eq:gadot} \nonumber \end{align} \tag{ 3.5 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\La}{\Lambda} \newcommand{\la}{\lambda} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\al}{\alpha} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \dot{K} = \frac{\al}{2} \left[ K^2 - 3 \La + \frac{3}{2}\, \ti{A}_{ij} \ti{A}^{ij} + {\rm e}^{-4\ph} \left(\frac{1}{2}\,B + \ti{\ga}^{ij} \mathcal{S}_{ij} \right) \right] + \be^i \pa_i K \,, \label{eq:Kdot} \nonumber \end{align} \tag{ 3.6 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\al}{\alpha} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \dot{\ti{A}}_{ij} = \al \left[ K \ti{A}_{ij} - 2 \ti{A}_{ik} \ti{A}_j^k + {\rm e}^{-4\ph} \left(B_{ij} - \mathcal{S}_{ij} \right)^{\rm T} \right] + \be^k \pa_k \ti{A}_{ij} + \ti{A}_{kj} \pa_i \be^k + \ti{A}_{ik} \pa_j \be^k - \frac{2}{3}\, \ti{A}_{ij} \pa_k \be^k \,, \label{eq:Adot} \nonumber \end{align} \tag{ 3.7 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\Ga}{\Gamma} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle B_{ij} := - \frac{1}{2}\, \ti{\pa}^2 \ti{\ga}_{ij} + \ti{\ga}^{kl} \ti{\ga}^{mn} \left(2 \ti{\Ga}_{km(i} \ti{\Ga}_{j)ln} + \ti{\Ga}_{kmi} \ti{\Ga}_{lnj} \right) - 8 \pa_i \ph \pa_j \ph \,, \nonumber \end{align} \tag{ 3.8 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\la}{\lambda} \newcommand{\Ga}{\Gamma} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle B := \ti{\ga}^{ij} B_{ij} = \ti{\ga}^{ij} \left(\ti{\Ga}^k_{li} \ti{\Ga}^l_{kj} - 8 \pa_i \ph \pa_j \ph \right) \,, \label{eq:B} \nonumber \end{align} \tag{ 3.9 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \displaystyle \mathcal{S}_{ij} \equiv {\rm e}^{4\ph} \ti{\ga}_{ij} p + \frac{\mathcal{P}_i \mathcal{P}_j}{\mathcal{E} + p} \,, \nonumber \end{align} \tag{ 3.10 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\La}{\Lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \mathcal{E} = - \frac{1}{2}\, \ti{A}_{ij} \ti{A}^{ij} + \frac{1}{3}\, K^2 - \La + \frac{1}{2}\, {\rm e}^{-4\ph} \left(B - 8 \ti{\pa}^2 \ph \right) \,, \nonumber \end{align} \tag{ 3.11 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\Ga}{\Gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \mathcal{P}_i = \pa_j \ti{A}_i^{\,j} - \ti{A}^{\,j}_k \ti{\Ga}^k_{ji} + 6 \ti{A}_i^{\,j} \pa_j \ph - \frac{2}{3}\, \pa_i K \,, \nonumber \end{align} \tag{ 3.12 }$$

'T' denotes the traceless part and $\newcommand{\ti}{\tilde} \newcommand{\ga}{\gamma} \newcommand{\Ga}{\Gamma} \newcommand{\e}{{\rm e}} \ti{\Ga}^i_{jk} \equiv \ti{\ga}^{il} \ti{\Ga}_{ljk}$ are the Christoffel symbols of $\newcommand{\ti}{\tilde} \newcommand{\ga}{\gamma} \ti{\ga}_{ij}$. Thanks to $\newcommand{\al}{\alpha} \newcommand{\si}{\sigma} \newcommand{\ph}{\phi} \al \sim {\rm e}^{-2\ph}$ and $\newcommand{\ti}{\tilde} \newcommand{\Ga}{\Gamma} \ti{\Ga}^i = 0$, the only double spatial derivatives that appear are conformal Laplacians $\newcommand{\pa}{\partial} \newcommand{\ti}{\tilde} \ti{\pa}^2$, a 'hyperbolic' characteristic which could be responsible for the good numerical performance.

Here it is important to stress that, although we do use the conformal/traceless decomposition, which is one of the defining features of the BSSN scheme, there is another crucial feature of BSSN which cannot be implemented here and therefore clearly distinguishes the two approaches. Indeed, in BSSN one evolves independently $\newcommand{\ti}{\tilde} \newcommand{\Ga}{\Gamma} \ti{\Ga}^i$, considering its definition as an extra constraint equation, and one uses the momentum constraint to alter its evolution equation. This last manipulation is 'essential for the numerical stability of the system' [55], because it is needed to make the system strongly hyperbolic [56, 57]. In our case, the fact that the constraint equations have been used to eliminate part of the fields means that we can no longer use them at all, since this would reintroduce the said fields in the equations. We are therefore far from the defining structure of the BSSN scheme on which many of its merits are based.

Let us now show that there are no growing linear modes, at the analytical level, in the evolved fields in this gauge for the simpler case of constant equation of state parameter $\newcommand{\ro}{\rho} w := p/\ro$. We thus linearize around FLRW and perform a SVT decomposition

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \label{eq13} \ph = \frac{1}{2} \left[ \log a + h \right] \,, \nonumber \end{align} \tag{ 3.13 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\lambda} \newcommand{\de}{\delta} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \label{eq14} \ti{\ga}_{ij} = \de_{ij} + 2 \left[ \left(\pa_i \pa_j - \frac{1}{3}\, \de_{ij} \pa^2 \right) \ti{h} + \pa_{(i} \ti{h}_{j)} + \ti{h}_{ij} \right] \,, ~~ \nonumber \end{align} \tag{ 3.14 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\lambda} \newcommand{\ka}{\kappa} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle K = - 3 H \left[ 1 + \ka \right] \,, \label{eq:Kpert} \nonumber \end{align} \tag{ 3.15 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\lambda} \newcommand{\ka}{\kappa} \newcommand{\de}{\delta} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \ti{A}_{ij} = - H \left[ \left(\pa_i \pa_j - \frac{1}{3}\, \de_{ij} \pa^2 \right) \ti{\ka} + \pa_{(i} \ti{\ka}_{j)} + \ti{\ka}_{ij} \right] \,, \label{eq:Aijpert} \nonumber \end{align} \tag{ 3.16 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\al}{\alpha} \newcommand{\bo}{\square} \displaystyle \label{eq17} \al = a \left(1 - h \right) \,, \nonumber \end{align} \tag{ 3.17 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\lambda} \newcommand{\ka}{\kappa} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \label{eq18} \be^i = - a H \left[ \left(\pa^i \ti{\ka} + \ti{\ka}^i \right) + \la \left(\pa^i \ti{h} + \ti{h}^i \right) \right] \,, \nonumber \end{align} \tag{ 3.18 }$$

where $a(t)$ is the scale factor and $\newcommand{\pa}{\partial} H(t) := (a^{-1} \pa_t) \log a$ is the physical Hubble parameter in conformal time. If we further set $\newcommand{\ti}{\tilde} \newcommand{\Ga}{\Gamma} \ti{\Ga}^i = 0$, then $\newcommand{\ti}{\tilde} \ti{h}$ and $\newcommand{\ti}{\tilde} \ti{h}_i$ vanish. However, we keep these components to verify that the system maintains their amplitude close to numerical error. Going to Fourier space and choosing the time variable $x := \log a$ the equations become

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \newcommand{\pa}{\partial} \displaystyle \pa_x S = M_s S,\quad \pa_x V = M_v V,\,\, \pa_x T = M_t T,\,\,\, \nonumber \end{align} \tag{ 3.19 }$$

for the scalar $\newcommand{\ti}{\tilde} \newcommand{\ka}{\kappa} S := (h, \ka, \ti{h}, \ti{\ka})$, vector $\newcommand{\ti}{\tilde} \newcommand{\ka}{\kappa} V := (\ti{h}_i, \ti{\ka}_i)$ and tensor $\newcommand{\ti}{\tilde} \newcommand{\ka}{\kappa} T := (\ti{h}_{ij}, \ti{\ka}_{ij})$ multiplets with

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\la}{\lambda} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \displaystyle M_s = \left(\begin{array}{@{}cccc@{}} -1 & 1 & \la k^2/3 & k^2/3 \nonumber \\ - w\ti{k}^2 & - X & 0 & 0 \nonumber \\ 0 & 0 & -\la & 0 \nonumber \\ 0 & 0 & -\ti{k}^2 & 3w - X \end{array} \right), \quad M_v = \left(\begin{array}{@{}cc@{}} -\la & 0 \nonumber \\ -\ti{k}^2 & 3w - X \end{array} \right) \,,\quad M_t = \left(\begin{array}{@{}cc@{}} 0 & 1 \nonumber \\ -\ti{k}^2 & 3w - X \end{array} \right) \,, \nonumber \end{align} \tag{ 3.20 }$$

where $\newcommand{\ti}{\tilde} \ti{k} := k/(a H)$, $\newcommand{\La}{\Lambda} \newcommand{\Om}{\Omega} X := (3/2) \left(1+w \right) \left(1 + \Om_{\La} \right)$ and $\newcommand{\La}{\Lambda} \newcommand{\Om}{\Omega} \Om_{\La} := \La/(3 H^2) < 1$. All eigenvalues have negative real part for the typical cases of interest $0 \leqslant w \leqslant 1/3$, $\newcommand{\La}{\Lambda} \newcommand{\Om}{\Omega} \Om_{\La} \geqslant 0$ and $\newcommand{\la}{\lambda} \la > 0$. Moreover, the determinant of the eigenvector matrices is generically non-zero, so that M_s, M_v and M_t are diagonalizable and one can thus trust the spectrum to deduce the behavior.

We first consider integrating the homogeneous and isotropic solution with zero spatial curvature. For the precise FLRW solution we will consider, the non-trivial fields are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \ph_{\rm FLRW}(t) = \frac{1}{2}\, \log a_{\rm FLRW}(t) = \log \left[ \frac{t}{t_0} + \sqrt{a_i} \left(1 - \frac{t}{t_0} \right) \right] \,, \nonumber \end{align} \tag{ 4.1 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle K_{\rm FLRW}(t) = - 3 H_{\rm FLRW}(t) = - \frac{6 \left(1 - \sqrt{a_i} \right)}{t + \sqrt{a_i} \left(t_0 - t \right)} \,, \nonumber \end{align} \tag{ 4.2 }$$

where the 'i' and '0' subscripts denote the initial and final times, respectively, and in particular we have chosen $t_i = 0$. Incidentally, the initial and final values of the scale factor are a_i and a₀, and here we have $a_0 = 1$. The value of t₀ is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle t_0 = 2 H_0^{-1} \left(1 - \sqrt{a_i} \right) \,, \nonumber \end{align} \tag{ 4.3 }$$

where H₀ is the final physical Hubble parameter and we set $H_0^{-1} = 3000$. Finally, we choose $a_i = 1/50$, meaning that we evolve from redshift 49 to 0, and leading to $t_0 \approx 5150$. The initial conditions of our integration being

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\de}{\delta} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \displaystyle \ph(0, \vec{x}) = \ph_{\rm FLRW}(0), \,\,K(0, \vec{x}) = K_{\rm FLRW}(0) \,, \quad \ti{\ga}_{ij}(0, \vec{x}) = \de_{ij} \,, \,\,\ti{A}_{ij}(0,\vec{x}) = 0 \,, \nonumber \end{align} \tag{ 4.4 }$$

the numerical evolution preserves both the lack of $\vec{x}$-dependence as well as the trivial values of $\newcommand{\ti}{\tilde} \newcommand{\ga}{\gamma} \ti{\ga}_{ij}$ and $\newcommand{\ti}{\tilde} \ti{A}_{ij}$. Thus, the non-trivial numerical information is solely $\newcommand{\ph}{\phi} \ph(t)$ and $K(t)$, whose values we pick at some arbitrary point $\vec{x}$ on the grid. More precisely, we are interested in the relative errors with respect to the analytical solution

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\de}{\delta} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \de_{\ph}(t) := \vert \ph(t) - \ph_{\rm FLRW}(t) \vert \,, \quad \de_K(t) := \left\vert \frac{K(t) - K_{\rm FLRW}(t)}{K_{\rm FLRW}(t)} \right\vert \,, \nonumber \end{align} \tag{ 4.5 }$$

where $\newcommand{\de}{\delta} \newcommand{\ph}{\phi} \de_{\ph}$ is indeed 'relative' since the actual gravitational field is $\newcommand{\ga}{\gamma} \newcommand{\si}{\sigma} \newcommand{\ph}{\phi} \ga_{ij} \sim {\rm e}^{4\ph}$. On the PF side, the only non-trivial quantity here is the energy density $\mathcal{E}(t)$, since $\newcommand{\e}{{\rm e}} p \equiv 0$, and through (2.18) we get $\mathcal{E} = K^2/3$. Therefore, the relative error on the energy density is of the same order

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\de}{\delta} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \de_\mathcal{E}(t) := \frac{\vert \mathcal{E}(t) - \mathcal{E}_{\rm FLRW}(t)\vert }{\mathcal{E}_{\rm FLRW}(t)} \approx 2 \de_K(t) \,, \nonumber \end{align} \tag{ 4.6 }$$

by construction. We perform the computation for three time resolutions $\newcommand{\De}{\Delta} \De t = 3.2 / N$, where $N \in \{1, 2, 4 \}$, and plot the normalized deviations $\newcommand{\de}{\delta} \newcommand{\ph}{\phi} N^4 \de_{\ph}$ and $\newcommand{\de}{\delta} N^4\de_K$ in figure 1. We first note that the relative errors are bounded in time and exhibit the same behavior for all resolutions. The fact that the curves decrease in magnitude as N increases means that we have achieved fourth-order global convergence, as expected with the RK4 method.

Zoom In Zoom Out Reset image size

In section 3.2 we demonstrated the absence of growing linear modes around the FLRW solution in theory, so the remaining non-trivial requirement is that the fully non-linear system is able to respect that condition numerically. We therefore consider the FLRW initial data plus random perturbation fields

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\vep}{\varepsilon} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \ph(0, \vec{x}) = \ph_{\rm FLRW}(0) + \vep^{\ph}(\vec{x}) \,, \nonumber \end{align} \tag{ 4.7 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\vep}{\varepsilon} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle K(0, \vec{x}) = K_{\rm FLRW}(0) \left[ 1 + \vep^K(\vec{x}) \right] \,, \nonumber \end{align} \tag{ 4.8 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\vep}{\varepsilon} \newcommand{\de}{\delta} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \displaystyle \ti{\ga}_{ij}(0, \vec{x}) = \frac{\de_{ij} + \vep^{\ga}_{ij}(\vec{x})}{\left[ \det \left(\de_{kl} + \vep^{\ga}_{kl}(\vec{x}) \right) \right]^{1/3}} \,, \nonumber \end{align} \tag{ 4.9 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\vep}{\varepsilon} \newcommand{\de}{\delta} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \displaystyle \ti{A}_{ij}(0,\vec{x}) = \left[ \de_i^k \de_j^l - \frac{1}{3}\, \ti{\ga}_{ij}(0,\vec{x})\, \ti{\ga}^{kl}(0,\vec{x}) \right] \vep^{A}_{kl}(\vec{x}) \,, \nonumber \end{align} \tag{ 4.10 }$$

where the entries of $\newcommand{\vep}{\varepsilon} \newcommand{\ph}{\phi} \vep^{\ph, K}$ and $\newcommand{\ga}{\gamma} \newcommand{\vep}{\varepsilon} \newcommand{\e}{{\rm e}} \vep^{\ga, A}_{ij} \equiv \vep^{\ga, A}_{ji}$ are drawn with a uniform distribution of amplitude $\newcommand{\vep}{\varepsilon} \vep$ independently at each space point and for each field component, using the GSL random generator with a different seed for each MPI process. We set periodic spatial boundary conditions and consider a fixed comoving box size of $L = 1024$, so that we start outside the horizon at early times and end up inside at late times. We perform the test for three resolutions, 64³, 128³ and 256³ grid points. We pick $\newcommand{\De}{\Delta} \newcommand{\vep}{\varepsilon} \vep = 10^{-12} \De x^2$, where $\newcommand{\De}{\Delta} \De x \in \{16, 8, 4 \}$ is the lattice spacing, so that the matter density contrast $\newcommand{\pa}{\partial} \newcommand{\de}{\delta} \newcommand{\vep}{\varepsilon} \newcommand{\si}{\sigma} \newcommand{\ph}{\phi} \de \sim \vep\, \pa^2 \ph$ is of the same order for all resolutions. It reaches a maximum value of $\newcommand{\de}{\delta} \newcommand{\si}{\sigma} \newcommand{\Ord}{\mathcal{O}} \newcommand{\Or}{{\rm O}} \de \sim \Ord(10^{-4})$ at t₀, so the linear regime is maintained. Finally the Courant factor is set to 10, i.e. $\newcommand{\De}{\Delta} \De t = \De x/10$, since the speed of light is unity. In terms of $\newcommand{\De}{\Delta} \De t$, our three resolutions therefore correspond to $\newcommand{\De}{\Delta} \De t \in \{1.6, 0.8, 0.4 \}$.

We next define the following measures of deviation from the analytical FLRW solution

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ph}{\phi} \newcommand{\de}{\delta} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\bo}{\square} \displaystyle \de_{\ph}(t) := \max_{\vec{x}} \vert \ph(t,\vec{x}) - \ph_{\rm FLRW}(t) \vert, \quad \de_K(t) := \max_{\vec{x}} \left\vert \frac{K(t,\vec{x}) - K_{\rm FLRW}(t)}{K_{\rm FLRW}(t)} \right\vert \,, \nonumber \end{align} \tag{ 4.11 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\De}{\Delta} \newcommand{\de}{\delta} \newcommand{\ga}{\gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \displaystyle \de_{\ga}(t) := \max_{\vec{x}} \sqrt{\De \ti{\ga}_{ij}(t,\vec{x})\, \De \ti{\ga}^{ij}(t,\vec{x})}, \quad \de_{A}(t) := \max_{\vec{x}} \sqrt{\ti{A}_{ij}(t,\vec{x})\, \ti{A}^{ij}(t,\vec{x})} \,, \nonumber \end{align} \tag{ 4.12 }$$

where $\newcommand{\ti}{\tilde} \newcommand{\ga}{\gamma} \newcommand{\de}{\delta} \newcommand{\De}{\Delta} \De \ti{\ga}_{ij} := \ti{\ga}_{ij} - \de_{ij}$, and the gauge violation

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\de}{\delta} \newcommand{\Ga}{\Gamma} \newcommand{\be}{\beta} \newcommand{\bet}{\boldsymbol\eta} \newcommand{\et}{\eta} \newcommand{\ti}{\tilde} \newcommand{\bo}{\square} \displaystyle \de_{\Ga}(t) := \max_{\vec{x}} \sqrt{\ti{\Ga}_i(t,\vec{x})\, \ti{\Ga}^i(t,\vec{x})} \,, \nonumber \end{align} \tag{ 4.13 }$$

respectively. Note that, in the case of $\newcommand{\ti}{\tilde} \newcommand{\ga}{\gamma} \ti{\ga}_{ij}$ and $\newcommand{\ti}{\tilde} \ti{A}_{ij}$, we can only consider absolute deviations since their FLRW values are trivial. In figures 2 and 3 we plot the evolution of these quantities for all three resolutions. We see that the perturbations are indeed bounded in time and, most importantly, that we have convergence under resolution increase. More precisely, here it is not the amplitude of the deviation that is relevant, because we are decreasing its initial value by hand as we increase the resolution $\newcommand{\De}{\Delta} \newcommand{\vep}{\varepsilon} \newcommand{\si}{\sigma} \vep \sim \De x^2$. Rather, what matters is that the behavior of the deviation remains qualitatively the same for all three resolutions. As for the violation of the gauge condition $\newcommand{\Ga}{\Gamma} \newcommand{\de}{\delta} \de_{\Ga}$, it is well controlled in that it decreases with time and its amplitude is several orders of magnitude below the one of $\newcommand{\ga}{\gamma} \newcommand{\de}{\delta} \de_{\ga}$. Finally, as in the previous subsection, here too the PF quantities have the same level of accuracy as the gravitational ones, since they are determined algebraically by the latter.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size