General Relativistic Cosmological N-body Simulations I: time integration

TL;DR

This paper introduces a stable numerical scheme for fully general-relativistic N-body simulations in cosmology, capable of accurately evolving inhomogeneous structures on small scales while maintaining physical constraints.

Contribution

It presents a novel combination of numerical relativity methods and time-integration algorithms tailored for cosmological N-body simulations, ensuring stability and accuracy in highly inhomogeneous regimes.

Findings

Successfully passes robustness tests

Accurately follows scalar linear modes in an expanding universe

Maintains constraint and energy-momentum conservation at small scales

Abstract

This is the first in a series of papers devoted to fully general-relativistic $N$-body simulations applied to late-time cosmology. The purpose of this paper is to present the combination of a numerical relativity scheme, discretization method and time-integration algorithm that provides satisfyingly stable evolution. More precisely, we show that it is able to pass a robustness test and to follow scalar linear modes around an expanding homogeneous and isotropic space-time. Most importantly, it is able to evolve typical cosmological initial conditions on comoving scales down to tenths of megaparsecs with controlled constraint and energy-momentum conservation violations all the way down to the regime of strong inhomogeneity.

Equations

$$8\pi G=c=\mathrm{Mpc}/h=1,$$

$${\mathrm{D}}_t\Theta$$

$${\mathrm{D}}_t{\hat{\Gamma }}^i$$

$${\mathrm{D}}_t\chi$$

$${\mathrm{D}}_{t}K$$

$${\mathrm{D}}_t{\tilde{\gamma }}_{ij}$$

$${\mathrm{D}}_t{\tilde{A}}_{ij}$$

$${\mathrm{D}}_t$$

$${\tilde{\Gamma }}^i$$

$${\tilde{Z}}^i$$

$${\tilde{\Gamma }}_{kij}$$

$${\tilde{\Gamma }}_{ij}^k$$

$${\hat{R}}_{ij}$$

$$D$$

$$D^{\prime }$$

$$\Theta$$

$${\tilde{Z}}^i$$

$$H$$

$$M_i$$

$$\mathrm{d}{s}^2=-{\alpha }^2\mathrm{d}{t}^2+{\gamma }_{ij}\left(\mathrm{d}{x}^i+{\beta }^i\mathrm{d}t\right)\left(\mathrm{d}{x}^j+{\beta }^j\mathrm{d}t\right),$$

$${\gamma }_{ij}:={\chi }^{-1}{\tilde{\gamma }}_{ij},K_{ij}:={\chi }^{-1}\left({\tilde{A}}_{ij}+\frac{1}{3}{\tilde{\gamma }}_{ij}K\right),$$

$${\mathrm{D}}_t\alpha =-\frac{1}{3}{\alpha }^2\left[K-2\Theta \right],$$

$${\beta }^i=0.$$

$$\alpha =Q{\chi }^{-1/2},$$

$$\dot{\Theta }=-\alpha \left[K+{\kappa }_1\left(2+{\kappa }_2\right)\right]\Theta +\dots ,{\dot{\tilde{Z}}}^i=-\alpha \left(\frac{2}{3}K+{\kappa }_1\right){\tilde{Z}}^i+\dots ,$$

$${\kappa }_1\ge 0,{\kappa }_2\ge -1.$$

$${\kappa }_1=-\frac{2}{3}K,{\kappa }_2=-\frac{1}{2},$$

$${\dot{x}}_a^i$$

$${\dot{p}}_i^a$$

$$E_a:=\sqrt{m^2+\chi {\tilde{\gamma }}^{ij}{p}_i^a{p}_j^a},$$

$${\partial }_{i}f(\vec{x})\to \frac{-2f(\vec{x}+2{\vec{\Delta }}_i)+8f(\vec{x}+{\vec{\Delta }}_i)-8f(\vec{x}-{\vec{\Delta }}_i)+2f(\vec{x}-2{\vec{\Delta }}_i)}{12\Delta x},$$

$$({\vec{\Delta }}_i{)}_j:={\delta }_{ij}\Delta x,$$

