GRAMSES: a new route to general relativistic $N$-body simulations in cosmology. Part I. Methodology and code description

TL;DR

GRAMSES is a novel relativistic N-body simulation code based on the ADM formalism, enabling fully nonlinear, background-independent cosmological modeling with adaptive mesh refinement, suitable for studying backreaction effects and GR impacts on cosmology.

Contribution

It introduces a new pipeline for relativistic cosmological N-body simulations using the ADM formalism with AMR, expanding the capabilities of existing cosmological simulation tools.

Findings

Successfully implemented the methodology and code tests.

Demonstrated initial results in a $\\Lambda$CDM universe.

Provides a framework for studying relativistic effects in cosmic structures.

Abstract

We present GRAMSES, a new pipeline for nonlinear cosmological $N$-body simulations in General Relativity (GR). This code adopts the Arnowitt-Deser-Misner (ADM) formalism of GR, with constant mean curvature and minimum distortion gauge fixings, which provides a fully nonlinear and background independent framework for relativistic cosmology. Employing a fully constrained formulation, the Einstein equations are reduced to a set of ten elliptical equations which are solved using multigrid relaxation with adaptive mesh refinements (AMR), and three hyperbolic equations for the evolution of tensor degrees of freedom. The current version of GRAMSES neglects the latter by using the conformal flatness approximation, which allows it to compute the two scalar and two vector degrees of freedom of the metric. In this paper we describe the methodology, implementation, code tests and first results for cosmological simulations in a $\Lambda$CDM universe, while the generation of initial conditions and physical results will be discussed elsewhere. Inheriting the efficient AMR and massive parallelisation infrastructure from the publicly-available $N$-body and hydrodynamic simulation code RAMSES, GRAMSES is ideal for studying the detailed behaviour of spacetime inside virialised cosmic structures and hence accurately quantifying the impact of backreaction effects on the cosmic expansion, as well as for investigating GR effects on cosmological observables using cosmic-volume simulations.

Figures

Equations

$$d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=-{\alpha }^{2}d{t}^2+{\gamma }_{ij}({\beta }^{i}dt+d{x}^i)({\beta }^{j}dt+d{x}^j),$$

$$R+K^2-K_{ij}{K}^{ij}=16\pi \rho ,$$

$$D_j(K^{ij}-{\gamma }^{ij}K)=8\pi S^i,$$

$$\rho$$

$$S_i$$

$$({\partial }_t-{\mathcal{L}}_{\beta }){\gamma }_{ij}=$$

$$({\partial }_t-{\mathcal{L}}_{\beta })K_{ij}=$$

$${\mathcal{L}}_{\beta }{\gamma }_{ij}=D_i{\beta }_j+D_j{\beta }_i,$$

$${\mathcal{L}}_{\beta }{K}_{ij}={\beta }^k{\partial }_k{K}_{ij}+K_{ik}{\partial }_j{\beta }^k+K_{kj}{\partial }_i{\beta }^k,$$

$${\gamma }_{ij}={\psi }^4{\bar{\gamma }}_{ij},$$

$$A^{ij}$$

$$8{\bar{D}}^2\psi -\psi \bar{R}-\frac{2}{3}{\psi }^5{K}^2+{\psi }^{-7}{\bar{A}}_{ij}{\bar{A}}^{ij}$$

$${\bar{D}}_j{\bar{A}}^{ij}-\frac{2}{3}{\psi }^6{\bar{\gamma }}^{ij}{\bar{D}}_{j}K$$

$${\bar{A}}^{ij}={\bar{A}}_{TT}^{ij}+{\bar{A}}_L^{ij},$$

$${\bar{A}}_L^{ij}$$

$$({\bar{\Delta }}_{L}W{)}_i-\frac{2}{3}{\psi }^6{\bar{D}}_{i}K=8\pi {\psi }^6{S}_i,$$

$$s_i\equiv {\psi }^6{S}_i.$$

$$K=-3H(t),$$

$${\bar{D}}^2(\alpha \psi )=\alpha \psi \left[\frac{7}{8}{\psi }^{-8}{\bar{A}}_{ij}{\bar{A}}^{ij}+\frac{5}{12}{\psi }^4{K}^2+\frac{1}{8}\bar{R}+2\pi {\psi }^4(\rho +2S)\right]-{\psi }^5\dot{K}.$$

$$H^2$$

$$H^2+\dot{H}$$

$${\bar{A}}^{ij}=\frac{{\psi }^6}{2\alpha }\left[(\bar{L}\beta {)}^{ij}+{\partial }_t{\bar{\gamma }}^{ij}\right],$$

$$D^i({\gamma }^{1/3}{\partial }_t{\bar{\gamma }}_{ij})=0.$$

$$({\bar{\Delta }}_L\beta {)}^i+(\bar{L}\beta {)}^{ij}{\bar{D}}_j\ln {\psi }^6=2{\psi }^{-6}{\bar{A}}^{ij}{\bar{D}}_j\alpha +16\pi {\psi }^4\alpha S^i,$$

$${\bar{\gamma }}_{ij}={\delta }_{ij},{\bar{A}}_{TT}^{ij}=0\forall t.$$

$$({\bar{\Delta }}_L\beta {)}^i=2{\partial }_j\left(\alpha {\psi }^{-6}{\bar{A}}_L^{ij}\right).$$

$$({\bar{\Delta }}_{L}W{)}_i=8\pi s_i,$$

$${\bar{A}}_L^{ij}=(\bar{L}W{)}^{ij}\equiv {\bar{A}}^{ij},$$

$${\bar{\nabla }}^2\psi =-2\pi {\psi }^{-1}{s}_0-\frac{1}{8}{\psi }^{-7}{\bar{A}}_{ij}{\bar{A}}^{ij}+2\pi {\psi }^5{\hat{\rho }}_m,$$

$${\bar{\nabla }}^2(\alpha \psi )=\alpha \left[2\pi {\psi }^{-1}(s_0+2s)+\frac{7}{8}{\psi }^{-7}{\bar{A}}_{ij}{\bar{A}}^{ij}+{\psi }^5\left(\frac{5K^2}{12}+2\pi {\rho }_{\Lambda }\right)\right]-{\psi }^5\dot{K},$$

$$d{s}^2$$

$$u^{\mu }{\nabla }_{\mu }{u}_{\nu }=0,$$

$$\frac{d{u}_i}{dt}=-\alpha u^0{\partial }_i\alpha +u_j{\partial }_i{\beta }^j-\frac{u_j{u}_k}{2u^0}{\partial }_i{\gamma }^{jk},$$

Full text

1 Introduction

The challenge of implementing cosmological simulations based on General Relativity (GR) for studying structure formation in our Universe has received increasing attention in the last few years, not only as a more accurate counterpart to Newtonian simulations in the era of precision cosmology, but also as a natural framework to study phenomena beyond the Newtonian approximation. Large upcoming surveys, such as LSST [1], Euclid [2] and DESI [3], aim at mapping billions of galaxies with unprecedented precision and will open up the possibility of detecting new general relativistic effects in cosmological observables, whose precise modelling might require a new generation of cosmological simulations based on GR. As an example, the long-standing issue of the back-reaction effect [4], i.e., the impact of matter inhomogeneities on the cosmic expansion rate, requires a unified GR treatment for spacetime and matter. While the Newtonian approximation offers an intuitive description for gravity, and simulations based on it enjoy a remarkable success in explaining the observable Universe so far [5, 6, 7], the assumption of a Friedmann-Lemaître-Robertson-Walker (FLRW) background turns the Hubble parameter into and external input that is dissociated from the structure formation process that might take place. Furthermore, attempts to include this back-reaction in Newtonian simulations result in null contributions by construction [8, 9, 10]. Even if back-reaction may not fulfill its original purpose of explaining the observed accelerated cosmic expansion, it can still contribute to our understanding of dark energy and possibly to alleviate the current $3\sigma$ tension between the Hubble constant determination from the local [11] and distant [12] Universe. Recent attempts to study back-reaction effects in a GR framework show that local deviations from the average expansion rate of the Universe can be substantial in underdense regions [13] and might impact our estimation of cosmological parameters in a $\Lambda$CDM model even at the percent level [14]. In any case, a fully non-linear, general relativistic framework is needed in order to reach conclusive answers to this problem in the late-time Universe [15, 16, 17, 18].

The rise of modern $N$-body cosmological simulations taking a leap beyond the Newtonian approximation was first achieved by the gevolution code [19, 20, 21, 22, 23, 24]. This code is based on a weak-field expansion around a FLRW background, which allows to linearise the equations for the gravitational sector while keeping the nonlinearities in the energy-momentum tensor, whose components become ‘dressed’ by the weakly perturbed metric. Following the standards set by state-of-the-art Newtonian simulations, in gevolution the dark matter content of the Universe is described as a particle ensemble, which allows to track them even after trajectory-crossing occurs. The latter is crucial for getting a realistic structure formation history as this phenomenon appears in the formation of virialised dark matter haloes.

Another interesting path for implementing general relativistic cosmological simulations is to resort to numerical relativity which condenses the theoretical and numerical tools needed for modelling relativistic systems. Even if the first applications of numerical relativity to study cosmological spacetimes were explored during the 80’s [25], modern developments in this direction arguably started with the cosmograph code [26, 27] and the investigation of Betivegna & Bruni [13], followed by that of Macpherson et al. [28]. Some of these GR codes based on numerical relativity share a common feature of being developed upon the Einstein Toolkit [29], an open-source community infrastructure for relativistic astrophysics. These works implement the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation for a numerically stable spacetime evolution [30, 31, 32] – one of the cornerstones in modern numerical relativity – and general relativistic hydrodynamics for the matter sector. This latter aspect, however, is different from standard cosmological simulations and the approximation of treating dark matter as a fluid means that particle trajectories after shell crossing in bound structures are not exactly followed during the simulation. As an example, a comparison of fully non-linear GR simulations with pressureless dark matter fluid and Newtonian $N$-body simulations shows partial discrepancies inside the (Newtonian) dark matter haloes where the weak field condition is violated [33], but the different treatments for matter make the results difficult to assess and hence demands a fully GR cosmological $N$-body simulation [34].

The first implementations of three-dimensional $N$-body simulations based on numerical relativity were carried out in the late 90’s for the study of black holes formation [35], and various groups have implemented these in general relativistic cosmological codes in recent years. In Ref. [36] the full Einstein equations are solved using the BSSN formalism with different gauge conditions, and featuring a tri-cubic spline (TCS) scheme for increasing the smoothness of the density field described by the particle ensemble. The code introduced in Ref. [37] adopts the damped CCZ4 formalism [38, 39, 40] for solving the Einstein equations coupled to a leapfrog (Kick-Drift-Kick) scheme for particles, and can achieve a stable evolution of cosmological simulations up to redshift $z\sim 10$. The numerical relativity code from [41, 34] is well-established and solves the Einstein equations using a generalised harmonic formalism [42] on adaptive mesh refinements (AMR). This has been extensively applied for running simulations of both compact objects and, more recently, to cosmological scenarios. While numerical relativity codes allow to run full GR simulations, their numerical implementation and long-time stability are very challenging due to the presence of dynamical degrees of freedom (DOF) in the metric. It is well known that the evolution equations for the gravitational sector in the ADM formalism are only weakly hyperbolic and hence numerically unstable [43], a problem that in the BSSN formalism is fixed by the introduction of additional dynamical variables. In addition, solving the evolution equations accurately requires fine timesteps, which in the case of cosmological simulations can become several orders of magnitude smaller than the characteristic time scale of the system. However, assuming that gravitational waves (GW, which are the only dynamical DOF in GR) do not play a significant role in the cosmological dynamics and its back-reaction on the spacetime, it is possible to follow a formulation in which we can neglect these so that time evolution is only due to the matter sector. This is very much in the similar spirit as the ‘waveless theories of gravity’ developed by Isenberg [44] and latter by Wilson & Mathews [45], who sought the natural generalisation of Newtonian gravity within GR.

In this paper we introduce gramses, a new code which combines a set of features of state-of-the-art general relativistic codes aimed at fully non-linear and background-independent cosmological structure formation simulations. In order to deal with the nonlinear general relativistic equations for the gravitational sector in an optimal way for cosmology, gramses adopts a so-called fully constrained formulation of GR [46, 47], in which only elliptic-type partial differential equations (PDEs) are solved to reconstruct the spacetime metric in the absence of GWs [48, 49]. gramses solves these PDEs using multigrid Gauss-Seidel relaxation, and offers an $N$-body description for non-relativistic dark matter particles that supports AMR to increase force resolution in high-density regions, so that the cosmic web formed in the simulations can be resolved to a high degree of details even after shell crossing. Our code is based on the publicly-available, free-licensed $N$-body and hydrodynamical simulation code ramses [50], which is efficiently parallelised using Message Passing Interface (MPI). Particles are evolved along geodesics using a leapfrog method, and the particle-to-mesh projection and force interpolation is performed in a Cloud-In-Cell (CIC) scheme.

The rest of this paper is organised as follows. In Section 2 we introduce the equations for the gravitational sector in the fully constrained ADM formulation of GR, while in Section 3 we describe the matter sector. In Section 4 we discuss the code structure and the numerical implementations of the GR equations in gramses, and in Section 5 we present some test results for the calculation of the relativistic source terms as well as for the multigrid and geodesic solvers. Finally, in Section 6 we present some first results for cosmological simulations in a $\Lambda$CDM universe, although the generation of initial conditions and physical implications will be addressed in subsequent papers.

Throughout this paper we adopt the $(-,+,+,+)$ signature for the spacetime metric. Geometric units, where $G=c=1$, are used in the theory part of the paper for brevity, but this are restored in Section 4 in order to introduce code units (which are used in the actual equations solved in gramses). Greek indices run from 0 to 3, whereas Latin ones from 1 to 3 only, with repeated indices implying summation.

2 Field equations for the gravitational sector

In this section we review some fundamental aspects of numerical relativity that are at the core of this work. For a comprehensive discussion on these topics we refer the reader to [51, 52, 53].

2.1 The Arnowitt-Deser-Misner formalism

In the ADM (3+1) formalism [54], the 4-dimensional spacetime is foliated into 3-dimensional hypersurfaces of constant times characterised by some unit vector $n^{\mu}$, with which we can write a set of evolution and constraint equations for the variables $(\gamma_{ij},K_{ij})$, representing the (induced) spatial metric of the 3-dimensional embedded manifold and its extrinsic curvature, respectively. The spacetime metric in the ADM formalism is given by

[TABLE]

in which the lapse function $\alpha$ and shift vector $\beta^{i}$ represent gauge (or coordinate) DOF. The projection of the Einstein equations into the 3-dimensional hypersurfaces yields the Hamiltonian constraint and the momentum constraint, respectively given by

[TABLE]

where $K=\gamma^{ij}K_{ij}$ is the trace of the extrinsic curvature, $D_{i}$ the covariant derivative associated with the spatial metric $\gamma_{ij}$ and $R$ the Ricci scalar. Here, we have introduced the energy density $\rho$ and the momentum density $S_{i}$ measured by a normal observer $n^{\mu}$, which are calculated by projecting the energy-momentum tensor as

[TABLE]

where $n_{\mu}=(-\alpha,0)$. In addition to the constraint equations (2.2) and (2.3), which offer no dynamics, the evolution equations for $(\gamma_{ij},K_{ij})$ are

[TABLE]

where

[TABLE]

and

[TABLE]

correspond to Lie derivatives along $\beta^{i}$. In addition to the matter source terms (2.4) and (2.5), in Eq. (2.7) we have defined the spatial stress $S_{ij}\equiv\gamma_{i\mu}\gamma_{j\nu}T^{\mu\nu}$, with $S=\gamma^{ij}S_{ij}$ its trace.

In order to disentangle the physical and gauge DOF at the nonlinear level we can resort to decompose $(\gamma_{ij},K_{ij})$. As a method to single out a particular degree of freedom in the spatial metric ${\gamma}_{ij}$, we use the conformal transformation [55]

[TABLE]

where $\psi^{4}=\gamma^{1/3}$ is the conformal factor, $\bar{\gamma}_{ij}$ the conformal metric, and $\gamma\equiv\det\left(\gamma_{ij}\right)$ the determinant of the metric $\gamma_{ij}$. We also introduce a conformal transformation for the traceless part of the extrinsic curvature $A_{ij}\equiv K_{ij}-\gamma_{ij}K/3$ as

[TABLE]

Notice that for raising and lowering indices of the conformal (overbarred) quantities we use $\bar{\gamma}^{ij}$ and $\bar{\gamma}_{ij}$, respectively. Then, the Hamiltonian (2.2) and momentum (2.3) constraints can be rewritten as

[TABLE]

2.2 The conformal transverse traceless decomposition

By applying the Conformal Transverse Traceless (CTT) decomposition [56] (also known as York-Lichnerowicz formulation), we can further isolate DOF in the extrinsic curvature by decomposing the traceless, symmetric tensor $\bar{A}^{ij}$ into a transverse-traceless (TT) part which is divergenceless, and a longitudinal part that is written in terms as the vector gradient of a vector potential, namely

[TABLE]

with $\bar{D}_{j}\bar{A}^{ij}_{TT}=0$ and

[TABLE]

where $W^{j}$ can be regarded as a vector potential and $\bar{L}$ is the longitudinal operator (also known as vector gradient or conformal Killing operator). Then, the momentum constraint (2.12) can be written in terms of the vector potential as

[TABLE]

where $(\bar{\Delta}_{L}W)^{i}\equiv\bar{D}_{j}\bar{A}^{ij}$ is the vector Laplacian, and we have used that $\bar{\gamma}_{ij}\bar{\gamma}^{jk}=\delta^{k}_{i}$. A convenient feature about the previous equation is that we can decouple it from the rest in some particular cases, e.g., if we take $K=0$ (a maximal slicing) and identify the conformal source term as

[TABLE]

Notice that the TT part of $\bar{A}^{ij}$ is not constrained by (2.15); in fact $\bar{A}^{ij}_{TT}$ are dynamical DOF connected with gravitational waves. Then, if we want to find a solution in the absence of GW for a given 3-dimensional hypersurface at $t=t_{0}$ we might take $\bar{A}^{ij}_{TT}=0$, but since this is a dynamical quantity, we do not have the freedom to fix this again at a $t>t_{0}$ time slice but we would need to solve the evolution equations to propagate them. We will come back to this point in the next section when we introduce the ‘waveless approximation’.

2.3 Gauge fixing

As discussed in the previous subsection, the conformal decomposition and the CTT approach recast the constraint equations in a convenient form without any assumption about the way that the system evolves, and then the lapse function $\alpha$ and shift vector $\beta^{i}$ remain completely unspecified. Even if we have complete freedom to choose them as they correspond to picking a coordinate system, in practice not all options are physically or numerically convenient. For instance, the simplest option $\alpha=1$ and $\beta^{i}=0$, known as geodesic slicing (or synchronous gauge) is not always suitable since coordinates can become ill-defined at some point during the evolution of the system, e.g., when shell crossing (or orbit crossing) occurs, as it is expected for collisionless particles.

In order to study cosmological (i.e., expanding/contracting) spacetimes, a convenient prescription to choose $\alpha$ is by applying the so-called Constant Mean Curvature (CMC) slicing condition [57], in which we can set

[TABLE]

where a fiducial Hubble parameter $H\equiv\dot{a}/a$ has been introduced (being $a$ a fiducial scale factor). Then, the lapse function can be found by solving the following constraint arising from (2.2) and the trace of (2.7) in terms of conformal variables,

[TABLE]

Notice that in this scheme $H$ (or $a$) is just a prescribed function for fixing the gauge, and in principle does not represent average (or background) properties of the universe. Nonetheless, we can still fix it by demanding that this satisfy the ‘reference’ (or ‘background’) Friedmann equations

[TABLE]

where $\hat{\rho}_{m}$ and $\rho_{\Lambda}$ are the homogeneous111Note that we use an overhat to denote the homogeneous matter quantities in order to avoid confusion with the overbars used above to denote geometric quantities constructed from the conformal metric $\bar{\gamma}_{ij}$. For $\rho_{\Lambda}$, on the other hand, we omit the overhat for brevity because it does not have an inhomogeneous part. matter and dark energy densities in the reference spacetime, respectively, and $\hat{P}=\hat{S}/3$. The advantage of introducing the fiducial Friedmann equations (2.19) and (2.20) is that we can subtract ‘background’ quantities from the full GR equations, which is more numerically convenient than to solve them directly. A similar idea is also exploited in [58], where a reference FLRW spacetime is conveniently subtracted (but under geodesic slicing), and the application of a fiducial Hubble parameter as part of the CMC slicing condition is later considered in [36]. We remark that, on itself, this ‘background’ subtraction does not constitute an approximation nor a perturbative approach, but rather it is simply a reformulation of the equations using a cosmologically motivated slicing condition.

On the other hand, to fix the remaining gauge freedom let us consider the Minimal Distortion (MD) condition, in which $\beta^{i}$ is chosen such that it minimises the time rate of change of $\bar{\gamma}_{ij}$ during the propagation of spatial coordinates from one hypersurface to the next one. From the traceless part of (2.6) in terms of conformal quantities, we find

[TABLE]

and the MD gauge condition corresponds to demand [57, 59]

[TABLE]

Let us remark here that, contrary to the decomposition (2.14) discussed in the CTT approach, in (2.21) $\bar{A}^{ij}$ has both longitudinal and transverse components even when the MD condition (2.22) is satisfied. Then, using the MD condition (2.22) the momentum constraint (2.12) translates into the following elliptic equation for the shift vector

[TABLE]

where we have also used the CMC condition to simplify the momentum constraint. In the rest of the paper the gauge is fixed by the MD condition and the CMC condition, so that (2.17), (2.18) and (2.23) are satisfied, although the latter is simplified by the approximation introduced in the next subsection. In addition, we use the reference FLRW equations (2.19) and (2.20) to determine $K$ and $\dot{K}$.