$$\beta^{i}\partial_{i}f(\vec{x})\to\frac{\beta^{i}}{12\Delta x}\times\left\{\begin{array}[]{ccc}f(\vec{x}+3\vec{\Delta}_{i})-6f(\vec{x}+2\vec{\Delta}_{i})+18f(\vec{x}+\vec{\Delta}_{i})-10f(\vec{x})-3f(\vec{x}-\vec{\Delta}_{i})&{\rm if}&\beta^{i}>0\\ -f(\vec{x}-3\vec{\Delta}_{i})+6f(\vec{x}-2\vec{\Delta}_{i})-18f(\vec{x}-\vec{\Delta}_{i})+10f(\vec{x})+3f(\vec{x}+\vec{\Delta}_{i})&{\rm if}&\beta^{i}<0\end{array}\right.\,.$$

$${\partial }_i{\partial }_{i}f(\vec{x})\to \frac{-f(\vec{x}+2{\vec{\Delta }}_i)+16f(\vec{x}+{\vec{\Delta }}_i)-30f(\vec{x})+16f(\vec{x}-{\vec{\Delta }}_i)-f(\vec{x}-2{\vec{\Delta }}_i)}{12\Delta x^2}.$$

$$E(\vec{x})$$

$$P_i(\vec{x})$$

$$S_{ij}(\vec{x})$$

$${\wedge }^{(3)}(\vec{x}):=\prod _{i=1}^3\wedge (x^i),$$

$$\wedge(x):=\left\{\begin{array}[]{cc}\frac{3}{4}+\left(\frac{x}{\Delta x}\right)^{2}&{\rm if}\,\,\,0\leq|x|\leq\frac{1}{2}\,\Delta x\\ \frac{1}{2}\left(\frac{3}{2}-\frac{x}{\Delta x}\right)^{2}&{\rm if}\,\,\,\frac{1}{2}\,\Delta x\leq x\leq\frac{3}{2}\,\Delta x\\ 0&{\rm otherwise}\end{array}\right.\,.$$

Full text

1 Introduction

The forthcoming advances in the observations of the cosmological large scale structure (LSS) [1, 2, 3, 4, 5] require a proportionate refinement of our theoretical predictions, not only to exploit the increased amount and precision of the data, but also in order to correctly interpret them. The standard numerical approach to study the non-linear LSS dynamics is the Newtonian $N$-body simulation [6, 7, 8], which essentially emulates the Boltzmann equation for “cold” collisionless matter in Newtonian gravity. Such simulations ignore the relativistic effects of General Relativity (GR) in the dynamics, but also in the reconstruction of observables, since they do not take into account the full geometrical information of space-time. The Newtonian approximation only applies to cosmological models where matter is non-relativistic and effectively decoupled from relativistic degrees of freedom, such as in $\Lambda$CDM, thus excluding several alternative descriptions of the dark sector. Moreover, it also fails at scales comparable to the Hubble radius, which the forthcoming missions will be able to probe. At such scales the causality imposed by relativity can no longer be ignored and relativistic effects are known to become important, at least in the observables [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30].

Nevertheless, these limitations can be circumvented to some extent with the help of analytical tools that were developed in the last decade. There are now refined perturbative expansions of the Einstein equations around the Friedmann-Lemaître-Robertson-Walker solution (FLRW) that are able to capture the non-linear matter dynamics [31, 32, 33, 34, 35, 36, 37, 38, 39, 40]. Along with mapping techniques or appropriate gauge choices, one can then use Newtonian $N$-body simulations to effectively solve the non-linear dynamics of the relativistic theory [41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63]. In [64, 65, 66, 67, 68] the authors went one step further by developing the first $N$-body code based on such a truncation of the Einstein equations, thus including all the information of the metric tensor and capturing the dominant relativistic effects.

Although the above methods are certainly very convenient, they are still defined within a perturbative approach. Some argue111See [69, 70, 71, 72, 73, 74, 75, 76] for reviews and discussions on the issue of “backreaction” of small scale inhomogeneities on the large scale dynamics, see [77, 31, 42, 78, 79, 80] for counter-arguments and [81, 82, 83, 84, 85, 86, 87, 88, 89] for related numerical investigations. that, in the presence of strong inhomogeneity, there could be important non-perturbative effects invalidating any perturbative treatment. This has motivated a recent interest in simulations that solve the fully non-linear Einstein equations [84, 85, 90, 91, 86, 92, 93, 94, 87, 88, 89], i.e. the application of numerical relativity (NR) to cosmology. In these cases, however, the matter sector has always been modeled as a pressureless perfect fluid field in a grid-based approach. Consequently, it cannot describe the correct (collisionless) dynamics at scales where shell-crossing occurs, which roughly coincides with the scales at which the dynamics become non-linear. One way of describing a cosmology with “granular” matter within NR, which has also received particular focus, are simulations of lattice black hole configurations [95, 96, 97, 98, 99, 100, 101, 102, 103], but the high degree of symmetry makes such solutions too idealized to describe realistic dark matter dynamics.

The status quo naturally leads us to consider the potential advantages of an $N$-body NR approach, i.e. taking into account all non-linear and relativistic effects, while solving the correct matter dynamics at small scales. On the one hand, such simulations would serve as a control reference for comparing with approximative methods, both analytical and numerical, thus testing their robustness and potentially settling issues if ambiguous results arise between different approaches. On the other hand, if any non-perturbative effects turn out to occur, be it those that have been speculated over or genuinely new ones, such codes would be the only way to capture them.

The first $N$-body NR simulations have been performed in the context of gravitational collapse dynamics in the mid-eighties [104, 105, 106, 107, 108, 109] for configurations of reduced dimensionality, while the first studies of three-dimensional configurations occurred in the late nineties [110, 111], with a recent revival in [112, 113, 114]. In [115] it is the case of massless particles that was considered to study the collision of plane-fronted gravitational waves. The only appearance of such simulations applied to cosmology is, to our knowledge, in [115, 116, 117]. However, the most complicated configuration considered in these papers is the triple-mode inhomogeneity around the FLRW space-time with equal comoving wavelengths, i.e. $\delta(\vec{x})=\sum_{i=1}^{3}A_{i}\sin(kx^{i})$. In [115, 116] the authors consider the case $A_{i}\sim 10^{-3}-10^{-2}$ and a comoving wavelength that is four times the initial comoving Hubble radius. In [117], where the aim is to study quantitatively the deviation from linear cosmological perturbation theory, the authors consider comoving wavelengths of 400 and 100 Mpc with an amplitude $A$ corresponding to the typical power at these scales. Finally, let us also mention another recent $N$-body code for cosmological simulations [118] which employs the so called “fully constrained formulation” of GR [119, 120, 121] with an approximation which essentially neglects its propagating degrees of freedom (gravitational waves). This reduces the problem to a set of non-linear elliptic field equations, as in the case of Newtonian $N$-body simulations.

Our aim is to develop a combination of numerical methods that will ultimately allow one to perform realistic three-dimensional cosmological $N$-body NR simulations. This endeavor presents new computational challenges compared to Newtonian $N$-body and grid-based NR codes, which we will address in a series of papers. The present paper focuses on the time evolution of the system and its stability, providing in particular the numerical relativity scheme, discretization method and time-integration algorithm. We show that this combination is able to solve the FLRW solution robustly, i.e. it is stable under the injection of white noise, and it can accurately follow a scalar linear mode fluctuation around that solution. The most important result is that it is able to evolve typical cosmological initial configurations for a pure dark matter universe, resolving comoving scales down to tenths of megaparsecs (see figures 1 and 2), with a controlled violation of the constraints and of the energy-momentum conservation. In particular, the Hamiltonian and momentum constraints are satisfied with an average relative precision of $\sim 10^{-6}$ and $\sim 10^{-2}$, respectively, while the energy-momentum conservation equations are satisfied at the level of $\sim 10^{-2}$, all the way down to the regime of strong inhomogeneity $z\sim 0$. At large enough scales, where the linear theory holds, we also follow the matter power spectrum with a relative precision of $\sim 10^{-2}$. The convergence is of the expected order, but not after redshift $z\sim 10$, and we believe that adaptive resolution, in both the mesh and the phase space sampling, will be able to resolve this issue.

In section 2 we present the involved numerical methods, in section 3 we present the results of the aforementioned tests and in section 4 we conclude. All the required equations related to our tests are derived in the appendices. We work in the following units

[TABLE]

and all the present tests are performed with zero cosmological constant. Our code is implemented on top of the latfield2 library [122].

2 Numerical methods

2.1 Evolution equations

On the gravitational side we consider the damped CCZ4 scheme [123, 124, 125], which can be seen as a generalization of the BSSNOK scheme [126, 127, 128] involving an additional “pure-constraint” field $\Theta$ that helps diluting constraint violation by propagating it away. We base our choice of scheme on a recent “Apples-with-Apples” comparison [129, 130] of CCZ4 with closely related ones that we performed in [131], and also the tests performed in this paper. The evolution equations are

[TABLE]

where

[TABLE]

and are subject to the constraint equations

[TABLE]

All indices are displaced using the conformal 3-metric $\tilde{\gamma}_{ij}$. The line-element reads

[TABLE]

where $\alpha$ is the lapse function, $\beta^{i}$ is the shift vector,

[TABLE]

are the 3-metric and extrinsic curvature of the $t={\rm const.}$ hypersurfaces, while $E$, $P_{i}$ and $S_{ij}$ are the energy, momentum and stress densities in the canonical frame $n=\alpha^{-1}\left(\partial_{t}-\beta^{i}\partial_{i}\right)$. By redefining $K\to K+2\Theta$ one obtains the equally well-performing Z4cc scheme considered in [131], which is then related to BSSNOK by simply setting $\Theta=0$.222By setting to zero some of the pure-constraint terms of Z4cc, one obtains the Z4c scheme proposed in [132, 133] and tested in [134], which we included in our comparison [131]. We refer the reader to [131] for a derivation of the above equations.

In this paper we will consider the following Bona-Masó slicing [135, 136]

[TABLE]

and zero shift vector

[TABLE]

This slicing choice corresponds to the conformal time parametrization on the FLRW solution, which we will keep denoting by “$t$” contrary to the usual convention in cosmology. For the analytical solution where $\Theta=0$, equation (2.21) can be solved analytically

[TABLE]

where $Q\equiv Q(\vec{x})$ is an arbitrary space field and corresponds to the residual gauge freedom of choosing the initial conditions of $\alpha$.

The parameters $\kappa_{1}$ and $\kappa_{2}$, introduced in [124], are free to choose and can be space-time dependent, since they multiply “pure-constraint” terms. The evolution equations of the Z4 fields $\Theta$ and $\tilde{Z}^{i}$ take the form

[TABLE]

where the ellipses denote either second-order terms in perturbations around FLRW or source terms, so these are the “linear” parts of the equations. We see that $\kappa_{1}$ and $\kappa_{2}$ are damping parameters, i.e. for appropriate values they push the system towards the constraint surface. For linear fluctuations around Minkowski space-time where $K=0$, demanding constraint stability leads to the bounds [124]

[TABLE]

Around FLRW space-time, however, we have that $K<0$ in our gauge, so the corresponding terms in (2.24) come with the wrong sign. Therefore, the “undamped” CCZ4 system $\kappa_{1}=0$ is not stable in the cosmological context, especially at early times where $|K|$ is large, because $\Theta$ and $\tilde{Z}^{i}$ diverge exponentially. Instead, the effectively undamped scheme in cosmology is the one corresponding to

[TABLE]

because this way the terms displayed in (2.24) cancel out. Greater values of $\kappa_{1}$ would then reintroduce a damping effect. For the tests performed in this paper, we will exclusively work with (2.26).

Finally, on the matter side we have a set of $N$ free-falling particles of mass $m$, with positions $x_{a}^{i}(t)$ and momenta $p_{i}^{a}(t)$, where $a=1,\dots,N$. The derivation of the corresponding evolution equations is given in appendix A, the result being the geodesic equation in first-order form

[TABLE]

where all gravitational fields appearing here are implicitly evaluated at $x^{i}_{a}$ and

[TABLE]

is the energy of the $a$-th particle. The corresponding energy-momentum tensor components are also derived in appendix A and are given explicitly in their discretized version in the following subsection.

2.2 Space discretization

We discretize the field equations on a Cartesian mesh using finite difference methods. In particular, for the spatial derivatives we use a centered five-point stencil, i.e.

[TABLE]

where

[TABLE]

with $\Delta x$ the lattice spacing. There are two exceptions to this. First, the $\partial_{i}$ appearing inside the convective derivative ${\rm D}_{t}:=\partial_{t}-\beta^{i}\partial_{i}$ is replaced with the up/down-wind five-point stencil, depending on the sign of $\beta^{i}$

[TABLE]

Second, whenever we have double derivatives $\partial_{i}\partial_{j}$, the diagonal terms are replaced with the second derivative centered five-point stencil

[TABLE]

As for the particle-mesh communication, the energy-momentum components are constructed by projecting the particle information according to

[TABLE]

where

[TABLE]

and $\wedge$ denotes the triangle-shaped cloud function

[TABLE]

For the interpolation of field values at particle positions we then use the inverse kernel.

2.3 Time integration

At the level of the FLRW space-time, our gauge corresponds to conformal time and comoving spatial coordinates, in terms of which light-like propagation corresponds to the same relation as in Minkowski space-time, i.e. $\Delta t=\Delta x$. For fluctuations around FLRW, the latter serves as a background space-time determining the causal structure of the dynamics, so it makes sense to consider a constant Courant factor in time

[TABLE]

which is the parameter relating the time step $\Delta t$ to the considered lattice spacing $\Delta x$. We have found that a satisfactory evolution is provided by a fourth-order Runge-Kutta (RK4) for the gravitational fields and a “drift-kick-drift” for the particles. However, the particles are evolved only every $N_{s}$ cycles with time step $N_{s}\Delta t$. We found that constraint violation is significantly reduced when $N_{s}$ is around 10 and in this paper we will consider for definiteness $N_{s}=8$.333Note that $N_{s}>1$ poses no problem for the resolution of the particle dynamics, because their typical velocities are smaller than the speed of light by several orders of magnitude. Given our particle integration method, the time step is therefore $N_{s}\Delta t$ for the particle positions and $2N_{s}\Delta t$ for their momenta. Denoting by $G$ the set of gravitational fields $\{\alpha,\partial_{i}\alpha,\beta^{i},\partial_{i}\beta^{i},\chi,\partial_{i}\chi,\tilde{\gamma}_{ij},\tilde{\Gamma}^{k}_{ij}\}$ that must be interpolated at the particle positions (2.27), the time-integration loop is described as follows:

The particles are displaced according to (drift)

[TABLE] 2. 2.

The matter fields $E$, $P_{i}$ and $S_{ij}$ are updated. 3. 3.

The gravitational fields are evolved by $\Delta t$ through RK4 $N_{s}$ times. 4. 4.

The particle momenta are updated according to (kick)

[TABLE] 5. 5.

The matter fields $E$, $P_{i}$ and $S_{ij}$ are updated. 6. 6.

The gravitational fields are evolved by $\Delta t$ through RK4 $N_{s}$ times. 7. 7.

The particles are displaced according to (drift)

[TABLE] 8. 8.

Send $t\to t+2N_{s}\Delta t$.

An illustration of this loop is found in figure 3. Note that we impose the constraint (2.14) by hand

[TABLE]

at each sub-step of the RK4. Moreover, we include Kreiss-Oliger dissipation [137] for the evolution of the gravitational fields, choosing the sixth-order one since we use a fourth-order time-integration. Thus, to the right-hand side of the evolution equation of some field $f$ we add

[TABLE]

The normalization of the parameter $\sigma$ is such that the stability bounds are [136]

[TABLE]

and here we will only consider the value $\sigma=0.05$. Note that, when updating the matter fields (2.33), it is the particle-dependent part that is updated every $N_{s}$ steps, while the factor $\chi^{3/2}$ is updated at every step. Let us also point out the difference between this integration method and a staggered leapfrog where step 1 and 7 would be glued together. With our approach the computation of $x^{i}_{a}(t+2N_{s}\Delta t)$ involves the gravitational fields evaluated at both $x^{i}_{a}(t)$ and $x^{i}_{a}(t+N_{s}\Delta t)$. This turns out to yield a better resolution of the particle dynamics and constraint violation control.

3 Tests

3.1 General definitions and specifications

As far as the constraints are concerned, we will only display the Hamiltonian and momentum ones $H$ and $M_{i}$, respectively, given in (2.17) and (2.18). We have monitored the rest of them $\{D,\Theta,\tilde{Z}^{i}\}$444Remember that $D^{\prime}$ is imposed algebraically as mentioned in the previous section. as well and found that they are controlled better than $H$ and $M_{i}$. We will consider the absolute values $|H|$ and $|M|:=\sqrt{\gamma^{ij}M_{i}M_{j}}$, but also the more relevant relative quantities suggested in [90]

[TABLE]

where $T^{H}_{n}$ and $T^{M}_{n,i}$ denote the $n$-th term appearing on the right-hand side of (2.17) and (2.18), respectively. These relative errors therefore capture the number of significant digits at which the cancellation in $H$ and $M_{i}$ occurs. Given some error field $\phi(\vec{x})$ on the lattice, the measures we will output are

[TABLE]

where $g$ indexes the grid points and $N_{g}^{3}:=(L/\Delta x)^{3}$ is their total number. In all three tests there is some relation to the FLRW solution which we consider in the case of zero spatial curvature and zero cosmological constant. On this solution the non-zero field components are

[TABLE]

where

[TABLE]

is the conformal Hubble parameter and we have chosen the normalization $Q=1$ (see (2.23)). Remember that, given our choice of lapse (2.21), $t$ is conformal time.555In cosmological perturbation theory the latter is usually denoted by “$\tau$” or “$\eta$” and its derivative is denoted by a prime, but here this is the variable with respect to which we solve our equations numerically, which is why we use $t$ and a dot, respectively. Indeed, from the perspective of numerical relativity, it makes no sense to change the symbol for the time variable only because we have specified the lapse function. For all simulations our initial time will always be $t=0$ and

[TABLE]

where $t_{f}$ is the final time. With these we have

[TABLE]

and we will express time-evolution either with respect to $t$, or the corresponding FLRW redshift

[TABLE]

In our units ${\rm Mpc}/h=c=1$ the Hubble constant is

[TABLE]

All of our runs start at $z(0)=1000$, so the corresponding final time $(z=0)$ is

[TABLE]

Note, however, that $H_{0}$ and the redshift parametrization cannot be given their realistic interpretations, because we are considering a pure-matter universe. Nevertheless, the corresponding cosmology has the correct orders of magnitude and we have access to strong inhomogeneity by going up to $z=0$.

Finally, we will provide no details on how we generate initial particle data $x_{a}^{i}(0)$ and $p_{i}^{a}(0)$ that reproduce the desired fields $E(0)$ and $P_{i}(0)$, as this will be addressed in another paper of this series. Let us just say that we consider regularly distributed particles with respect to the lattice, before displacing them to obtain the initial positions. The corresponding mass is then determined by

[TABLE]

where $N_{\rm ppc}$ denotes the number of particles per grid cell before displacement. The number $N_{\rm ppc}$ will always be considered constant, meaning that the total number of particles scales as $N\sim\Delta x^{-3}$.

3.2 FLRW robustness test

Here we conduct a robustness test, as defined in [129, 130], but adapted to the FLRW solution instead of Minkowski and also to the inclusion of particles. We thus consider the following perturbation of the FLRW initial conditions

[TABLE]

while for the particles we have

[TABLE]

where $\delta x_{a}^{i}$ is the displacement from the regular configuration. The $\epsilon_{\dots}$ numbers are drawn randomly out of a uniform distribution independently for each component, for each point $\vec{x}$ for the fields and for each particle. The amplitude of the distribution is $10^{-7}N_{g}^{-2}$ for the fields and $10^{-7}N_{g}^{-3}$ for the particles.666For the fields the $\sim N_{g}^{-2}$ dependence is required because of the presence of second-order spatial derivatives, which grow like $\sim N_{g}^{2}$ on random noise. For the particles, we have that the initial displacement field $\delta x^{i}(0)$ is related to the corresponding density contrast through $\partial^{2}\delta x^{i}\sim\partial_{i}\delta$. Thus, since we keep $N_{\rm ppc}$ fixed, the amplitude of $\delta$ grows with resolution increase like $\sim N_{g}$, so the $\delta x^{i}$ amplitude must scale as $\sim N_{g}^{-3}$ for $\delta$ to follow the same trend as the other fields $\sim N^{-2}_{g}$. Indeed, we have checked that there is no convergence for the error on $E$ if the particle amplitude follows $\sim N_{g}^{-2}$ instead of $\sim N_{g}^{-3}$. The runs are performed in a box with comoving size $L=64$ at three spatial resolutions $\Delta x\in\{8,4,2\}$, meaning that $N_{g}\in\{8,16,32\}$. The Courant factor is $C=0.1$ and we use one particle per grid cell $N_{\rm ppc}=1$. Figures 4 and 5 show the relative errors of $\chi$ and $E$ with respect to their respective analytical solutions

[TABLE]

while figure 6 shows the absolute constraint violations $H$ and $M$. In all cases we plot the $L_{\infty}$ measure. We see that $\delta_{E}$ grows in the presence of noise, which can be understood by the fact that this quantity has a growing mode $\sim a$ already at the analytical level, i.e. any inhomogeneity must grow. In the bottom panel of figure 5, we plot the evolution of $\delta_{E}/a$ and see that it is bounded in time, so that the noise is under control in this particular context. The overall verdict is that we are able to follow the analytical solution with good stability and convergence.

3.3 Scalar linear mode test

In this test we check whether the code can accurately evolve a single scalar mode of inhomogeneity in the linear regime around the FLRW solution. We use the definitions, residual gauge choices and the Zel’dovich condition described in appendix B. The initial conditions are therefore completely determined by the gauge-invariant density contrast and here we consider a single mode profile

[TABLE]

where $L$ is the comoving box size. We choose $A=10^{-7}$, which will lead to $\delta_{\star}\sim 10^{-4}$ at redshift zero, thus remaining inside the regime of validity of linear perturbation theory at all times. Moreover, we choose $L=128$, meaning that the mode starts outside the initial comoving Hubble radius ${\cal H}^{-1}(0)\approx 95$ and finishes inside. The rest of the parameters are $C=0.1$, $N_{\rm ppc}=27$ and $\Delta x\in\{4,2,1\}$.

Note that the gravitational potential $\varphi$ must remain constant in time, which is what we see in figure 7, where we have plotted its profile for all three resolutions at both the initial and final times. The right panel is a magnification of the region around the maximum and shows that the initial and final profiles converge towards each other with increasing resolution, although from opposite sides. On the left panel of figure 8 we plot the relative error of $\delta$ with respect to the amplitude of the analytical solution

[TABLE]

which is controlled and converges with resolution. On the right panel we plot the $L_{\infty}$ measure of the relative Hamiltonian constraint $H_{r}$. We see that it diverges with resolution and, in fact, this is the case for all the constraints, for both the absolute and relative cases and for both the $L_{\infty}$ and $L_{1}$ measures. Note that this divergence under resolution increase also occurs in the linearized (gravitational) wave test around Minkowski space-time [131] for the BSSNOK, CCZ4, Z4cc and Z4c schemes (with and without constraint violation damping) for three different gauge choices (see [131] for details). However, this is not observed in the linear regime of the typical cosmology test (multi-mode) of the next subsection, which is the relevant one for cosmology. It therefore seems that the present divergence is an artefact of the plane symmetry of this special configuration.777This divergence is also not observed in the triple mode tests performed in [115, 116, 117], which supports this explanation of the problem. However, there might be also other factors involved in this issue, because in [117] the authors also perform planar-symmetric single mode tests and obtain convergence. They use a BSSNOK scheme, but with different amplitude, wave-length and redshift range, and also a more sophisticated particle-mesh communication method. Moreover, despite this divergence of $L_{\infty}(H_{r})$, its magnitude is still several orders of magnitude smaller that the relative error on $\delta$.

3.4 Typical cosmology test

In this subsection we test the behavior of our code for initial conditions that exhibit the typical inhomogeneities one encounters in cosmology. As in the previous test, we consider again the scalar linear perturbation theory around the FLRW solution with the residual gauge choices and the Zel’dovich condition described in appendix B, but now only for our initial conditions at redshift 1000. For the initial gauge-invariant density contrast $\delta_{\star}$, we use the power spectrum provided by the linear Boltzmann code CLASS [138, 139] to generate a corresponding random field, which is then used to determine the initial field and particle data.888We have checked that our initial conditions respect the Zel’dovich approximation well enough $\varphi+h\sim 10^{-3}\varphi$ and also that the divergence of the momentum dominates its curl by $|\vec{k}\times\vec{P}|\sim 10^{-3}\vec{k}\cdot\vec{P}$ on the resolved scales.

We will consider two simulations: one with box size $L=704$ and one with $L=256$, cutting-off the power at wavelengths $L_{\rm cut}=220$ and $L_{\rm cut}=40$, respectively, and with resolutions $\Delta x\in\{\frac{L}{32},\frac{L}{64},\frac{L}{128}\}$ for the former and $\Delta x\in\{\frac{L}{64},\frac{L}{128},\frac{L}{256}\}$ for the latter. Thus, the cut-off scales $L_{\rm cut}$ correspond to ten times the lattice spacing of the poorest resolution. Note that these cut-offs are imposed at the initial conditions, but evolution will generate structure at smaller scales. The figures 1 and 2 correspond to the smaller box simulation $L=256$. The rest of the parameters are a Courant factor of $C=0.05$ and twenty-seven particles per grid cell $N_{\rm ppc}=27$. We will also denote by

[TABLE]

the convergence ratio of a given quantity $X$, where $X_{\Delta x}$ denotes the value computed with resolution $\Delta x$. Since we have three resolutions for each run, we will have two $C_{\Delta x}[X]$ values for each quantity $X$. Note that, although the field derivatives are computed with fourth-order precision, the particle time integration is of second order, so a successful convergence corresponds to $C_{\Delta x}[X]\geq 2$.

In the $L=704$ case we are mainly considering scales at which shell-crossing is negligible, so the particle dynamics should coincide with the ones of a pressureless perfect fluid. We have chosen the particular numbers $L=704$ and $L_{\rm cut}=220$ in order to compare with the perfect fluid NR simulation performed in [88].999Note that in [88] the units are ${\rm Mpc}=1$ instead of Mpc$/h=1$, as we have here, and the authors use $h=0.704$. More precisely, this matches the simulations of [88] with “controlled number of modes” for which the constraints are plotted as a function of redshift. Moreover, we work with the same gauge as [88]. In figure 9 we plot the $L_{1}$ measure of the relative Hamiltonian and momentum constraints (3.1) as a function of redshift. Note that the $L_{1}$ measure of [88], given in their equation (C4), is the ratio of averages, instead of the average of ratios which we use (3.2), and these two do not obey a definite order relation. Nevertheless, we verify that they are of the same order of magnitude by displaying both. In figures 10 and 11 we plot the $L_{1}$ measure of the absolute constraints $H$ and $M$, respectively, along with the corresponding convergence ratios.

Let us start with the magnitude of constraint violation, thus focusing on figure 9 and the right panels of figure 11 of [88]. We first note that our relative constraints are quite stable in time (except for $H_{r}$ of the poorest resolution run), as opposed to the relative Hamiltonian constraint of [88] which grows until it reaches a plateau value. Moreover, we have three orders of magnitude less error for $L_{1}(H_{r})$ and one order of magnitude less error for $L_{1}(M_{r})$ at redshift zero. As for the convergence of the constraint violation, we compare our figures 10 and 11 to the left panels of figure 11 and to figure 12 of [88]. We find that our Hamiltonian constraint behaves less well, in that it is not really converging at second order after $z\sim 1$. For the momentum, we obtain a clearer separation of the curves at all times, but end up with a convergence of first order only. Our verdict is therefore that we are able to control constraint violation a lot better than in [88] and that our convergence over the full evolution is of comparable quality. It must be stressed, however, that [88] employ the BSSNOK scheme without constraint-damping mechanisms.

Let us now consider the simulation of smaller size $L=256$. The analogues of figures 9, 10 and 11 are now figures 12, 13 and 14, respectively. At high redshift the behavior of the constraints is similar to the previous simulation, i.e. their amplitude and convergence is controlled. However, we observe two significant “jumps”, first at $z\sim 10$, after which convergence is no longer achieved at all, and then at $z\sim 1$. The latter is clearly due to the strong inhomogeneity that develops at these times and should be avoided by using adaptive mesh techniques. As for the first jump at $z\sim 10$, we observe that it coincides with the moment at which the particle number per cell develops spurious sharp variations in space. These are smoothed out when projected on the mesh to build the corresponding fields, but still strong enough to affect constraint violation. These features are subsequently washed away by the formation of structure. We believe that this effect is due to our initial over-sampling of phase space, i.e. $N_{\rm ppc}=27$. The problem is that we cannot lower this parameter, because then we under-sample the voids at late times and this leads to important constraint violation. It therefore seems that adaptive phase space resolution methods will cure this problem. A detailed investigation of this issue will be presented in another paper of this series. Note, also, that despite the aforementioned issue, we are able to control the relative constraints with an $L_{1}$ measure of at most $\sim 10^{-6}$ for the Hamiltonian and $\sim 10^{-2}$ for the momentum, all the way down to redshift zero, for all three resolutions and for both box sizes. Finally, note that the present test is a scalar linear multi-mode test at high enough redshits where evolution is linear. It is therefore interesting to compare with the results of the single linear mode test of the previous subsection. As already mentioned there, we see that here the relative constraints $H_{r}$ and $M_{r}$ converge with resolution (figures 9 and 12), contrary to the single mode case (figure 8).

Although controlling constraint violation is a necessary condition for an accurate resolution of the dynamics, it is not sufficient, and it can even be a misleading check when one uses a scheme that is precisely designed to dissipate that violation. We therefore now consider two more types of error for the $L=256$ simulation. First, we verify that we are accurately resolving the matter power spectrum in the linear regime, for which we have the analytic solution (B.29). Numerically, it is computed by Fourier transforming the field $\delta(t,\vec{x})$ and averaging its modulus squared over the angles

[TABLE]

that is approximated by averaging in $\Delta k=1.28$ shells in practice. In figure 15 we plot this function of $k$ for the three redshift values $z=\{100,10,1\}$ and also the relative difference with respect to the linear analytical solution

[TABLE]

Here $P_{\rm an.}$ is constructed using (B.29) where the initial $P_{\rm an.}(0,k)$ is the one of the initial state of the simulation. Each plot contains all three resolutions and also the analytical one (dashed) and we observe a deviation of the order of $\sim 10^{-2}$ at worst for the best resolution. The growth of error with increasing $k$ is to be expected, since the linear approximation becomes less and less valid.

Second, we also consider whether the violation of energy-momentum conservation is under control. Contrary to the previous check, here we can measure an error at the fully non-linear level, so we do not need to restrict our attention to the large scale Fourier modes only. In terms of the present variables, the conservation equations $\nabla_{\mu}T^{\mu\nu}=0$ read

[TABLE]

and the time-derivative is computed using a three-point stencil, involving three successive loop time-steps, i.e. separated by $N_{s}\Delta t$ (see subsection 2.3). As with the constraint equations, here too the relevant quantities are the relative ones with respect to the typical magnitude of the involved terms

[TABLE]

where $T^{E}_{n}$ and $T^{P}_{n,i}$ denote the $n$-th term appearing on the right-hand side of (3.25) and (3.4), respectively. In figure 16 we show a rough plot of the time evolution of the $L_{1}$ norm of these two quantities, by displaying in particular the value for the redshifts $z=\{1000,800,400,100,50,10,1\}$. We observe an average relative error of $\sim 10^{-2}$ for both $C^{E}$ and $C_{i}^{P}$ for all considered resolutions and redshifts. Finally, as in the case of the constraints, convergence is generically lost at small redshifts.

4 Conclusion

In this paper we have proposed a numerical prescription for the time-evolution of $N$-body NR simulations in cosmology. We have shown that it passes the robustness and scalar linear mode tests around the FLRW solution. We then considered the evolution of typical cosmological initial conditions and showed that our code follows the linear part of the power spectrum accurately and controls well the violation of the constraints and energy-momentum conservation. However, convergence is not achieved for low redshift and small scale simulations. Nevertheless, we have argued that this problem is related to the fact that we work with a Cartesian mesh and with fixed number of particles. It could therefore probably be resolved by considering adaptive resolution methods both for the mesh and the phase space samplers.

Appendix A Derivation of the particle equations

The action for a set of $N$ minimally coupled particles of mass $m$ is given by

[TABLE]

where $\lambda$ is an arbitrary parameter and here the dot denotes the derivative with respect to it. We fix the $\lambda$-reparametrization gauge by requiring that $\lambda$ coincides with the space-time coordinate $x^{0}\equiv t$

[TABLE]

so now the dot coincides with $\partial_{t}$. Note that this is not the proper time of the particle, unless we also choose to fix $\alpha=1$, which is not the case here. We next express the action in terms of the ADM variables

[TABLE]

compute the conjugate momenta of $x^{i}_{a}$

[TABLE]

whose inverse relation reads

[TABLE]

and Legendre transform with respect to $\dot{x}_{a}^{i}$ to get the canonical action

[TABLE]

The canonical energy-momentum components are therefore

[TABLE]

and the equations of motion are

[TABLE]

where it is understood that the above fields are evaluated at $\vec{x}_{a}$ and $\Gamma^{k}_{ij}$ are the Christoffel symbols of $\gamma_{ij}$. Expressing these equations in terms of the conformally decomposed variables one obtains (2.27), while the “discretization” of the Dirac delta in (A.7) yields (2.33).

Appendix B Linear perturbation theory equations

Here we work with the analytical solution (2.23) and model matter as a pressureless perfect fluid, which is a valid description in the regime where linear perturbation theory applies since the velocity field is smooth. We introduce perturbations around the FLRW solution given in (3.3)

[TABLE]

where $\psi$ is time-independent and its profile amounts to the residual gauge choice of initial conditions for $\alpha$. Note also that $v$ is a comoving velocity. Going to Fourier space, the linearized constraint equations read

[TABLE]

while the evolution equations yield

[TABLE]

Using (B.9) to eliminate $\delta$, the scalar evolution equations can be written in closed second-order form

[TABLE]

In terms of these variables the scalar Bardeen potentials are

[TABLE]

while the gauge-invariant matter quantities are

[TABLE]

The scalar constraint equations can now be expressed as

[TABLE]

while (B.15) is nothing but the absence of anisotropic stress $\Psi=\Phi$. Let us now fix the residual gauge freedom, starting with

[TABLE]

so that the linearized gauge conditions read

[TABLE]

These are preserved under a gauge transformation with generating vector $\xi^{\mu}=\left(\xi^{t},\delta^{ij}\partial_{j}\xi\right)$ obeying

[TABLE]

This is a first-order system for $\xi^{t}$ and $\xi$, so their initial conditions are free to choose. To use this freedom we note that, under a general gauge transformation

[TABLE]

where we used (B.11) for the latter. Thus, under a residual gauge transformation at initial time $t=0$, using (B.21), we find

[TABLE]

which allows us to set

[TABLE]

The advantage of this gauge is that, under the assumption of the Zel’dovich condition at $t=0$

[TABLE]

the gravitational fields remain constant in time

[TABLE]

and the corresponding line-element takes the conformal Newtonian form

[TABLE]

We are therefore in the conformal Newtonian gauge in the scalar sector, but only with our choice of initial conditions (B.25) and evolution equations, i.e. the fact that the considered theory is GR. The only evolving quantities are the matter density and velocity, because of the ${\cal H}$ factors in (B.9), and now read

[TABLE]

With the first of these equations we can then infer the evolution of the matter power spectrum

[TABLE]

Acknowledgments

We are grateful to Hayley Macpherson, Joachim Stadel and Romain Teyssier for useful discussions and especially to Martin Kunz for helpful support. This work has been supported by an Advanced Postdoc.Mobility grant of the Swiss National Science Foundation for DD, by a Consolidator Grant of the European Research Council (ERC-2015-CoG grant 680886) for YD and EM and by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID s751. The simulations have been performed on the Baobab cluster of the University of Geneva and on Piz Daint of the CSCS.