2.4 The fully constrained formulation of GR

One of the fingerprints of GR is that gravity is no longer a static field as in Newton’s theory, where the gravitational potential $\Phi_{N}$ is completely ‘slaved’ by the matter distribution, but it hosts two dynamical DOF representing ripples in the spacetime. Due to their faint nature, the existence of these GW has been only recently confirmed by LIGO [60, 61] – about a century after its theoretical prediction [62] – and has opened up a plethora of new possibilities for exploring our Universe. However, in the context of cosmological structure formation and its back-reaction these play a subdominant role, and then one might try to reconstruct the spacetime in absence of GW while retaining all other virtues of GR. Therefore, as a natural extension of Newtonian $N$-body simulations, we propose to use a formulation of GR featuring only elliptic equations as a first step.

In order to achieve this within the CTT approach, the first step would be to construct the initial data by choosing a conformally flat metric $\bar{\gamma}_{ij}=\delta_{ij}$ as well as $\bar{A}^{ij}_{TT}=0$. However, after we fix the gauge assuming the MD condition (2.22), there is no remaining freedom to enforce these conditions for $t>0$ since $h_{ij}\equiv\bar{\gamma}_{ij}-\delta_{ij}$ and $\bar{A}^{ij}_{TT}$ actually satisfy evolution equations. Then, as a possible way to have a fully constrained system of GR equations where the effect of GW in the cosmological dynamics is neglected, following [46, 47] we make the approximations

[TABLE]

This approach follows the similar spirit as in the ‘waveless theories of gravity’ developed originally by Isenberg [44] and later by Wilson & Mathews [45], and its application is supported by previous works from both theoretical and numerical standpoints. In [48] it has been explicitly shown, by using post-Newtonian expansions, that the conformal flatness approximation and the neglect of the TT term (2.24) are accurate even in highly relativistic regimes (further details on this point are discussed the in Appendix of [48]). Furthermore, in [49] the authors discuss a ‘passive’ method to compute the GW emission within this formulation by solving the hyperbolic evolution equations, but without including its back-reaction onto the dynamics of the system, as well as the so-called fully constrained formulation [FCF, 46, 47, 63] in which the latter is properly included. These approaches to go beyond the simple approximation  (2.24) can be potentially implemented in gramses, but this will be left for a future version. Then, due to the conformal flatness approximation intrinsic to this scheme, using that $\partial_{t}\bar{\gamma}^{ij}=0$ in (2.21) it can be shown that the momentum constraint (2.23) reduces to [47]

[TABLE]

Let us now briefly discuss the impact of the MD gauge condition and the conformal flatness approximation on the spacetime metric components. We know that (in 4 dimensions) $g_{\mu\nu}$ has $4$ gauge (redundant) DOF which are fixed by the gauge choices, as well as $6=2+2+2$ independent physical DOF which correspond to the scalar, vector and tensor modes, respectively. The CMC and MD gauge conditions, Eq. (2.17) and Eq. (2.22), place conditions on the metric variables $\alpha$ and $\bar{\gamma}_{ij}$ and fix the 4 gauge DOF. It would then appear that the approximations Eq. (2.24) place further conditions on the metric variables and thus over-constrain the system. This is because $\bar{\gamma}_{ij}$ effectively contains 5 DOF (notice that one scalar dof, $\psi$, has already been factored out of $\gamma_{ij}$ when defining $\bar{\gamma}_{ij}$), all of which are set to zero by demanding that $\bar{\gamma}_{ij}=\delta_{ij}$. However, in the linear perturbation regime, it can be shown that the MD conditions (2.22) are a subset of $\bar{\gamma}_{ij}=\delta_{ij}$, while the latter additionally sets the tensor modes (gravitational waves) to zero. In the nonlinear regime, it is no longer possible to cleanly separate the scalar, vector and tensor modes: in this case, the condition (2.24) is effectively removing all GW content and its back-reactions on spacetime from the system, and in so far as these back-reactions have a negligible effect on the structure formation in cosmology, the approximation (2.24) is good. The conformal flatness approximation has been shown to be accurate for astrophysical systems such as the rotational collapse of cores of (super)massive stars and merger of binary neutron stars [48], where gravity is not exceedingly strong. Under this approximation, there are four physical non-dynamical DOF in the system, the scalar variable $\psi$ plus a scalar and two vector models contained in the shift vector $\beta^{i}$ (this is different from the commonly-used Poisson gauge, e.g., in gevolution), which respectively satisfy the Hamiltonian and momentum constraint equations. Below in this paper we shall demonstrate that GR simulations using this formula are able to predict the scalar and vector modes accurately. A detailed discussion on how the MD gauge singles out the dynamical wave modes of $\bar{\gamma}_{ij}$ can be found in Appendix A of [59].

Under this fully constrained formulation of GR, the gravitational sector equations can be conveniently solved following the next sequence of steps:

Using the CTT decomposition, we solve the momentum constraint (2.15) as an elliptic equation for the longitudinal part of $\bar{A}^{ij}$,

[TABLE]

where $s_{i}\equiv\psi^{6}S_{i}=\gamma^{1/2}S_{i}$ is the conformal momentum density. After this, we construct the traceless part of the extrinsic curvature as

[TABLE]

where in the last step we have neglected $\bar{A}^{ij}_{TT}$ due to the approximation (2.24). 2. 2.

We solve the Hamiltonian constraint (2.11), with (2.19) subtracted, as an elliptic equation for the conformal factor $\psi$

[TABLE]

where $s_{0}\equiv\psi^{6}\rho_{m}=\gamma^{1/2}\rho_{m}$ is the conformal matter density. Here $\bar{\nabla}$ is $\bar{D}$ with $\bar{\gamma}_{ij}=\delta_{ij}$ due to the conformal flatness approximation. 3. 3.

We determine the lapse function from the CMC slicing condition (2.17)

[TABLE]

where we have used (2.19)-(2.20) as well as defined $s\equiv\psi^{6}S=\sqrt{\gamma}S^{i}_{i}$. 4. 4.

Finally, we determine the shift vector from the momentum constraint (2.25). 5. 5.

Then, the current state of the gravitational field is characterised by the spacetime metric

[TABLE] 6. 6.

The metric (2.30) is used to solve the equations of motion (EOM) for particles, and evolve the system.

Notice that in this scheme (2.29) plays the role analogous to the Poisson equation in Newtonian gravity as it determines the $g_{00}$ metric component $\alpha\sim 1+\Phi$. However, (2.29) includes a term that is linearly proportional to $\alpha$ and therefore it formally resembles a finite-range (Yukawa-like) potential, which means that the gravitational force appears ‘screened’ in the CMC gauge. This seems to indicate that it would translate in a suppression of the matter power spectrum at large scales. However, as we will discuss later in Section 6, this ‘screening’ simply represents the effect by the choice of gauge on the simulated matter or velocity power spectrum.

3 The matter sector

In gramses we implement a fully general relativistic $N$-body system for dark matter following the ADM formalism, with which we can describe the matter content in the 3-dimensional hypersurface and its time evolution, rather than dealing directly with 4-dimensional quantities. In GR the equation of motion for collisionless particles is the geodesic equation

[TABLE]

which for the spatial components reads

[TABLE]

where the time-component $u^{0}$ is determined from the normalisation condition $u^{\mu}u_{\mu}=-1$ as

[TABLE]

The relation between $u_{j}$ and the 3-velocity (coordinate velocity) $v^{i}\equiv dx^{i}/dt=u^{i}/u^{0}$ is

[TABLE]

Naturally, the Newtonian limit is recovered with $u_{i}\ll 1$, $\gamma^{ij}\to\delta^{ij}$, $\alpha\to(1+\Phi_{N})$ and $\beta^{i}\to 0$, where $\Phi_{N}$ is the Newtonian gravitational potential, in which case (3.2) reduces to Newton’s second law (in a comoving coordinate system) and $dx^{i}/dt=u^{i}$. Eq. (3.4) shows the nontrivial relationship between the velocities $u_{i}$ and $u^{i}$. To be clear, in the code implementation below we shall only use $u_{i}$ with a lower index, and similarly we always use $\beta^{i}$ with an upper index.

After we evolve particles with (3.2) and (3.4) we can calculate the matter source terms appearing in the equations for the gravitational sector (2.26)-(2.29). For this purpose, it is convenient to introduce a set of conformal source terms defined as

[TABLE]

We will show below that these quantities are analogous to the usual ‘comoving’ ones and correspond to those actually determined numerically in the CIC scheme. For example, the ‘density’ contrast for $s_{0}$, defined as $\delta s_{0}/s_{0}$, corresponds to the local fluctuation of the particle number count rather than the relativistic energy density. This is more convenient in practice because we naturally would like to follow ‘particles’ rather than the ‘density field’ in simulations: a same particle can contribute different energy densities at different places, and the relativistic correction effect can be calculated according to the local metric $\gamma_{ij}$. We will discuss the implication of this on the generation of initial conditions in a forthcoming paper.

4 Numerical implementation

Let us now describe the actual implementation of the fully constrained formulation of GR and the matter evolution equations in gramses. For this purpose, we recast the GR equations using the code units detailed next. Notice that these are derived from the GR equations in physical units rather in geometric units, i.e., accounting for all $G$ and $c$ factors.

4.1 Code units

4.1.1 The gravitational sector

In order to implement the GR equations in gramses we introduce a set of dimensionless quantities that are based on $H^{-1}_{0}$ for measuring time, the box size $L$ for spatial coordinates, the critical density $\rho_{c}={3H^{2}_{0}}/{8\pi G}$ and the fractional matter density $\Omega_{m}\equiv\hat{\rho}_{m,0}/\rho_{c}$ (today’s values satisfying the fiducial FLRW equations (2.19) and (2.20)):

[TABLE]

Notice that, in order to simplify the equations in code units, we have introduced the supercomoving coordinate time $\tilde{t}$ [64], and for $\tilde{s}_{i}$ and $\tilde{s}$ we have introduced a ${c}^{-1}$ factor that’s not present in $\tilde{s}_{0}$. In addition, rather than solving the nonlinear equations for the conformal factor and the lapse function, it is more convenient to reparametrise them by definining new variables $\tilde{\Phi}$ and $\tilde{\Psi}$ as

[TABLE]

where $\tilde{c}$ is the speed of light in code units. This way, both $\tilde{\Phi}$ and $\tilde{\Psi}$ are quantities measuring deviations from their reference FLRW values (but no linearisation on $\tilde{\Phi}$ or $\tilde{\Psi}$ is carried out). Using this scheme, we can write the momentum constraint (2.26) as

[TABLE]

The Hamiltonian constraint (2.28) can be written as

[TABLE]

Next, using the the $2^{\text{nd}}$ Friedmann equation (4.12) and the Hamiltonian constraint (4.7), the CMC condition (2.29) in code units becomes

[TABLE]

where

[TABLE]

Finally, the MD condition (2.25) in code units is

[TABLE]

In addition, the Friedmann equations (2.19) and (2.20) in code units are

[TABLE]

where $\hat{\tilde{\rho}}_{m,0}\equiv\hat{\tilde{S}}_{0}=1$ is the homogeneous (comoving) density field in code units, $\Omega_{\Lambda}\equiv 8\pi G\rho_{\Lambda}/3H^{2}_{0}$, and $\hat{\tilde{S}}=3\hat{\tilde{P}}$ (with $\tilde{s}=\psi^{6}\tilde{S}$).

4.1.2 The matter sector

For writing the EOM for particles in code units, we introduce the following dimensionless quantities for the particle’s 4-velocity and mass:

[TABLE]

Then, the system consisting of Eqs. (3.2)-(3.4) becomes

[TABLE]

where

[TABLE]

is the Lorentz factor in code units.

Finally, for the numerical implementation of the matter source terms in Eqs. (3.5)-(3.7) we need to specify a prescription to calculate them. For this we consider an ensemble of $N$ identical particles of rest mass $m$ treated in a Cloud-in-cell (CIC) scheme. Then, the contributions to each matter source term in Eqs. (3.5)-(3.7) due to a particle at position $(\tilde{x}_{p},\tilde{y}_{p},\tilde{z}_{p})$ can be calculated [35] as

[TABLE]

where

[TABLE]

represent the relative position of the particle inside the $i$-th cell in the $x$, $y$ and $z$ direction, respectively, with $\tilde{f}^{-}_{y}=1-\tilde{f}^{+}_{y}$, $\tilde{f}^{-}_{y}=1-\tilde{f}^{+}_{y}$, $\tilde{f}^{-}_{z}=1-\tilde{f}^{+}_{z}$, and $\Delta\tilde{V}_{ijk}=\Delta\tilde{x}\Delta\tilde{y}\Delta\tilde{z}$ is the cell volume.

In the remainder of this paper we will only deal with quantities in code units and the tilde will be dropped to avoid cluttered notation.

4.2 Code structure

Let us now describe the logic flow of the global solution scheme implemented in gramses. We recall that in this code time evolution is only due to particles, and these are evolved in a leapfrog scheme detailed in Appendix A. At a given timestep $n$, using the positions and velocities $\{x^{n},u^{n-1/2}\}$ for particles and the GR fields values $\{\Psi^{n},\Phi^{n},{\beta^{i}}^{n}\}$ the code takes the following main steps:

Optional (gr_newtonian): Solve Newtonian gravity to get $\Phi^{n+1}_{N}$. 2. 2.

Optional (gr_newtonian): Temporarily synchronise velocities with $\Phi^{n+1}_{N}$: $u^{n-1/2}\longmapsto{u^{n}_{N}}$. 3. 3.

Calculate the source terms of the GR equations using the currently-available velocity and GR potential fields. 4. 4.

Optional (gr_newtonian): Revert the temporary synchronisation done in step 5: $u^{n}_{N}\longmapsto{u^{n-1/2}}$. 5. 5.

Solve the ten GR field equations: $\{\Psi^{n},\Phi^{n},{\beta^{i}}^{n}\}\longmapsto{\{\Psi^{n+1},\Phi^{n+1},{\beta^{i}}^{n+1}\}}$. 6. 6.

Synchronise velocities (last ‘Kick’ from previous timestep: $u^{n-1/2}\longmapsto{u^{n}}$). 7. 7.

Update velocities (first ‘Kick’ of the current timestep: $u^{n}\longmapsto{u^{n+1/2}}$). 8. 8.

Update positions using $u^{n+1/2}$ (‘Drift’: $x^{n}\longmapsto{x^{n+1}}$). 9. 9.

Go to step 1 with the replacement $n+1\longmapsto n$ and repeat the process for the next timestep.

We address these points in detail in the next subsections.

4.3 Multigrid solver for the gravitational sector

Let us now discuss the main points of the multigrid [65] implementation of the GR equations for the gravitational sector. For this purpose, it is convenient to first split the equations into the Poisson-like ones and the non-Poisson-like ones.

4.3.1 Poisson-like equations

In the solution scheme for the gravitational sector we have eight linear Poisson-type equations arising from the vector Laplacians (4.6) and (4.10) respectively, i.e.,

[TABLE]

where ${W}_{i}\equiv{V}_{i}+{\partial}_{i}{U}$ and ${\beta}^{i}\equiv{B}^{i}+{\partial}^{i}{b}$ have been used to cast these equations in the form of a standard Poisson equation.

In order to solve the above equations numerically, we need to represent them in a discrete form. If we consider a uniform grid with cubic cell size $h$, these equations are formally equivalent to

[TABLE]

where $\mathcal{L}^{h}$ is the Laplacian operator and $f^{h}$ a source term. The former is discretised using the standard second order formula

[TABLE]

where $\varphi_{l,m,n}$ is the value of the field $\varphi$ in the grid cell with index $(l,m,n)$.

The actual discretisation method for $f^{h}$ depends on the particular source term under consideration. In our implementation, the source $f^{h}$ for Eq. (4.21) is calculated using the CIC scheme, while the sources for the rest of these equations correspond to divergences which are calculated using finite differences with a 3-point stencil, e.g.,

[TABLE]

where $\bar{A}^{\prime lm}\equiv\left(1+\frac{{\Phi}}{a^{2}{c}^{2}}\right)\left(1-\frac{{\Psi}}{2a^{2}{c}^{2}}\right)^{-6}{\bar{A}}^{lm}$. The extrinsic curvature term itself is calculated as

[TABLE]

Notice that for the non-diagonal components of $\bar{A}_{ij}$ we need to calculate cross-derivatives of $U$ that depend on the diagonal neighbours of the central cell, for which we use a second-order-accuracy formula. Then, for a cross-derivative in the $xy$ plane (contributing to $\bar{A}_{xy}$) we use the discrete expression

[TABLE]

and equivalent ones for cross-derivatives in the $xz$ and $yz$ planes.

As a side point, we note that the left-hand sides of Eqs. (4.21)-(4.24) have the ‘cyclic’ property, e.g., if one sums up the values of $\nabla^{2}U$ from all cells with the same $y,z$ (but different $x$) coordinates, the result is zero, which is guaranteed by the nature of finite difference. The same periodic property is also held by the right-hand sides of Eqs. (4.22)-(4.24) by default, which means that these equations are numerically ‘self-consistent’. For Eq. (4.21), however, its right-hand side does not automatically satisfy the ‘cyclic’ property, which could potentially lead to the situation where the two sides of the equation cannot be numerically identical. To prevent this issue, for Eq. (4.21) we have redefined its source term $f^{h}$ as $f^{h}-\langle{f}^{h}\rangle$, where $\langle{f}^{h}\rangle$ is the mean value of $f^{h}$ on the whole domain grid (the finest grid with uniform resolution that covers the whole simulation box).

The discrete version of Eq. (4.25) is then solved using a Gauss-Seidel (pseudo-time) relaxation method, which after a given number of iterations converges to a solution $\hat{\varphi}^{h}$. Then, the residual at the fine level is defined as

[TABLE]

In order to improve the accuracy of the solution found using relaxation on the fine grid $h$ (and hence reduce $\epsilon^{h}$), instead of keeping iterating on the same grid we can accelerate the process by moving to a coarser grid with cell size $H=2h$, which is then associated with the multigrid level $\ell-1$ (with $\ell$ representing the finest grid). Since in Eq. (4.25) the operator $\mathcal{L}^{h}(\varphi^{h})$ is linear, the coarsified equation is

[TABLE]

where $\mathcal{R}$ is the restriction operator and $\delta\varphi^{H}\equiv\varphi^{H}-\mathcal{R}\hat{\varphi}^{h}$. The coarse-level equation (4.32) is then solved using Gauss-Seidel relaxation and the solution $\hat{\delta\varphi^{H}}$ found after a given number of iterations is used to correct the fine solution $\hat{\varphi}^{h}$ by means of a prolongation step

[TABLE]

where $\mathcal{P}$ is the prolongation operator that maps the information from coarse to fine grids. In practice, we use all coarse levels from $(\ell-1)$ down to $1$ (the level at which the whole simulation box is divided into 8 cells) for our multigrid implementation, and arrange them in a ‘V-cycle’: instead of correcting the level-$\ell$ solution immediately after solving the equation on level-$(\ell-1)$, the process of solving coarser-level equations goes all the way down to level $1$, and the correction process then goes all the way up to level $\ell$. This has greatly improved the speed at which a solution with acceptable accuracy is achieved.

4.3.2 Non-Poisson-like equations

Let us now discuss the implementation of the Hamiltonian constraint and the CMC condition, which correspond to the nonlinear PDEs solved in the code.

4.3.2.1   Hamiltonian constraint

In the multigrid language, the Hamiltonian constraint (4.7) at the fine level can be formally written as

[TABLE]

where

[TABLE]

Notice that by keeping the density term ${s}_{0}$ in the source $f^{h}$ (rather than absorbing it in the definition of $\mathcal{L}^{h}$) we avoid restricting $s_{0}$ from fine to coarse level, and we only need to restrict the operator coefficient $\bar{A}_{ij}\bar{A}^{ij}\subset\mathcal{L}^{h}(\Psi^{h})$. However, as we discussed above for the case of Poisson-like equations, since our solution needs to satisfy periodic boundary condition, we need to regularise this equation to ensure that the numerical self-consistency condition $\langle\mathcal{L}^{h}\left(\Psi^{h}\right)\rangle=0$ is satisfied222Note that in this case $\langle{f}^{h}\rangle$ is guaranteed to be zero by the periodic boundary condition of $s_{0}$ and the subtraction of $1$ in Eq. (4.35).. Then, dividing both sides of Eq. (4.34) by $\left(1-\frac{\Psi^{h}}{2a^{2}{c}^{2}}\right)$ and taking the mean on both sides, we identify the regularisation term

[TABLE]

and the regularised version of (4.7) solved in the code reads

[TABLE]

which corresponds to a redefinition of our differential operator at the fine level. Then, the residual at the fine level is defined as

[TABLE]

Next, at level $\ell-1$ the coarsified version of this nonlinear equation has the generic form

[TABLE]

where $\hat{\Psi}^{h}$ is the solution obtained at the fine level, $\mathcal{R}$ is the restriction operator. In principle we could also regularise the equation at the coarse level, but in practice this is not needed if the fine level is already regularised.

4.3.2.2   The CMC condition

In order to optimally solve the CMC condition (4.8) in the multigrid scheme we introduce the combination

[TABLE]

which allows to avoid the usage of additional arrays for storing the restricted field ${\Psi}$ while solving for ${\Phi}$ on coarse levels. Then, in terms of $\xi$ the CMC equation (4.8) at the fine level can be written formally as

[TABLE]

where

[TABLE]

and $Q$ is given by (4.9). After subtracting the mean from both sides, the regularised equation becomes

[TABLE]

This means that for solving this equation at the coarse level it is sufficient to restrict the operator coefficient $\left(1-\frac{{\Psi}}{2a^{2}{c}^{2}}\right)^{-1}\frac{Q}{a^{2}{c}^{2}}\subset\mathcal{L}^{h}(\xi^{h})$. Further, similar to the Poisson-like equations, the regularisation term in (4.43) is absorbed in the source term, so there is no redefinition of the differential operator but of $f^{h}$. Then, the residual at the fine level in this case is

[TABLE]

On the other hand, at $\ell-1$ the coarsified equation is

[TABLE]

where $\hat{\xi}^{h}$ is the solution obtained at the fine level. As in the case of the Hamiltonian constraint, there is no need to regularise this equation at the coarse level.

4.4 Particles evolution

After having reconstructed the spacetime by computing the metric components $(\gamma_{ij},\alpha,\beta^{i})$ in a fully nonlinear fashion, we can then solve the EOM for particles. Drawing the analogy with the Newtonian case, the geodesic equation can be rewritten effectively as

[TABLE]

where

[TABLE]

In the Newtonian limit the ‘force term’ given by (4.48) reduces to $F_{i}\to-\partial_{i}\Phi$, but in the GR case this depends not only on gradients of the various gravitational fields (and on the fields themselves), but also on $u_{i}$. In practice, this means that we cannot compute all the contributions to (4.48) in the same way as the default ramses code (or any standard Newtonian code) does. Therefore, in gramses we divide the Kick step (4.46) into a sequence of 5 substeps, each one updating the particle velocity using one force contribution on the right hand side of (4.48), which is decomposed as $F_{i}=\sum_{j}c_{j}(f_{j})_{i}$ with

[TABLE]

Notice that in doing this we are using the (partially) updated velocity during each substep. Of these five substeps, (4.52) is the dominant one in most realistic situations as it corresponds to the standard Newtonian force in the Newtonian limit.

As an attempt to preserve the Stormer-Verlet scheme as best as we can, the last Kick in the KDK scheme is done following (4.50) to (4.52), i.e., with the largest contribution to $F_{i}$ included at last (see Appendix A). In contrast, for the first Kick step (before updating the particles’ positions) we use the reverse order, i.e., with the largest contribution included first. Again, this is because that, according to the Stormer-Verlet scheme, during the first Kick the ‘force’ should be evaluated using $u^{n+1/2}$, which implies an implicit equation for the latter (since this is the very velocity that we want to update to). Hence, as an approximation we use the synchronised velocity $u^{n}$ and the largest contribution during the first substep of the Kick, which then yields a velocity $u^{n+1/2}_{*}$ that will be the close to $u^{n+1/2}$, and this is then used in the next substeps to calculate the corrections. Finally, the positions are updated in the Drift step (4.47) in a single calculation once the velocities have been fully updated by the previous prescription.

4.5 Calculation of matter sources

A key difference in the calculation of the general relativistic matter source terms in Eqs. (3.5)-(3.7) with respect to the Newtonian case is that the former depend not only on $u_{i}$ but also on $\Psi$ through the Lorentz factor Eq. (4.49). For calculating these we use the already-known values $\Psi^{n}$, with which the GR equations are solved to get the updated metric components ${\{\Psi^{n+1},\Phi^{n+1},{\beta^{i}}^{n+1}\}}$. This is equivalent to the numerical implementation in the gevolution code, where the geometric corrections in the energy-momentum tensor at a given timestep are calculated using the values from the previous timestep [23]. More explicitly, the CIC quantities depend on these as

[TABLE]

Note that as the second ‘Kick’ step – which takes particle velocities from $u^{n+1/2}$ to $u^{n+1}$ and therefore ‘synchronises’ particle velocities – is done after we solve the GR equations, at the time when the code calculates the matter source terms for the GR equations what are available are the fully updated positions at timestep $(n+1)$, $x^{n+1}$, and the partially-updated velocities $u^{n+1/2}$ (which are still delayed by half a timestep). This issue is not present in the Newtonian case since the gravitational potential (and hence the force) is independent of the particle velocities, and we expect that the use of $u^{n+1/2}$ instead of $u^{n+1}$ in the matter source calculation should be a good approximation given the generally small timesteps for simulations.

Nevertheless, gramses has an option to remedy the fact that we only have $u^{n+1/2}$ to calculate the GR source terms, by using a temporary ‘Newtonian’ synchronisation from solving the standard Newtonian gravity. With this option switched on, the code uses the Newtonian gravitational potential $\Phi_{N}$ to temporarily update the velocities $u^{n+1/2}\longmapsto{u^{n+1}_{N}}$, which are then used to calculate the source terms as a better approximation than using $u^{n+1/2}$ directly in Eqs. (4.53)-(4.55). Then, after the GR equations are solved, we can exactly revert the velocities back to $u^{n+1/2}$, before carrying the Kick step normally (see appendix A).

5 Code tests

We have performed several code tests for gramses, particularly aimed to test the implementations of the linear and nonlinear solvers of the ten GR potentials, the subroutines that calculate new GR quantities, as well as the geodesic solver.

5.1 Static tests

Let us first discuss the tests that require no cosmological evolution, for which we set the fiducial scale factor to $a=1$. The results shown in this section correspond to simulations with a box size $L=256$ Mpc$/h$ and $N_{p}=256^{3}$ particles, and they are used to check the subroutines in gramses to calculate the matter and geometric source terms, and to solve the relevant PDEs for the gravitational sector.

5.1.1 Matter and geometric source terms

While in Newtonian $N$-body simulations the matter density field is the only source term feeding the gravitational potential, in gramses the picture is more intricate: the Poisson-like equations (4.21)-(4.24) feature the momentum density $s_{i}$, the divergence of GR potentials and also that of the traceless part of the extrinsic curvature, i.e., $\partial_{j}\bar{A}^{ij}$. Furthermore, for the non-Poisson-like equations, (4.7) and (4.8), we also need to calculate terms such as the contraction $\bar{A}_{ij}\bar{A}^{ij}$, the density $s_{0}$ and the trace $s=\gamma^{ij}s_{ij}$. Clearly, the calculation of the last two quantities is more complicated than that of $s_{i}$ since they depend nonlinearly on $u_{i}$ through the Lorentz factor $W$ as can be seen from (4.17)-(4.19). Since in the CIC scheme particles may contribute to different cells depending on their positions, in order to assess the calculation of these quantities, in the tests we fix the particle positions and velocity values by hand, as well as $\Psi$, so we can compare against analytical expressions for the matter source terms.

Figure 1 shows the results for the CIC calculation of $s_{0}$ and two $s_{i}$ components when using a uniform particle distribution, a velocity field $u_{i}=0.3c\sin{2\pi{x_{i}}}$ and $\Psi=c^{2}\sin{2\pi{x}}\sin{2\pi{y}}\sin{2\pi{z}}$, where $c$ is the speed of light in code units. We find good agreement since the difference with respect to their analytical counterparts remains below $10^{-5}$, and the structures observed for this in the bottom panels of Figure 1 stem from the shape of the functions on the top panels.

Regarding the geometric source terms, in Figure 2 we show the calculation of $\partial_{i}\bar{A}^{yi}$ and $\bar{A}_{ij}\bar{A}^{ij}$ which involve the discrete formulae, (4.28), (4.29) and (4.30). In order to compare with analytical expressions for this quantities we use $V_{i}=U=\sin{2\pi{x}}\sin{2\pi{y}}\sin{2\pi{z}}$ as an input for the relevant code subroutines. The agreement found from Fig 2 is roughly one order of magnitude worse than for the matter source terms in Fig 1 but is still nonetheless stays around $10^{-4}$ for both $\partial_{i}\bar{A}^{yi}$ and $\bar{A}_{ij}\bar{A}^{ij}$. Naturally, the accuracy in the calculation of these quantities is expected to be different as they depend on how fine is the mesh used to perform the finite differences involved in their calculations, and the agreement shown in Fig. 2 is consistent with the grid used for these tests. Similar to Fig. 1, the structure of the difference observed in the bottom panels of Fig. 2 follows from the shape of the testing functions.

We have done similar tests as in Figures 1 and 2, using different choices of particle velocities, and different functional forms for $\Psi$, $U$ and $V_{i}$. For all these tests we found similar agreement between the code results and the analytical predictions. For simplicity, the extra tests are not shown here.

5.1.2 Homogeneous density field

We next show some test cases for the linear and nonlinear equation solvers implemented in gramses. The simplest possibility of such tests for the gravitational sector is to reconstruct the spacetime due to a homogeneous density configuration which corresponds to an FRLW (Einstein-de Sitter) solution. In our test, we set $s_{0}=1$ and $s_{i}=s=0$.

Recall that the CMC gauge condition is fixed with the aid of the fiducial Friedmann equations (4.11) and (4.12) and hence the Hamiltonian constraint (4.7) and CMC condition (4.8) have the ‘background’ solutions $\Psi=\Phi=0$. Likewise, it is also straightforward to find that the linear GR equations (4.21)-(4.24) are trivially solved, i.e., $V_{i}=U=B^{i}=b=0$, in this case. In Figure 3 we show the numerical solutions for $\Psi$ and $\Phi$ that are obtained by the nonlinear Gauss-Seidel relaxation solver333While technically speaking the CMC gauge condition (4.41) is a linear equation, in gramses it is solved using the same nonlinear relaxation solver as in the Halmiltonian constraint (4.34). implemented in gramses after performing the relaxation starting from two sets of random initial guesses (green triangles and red circles). As the figure shows, there is a very good agreement in the numerical solutions to both GR potentials regardless of the initial guess from where the relaxation is started.

We have tested the solutions to the other GR potentials finding similarly good agreements. To save space, here we only show the results for the two most complicated equations.

5.1.3 Point-like, sinusoidal and spherically symmetric sources

In order to add some nontrivial features to the test (as opposed to the case with a homogeneous matter field), we have also tested gramses’s relaxation solvers for the Poisson-type and non-Poisson-type equations individually by using various configurations for the source terms. For instance, for the Poisson-type equations (4.22)-(4.24) we know that the solution for a point-like source located at $x_{p}$ is simply given by Green’s function

[TABLE]

We test this point-like source scenario by initialising the value of the source in a single cell of the simulation box (at position $x_{p}$) to a constant: in the case of (4.21) this is equivalent to having a single particle in the center of a cell so that it contributes only to that same cell in the CIC scheme, while for the other linear equations this corresponds to having a non-zero value for their geometric source terms in a single cell of the domain grid. The left panel of Figure 4 shows the numerical solution to (4.22) in such case and its comparison with the exact analytical solution, (5.1). We note that the numerical solution deviates from the exact one towards the center of the simulation box as well as towards the box boundary. Both discrepancies can be understood in terms of the discrete nature of the numerical simulation: on the one hand, when approaching the source position resolution effects become important and the isotropy respected by the exact solution is broken. On the other hand, towards the box boundary the effect of the finite box size and the periodic boundary conditions imposed on the numerical solution causes deviations from the exact solution (5.1) which decays at infinity. Analogue results for the case of a Schwarzschild spacetime in isotropic coordinates are found in [23] and for point-like mass tests of modified gravity codes, e.g., [66].

The non-Poisson-type equations (4.7) and (4.8) do not have exact analytical solutions for configurations such as point mass to compare against. As an alternative, we have used an inverse approach in which we choose some functional form for $\Psi$ or $\Phi$ by hand and solve (4.7) and (4.8) for some of the source terms under certain simplifications, which can be then used as input to gramses, to check the agreement of its resulting numerical solutions to $\Psi,\Phi$ with the above choices.

In Figure 4 we show the results of such procedure using $\Psi$ with two different functional forms: $\Psi(x)=10^{-2}\sin{2\pi x}$ and $\Psi(r)=\exp{[-(r-0.5)^{2}/10^{2}]}$ which are then used in Eq. (4.7) to analytically solve for $s_{0}=s_{0}(x)$ and $\bar{A}_{ij}\bar{A}^{ij}=\bar{A}_{ij}\bar{A}^{ij}(r)$, respectively, assuming all other source terms to be zero (here $r$ is measured from the center of the box). These are used as inputs for the gramses nonlinear solver. For the sinusoidal test (middle panel of Figure 4) we find that the numerical solution deviates less than $10^{-5}$ from the exact solution (we also include the initial random guess for the relaxation in that plot), while deviations for the spherically symmetric test (right panel of Figure 4) are larger but still better than those for the point-like test. Similar to the latter case, the numerical solution in the spherically symmetric test also suffers from the effects of periodic boundary conditions which depends on the rate at which the tail of the exponential function decays. Like before, we have tried more test settings and carried out similar tests for other equations as well. As all the tests result in similar agreement between numerical and exact solutions, we shall not show those tests here to save space.

Note that, in gramses, the numerical errors of these solutions are controlled by using a measure of the magnitude of the residual $\epsilon$ defined in Eq. (4.31) for the Poisson-like equations and defined in Eq. (4.38) for the non-Poisson-like equations, which is fundamentally different from standard numerical relativity codes in which the full evolution equations for the gravitational sector are solved. In the latter, the accuracy of the numerical solutions is monitored by substituting them into the momentum and Hamiltonian constraint equations, thereby quantifying the measure of ‘constraint violation’. However, since in the GR formulation implemented in gramses the constraint equations are actually used as part of the the solution scheme, the numerical errors are entirely controlled by the residual $\epsilon$. Some threshold value for $\epsilon$ is specified in order to achieve the desired accuracy through the multigrid Gauss-Seidel relaxation solvers, which can be made smaller to increase the solution accuracy at the cost of using more iterations. For example, for our chosen threshold value for $\epsilon$, the $\mathcal{O}(10^{-5})$ error in the numerical solution for $\Psi$ shown in Fig. 4 (middle panel) is comparable to the level of Hamiltonian constraint violations typically found in other numerical relativity codes [26, 37].

5.2 Dynamical tests

To test the time integration part of the code, which is determined by the evolution of particles along geodesics as discussed in Section 4.4, we consider two typical scenarios, including an FLRW robustness test and the evolution of a linear density perturbation. For these tests we use a box size $L=256$ Mpc$/h$ with $N_{p}=128^{3}$ particles. The discussion on the cosmological runs using higher-resolution specifications is presented in Section 6.

5.2.1 FLRW robustness test

We now present the test for the robustness of the evolution of a FLRW universe. For this, we initially set a completely homogeneous and static particle distribution at $z=99$ as determined from the reference Friedmann equations (2.19) and (2.20). This scenario represents the simplest test of the geodesic solver as ideally particles should remain at fixed positions. In order to test the robustness of the numerical FLRW evolution we introduce tiny initial inhomogeneities by using a space-dependent initial random guess of order $\mathcal{O}(10^{-8})$ for the solutions of $\{V_{i},U,B^{i},b,\Psi,\Phi\}$, with which the multigrid relaxation solvers of the respective PDEs converge to solutions that are not exactly zero. After this, the geodesic equations (4.46) and (4.47) are solved using these solutions and, because the solutions to the potentials are not exactly homogeneous due to the numerical errors, non-zero gradients of them and therefore forces result and the particles are driven away from their original positions, which effectively introduces an initial deviation from homogeneity in the density field of $\mathcal{O}(10^{-14})$. Then, in the new timestep the CIC source terms use the updated particle information and the non-exact solutions from the previous timestep are used as initial guesses for the gravitational solver. In this way, the cycle is repeated in such a way that the numerical error becomes accumulated over time. In Fig. 5 we show the resulting matter density field obtained at $z=1$ along three different directions across the simulation box. The error found in Fig. 5 is of around the same order of magnitude as the typical values obtained for the FLRW robustness test in other cosmological GR codes [26].

As a sanity check, we have also carried out the same test but by using zeros as the initial guesses for the relaxation solvers for the different gravitational potentials. In this case, the initial guess (zero) is the exact solution for all the potentials, and no force is exerted on the particles, which in turn stay at their initial positions during the whole simulation (exact FLRW solution). We have also checked that the long wavelength feature of $s_{0}$ along each direction shown in Fig. 5 is an artefact of the Gauss-Seidel relaxation solver. This type of solver efficiently eliminates the short wavelength modes of the error in the initial guess, while long wavelength modes require more iterations and can survive the convergence criteria. Then, these long wavelength modes in the gravitational potentials get imprinted in the matter sector throughout the evolution of the system, although they do not affect the good agreement with the exact solution as shown in Fig. 5.

5.2.2 Linear density perturbation

As a second test for the geodesic solver, here we present the evolution of a perturbed density field that is initialised at $z_{\rm ini}=99$ to $\delta s_{0}=A\sin{2\pi nx}$, with $A=10^{-2}$ and $n=4$. This can be considered as a single mode of a density field. The linear velocity field is inferred from the Zel’dovich approximation [67],

[TABLE]

This density perturbation is then evolved by the code on a $\Lambda$CDM FLRW background (the details of this $\Lambda$CDM model are the same as for the cosmological tests below and, as they are not important for the test here, we postpone the mentioning of them to the next section).

In this test, the density perturbation is expected to scale at linear order as $\propto a$, and the momentum (or velocity) field as $\propto a^{-3/2}$. The density field profile at different output redshifts (circles) and the corresponding analytical predictions (solid line) is shown in the left panel of Fig. 6. We find good agreement at all considered redshifts, though it becomes worse towards low $z$ as expected. At $z=9$, the relative deviation with respect to the linear theory solution is $4.0\%$ at the maxima.

For the momentum density (or equivalently velocity) field, we show the results for the $V_{1}$ potential, which is related to it through the momentum constraint (4.21), in the right panel of panel of Fig. 6. We find that in this case the solution deviates slightly more rapidly from the analytical prediction toward late times compared to the density case, but still follows expectation very well. A somewhat similar situation is also encountered in other numerical relativity codes, where the momentum constraint violations are usually more severe than those found for the Hamiltonian constraint [26, 14, 37]. In this case, at $z=9$ we find a $9.6\%$ of relative deviation with respect to the linear theory solution at the maxima.

6 Cosmological simulations

In this section we present some of the first results of gramses cosmological simulations in a $\Lambda$CDM universe using three different setups. We use a simulation with comoving box size $L=512h^{-1}$Mpc and $N_{p}=512^{3}$ dark matter particles, with the AMR option switched off, to generate the different maps shown in Fig. 7 for visualisation. For the analyses of the matter and velocity divergence power spectra we use big box size simulation with $L=4h^{-1}$Gpc and $N_{p}=1024^{3}$, with AMR switched on, and we compare against a Newtonian simulation run using the default ramses code with identical specifications. Finally, we use a high-resolution simulation with $L=256h^{-1}$Mpc and $N_{p}=512^{3}$ for the power spectra of the scalar and vector modes of the shift vector (again with AMR switched off). The initial conditions (IC) for both the GR and the Newtonian simulations are generated at $z_{\rm ini}=49$ using the same random number sequence as seed in order to suppress the effect of realisation scatter in our results. For all the simulations, the cosmological parameters used for the $\Lambda$CDM model are $\{\Omega_{\Lambda}=0.6928,\Omega_{m}=0.3072,\Omega_{K}=0,h=0.68\}$. Since gramses works in a different gauge than standard $N$-body simulations, the generation of IC is nontrivial; here we simply mention that the ICs were generated using a new technique specifically developed for gramses, and will defer a detailed description of it to a forthcoming paper. Throughout the analysis, the redshift values quoted for the relativistic simulation correspond to those determined from the reference cosmology Eq. (2.19) and Eq. (2.20) (fixed by the CMC slicing), which coincides with the background of the Newtonian counterpart.

Fig. 7 is a visual illustration of the maps of three GR matter source terms, $s_{0}$ (top row), $\theta=\nabla\cdot u$ (middle row) and $s=\gamma^{ij}s_{ij}$ (bottom row), in a slice of constant $z$ coordinate randomly selected from the smaller gramses simulation. The three columns correspond to three different redshifts, $z=9$, $4$, $1$ (in the reference cosmology) from left to right. One can see that as time advances and structure formation progresses, finer features start to appear in all three quantities, and their amplitudes also increase; these results are as expected. In addition, the main features in all three quantities have good correspondences, with high-density regions having larger values of $s_{0},\theta$ and $s$, and vice versa. Note the amplitude of $s$, which has a maximum of order $10^{-4}$ at $z=1$ that is much smaller than the perturbation in $s_{0}$; according to Eq. (4.19), this indicates that $(u/c)^{2}=(\tilde{u}/\tilde{c})^{2}\lesssim\mathcal{O}(10^{-4})$ because in code units $\tilde{m}=1$.

Figure 8 presents maps from the same slice as Figure 5, but for various GR quantities or their scalar combinations. From top to bottom the rows show respectively $\alpha$, $\psi$, $|\beta|$ (the amplitude of $\beta_{i}$) and $\bar{A}_{ij}\bar{A}^{ij}$. A logarithmic colour scale is used for $\alpha$, $\psi$ and $\bar{A}_{ij}\bar{A}^{ij}$, while a linear colour scale is used for $|\beta|$, in order to make the features in the maps clearer. To avoid the plot getting too cluttered, we have not shown the colour bars. In all cases, the same colour scheme in Python (the ‘jet’ scheme) is used, but is ‘regulated’ such that the reddest (bluest) colour represents the maximum (minimum) of the field values in all pixels of a given map. This is done deliberately: had we used a fixed maximum and minimum value for a quantity at all redshifts, the $z=9$ maps would appear uniform and almost completely erased of details. In contrast, the ‘regulation’ not only makes the features at high redshifts clearer, but also demonstrates that, apart from amplitudes, the qualitative patterns of these features barely evolve in time: indeed, there is hardly any visual difference between the left and the middle columns. This is because the GR potentials usually satisfy a Poisson-like equation444The deviations from Poisson, such as the inclusion of nonlinear source terms or, like the case of $\bar{\alpha}$, the addition of an extra Yukawa term (i.e., the term linear in $\xi$ in Eq. (4.42)), do not affect the qualitative discussion here., which in Fourier space takes the form $k^{2}\times{\rm field}={\rm source}$, so that the field value scales as $k^{-2}$ and so is dominated by the large-scale modes (the reason why the maps in Figure 8 generally lack the much finer details present in the maps of Figure 7). These modes remain linear over time, so that their amplitudes grow quickly while the qualitative patterns change much more slowly. The maps for $\alpha$ and $\psi$ are almost identical, with their colours flipped, which is what is expected from Eqs. (2.28) and (2.29).

6.1 Matter and velocity power spectra

Having displayed some visualisation of the simulation outputs from gramses in the previous subsection, we next show a few more quantitative results to illustrate that the code works properly. Figure 9 shows the matter power spectra from the $4h^{-1}$Gpc simulations mentioned above, for six outputting redshifts from $z=49$ (initial redshift, top left panel) to $z=1$ (bottom right panel). Within each panel, the top subpanel shows the absolute power spectra and the bottom subpanel shows the relative difference between gramses (CMC-MD-gauge spectrum for the particle number count perturbations, $\delta s_{0}/s_{0}$) and ramses (synchronous-gauge spectrum for the energy density perturbations, $\delta\rho/\rho$). For comparison, we also show the linear theory predictions for these two spectra using a modified version of the camb code [68]: these are respectively the blue and black solid lines in the top subpanels, while the solid lines in the bottom subpanels denote their relative differences. The magenta symbols in the top subpanels show the gramses results as measured using the power spectrum code powmes [69]: while they agree with linear theory rather well on intermediate and small scales, for the largest scale probed at $z=49$ powmes fails to recover the upturn predicted by linear theory (the blue solid line). The orange and red symbols, on the other hand, are the power spectra for the gramses and ramses simulations measured using the dtfe code555The dtfe code tessellates the simulation volume following the Delaunay triangulation scheme, where the 3D space is decomposed into tetrahedrons whose vertices are simulation particles. The (density or velocity divergence) field value in each tetrahedron is determined by the corresponding particle quantities (mass and velocity) on its four vertices. The values in the tetrahedrons are then interpolated to a regular grid to give the field values on the latter, from which the corresponding power spectra can be measured using normal fast Fourier transform. This is particularly useful for the velocity field, since it can help to avoid the numerical problem of directly interpolating to a regular grid, which often leads to certain grid cells having zero velocity. See [70] for more details. [70], which does capture the upturn. Therefore, in the bottom panels we show the relative differences between the spectra measured using dtfe (red symbols). This comparison poses an interesting question regarding the applicabilities of the different methodologies to calculate the power spectrum when it deviates significantly from the usual behaviour observed in the synchronous gauge666Note that at the largest scale probed by powmes the CMC-MD-gauge power spectrum is nearly two orders of magnitude larger than that in the synchronous gauge. (e.g., on large scales due to the gauge difference). Furthermore, the prediction of actual observable quantities in GR is not straightforward and requires the application of ray tracing algorithms acting either in real time (on-the-fly) or in a post-processing step. The former has been implemented for ramses in [71] and could serve as a starting point to implement a general relativistic version in gramses in the future.

From Figure 9 we see that on large scales the simulation result agrees with linear-theory prediction very well for all redshifts shown. Note that the relative difference between the two power spectra at the largest scale probed by dtfe starts from $\sim 1500\%$ at $z=49$ and decreases to $\lesssim 20\%$ at $z=1$, which is a very large range of change that is properly reproduced by gramses. Towards low redshift a small discrepancy from linear theory appears: this is partly because of the smaller differences between the power spectra in the CMC-MD and synchronous gauges, and partly because of the coarse time resolution in our simulation (which has fewer than 50 time steps between $z=49$ and $z=1$). Due to the low resolution of the simulation we shall not focus on the results at small scales.

Figure 10 is similar to Fig. 9, but shows the power spectrum of the velocity divergence ($\theta$ defined above). The blue and black solid lines in the upper subpanels are respectively the velocity convergence power spectra in the CMC and Newtonian gauge calculated by a modified version of the camb code, and the solid lines in the lower subpanels are their relative differences. Here all the simulation spectra have been measured by the dtfe code. A similar level of agreement between the linear theory predictions and simulations to what is shown in Figure 9 can be found here, where at the largest scale probed by dtfe the relative difference between the gramses and ramses results changes strongly from $\sim 1000\%$ at $z=49$ to $\sim 20\%$ at $z=1$. While the matter power spectra from the gramses GR simulations are consistently higher than those from the ramses Newtonian simulations, the velocity power spectra show the opposite trend (this is in agreement with the findings of [72], whose uniform-expansion gauge corresponds to our CMC gauge) – this is partially why the difference between the matter power spectra in the two simulations decreases over time.

Let us remark here that neither the matter power spectrum $P_{\delta\delta}(k)$ nor the velocity divergence power spectrum $P_{\theta\theta}$ are gauge invariant quantities and the differences on large scales seen in Figures 9 and 10 are therefore not physical effects. However, being able to reproduce the expected gauge effects in our simulations is a useful test of the numerical implementation itself. We have stated briefly above that the $P_{\delta\delta}$ measured from the gramses simulations are the power spectra for the particle number count perturbations in the CMC-MD gauge, while the $P_{\delta\delta}$ measured from the ramses simulations are for the energy density field in the synchronous gauge. We have also compared $P_{\theta\theta}$ from ramses with the Newtonian gauge results from linear theory. These issues are indeed intricately related to the fact that the initial conditions in our GR simulations have to be generated in a way compatible with our gauge choice (see, e.g., [73, 74, 75]), and as mentioned above these topics will be discussed in depth in a forthcoming paper [76].

6.2 The shift vector power spectrum

Let us now discuss some results on the shift vector power spectrum from gramses. This is a particularly interesting quantity since it is related to ‘frame-dragging’, a GR effect which has been measured in the Solar System [77]. Contrary to the longitudinal gauge commonly adopted for late-time cosmology, in the MD gauge the shift vector has both scalar and vector components, i.e. $\beta^{i}=\beta^{i}_{s}+\beta^{i}_{V}$. In perturbation theory, $\beta^{i}_{V}$ appears as a gauge-invariant, second order quantity (at lowest order) that is only sourced by the product of the first order scalar perturbations $\delta$ and $v$ [78, 79] in the perfect fluid case. However, in $N$-body simulations the shell-crossing of dark matter particles induces velocity dispersion and sources the rotational velocity. Expanding Eq. (4.21) and (4.23) up to second order in perturbation theory, it can be shown that the dimensionless power spectrum $\Delta(k)\equiv k^{3}P(k)/(2\pi^{2})$ for $\beta^{i}_{V}$ is given by [79]

[TABLE]

where

[TABLE]

and $v=k^{\prime}/k$, $u=\sqrt{1+v^{2}-2v\cos{\theta}}$, with $\cos{\theta}=k^{\prime i}k_{i}/\sqrt{k^{\prime i}k^{\prime}_{i}k^{j}k_{j}}$. Due to the convolution in (6.1) a given $k$-mode of $\Delta_{\beta^{V}}$ can receive contributions from arbitrarily short and long-wavelength $k^{\prime}$-modes of $\Delta_{\delta\delta}$, $\Delta_{vv}$ and $\Delta_{\delta v}$. This makes its comparison against simulation results intrincate as the latter has intrinsic cuttoff scales due to its discrete nature; on the one hand, the simulation cannot access long-wavelength modes beyond the fundamental mode $k=2\pi/L$, while on the other hand the contributions coming from short-wavelength modes can be contaminated by modes beyond the Nyquist frequency $k=\pi N_{p}^{{1/3}}/L$. In principle, there is no clear correspondence between the aforementioned modes and the cutoff scales needed in (6.1) in order to compare faithfully against the simulation results. For our current comparison we adopt such modes as hard cuttoffs and show results only up to $25\%$ of the Nyquist frequency [20].

We note that, as a consequence of the method for solving the vector Laplacian equations (4.6) and (4.10), there is no complete separation of scalar and vector modes of $\beta^{i}$ and $W_{i}$ in (4.21)-(4.23). While (4.22) and (4.24) guarantee that $U$ and $b$ are scalars, $B^{i}$ and $V_{i}$ contain both scalar and vector modes. Hence, in order to safely extract $\beta^{i}_{V}=B^{i}_{V}$ from $B^{i}$ as a post-processing we apply a discrete curl operator ($\nabla\times$) to remove any scalar component in the latter, after which we can use the relation $k^{2}P_{\beta^{V}}(k)=P_{\nabla\times\beta}(k)$ to get $\Delta_{\beta^{V}}$. As a consistency check we have also calculated the curl of the full momentum density $s_{i}$ which sources $\beta^{i}$ through the momentum constraint Eq. (4.22), a procedure that has been previously applied to extract the shift vector from Newtonian simulations under a post-Friedmann approach in [80, 81].

Figure 11 shows the power spectra for $\beta^{i}_{s}$ (left panel) and $\beta^{i}_{V}$ (right panel) extracted from the high-resolution gramses simulation at three redshifts: $z=9$ (red), $4$ (blue) and $1$ (magenta). The solid curves denote the predictions of linear perturbation theory (for $\beta^{i}_{s}$) and second order perturbation theory (for $\beta_{V}^{i}$), while the symbols are simulation results; for $\beta_{V}^{i}$ the open triangles and filled circles represent respectively the power spectra for $s_{i}^{V}$ and $\beta_{V}^{i}$. In both cases, the solid curves and symbols agree well, especially at higher redshifts.

We find that on large scales the scalar mode of the shift vector (left panel) can be many orders of magnitude larger than the vector mode (right panel), and thus the curl method described above is necessary to isolate the latter. As we discussed below Eq. (2.25), the existence of this scalar mode within the shift vector is a result of the MD gauge condition (2.22), and shows that the fully constrained formulation adopted in gramses successfully accounts for all the scalar and vector modes of the metric. The specification of cutoffs for the second order perturbation theory prediction (6.1) leads to a power loss in $\Delta_{\beta^{V}}$ and can affect the matching to simulation results specially on large scales. However, the good agreement between perturbative and numerical results shown in the right panel of Fig. 11 for $z=9$ and $z=4$ is consistent with previous studies showing that the vorticity component of the velocity field (which is absent in the PT calculation) is subdominant with respect to the product of the first order scalar perturbations $\delta$ and $v$ [80, 82].

7 Discussion and conclusions

In this paper we have presented gramses code, a new implementation of general relativistic $N$-body AMR simulations in cosmology. This code brings together several advantages of current GR codes under a fully constrained formulation in which dynamical DOF of the metric can be isolated in a consistent and convenient way for cosmological simulations of structure formation. The current version of gramses does not include the evolution equations for tensor modes, the omission of which is expected to have negligible impact on this kind of simulations and whose implementation is therefore left as a future project. Combined with the $N$-body particle methods for standard cosmological simulations, it can accurately solve the nonlinear Einstein equations all the way down to deeply nonlinear scales on which particle orbit crossing is frequent. We have discussed the relevant theoretical background upon which this is based as well as its numerical implementation in detail.

Multigrid relaxation, adaptive mesh refinement and efficient MPI parallelisation are three of the key features of gramses, which are inherited from its parent code – the publicly-available $N$-body and hydrodynamical simulation solver ramses [50]. These are what will enable gramses to run large simulations for cosmological volumes while at the same time resolving scales well within virialised structures. In this paper we have opted to not devote much space to the description of these features, as they are well-established techniques nowadays and a lot of details can be found from the ramses code paper, references therein and its derived literature. However, let us briefly mention that the computation of geometric source terms for the Einstein equations in the ADM formalism – such as $\bar{A}_{ij}\bar{A}^{ij}$ and $\partial_{i}\bar{A}^{ij}$ – are actually nontrivial at the boundaries of refined meshes, which is largely due to the fact that $\bar{A}_{ij}$ itself already involves up to second-order derivatives of the GR potentials $V_{i}$ and $U$. Furthermore, the use of AMR increases the memory requirement of the code: without AMR, the different GR potentials, for a given time step, can be calculated sequentially so that a same array can be recycled for storing them; with AMR and the associated adaptive timestepping, the GR potentials must be kept for longer because they are needed to set up the boundary conditions for the refinement in all subtimesteps, and so we opt to create a separate array for each GR potential.

We have conducted various code tests to verify that all the source terms are computed correctly and that the implemented multigrid relaxation and geodesic solvers are reproducing the correct solutions for test cases where these can be calculated analytically or can be derived in alternative ways. We also successfully run a cosmological simulation for a $\Lambda$CDM universe in a general relativistic setting. The maps produced in this GR simulation show expected features, and a more detailed quantitative analysis shows that the matter and velocity power spectra from the GR simulation display the expected behaviour for the CMC-MD gauge on large scales. The $L=4h^{-1}$Gpc test simulations used to compare the relativistic and Newtonian codes also show that gramses is roughly ten times more computationally expensive than its parent code – ramses. This number can be estimated from the ratio between the numbers of PDEs that are needed to solve the gravitational sector in each code, taking into account the fact that the non-Poisson-like equations solved in gramses are only mildly nonlinear so that they are not significantly more expensive to solve than the standard Poisson equation. Nonetheless, compared to standard numerical relativity codes, an advantage of gramses is that no evolution equations for gravity need to be solved in the formalism implemented so that relatively large time steps can be used. When GW are also implemented, their finite propagation speed impose additional requirements for the time stepping of the simulations, which could substantially increase the computational costs (for cosmological simulations, this should not be a serious issue).

A few words to compare gramses and gevolution are in order here, considering that they appear to share some common features. While the current implementation of gramses solves all the used metric DOF from nonlinear elliptic equations, gevolution employs a mixture of linearised elliptic, parabolic and hyperbolic equations. As gevolution adopts the Poisson gauge, the vector modes are cleanly isolated within the shift vector $\beta^{i}$. This is not the case in gramses, in which the elliptic equations mix one scalar mode in $\beta^{i}$, and a post-processing step is required if their separation is needed, as discussed in Section 6.2. For tensor modes, in gevolution they are solved from a linearised hyperbolic equation, while these are not taken into account in the current gramses version, although the framework allows to include them by solving additional nonlinear time evolution equations. In any case, the facts that gevolution is based on the weak field expansion and utilises Fourier transform in the gravitational solver may boost its overall performance, while the AMR nature of gramses may make it more suitable in situations where one wishes to focus on high-density regions.

The generation of initial condition for gramses simulations will be addressed in a forthcoming paper. Nonetheless, we can mention here that in order to generate initial conditions in the CMC-MD gauge we have modified the camb and 2lptic [83] codes. The former is used to generate the matter power spectrum implementing the gauge transformations at the initial redshift $z_{\rm ini}$ of the simulation, but also on two neighbouring ones (one at a slightly higher while the other a slightly lower redshift). Then, our modified 2lptic code realizes this power spectrum to calculate the particle positions in the standard way, but the velocity is calculated by finite-differentiating the particle displacements of the two neighbouring snapshots. This has the advantage of being independent of an explicit parameterisation of the growth factor (and growth rate) which can become scale-dependent in certain gauges and dark energy/modified gravity models and thus violate a basic assumption (scale-independent linear growth rate) of the default 2lptic code.

The code can be particularly useful for capturing relativistic effects in large-scale simulations as well as for studying problems beyond the Newtonian approximation such as the effect of large density contrasts on the surrounding spacetime and its impact on the estimation of cosmological parameters. Due to perturbed photon geodesics, distant objects can have observed redshifts, angular positions and fluxes which deviate from the true values. This can lead to various effects in observations such as in the galaxy density field [84, 85, 86]. While subtle, these effects are detectable using suitable estimators, and their detectability varies with redshift, scale and estimator [e.g., 87, 88, 89, 90, 91, 92, 93]. To fully accurately study such effects in future galaxy surveys and how they could impact on the estimation of cosmological parameters and constraints on models, realistic mock galaxy catalogues based on cosmological simulations are needed [90]; such simulations should ideally have large volumes to capture the very large scale effects, but also resolve nonlinear scales and even virialised objects to more accurately predict the nonlinear effects [e.g., 94, 95]. Another interesting topic is the back-reaction effect of space-time averaging on the expansion rate [4], an effect the exact size of which is still being debated and which can have important implications on the understanding of the cosmic acceleration or cosmological parameter estimations, which is a nonlinear effect that cannot be fully captured by linear perturbation theory [13] and is therefore better to be quantified using $N$-body simulations that solve the nonlinear Einstein equations inside the most nonlinear virialised objects. In addition, the implementation of scalar field dark energy or modified gravity models in gramses is potentially interesting as it would allow to study the evolution of the new DOF inside virialised dark matter haloes and their impact on the GR effects [96, 97].

Acknowledgments

We are grateful to Marius Cautun for the help with the dtfe code and to Alexander Eggemeier for help on the shift vector power spectrum calculations. We thank Julian Adamek, Marco Bruni, Kazuya Koyama, James Mertens and Boudeiwjn Roukema for useful discussions and comments on this work. We also thank the anonymous referee for their useful comments which helped us to improve the manuscript. CB-H is supported by the Chilean National Commission for Scientific and Technological Research through grant CONICYT/Becas-Chile (No. 72180214). BL is supported by STFC Consolidated Grant (Nos. ST/I00162X/1, ST/P000541/1) and a European Research Council Starting Grant (ERC-StG-716532-PUNCA). This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operation grant ST/R000832/1. DiRAC is part of the UK National e-Infrastructure.

Appendix A Leapfrog scheme for time evolution of particles

gramses uses the usual leapfrog or Stormer-Verlet scheme for particle movements. In this scheme, the position and momentum (or velocity) of a given particle from step $n$ to step $(n+1)$, with a time interval $\Delta t$, are updated using the following prescription,

[TABLE]

where $n+1/2$ is the middle between the two neighbouring timesteps $t^{n}$ and $t^{n+1}$, $p$ is the conjugated momenta to the canonical variable $x$, $H(x,p)$ is the Hamiltonian of the system, and $H_{x},H_{p}$ are the partial derivatives of $H(x,p)$ with respect to $x$ and $p$ respectively. In the case of Newtonian gravity $H=p^{2}/2m+\Phi_{N}(x)$, where $p=mv$ and $\Phi_{N}(x)$ is the Newtonian gravitational potential, the Hamiltonian $H$ is completely separable for $x$ and $v$, and the above operations reduce to the standard Kick-Drift-Kick (KDK) scheme

[TABLE]

Apparently, this makes the system explicit, i.e., the right-hand sides of Eqs. (A.4)-(A.6) do not depend on the quantities (which are at step $(n+1)$) on the left-hand sides. Although from Eqs. (A.4) and (A.6) it would seem that we need to do two force calculations per time step, this is not actually the case since $\partial_{x}\Phi_{N}(x^{n+1})$ in (A.6) is the same as the force in (A.4) in the next step, so the second Kick (A.6) operation can wait until the following (i.e., the $(n+1)$th) timestep when $\Phi_{N}(x^{n+1})$ has been solved. In other words, in practice the second Kick operation (A.6) for the $n$the timestep is done after the Newtonian potential is solved in the $(n+1)$th timestep; in ramses this is called synchronisation as it finally ‘synchronises’ all particle velocities to the correct time before these velocities can be used to move (Drift) the particles.

In the case of GR, the Hamiltonian of a free particle,in the 3+1 formalism is

[TABLE]

where $p_{i}=mu_{i}$ is the (spatial) momentum. Using Hamilton’s equations

[TABLE]

we can derive the equation of motion for this system, i.e., the geodesic equations (4.46)-(4.48). We note that in this case the Hamiltonian (A.7) is not separable, because there is the multiplication of $\gamma^{ij}$ (which depends on $x$) and $p_{i}p_{j}$ under the square root. Therefore, the leapfrog system (A.1)-(A.3) is implicit and not straightforward to implement as in the Newtonian case. The simplest approximation to make the system explicit is to evaluate the Hamiltonian derivatives at the wrong phases, i.e.,

[TABLE]

where

[TABLE]

Notice that here, for evaluating $\mathcal{V}(x^{n},u^{n+1/2})$ used in (A.11) according to (A.14), we use the current value of the gravitational fields ($\Phi,\Psi,\alpha$ and $\beta_{i}$) and only the quantities depending on $u_{j}$ (including $W$) are updated (to $t^{n+1/2}$). Likewise, in (A.12) the force term $F(x^{n+1},u^{n+1/2})$ uses the updated velocities for the explicit dependence on ${u}^{n+1/2}$ as well as for the source terms for the fields at the new timestep. Finally, for repeating the process, in (A.10) we use the updated velocities for the explicit dependences on ${u}^{n}$ and geometric fields based on the updated particles positions (with sources at ${u}^{n+1/2})$.

Let us remark that even if in the Hamiltonian formalism the variables $(x,p)$ are independent (conjugated) variables, and in the second Kick step (A.12) the various gravitational fields appearing in $F(x^{n+1},u^{n+1/2})$ are solved at the final positions, e.g. $\Phi(x^{n+1})$, the source terms for their equations have used the velocities $u^{n+1/2}$ because we have not yet synchronised by the time we evaluate these sources at timesteim $(n+1)$, and thus the fields carry a delayed information about the velocities by half a timestep. Again, this issue is not present in the Newtonian case since the gravitational field $\Phi_{N}$ is only sourced by the mass density field which depends only on the particles position but not on their velocities. A possible way to get around this is to temporarily update the velocity before carrying (A.12) using Poisson equation for the Newtonian gravitational potential $\Phi_{N}$, giving us an estimated updated velocity, namely

[TABLE]

which can be used (as an approximation) in the source terms for the GR potentials. After solving the field equations for these the velocity is then reverted back to $u^{n+1/2}$, after which ‘true’ synchronisation (A.12) is performed.

In principle, the above scheme could be further improved by introducing an extra step to update the position in such a way that the symplecticity of the scheme is restored (although time-reversal invariance is still broken). However, since for simulations with AMR the adaptive timesteps render the KDK scheme non-symplectic even in the Newtonian case, we shall not explain these alternatives, which are more complicated, in detail here.