#  GrGadget: an N-body TreePM relativistic code for cosmological simulations

Eduardo Quintana-Miranda,1,2,3 Pierluigi Monaco1,2,3,4 and Luca Tornatore1,2,3
1 Dipartimento di Fisica, Sezione di Astronomia, via G.B. Tiepolo 11, I-34143
Trieste, Italy
2 INAF – Istituto Nazionale di Astrofisica, Osservatorio Astronomico di
Trieste, via G.B. Tiepolo 11, I-34143 Trieste, Italy
3 IFPU – Institute for the Fundamental Physics of the Universe, via Beirut 2,
I-34100 Trieste, Italy
4 INFN – Istituto Nazionale di Fisica Nucleare, Via Valerio 2, I-34127
Trieste, Italy E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present the merging of the Particle-Mesh (PM) relativistic Gevolution code
with the TreePM Gadget-4 code, with the aim of studying general relativity
effects in cosmology. Our code, called GrGadget, is able to track the
evolution of metric perturbations in the weak field limit by using
Gevolution’s implementation of a relativistic PM in the Poisson gauge. To
achieve this, starting from Gevolution we have written a C++ library called
Libgevolution, that allows a code to access and use the same abstractions and
resources that Gevolution uses for its PM-only N-body simulations. The code
works under the assumption that particle interactions at short distances can
be approximated as Newtonian, so that we can combine the forces computed with
a Newtonian Tree with those computed with a relativistic PM. The result is a
TreePM simulation code that represents metric perturbations at the scales
where they are relevant, while resolving non-linear structures. We validate
our code by closely matching Gadget-4 forces, computed with the Tree switched
off, with those computed with Libgevolution in the Newtonian limit. With
GrGadget we obtain a matter power spectrum that is compatible with Newtonian
Gadget-4 at small scales and contains GR features at large scales that are
consistent with results obtained with Gevolution. We demonstrate that, due to
the better resolution of the highly non-linear regime, the representation of
the relativistic fields sampled on the mesh improves with respect to the PM-
only simulations.

###### keywords:

cosmology: theory – large-scale structure of the Universe

††pubyear: 2022††pagerange: GrGadget: an N-body TreePM relativistic code for
cosmological simulations–References

##  1 Introduction

The state of the art of precision cosmology provides a standard cosmological
model, $\Lambda$CDM, that is consistent with most observational evidence on
large scales, but relies on the existence of a dark sector populated by Dark
Matter (DM) and Dark Energy (DE). The first is responsible for the formation
of cosmological structures such as galaxies and their large-scale density
field, while the second causes the observed accelerated expansion of the
universe in the present epoch. Their physical nature is an open problem, since
the only evidence of their existence comes from their gravitational
interaction with visible matter. A possible explanation is that the dark
sector is due to a misrepresentation of gravity, that on large scales does not
follow Einstein’s General Relativity (GR), at the basis of the $\Lambda$CDM
model.

This fact has triggered a wave of interest in modifications of GR, that can
lead to extra terms that explain dark energy or dark matter (see, e.g.,
Silvestri & Trodden, 2009; Capozziello & De Laurentis, 2012, and references
therein). Such modifications must be significant only on large scales or low
density, because GR is very accurate in predicting planetary orbits, light
deflection and Doppler effects in solar system tests and has more recently
been successfully tested with the detection of gravitational waves (Abbott et
al., 2016) and the direct imaging of black hole event horizons (Event Horizon
Telescope Collaboration et al., 2019).

In order to characterize dark energy in the age of its dominance, many
projects have been planned to survey large parts of the sky and probe the
large-scale distribution of matter using galaxy clustering and galaxy lensing,
both from the ground (DES111www.darkenergysurvey.org, Krause et al. 2017;
DESI222www.desi.lbl.gov, DESI Collaboration et al. 2016; Rubin’s
LSST333www.lsst.org, Ivezić et al. 2019; SKAO444www.skao.int surveys) and from
space (Euclid555sci.esa.int/web/euclid, Laureijs et al. 2011;
Roman666roman.gsfc.nasa.gov, Spergel et al. 2015;
SphereX777spherex.caltech.edu, Doré et al. 2014). Some of these surveys have
already started to produce a flood of data that will soon lead to a precise
characterization of the galaxy and matter density fields. A comparison of
these observations to model predictions, either using summary statistics or
field-level inference, will lead to unprecedented tests not only of the
cosmological model but also of the gravity theory behind it. With precision
being guaranteed by the amount of available high-quality data, accuracy will
be achieved only by rigorous control of systematics, both in the data and in
theory predictions.

The highly non-linear nature of the observed density field and the non-
locality of gravity make cosmological simulations necessary to compare the
predictions of current theories with the observations at an increasing level
of accuracy. Yet, most of the widely adopted simulation codes, like e.g.
Gadget-4 (Springel et al., 2021), use Newtonian dynamics for the evolution of
matter perturbations. This is not the ideal configuration to pass from the
unobservable distribution of matter in a periodic comoving box to the
observable distribution of light in the past light cone. Relativistic
corrections can be added _a posteriori_ by post-processing Newtonian
simulation outputs; one specific example of this approach is the modeling of
lensing due to the distortion of null geodesics (Bartelmann & Schneider,
2001), while a more comprehensive approach to adding relativistic effects is
presented by Borzyszkowski et al. (2017). However, even though the biases
introduced by this approach are expected to be small, a fully self-consistent
approach is necessary to convincingly demonstrate our ability of controlling
theory systematics. For instance, galaxy clustering is affected by
magnification bias due to lensing, and neglecting this effect induces a non-
negligible bias in parameter estimation (Lepori et al., 2020; Alam et al.,
2021). This is even more true when modified gravity theories are used:
extensions of gravity are typically derived in a full relativistic context,
and while they influence the Newtonian limit of gravity, the small but
measurable relativistic effects may provide smoking-gun signals of a specific
class of gravity theories. In this sense, restricting to the treatment of the
Newtonian limit of modified gravity theories (as, e.g., in Puchwein et al.,
2013) may leave out crucial observable signatures.

Two examples of fully relativistic N-body codes for the evolution of cosmic
perturbations, that integrate Einstein’s equations to follow the motion of
massive particles along their geodesics, are the Adaptive Mesh Refinement
(AMR) code Gramses (Barrera-Hinojosa & Li, 2020) and the Particle-Mesh (PM)
code Gevolution (Adamek et al., 2016). These have proven to be precious tools
to produce accurate cosmological predictions, like a self-consistent treatment
of massive neutrinos (Adamek et al., 2022), and to explore phenomena that were
previously overlooked, like the strength of the frame dragging field acting on
dark matter haloes (Barrera-Hinojosa et al., 2020). These codes sample the
fields in a mesh that fills the simulated volume, but while Gramses uses an
AMR scheme to increase resolution only where it is needed, PM schemes working
on a single non-adaptive mesh are well known to be limited by memory, so they
are unable to achieve the large dynamic range required, e.g., to resolve DM
halos in large cosmological volumes. The integration of Newtonian particle
trajectories has historically been addressed with the introduction of an oct-
tree data structure (Barnes & Hut, 1986), that provides a $N\log N$ scaling
for the computation of gravity without compromising its accuracy. Because the
integration of large-scale perturbations is very slow in this scheme, such an
oct-tree is used to compute short-range forces, and is complemented by a
Particle-Mesh (PM) code on large scales. The resulting algorithm is commonly
called TreePM, and it is the standard gravity solver for Gadget-4.

As we will show in next Section, deviations from a pure Newtonian approach
become significant on scales that are comparable with the Hubble horizon, so a
Newtonian treatment of small-scale clustering, performed by the Tree
algorithm, would introduce a negligible error if large scales are treated by a
fully relativistic gravity solver. This can be achieved, in a TreePM scheme,
by using a relativistic PM code for large-scale gravity, where relativistic
potentials are sampled on a small enough mesh so as to be effectively
Newtonian on the scales where the Tree code gets in.

In this paper we present an implementation of Gadget-4 that uses a PM library,
based on Gevolution relativistic code, as the PM part of the TreePM solver.
This is a step toward the construction of an ecosystem of codes and post-
processing tools to perform end-to-end simulations of future surveys, with the
aim of achieving optimal control of all systematics, including theoretical
ones. The paper is organized as follows: Section 2 gives an overview of the
theory of relativistic perturbations, with a focus on the approach used in
Gevolution. Section 3 gives a description of the Gadget-4 and Gevolution
codes, and describes the implementation of Libgevolution and GrGadget. Section
4 presents the tests performed to validate GrGadget, while Section 5 gives our
conclusions.

##  2 Theory of relativistic perturbations

The success of Newtonian simulations in describing the large-scale structure
of the universe follows from the fact that, for an observer at rest with
respect to the CMB, the metric of spacetime is very close to Friedmann-
Lemaitre-Robertson-Walker’s (FLRW). Deviations from the Newtonian approach are
expected to be significant, albeit small, on scales near the Hubble horizon,
or when the energy-momentum tensor has relativistic components like radiation
or fast massive neutrinos. Deviations from FLRW metric are expected to be
strong in the proximity of compact objects, but this happens on scales that
are far smaller than the resolution that can be afforded in simulations of
large comoving volumes. It is thus fair to assume that the perturbations to
the metric are small and can be described in a weak-field regime. This does
not imply that deviations of the components of the energy-momentum tensor from
homogeneity are assumed to be small, density perturbations can be highly non-
linear: what we require is that the size of self-gravitating objects is much
larger than their gravitational radius.

The Gevolution code (see Adamek et al., 2016) models the spacetime metric with
a perturbed FLRW metric in the weak field regime. In the _Poisson gauge_ the
metric can be written as:

$\begin{split}ds^{2}=&a^{2}\Big{(}-c^{2}\,d\tau^{2}(1+2\Psi)-2c\,d\tau
dx^{i}B_{i}+\\\
&+dx^{i}dx^{j}\big{(}\gamma_{ij}(1-2\Phi)+h_{ij}\big{)}\Big{)},\end{split}$
(1)

where $a(\tau)$ is the scale factor of the FLRW background, $\tau$ is the
conformal time and $x^{i}$ are the space coordinates. It is possible to
exploit the residual degrees of freedom of the metric to impose the conditions
$B_{i|}{}^{i}=0$, $h_{i}{}^{i}=0$ and $h_{ij|}{}^{j}=0$. In our notation,
repeated latin indexes denote Einstein’s summation over the spatial
coordinates $1,2,3$ and the vertical bar subscript, e.g. $B_{i|j}$, denotes a
covariant derivative with respect to the affine connection that emerges from
the background spatial metric $\gamma_{ij}$.

The choice of the Poisson gauge is convenient because the two potentials
$\Psi$ and $\Phi$ are the gauge-invariant Bardeen potentials, and in the
Newtonian limit the the field $\Psi$ can be interpreted as the gravitational
potential. In other words, this is the gauge in which the standard N-body
solver is integrating the right equations of motion in the Newtonian limit
(Chisari & Zaldarriaga, 2011).

###  2.1 Field equations

The background, characterized by $a(\tau)$, is by construction a solution of
the Einstein’s equations in the presence of a homogeneous and isotropic
energy-momentum tensor $\bar{T}^{\mu}{}_{\nu}$:

$\bar{G}^{\mu}{}_{\nu}=-\frac{8\pi G}{c^{4}}\bar{T}^{\mu}{}_{\nu}\,,$ (2)

where $\bar{G}^{\mu}{}_{\nu}$ is Einstein’s tensor constructed from the metric
(1) with the perturbations $\Psi,\Phi,B_{i},h_{ij}$ set to zero. Applying
equation (2) to the FLRW metric one obtains Friedmann’s equations.

To solve for the perturbations of the metric, the usual procedure consist in
subtracting (2) from the full Einstein’s equations:

$G^{\mu}{}_{\nu}-\bar{G}^{\mu}{}_{\nu}=-\frac{8\pi
G}{c^{4}}(T^{\mu}{}_{\nu}-\bar{T}^{\mu}{}_{\nu})\,.$ (3)

The right hand side now contains the perturbation of the energy-momentum
tensor due to inhomogeneities in mass and energy distributions, while the left
hand side is a very complicated non-linear expression containing the
potentials $\Psi,\Phi,B_{i},h_{ij}$ and their space-time derivatives up to
second order.

To reach a tractable set of equations that we can interpret and solve
numerically, we apply the weak field assumption. The perturbations
$\Psi,\Phi,B_{i},h_{ij}$ are assumed to be of order $\epsilon\ll 1$. Spatial
derivatives are known to increase their amplitude by a factor of
$\epsilon^{-1/2}$, accounting for the presence of shortwave fluctuations
induced by the non-linear structure in the energy-momentum tensor, while time
derivatives are assumed to preserve the perturbation order. Then one can
expand $G^{\mu}{}_{\nu}-\bar{G}^{\mu}{}_{\nu}$ in terms of the metric
perturbations, neglecting contributions with order higher than $\epsilon$. For
example: $\Phi$ is a term of order $\mathcal{O}(\epsilon)$, $\Phi_{,i}$ has
order $\mathcal{O}(\epsilon^{1/2})$, $\Phi_{|n}{}^{n}$ is a leading term
(order 1, because of the second derivative), quadratic terms like
$\Phi_{,n}\Phi{,}{}^{n}$ are $\mathcal{O}(\epsilon)$, and a term like
$\Phi_{,00}$ is considered as $\mathcal{O}(\epsilon)$. This type of expansion
is known as the _shortwave correction_ (Adamek et al., 2014).

Furthermore, experience has shown that the scalar perturbations $\Phi$ and
$\Psi$ are generally larger than the vector and tensor perturbations $B_{i}$
and $h_{ij}$. Indeed, the scalar potentials, that are sourced by the density
perturbation $\Delta T^{00}$, become the Newtonian potential in the Newtonian
limit, while the vector perturbation $B_{i}$ is sourced by $\Delta T^{0i}$,
that is small by a factor of $v/c$ for non-relativistic matter perturbations,
and $h_{ij}$ by $\Delta T^{ij}$, that is suppressed by a $(v/c)^{2}$ factor.
Hence, it is fair to drop quadratic terms of $B_{i}$ and $h_{ij}$ in this weak
field limit approximation.

In this approximation, from Eq. (3) it descends that its time-time component
yields a Poisson-like equation for the scalar $\Phi$:

$\begin{split}\Phi_{|n}{}^{n}(1+4\Phi)&-3\frac{\mathcal{H}}{c^{2}}\Phi_{,0}+3\frac{\mathcal{H}^{2}}{c^{2}}(\chi-\Phi)+\frac{3}{2}\Phi_{|n}\Phi_{|}{}^{n}\\\
&=\frac{4\pi Ga^{2}}{c^{4}}\Delta T^{0}{}_{0}\,,\end{split}$ (4)

where $\mathcal{H}=a^{-1}\frac{da}{d\tau}$ and $\chi=\Phi-\Psi$. From the
time-space section of eq. (3) we obtain:

$-\frac{B_{i|n}{}^{n}}{4c}-\frac{\Phi_{,i0}}{c^{2}}-\frac{\mathcal{H}}{c^{2}}(\Phi_{,i}-\chi_{,i})=-\frac{4\pi
Ga^{2}}{c^{4}}\Delta T^{0}{}_{i}\,,$ (5)

that, taking advantage of the condition $B_{n|}{}^{n}=0$, can be reduced to:

$-\frac{B_{i|n}{}^{n}}{4c}=-\frac{4\pi Ga^{2}}{c^{4}}P_{\perp}\Delta
T^{0}{}_{i}\,,$ (6)

where $P_{\perp}$ is a linear operator that selects from a vector field its
divergenceless component.

The traceless part of the spatial section of eq. 3 leads to:

$\begin{split}&\left(\delta^{j}{}_{b}\delta^{a}{}_{i}-\frac{1}{3}\delta^{a}{}_{b}\delta^{j}{}_{i}\right)\Bigg{[}\chi_{|j}{}^{i}-2\Phi_{|j}{}^{i}\chi+4\Phi\Phi_{|j}{}^{i}+2\Phi_{|j}\Phi_{|}{}^{i}\\\
&\qquad+\frac{1}{2c^{2}}h^{i}{}_{j,00}+\frac{\mathcal{H}}{c^{2}}h^{i}{}_{j,0}-\frac{1}{2}h^{i}{}_{j|n}{}^{n}\\\
&\qquad+\frac{1}{2c}\left(\frac{\partial}{\partial\tau}+2\mathcal{H}\right)\left(B^{i}{}_{|j}+B_{j|}{}^{i}\right)\Bigg{]}\\\
&=\left(\delta^{j}{}_{b}\delta^{a}{}_{i}-\frac{1}{3}\delta^{a}{}_{b}\delta^{j}{}_{i}\right)\left(-\frac{8\pi
G}{c^{4}}\Delta T^{i}{}_{j}\right)\,,\end{split}$ (7)

from which we can determine the rest of the metric degrees of freedom $\chi$
and $h_{ij}$. Since the source of $\chi$ and $h_{ij}$ are the perturbation of
the of the energy-momentum tensor $\Delta T^{i}{}_{j}$, their amplitude in a
matter dominated universe is suppressed by a factor $(v/c)^{2}$. That is
equivalent to say: since dark matter is non-relativistic, $\chi$ and $h_{ij}$
must be very small with respect to $\Phi$ or even $B_{i}$.

As a matter of fact, V1.2 of Gevolution implements an improved expansion of
the metric perturbations, that has been presented in Adamek et al. (2017). For
our tests we used the implementation of the original expansion, the one
presented above. However, the improved expansion has been ported to
Libgevolution and will be used when analysing result on the past light cone.
We do not expect the results presented in this paper to depend on the specific
expansion used.

###  2.2 Motion of particles along geodesics

Massive particles move along geodesics, whose equation can be expressed as:

$\begin{split}\frac{dx^{i}}{d\tau}=&\frac{cp^{i}}{\sqrt{(mca)^{2}+p^{2}}}+cB^{i}\\\
&+\frac{cp^{i}}{\sqrt{(mca)^{2}+p^{2}}}\left(\Psi+\Phi\frac{2(mac)^{2}+p^{2}}{(mac)^{2}+p^{2}}\right)\,,\end{split}$
(8)
$\begin{split}\frac{dp_{i}}{d\tau}=&-c\Big{(}p^{n}B_{n|i}+\Psi_{,i}\sqrt{(mca)^{2}+p^{2}}+\frac{p^{2}\Phi_{,i}}{\sqrt{(mca)^{2}+p^{2}}}\Big{)}\,,\end{split}$
(9)

where $p^{i}$ is the space part of the particle momentum and $p$ its norm. The
right hand side in the last equation is the generalized force acting on the
particles. The term proportional to $\Psi_{,i}$ becomes the Newtonian force in
the limit of small velocities, while $p^{n}B_{n|i}$ represent the corrections
due to _frame dragging_ and the third term in parenthesis is a further
relativistic correction.

The energy-momentum tensor is constructed from the knowledge of particle
positions and momenta, but its computation depends on the perturbed metric.
This means that, in Eqs. (4), (6), and (7), the source terms on the right hand
sides depend on the potentials themselves. These implicit equations may be
solved starting from the potentials at the previous time step and solving the
equations iteratively until convergence. The integration scheme that
Gevolution implements is simpler: at each time step the energy momentum tensor
is computed using the potentials from the previous step, then the Poisson
equations are solved to find the updated potentials, that will be used in the
next time step to compute the energy-momentum tensor.

The Newtonian limit is recovered when we consider Fourier modes larger than
$\mathcal{H}/c$ and we further neglect $B_{i}$ and consider $\Phi\ll 1$; then
equation (4) becomes:

$\begin{split}\Phi_{|n}{}^{n}=\frac{4\pi Ga^{2}}{c^{4}}\Delta
T^{0}{}_{0}\end{split}$ (10)

while (8) and (9) become:

$\begin{split}\frac{dx^{i}}{d\tau}=&\frac{p^{i}}{ma}\,,\end{split}$ (11)
$\begin{split}\frac{dp_{i}}{d\tau}=-\Phi_{,i}mc^{2}a\,.\end{split}$ (12)

##  3 Algorithms and code infrastructure

###  3.1 Gevolution

Gevolution888https://github.com/gevolution-code (Adamek et al., 2016) is an
N-body relativistic cosmological code, written in C++ and parallelized with
the MPI paradigm. The physical theory behind this code has been described at
length in Section 2. Numerically, this code implements a PM scheme to follow
the evolution of energy-momentum tensor perturbations. As in PM codes, the
advantage of working with a single grid and using Fast Fourier Transforms
(FFTs) to solve the Poisson-like equations for the fields is paid with a high
cost in memory, of $\mathcal{O}(N^{3})$ where $N$ is the number of grid points
per dimension.

Gevolution, can run in either _Newton_ or _General Relativity_ modes. The
Newtonian gravity solver inverts the Laplace operator in the Poisson equation
for the Newtonian potential, Eq. 10. When running the General Relativity mode,
the code solves Eqs. 4, 6 and 7, that require the computation of the perturbed
energy-momentum tensor. This is performed using a Cloud-In-Cell (CIC) scheme
both for the density and for particle velocities; details are given in the
presentation paper. Then the Hamiltonian forces to which particles are
subjected are computed from Eqs. 8 and 9.

Gevolution solves the field equations in Fourier space, using a C++ library
called LATfield2 to operate FFTs on classical fields in massively parallel
applications with distributed memory. LATfield2 provides a programming
interface to perform operations on the fields, either in their real or Fourier
space representations. This library implements FFTs of 3-dimensional fields
whose memory is distributed among parallel processes following a 2-dimensional
uniform decomposition of space, in which each process owns in memory a portion
of the grid with a _rod_ shape (Daverio et al., 2015). In this way LATfield2
overcomes the scaling limitations of a simpler 1-dimensional domain (_slab_)
decomposition provided by the mainstream FFTW3 library999http://fftw.org/.
FFTW3 is used, however, to compute 1D FFTs.

###  3.2 Gadget-4

Gadget-4101010https://wwwmpa.mpa-garching.mpg.de/gadget4 is a state-of-the-art
TreePM N-body hydrodynamical cosmological code written in C++ (see Springel et
al., 2021); it is massively parallelized in a distributed-memory paradigm
using MPI.

As in most N-body codes, gravity in Gadget-4 is represented in the Newtonian
limit, but the equations of motion are modified to take into account the
Universe expansion, obtained by integrating the Friedmann equations
separately. As mentioned above, this approach is consistent with General
Relativity in the Poisson gauge, and gives the leading-order term of weak
field expansion. This amounts to neglecting the metric degrees of freedom
$B_{i}$, $\chi$ and $h_{ij}$, and is valid on scales much smaller than the
Hubble horizon. In a typical configuration that is convenient for large
cosmological volumes, the code solves for the forces acting on each particle,
representing them as the sum of two contributions, one due to the interactions
with nearby particles, computed with a Tree algorithm, and one due to long-
range interactions, computed with a PM algorithm.111111The code can work in
other configurations (a non-cosmological volume, switching off the PM,
enhancing the Tree part using multipole expansion) that are however not
relevant for this paper.

The Tree algorithm works by partitioning the space into cubic cells, called
nodes; in turn, each node is recursively partitioned into 8 children nodes
down to a pre-determined maximum refinement level. A tree structure tracks the
list of particles that are located within each node. This structure is used to
speed up the computation of gravitational force on a particle: in a particle-
particle integration scheme, this force is computed by adding up a series of
$\vec{r}\,m/r^{3}$ terms, one for each particle pair, but we know that the
accuracy of force evaluation does not depend strongly on the small-scale
distribution of distant particles, so in the Tree scheme the evaluation of
gravity is performed by grouping particles that belong to the same node, under
the condition that the node subtends a given aperture angle $\theta$.
Particle-particle computation is then used only for the nearest neighbours.
This is equivalent to considering the leading order in a multipole expansion
of the gravity force from particles belonging to a distant cell. While the
construction of the Tree is expensive in terms of computing time, it allows to
achieve $\mathcal{O}(N_{p}\log N_{p})$ scaling for the force computation,
where $N_{p}$ is the total number of particles in the simulation. Thus the
Tree is able to compute with high accuracy the short wavelength modes of the
gravitational interaction, while keeping the computational time low for large
simulations. However, the Tree code is slow in integrating particle motions
near the initial conditions, when the departures from homogeneity are small.
This is why it is often coupled with a PM code to speed up the first time
steps of a cosmological box.

The PM algorithm represents gravity through the gravitational potential field
$\Phi$, evaluated on a Cartesian cubic mesh of fixed size. The potential is
found from the density field by solving the Poisson equation in Fourier space,
while the force is computed from the gradient of the potential, obtained with
a finite differences scheme. According to the Nyquist-Shannon theory, this
implies that the information handled by the PM is limited to the long modes,
up to the Nyquist frequency.

To combine the forces provided by the PM and Tree codes, the gravitational
potential is split into the sum of two fields:

$\Phi=\Phi^{(L)}+\Phi^{(S)}\,,$ (13)

where $\Phi^{(L)}$ represents long-range modes from the PM, and $\Phi^{(S)}$
represents short-range modes from the Tree. Written in Fourier space (tilde on
top of symbols denotes a Fourier transform), the Poisson equation reads:

$\tilde{\Phi}_{k}=-\frac{4\pi}{k^{2}}\tilde{\rho}_{k}\,,$ (14)

where $\rho$ denotes the mass density. We can split the density as a sum of
short-range and long-range terms, using Gaussian filters:

$\tilde{\Phi}_{k}=-\frac{4\pi}{k^{2}}\tilde{\rho}_{k}\left(1-\exp(-k^{2}r_{a}^{2})\right)-\frac{4\pi}{k^{2}}\tilde{\rho}_{k}\exp(-k^{2}r_{a}^{2})\,.$
(15)

The scale $r_{a}$ is the one at which we split long- and short-range modes. We
can obtain $\Phi^{(S)}$ by solving the modified Poisson equation for short
modes:

$\tilde{\Phi}^{(S)}_{k}=-\frac{4\pi}{k^{2}}\tilde{\rho}_{k}\left(1-\exp(-k^{2}r_{a}^{2})\right)\,,$
(16)

and $\Phi^{(L)}$ by solving the modified Poisson equation for long modes

$\tilde{\Phi}^{(L)}_{k}=-\frac{4\pi}{k^{2}}\tilde{\rho}_{k}\exp(-k^{2}r_{a}^{2})\,.$
(17)

The long-mode Poisson equation (17) is solved by the PM in Fourier space, so
the convolution with the kernel is a simple multiplication. The Tree on the
other hand works in real space, hence equation (16) has to be transformed;
this can be done analytically, yielding:

$\Phi^{(S)}(\vec{x})=-G\sum_{i}\frac{m_{i}}{|\vec{x}-\vec{r}_{i}|}\operatorname{erfc}\left(\frac{|\vec{x}-\vec{r}_{i}|}{2r_{a}^{2}}\right)\,.$
(18)

###  3.3 GrGadget

####  3.3.1 Libgevolution library

In order to have a relativistic PM code working in Gadget-4, we developed a
library that implements both the Newtonian and the relativistic PM algorithms
of the monolithic Gevolution code. This was done by forking the Gevolution
github repository into Libgevolution, a library that is publicly available on
github121212https://github.com/GrGadget/gevolution-1.2 under MIT license.

The rationale behind the development of Libgevolution is to encapsulate
Gevolution’s resources and methods into abstract objects. This yields several
benefits. Firstly, Gevolution maintenance is eased by the logical
modularization of the code, i.e. instead of a monolitic code with a unique
workflow we can divide Gevolution into components (C++ classes and/or
namespaces) with well defined purposes. Secondly, we are allowed to re-use
Gevolution components within other applications, such as we do within Gadget-4
in the present paper.

We give here an overview of the library; the precise signature of all the
defined functions, methods and data structures is described in the technical
documentation of the code. Libgevolution is based on three cornerstones: (i) a
particle container implemented through the class `Particles_gevolution`; (ii)
a PM data structure named `particle_mesh`, templated on the particle container
type, that can be used either as a `relativistic_pm` or a `newtonian_pm`;
(iii) an executable application that uses the previous components to produce
N-body simulations as the original code does. `particle_mesh` has to be
understood as a container that is aware of the parallelization of the tasks
and distribution of memory; it holds the gravitational fields and it allows
the user to compute the forces acting on the simulation particles. The user
interface declared in `particle_mesh` consists of the following functions:

  * •

sample(...), that builds the sources (density field or energy-momentum tensor)
by sampling particle properties in the mesh;

  * •

compute_potential(...), that solves Poisson equations to compute the potential
fields;

  * •

compute_forces(...), that computes the forces acting on particles.

`particle_mesh` is specialized to solve the Newtonian problem or the General
Relativistic problem using class inheritance; Figure 1 illustrates the class
hierarchy of Libgevolution’s `particle_mesh`. The expert user will be able to
specialize `particle_mesh` to his/her own needs, for example by deriving a PM
that solves a modified gravity problem.

particle_meshrelativistic_pmnewtonian_pm Figure 1: PM class hierarchy in
Libgevolution.

`newtonian_pm` is the specialization of `particle_mesh` that contains a real
`LATfield2::Field` scalar field $\Phi_{\text{Newton}}$ and its complex
`LATField2::Field` Fourier transform $\tilde{\Phi}_{\text{Newton}}$, plus a
`LATField2::PlanFFT` that connects $\Phi_{\text{Newton}}$ with
$\tilde{\Phi}_{\text{Newton}}$ through discrete Fourier transform.
`relativistic_pm` is the specialization of `particle_mesh` that contains the
above quoted degrees of freedom of the perturbed FLRW metric, $\Phi$, $B_{i}$
and $\chi$. These are represented as real `LATfield2::Field`, with complex
`LATField2::Field` counterparts to represent their Fourier transforms and a
`LATField2::PlanFFT` for each field.

As a first testing phase, we run Libgevolution, called with a simple wrapper,
and the native Gevolution code, applying them to the same set of initial
conditions, checking that the results were identical both in the Newtonian and
relativistic cases. Then we stripped down Gadget-4 by switching off the Tree
code, and compared its results to the Newtonian results of Libgevolution. It
is necessary that this comparison gives nearly identical results if we want
Libgevolution to substitute the native PM code of Gadget-4 without loss of
accuracy. To achieve a satisfactory match of the two PM codes we had to change
the Gevolution scheme in a few points.

We started from V1.2 of Gevolution, that implemented a first-order version of
finite differences instead of the fourth-order scheme of Gadget-4. This
resulted in a difference with Gadget-4 run on the same initial conditions, and
in a percent-level offset of the matter power spectrum on large scales at low
redshift. We upgraded the computation of spatial derivatives to fourth order,
in parallel with the Gevolution developers that had noticed the same problem;
our implementation is equivalent the most recent issue of Gevolution (used,
e.g., in Adamek et al., 2022). The upgrade is the following: let’s consider
the gravitational potential along one direction of the mesh, and let’s call
its values $\Phi_{i}$, where the index $i$ denotes its position along that
direction. Its first derivative is computed with finite differences at the
first order as:

$\frac{\partial\Phi_{i}}{\partial
x}=\frac{\Phi_{i+1}-\Phi_{i}}{h}+\mathcal{O}(h),$ (19)

where $h$ is the size of the mesh cell. Fourth-order Taylor expansion gives:

$\frac{\partial\Phi_{i}}{\partial
x}=8\frac{\Phi_{i+1}-\Phi_{i-1}}{12h}-\frac{\Phi_{i+2}-\Phi_{i-2}}{12h}+\mathcal{O}(h^{4})\,.$
(20)

This has a smaller error of order $\mathcal{O}(h^{4})$, so it achieves higher
precision than (19) with the little cost of knowing the potential value at the
second-nearest cell, that implies a negligible communication overhead.

Another improvement with respect to V1.2 of Gevolution, that follows an
implementation of Gadget-4, was the application of correcting filters to the
density in Fourier space to compensate for cloud-in-cell (CIC) interpolation.
Indeed, as discussed e.g. in Springel (2005) or Sefusatti et al. (2016), CIC
interpolation at some finite order leads to some loss of power that can be
compensated for in Fourier space using suitable kernels. This was applied both
to the computation of the density and to the computation of energy-momentum
tensor components in the relativistic case.

Lastly, to make the Newtonian PM scheme equivalent to that of Gadget-4 we
changed the form of the discrete Laplacian operator in the Poisson equation
solver from its original form

$\nabla^{2}\to-\frac{4N^{2}}{L^{2}}\Big{(}\sin^{2}\frac{\pi
k_{x}}{N}+\sin^{2}\frac{\pi k_{y}}{N}+\sin^{2}\frac{\pi k_{z}}{N}\Big{)}\,,$
(21)

described in Adamek et al. (2016), equation (C.5), to the form used in
Gadget-4:

$\nabla^{2}\to-\frac{4\pi^{2}}{L^{2}}\Big{(}k_{x}^{2}+k_{y}^{2}+k_{z}^{2}\Big{)}.$
(22)

####  3.3.2 Calling Libgevolution from Gadget-4

The implementation of Libgevolution in Gadget-4 was performed as follows. We
created a new PM class with a similar interface as the original one in
Gadget-4, so that it is initialized and executed with the same functions as
Gadget-4, i.e. init_periodic() and pmforce_periodic(). A new class
relativistic_pm was implemented within an gadget::gevolution_api namespace,
avoiding to use the wider gadget namespace to make a clear distinction of
purpose between the original Gadget-4 code and our additional features. This
relativistic_pm class acts much like a mediator taking information in and out
of gadget simulation particles, processing the correct units conversion and
calling the methods on gevolution namespace. Figure 2 shows a diagram that
summarizes the contents of this PM class, its relation with Gadget-4’s
resources and the entry points for gevolution’s api.

relativistic_pm consists of:

  * •

A variable of type `simparticle_handler` that acts as a wrapper for providing
particle information from Gadget-4’s simparticles global variable and writing
back the data produced by gevolution’s PM.

  * •

A variable of type `latfield_handler` that takes care of correctly
initializing LATfield global state. Indeed, while Gadget-4 can run with any
number of MPI processes, LATfield has limitations that depend on the number of
grid points in the PM. `latfield_handler` also takes care of creating a sub-
communicator from Gadget-4’s MPI global communicator that satisfies the
constraints set by LATfield.

  * •

A variable of type gevolution::cosmology that contains the parameters for the
background evolution.

  * •

A container of type gevolution::Particles_gevolution that holds particle
information, stored according to their location on the PM grid.

  * •

Variables of type gevolution::relativistic_pm and gevolution::newtonian_pm
that perform the actual PM computations, i.e. construct the sources, either
density or the components of the energy-momentum tensor, compute the
gravitational potential or the metric perturbation fields and the forces that
act upon the particles.

  * •

The methods `pm_init_periodic` and `pmforce_periodic`, for initialization and
execution of the PM, respectively.

Figure 2: Diagram of resource ownership and relations for Libgevolution
integrated into Gadget-4’s workflow. Each solid box represent a memory
resource (an instantiation of a variable type) while the dashed boxes indicate
ownership. The newly developed code, represented in the right part of the
diagram denoted with the namespace gadget::gevolution_api, consists in a class
named relativistic_pm that owns a particle_handler object that reads and
writes directly into gadget::simparticles, a latfield_handler that takes care
of setting up and inspect the state of LATfield2::parallel, and some types
defined in Libgevolution, that are defined in gevolution namespace, like
cosmology, Particles_gevolution and relativistic_pm. The methods
sim::begrun1() and sim::gravity_long_range_force() in gadget:: interact with
the relativistic_pm through their interface init_periodic() and
pmforce_periodic().

####  3.3.3 Kick and drift operators

In order to keep the Hamiltonian character of the equations of motion in
Gadget-4, we have to describe the state of each particle through its position
and momentum, not velocity. Following a leap-frog scheme, the momentum should
be updated with a _kick_ operation using the full relativistic Eqs. (8) and
(9). However, velocities in Gadget-4 are to be interpreted as momenta (per
unit mass) of non-relativistic particles in the Newtonian limit. Then we
redefine the Gadget-4 _kick_ and _drift_ operators assuming non-relativistic
matter, $p\ll mca$, and further neglecting the very small contribution coming
from $\chi$:

$\begin{split}\frac{dx^{i}}{d\tau}=&\frac{p^{i}}{ma}\left(1+3\Phi\right)+cB^{i}\,,\end{split}$
(23)

$\begin{split}\frac{dp_{i}}{d\tau}=&-cp^{n}B_{n|i}-\Phi_{,i}mc^{2}a\,.\end{split}$
(24)

The right hand side of (24) is what we call _force_.

####  3.3.4 Adding long-range and short-range forces

To combine the forces computed with the relativistic PM and Gadget-4’s
Newtonian Tree we have extended the idea of the TreePM coupling. From equation
(13) one obtains that the force acting on a particle in a TreePM scheme
consists of two terms:

$\vec{F}=S_{r_{a}}[\vec{F}^{\mathrm{Tree}}_{\mathrm{Newton}}]+L_{r_{a}}[\vec{F}^{\mathrm{PM}}_{\mathrm{Newton}}].$
(25)

The first term is the force computed using the Tree on which an exponential
high-pass filter $S_{r_{a}}$ is applied, leaving short-wavelength modes. The
second term corresponds to the PM force on which the complementary low-pass
filter $L_{r_{a}}$ is applied to leave long-wavelength modes. The symbols
$S_{a}$ and $L_{a}$ formally denote these linear operators:

$S_{r_{a}}[f](\vec{r})=\frac{1}{N}\sum_{\vec{k}}\tilde{f}_{\vec{k}}(1-\exp(-k^{2}{r_{a}}^{2}))\exp(-i\vec{k}\cdot\vec{r})\,,$
(26)

and

$L_{r_{a}}[f](\vec{r})=\frac{1}{N}\sum_{\vec{k}}\tilde{f}_{\vec{k}}\exp(-k^{2}{r_{a}}^{2})\exp(-i\vec{k}\cdot\vec{r})\,.$
(27)

The _grid smoothing scale_ $r_{a}$ scales with the PM mesh size, and its value
is optimized in Gadget-4, in a way that will be tested below, to minimize the
impact of the two different treatments of the gravitational force.

In order to account for the relativistic dynamics while preserving the match
between tree and PM contributions that is valid in the Newtonian case, we
choose the following strategy: Gadget-4 calls both newtonian_pm and
relativistic_pm, the Newtonian value of the force is added to the Tree force
as in a standard Newtonian simulation, while the difference between the
Newtonian and the relativistic forces is added on top as a correction, but
filtered on a different scale $r_{b}$, that we call _gr-smoothing scale_. Eq.
(25) then becomes:

$\vec{F}=S_{r_{a}}[\vec{F}^{\mathrm{Tree}}_{\mathrm{Newton}}]+L_{r_{a}}[\vec{F}^{\mathrm{PM}}_{\mathrm{Newton}}]+L_{r_{b}}[\vec{F}^{\mathrm{PM}}_{\mathrm{GR}}-\vec{F}^{\mathrm{PM}}_{\mathrm{Newton}}]\,.$
(28)

The case $r_{a}=r_{b}$ would correspond to simply adding the relativistic
force to the Tree:

$\vec{F}=S_{r_{a}}[\vec{F}^{\mathrm{Tree}}_{\mathrm{Newton}}]+L_{r_{b}}[\vec{F}^{\mathrm{PM}}_{\mathrm{GR}}]\,.$
(29)

However, while the size of $r_{a}$, that regulates the match between Newtonian
Tree and PM forces, is very well tested within Gadget-4, the optimal value of
$r_{b}$ is to be found; we will show in the next Section that using $r_{b}$
larger than $r_{a}$ allows us to achieve percent accuracy at small scales.

##  4 Validation

The GrGadget code has been validated by running it on a few realizations of
initial conditions, listed in table 1. These were generated with Gadget-4’s
ngenic code at $z=19$, starting from a linear power spectrum generated with
CAMB131313https://camb.info/ and with cosmological parameters consistent with
Planck 2018 result (Planck Collaboration et al., 2020):
$\Omega_{b}h^{2}=0.0223$, $\Omega_{c}h^{2}=0.120$,
$H_{0}=67.3\,$\text{\,}\mathrm{k}\mathrm{m}\,\mathrm{s}^{-1}\,\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$$,
$A_{s}=2.097\times 10^{-9}$ and $n_{s}=0.965$.

name |  $N_{p}$ (particles) |  $N$ (PM grid points) |  $L$ (box size)
---|---|---|---
N64 | $64^{3}$ | 64 | $1\,$\text{\,}\mathrm{G}\mathrm{p}\mathrm{c}$/h$
N256 | $256^{3}$ | 256 | $1\,$\text{\,}\mathrm{G}\mathrm{p}\mathrm{c}$/h$
high_res | $512^{3}$ | 512 | $500\,$\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$
Table 1: Cosmological simulation configurations used to validate GrGadget.

###  4.1 Gevolution and Gadget-4 original codes

As already discussed in Section 3.3.1, the newtonian_pm implementation in V1.2
of Gevolution computes the Newtonian forces differently from those obtained
with Gadget-4’s PM. Before implementing Libgevolution as the PM engine of
Gadget-4, we need to make the two algorithms work in the same way.

To this aim, we have run a set of simulations with the configuration N64
(described in table 1) with a small number of particles $N_{p}=64^{3}$ to be
able to compute forces using a straightforward particle-particle (PP) scheme,
that can be taken as the true force that we are trying to approximate. The
same initial conditions at $z=19$ have been fed to both Gadget-4 (with Tree
either on or switched off to have a pure PM run) and Gevolution (in Newtonian
mode) codes. At later times, $z=8$ and $z=0$, we have written snapshots of the
forces that the simulation particles experience, separating the PM and the
TreePM components; we have then compared those to the true Newtonian force
computed with the PP scheme. The data we have obtained are summarized in the
plots shown in figure 3. We have binned particles according to the value of
the true force, then for each bin we have computed the mean (colored lines)
and standard deviation (shaded regions) of the difference between the force
computed with approximate methods (PM or TreePM) and the true value. Forces
are given in Gadget-4’s default units, which is actually acceleration,
measured in units of $10H_{0}\
$\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$=h\,$\text{\,}\mathrm{k}\mathrm{m}^{2}\,\mathrm{s}^{-2}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-1}$$.
The green line shows the PM result using the original Gevolution code (the
true force is anyway computed with Gadget-4 and matched particle by particle)
while the red line is obtained from a pure PM using Gadget-4’s original code.
The black line gives the TreePM method precision, obtained using Gadget-4.

Figure 3: Difference of gravity force with respect to the true PP value,
binned according to the true force, for N64 initial conditions, at $z=8$ (left
panel) and $z=0$ (right panel). Lines represent the mean value of force
difference in the bins, with colours explained in the legend; the shaded
regions give the standard deviation of the corresponding force difference.

Looking at the red and green lines (and their shaded areas) we find two known
results. Firstly, the TreePM method produces far less bias and dispersion when
estimating forces; for instance, in the left panel of Fig. 3 the error is of
the order141414This quantification is in code units, we can take this value as
a reference for a high accuracy gravity solver. of
$0.1\,h\,$\text{\,}\mathrm{k}\mathrm{m}^{2}\,\mathrm{s}^{-2}\,\mathrm{k}\mathrm{p}\mathrm{c}^{-1}$$,
while in the right panel it is larger but barely visible when compared with
the other curves. Secondly, while the PM force has low bias but a much larger
variance than the TreePM one at high redshift, at low redshift, i.e. at higher
level of non-linearity, it underestimates the value of the Newtonian force as
its magnitude increases. This underestimation is due to the failure of PM in
resolving interaction at scales smaller than the grid resolution.

When comparing Gevolution PM and true forces, we notice an _S_ -shaped feature
in the plot, much more visible at high redshift. As anticipated in Section
3.3.1, this is mostly due to the first-order interpolation used to find the
gradient of the potential in the code version that we tested.

In Fig. 4 we show the matter power spectra151515In this paper all particle
power spectra were computed using PowerI4 code presented in Sefusatti et al.
(2016). Unless otherwise stated, all power spectra are computed up to the the
Nyquist frequency of the PM mesh. obtained at $z=0$ from a set of larger
simulations with the configuration N256 (see table 1). The red solid line
shows the result obtained with the original Gadget-4 code with its TreePM
method, while the red dotted line shows the results obtained by switching off
the Tree so that the forces are computed using the PM alone. The green lines
show results obtained with the latest develop version of Gevolution that
implements higher order schemes for finite differences; the dotted line gives
results obtained with GRADIENT_ORDER=1 and is identical to the result obtained
with V1.2 of Gevolution, the green solid line uses GRADIENT_ORDER=2, that
corresponds to a second-order scheme. These power spectra show that the matter
distribution in Gevolution using first-order gradients loses power in what
seems to be a uniform trend for large-scale modes. This is a behaviour which
is not inherent to the PM nature of the code, since that type of numerical
approximation should predict very well the linear evolution at large scales;
indeed, the higher-order scheme recovers power on large scales to sub-percent
accuracy. Conversely, Gadget-4’s PM and TreePM agree very well at wavenumbers
below $k\sim 0.1h/{\rm Mpc}$ scale,

The higher-order differentiation worsens the loss of power of Gevolution for
high values of $k$, that is not present in Gadget-4. This can be explained as
a consequence of the particle-to-mesh sampling and mesh-to-particle
interpolation described in section 3.3.1. As discussed there, Gadget-4’s PM
corrects for these effects, resulting in a power spectrum that degrades only
at very high values of $k$ as we approach the Nyquist frequency, while
producing a $\sim 2$ percent overcorrection at $k\sim 0.4\ h/{\rm Mpc}$.

After implementing the higher-order differentiation scheme, the correction for
the loss of power discussed above and the change in the discrete Laplacian
operator (Section 3.3.1), the results of native Gadget-4 and Libgevolution PMs
become indistinguishable.

Figure 4: Matter power spectrum of N256 cosmological simulations. The lower
panel shows residuals with respect to Gadget-4’s original code (in red), used
as baseline. The black line shows the linear power spectrum obtained with
CAMB. Red lines show results obtained with Gadget-4, with the Tree part on
(solid line) or switched off (dotted line). Green lines show results obtained
with Gevolution in Newtonian configuration, with finite differences at first
order (dotted line) or second order (solid line).

###  4.2 Newtonian forces

We have tested our implementation of the GrGadget code by running a standard
test in Gadget-4: we create an N-body configuration in which there is a single
massive particle in the entire simulation box, while other massless test
particles are placed at different distances from the first. In this setting
the exact value of the force on each particle is known, hence one can compare
the numerical results coming from the TreePM algorithm to the analytical
solution.

The results are shown in figure 5, where each dot represents a test particle.
The x-axis gives the distance to the massive particle that sources the
gravitational field, in units of the PM resolution ($L/N$), while the y-axis
gives the corresponding absolute value of the relative difference of the true
and estimated forces acting on the test particle. The red and blue lines
correspond to the mean value of force residuals, for particles binned into
distance bins; the red line denotes the statistics obtained from a simulation
using Gadget-4’s original TreePM implementation and the blue line was produced
using GrGadget, in this case with the Newtonian gravity engine.

This figure shows that the accuracy with which the TreePM code reproduces the
gravitational force is at worst at percent level on scales of a few mesh
cells, corresponding to the scale where the PM and Tree contributions are
matched, and gets very accurate in the limits where either the Tree (small
scales) or the PM (large scales) dominates.

Gadget-4’s and GrGadget’s Newtonian PMs show basically the same accuracy, even
though their PM implementations are very different.

Figure 5: Forces due to a point source: the points are test particles located
at different distances (in units of the mesh resolution $L/N$) from the source
and the lines represent the RMS of the difference between real and TreePM
forces in different distance bins. The red line corresponds to Gadget-4
original TreePM while the blue line was obtained with GrGadget in Newtonian
mode. As for the the grid smoothing scale, the default value was used:
$r_{a}=1.25L/N$. For this test we have used $N=256$ and
$L=1\,$\text{\,}\mathrm{G}\mathrm{p}\mathrm{c}$/h$.

In Fig. 6 we show the matter power spectra of a set of N256 simulations (see
table 1). In this case we are comparing the matter clustering of GrGadget, in
blue (with Newtonian forces for testing purposes), against Gadget-4, in red.
In agreement with the previous test of force differences, we find that both
codes produce the same matter power spectrum up to floating point errors. This
is verified both in the case of simulations computing forces using a pure PM
and in the case of TreePM.

Figure 6: Matter power spectrum of four simulations starting from the same
initial conditions high_res: blue lines give results for Gadget-4 original
code, red lines give results for GrGadget. In both cases dotted lines refer to
runs with PM-only, solid lines refer to runs with full TreePM.

###  4.3 Relativistic simulations with GrGadget.

We present here results obtained by running GrGadget with relativistic_pm,
comparing them with the corresponding relativistic version of Gevolution. We
expect that the power spectrum of the matter density displays some
relativistic features at large scales due to terms preceded by $\mathcal{H}$
in the field equation (4), while at small scales results should be compatible
with Gadget-4’s Newtonian simulations. However, the matter power spectrum
shown here is not an observable quantity, so this comparison is just meant to
give a first validation of the results. A more thorough comparison of
observables reconstructed on the past light cone will be presented in a future
paper.

Figure 7 shows the matter power spectra for a series of N256 simulations (see
table 1). In this case Gevolution and GrGadget are run in GR mode. The
parameter that regulates the scale of the relativistic correction (Eqs. 28 and
29) is set to $r_{b}=6\,L/N\approx
23$\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$, i.e. the relativistic
corrections of the PM method are smoothed at a distances below $6$ grid cells.
The plot shows that relativistic PM-only simulations, GrGadget (blue dotted
line) and Gevolution (green lines) are compatible on large scales
($k<0.03h/$\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$$) up to a small percent-
level difference that it is likely caused by the use of different orders for
finite difference gradient; indeed, going from first- to second-order
differences (from dotted to solid green line) the power spectrum gets nearer
to GrGadget’s fourth-order one. The plot also confirms that our combination of
Tree and PM forces in the relativistic weak field limit with GrGadget (blue
solid line) reproduces the Newtonian non-linear features to sub-percent level
at small scales, that is for
$k>0.1h/$\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$$; here Gadget-4 (red solid
line) is again our reference for the non-linear clustering.

Figure 7: Matter power spectrum of Gadget-4, Gevolution and GrGadget runs, the
last code being run in relativistic mode. The upper panel shows the absolute
value and the lower panel the relative difference with respect to Gadget-4’s
TreePM. The black line gives the linear matter power spectrum; red and blue
lines give Gadget-4 and GrGadget results, with full TreePM forces (solid
lines) or with the Tree switched off (dotted lines). Green lines give
Gevolution results, dotted line referring to first-order finite differences
(GRADIENT_ORDER=1) and solid line referring to second-order calculation
(GRADIENT_ORDER=2).

Being designed for the use of Fourier methods from the beginning,
Libgevolution offers an interface for the computation of the power spectrum of
the fields defined through the library’s interface. Thus we can also extract
and analyse the power spectra of the individual components of the metric
perturbations from the relativistic simulations. Figures 8 and 9 show the
power spectra of the relativistic potentials, $\Phi$, $B_{i}$ and $\chi$, for
a high resolution configuration high_res (see table 1). These plots show a
comparison of PM (blue lines) and TreePM (red lines) simulations. The power
spectrum of the gravitational potentials converge for both methods on large
scales. However, below $$1\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$ the PM-
only simulation loses power with respect to the TreePM one; the differences
can reach up to $40\%$ as we approach the Nyquist frequency. This pattern is
equally found for the scalar fields $\Phi$ and $\chi$, as well as for the
individual components of $B_{i}$.

The right plot in Fig. 8 helps to understand the reason behind this result.
Generally speaking, energy density, momentum density and their respective
density current (the components of the Energy-Momentum tensor) are sources of
the metric perturbations. Even though those quantities, as fields, are found
at discrete positions of space defined by the mesh, their values are computed
by sampling the energy and momentum carried by the particle distribution,
which contain information on the clustering due to the short range
interactions (through the Tree) that goes well below the mesh resolution
$L/N$. Therefore, TreePM simulations, having power on scales well smaller than
the PM mesh, give a better representation of the source of metric
perturbation, and thus allow to recover power at frequency modes right below
Nyquist. Fig. 8 highlights the particular case of $T^{0}{}_{0}$ (the matter
density) as a source for $\Phi$; by comparing $T^{0}_{0}$ with $k^{2}\Phi$, we
are verifying the Poisson equation
$k^{2}\tilde{\Phi}\approx\tilde{T}^{0}{}_{0}$ that is valid for wavelengths
below the Hubble horizon. This confirms that the presence of small-scale
clustering in the particle distribution propagates to the gravitational fields
up to the maximum resolution that the PM allows. The same thing is visible in
the vector modes $B_{i}$ and in $\chi$ (Figure 9), where we also notice a
small, few-percent mismatch on large scales. These fields are known to give
sub-percent effects on observables, so this difference, that is likely due to
some degree of numerical mode coupling, is non considered as a problem.

In figure 10 we show how the matter power spectrum obtained using GrGadget is
affected by the choice of the gr-smoothing scale parameter $r_{b}$. We have
used an N256 box configuration to perform this test, and tested values of
$r_{b}=1.5,3,6$ in units of
$L/N\approx$4\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$. We find that large-
scales power is independent of the value of $r_{b}$ parameter; structures one
scales below the PM resolution are resolved by the Tree algorithm, hence for
$k>k_{\mathrm{Nyquist}}$ there is a convergence of all simulations to a common
non-linear power spectrum tail. It is in the medium to small scales
$k_{\mathrm{Nyquist}}>k>$0.2\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$h$
that we notice differences in the power spectrum above the $\sim 1\%$ (dashed
grey line). For small values of $r_{b}$ ($\sim 1.5\,L/N$), we obtain
discrepancies in the power spectrum at
$k\sim$0.5\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$h$ that can be as large
as 5 percent and indicate the limitations of our force summation scheme, Eq.
(28). A value of $r_{b}=3\,L/N$ or possibly higher is needed to obtain a good
compatibility of GrGadget and Gadget-4 for all modes greater than
$$0.1\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$h$, where relativistic
features in the matter clustering is negligible.

The last test we present here regards the convergence of the numerical results
for increasing resolution. Figure 11 shows the matter power spectrum obtained
from running Gadget-4’s TreePM (red lines), GrGadget with PM-only (blue dotted
lines) and GrGadget with TreePM (blue continuous line). These various code
configurations were run with different combinations of the number of grid
points per dimension $N=256$, $N=512$ and box length $L=250$, 500, 1000,
$$2000\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$; the number of particles was
fixed as $N_{p}=N^{3}$. In all cases we have set the PM smoothing scale to
$r_{a}=1.5\,L/N$ and the gr-smoothing scale to $r_{b}=3\,L/N$. It can be
observed with the finest resolution, in the top plots, that there is a
matching between General Relativity and Newtonian dynamics in the small
scales. Then as the mesh size becomes coarser, in the middle plots, some
discrepancies in the power spectrum start to appear which become more evident
for even coarser meshes, in the bottom plots. This mismatch may be caused by
$r_{b}=3\,L/N$ moving towards larger scales, so that the assumption that PM
forces are Newtonian on the small scales breaks. Indeed, while with $L/N=1\
h^{-1}$ Mpc ($r_{b}=3\ h^{-1}$ Mpc) the scales where relativistic effects
become evident in the matter power spectrum and the scales where the pure PM
prediction starts to deviate from TreePM are well separated, for larger $L/N$
values the two scales get nearer, indicating that the assumption of pure
Newtonian forces on the mesh scale may not be very good. This conclusion is
apparently at variance with the discussion of Figure 10, where a larger value
of $r_{b}$ was preferred; however, that figure refers to $L/N=1$ and is shown
at $z=0.5$, where clustering is a bit weaker. swe thus recommend to work with
mesh sizes of $L/N\sim$1\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$.

Figure 8: In the left plot: power spectrum of the metric perturbation $\Phi$
in a high_res simulation obtained with GrGadget. In the right plot: power
spectrum of $k^{2}\Phi$ and $T^{0}{}_{0}$. For modes well below the Hubble
horizon and small perturbations it should be verified that
$k^{2}\tilde{\Phi}\approx\tilde{T}^{0}{}_{0}$.

Figure 9: In the left plot: power spectrum of the metric perturbation $B_{i}$
(the $x$ component) in a high_res simulation obtained with GrGadget. In the
right plot: power spectrum of $\chi$. Figure 10: Power spectrum of matter
density for Gadget-4 and GrGadget, on a N256 simulation configuration. The
upper panel shows the absolute value and the lower panel the relative
difference with respect to Gadget-4’s TreePM. Different shades of blue
indicate different values of the gr-smoothing scale parameter $r_{b}=1.5,3,6$
in units of $L/N$. The PM smoothing scale is $r_{a}=1.5\,L/N$. The power
spectra in this plot are computed beyond the Nyquist frequency to show the
convergence of the matter distribution correlations for distances below the
grid resolution, the Tree regime. Figure 11: Matter power spectrum from
cosmological simulations at $z=0$ using GrGadget (the blue lines) and compared
to Gadget-4 (the red line) at $z=0$. The dotted line is obtained with a
simulation in which only the PM is used to compute forces. The plots show the
relative difference with respect to the power spectrum obtained with Gadget-4.
The left column corresponds to simulations with $N=256$ grid points per
dimension while for the right column $N=512$. The boxsize changes along the
ranks so that for the top plots the resolution is the highest
$L/N\approx$1\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$, in the middle
$L/N\approx$2\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$ and the bottom plots
correspond to $L/N\approx$4\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$. In all
cases $r_{a}=1.5\,L/N$ and $r_{b}=3\,L/N$. The grey dashed line indicate a
$1\%$ error.

##  5 Conclusions

We have constructed a relativistic TreePM code, that we call GrGadget, where
the large-scale contribution to the gravitational force is computed using the
relativistic C++ PM library Libgevolution, based on Gevolution code, while
gravity coming from small scales is computed by the Tree code of Gadget-4. The
code works under the assumption that, in the context of cosmological
simulations, dark matter can be treated non-relativistically and then the
equations of motion of tracer particles tend to the Newtonian limit at scales
well below the Hubble horizon. Following the Gevolution approach, we use a
weak field approximation of GR, where the perturbations of the space-time
metric with respect to FLRW background are encoded as fields and simulated by
the PM. Comparing the matter power spectrum from GrGadget simulations with
that of original Gadget-4 and Gevolution codes, we conclude that the code
produces consistent results as long as the PM cell size $L/N$ is smaller than
$$2\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$ and the gr-smoothing parameter
is $r_{b}\approx 3\,L/N$.

With respect to the pure PM implementation of Gevolution, the predictive power
of GrGadget gives an improvement even on the scales sampled by the mesh. This
is due to the fact that the energy-momentum tensor, that sources the equations
of the fields that represent the perturbations of the metric, is computed from
a fully non-linear distribution of particles, with gravity being resolved down
to a much smaller softening length and not down to the mesh size. This may be
very useful, e.g., when assessing the possibility of detecting the frame-
dragging effect of a rotating dark-matter halo, if not of a spiral galaxy
(Bruni et al., 2014). Furthermore, this code is a development of the widely
used Gadget-4 code, and because the PM sector of the code is called only by
the computation of the gravity force, our code can be easily extended to
simulations of galaxies or galaxy clusters by switching on the hydrodynamics,
star formation and feedback sectors. All the physics described by these
sectors can safely be treated in the Newtownian limit; one should in principle
add thermal energy of gas particles to the energy-momentum tensor, but while
this extension is straightforward, it is likely to provide a negligible
contribution.

This is, for our group, a further step in the construction of an ecosystem of
simulation codes and post-processing tools for modeling the evolution of
structure in the Universe, with the aim of making predictions for precision
cosmology. Sub-percent accuracy in cosmological predictions, that matches the
smallness of the statistical error that will be obtained with forthcoming
galaxy surveys mentioned in the Introduction, can only be obtained taking into
account relativistic effect (e.g. Lepori et al., 2020), and we can foresee
that a self-consistent treatment of these effects (to within the required
accuracy) will soon become the standard in cosmological simulations. These
effects can also be added by post-processing Newtonian simulations, but a
validation of these procedures requires validation against a more self-
consistent approach. Conversely, a large community is developing Gevolution in
the direction of adding modifications of gravity, whose formulation is
typically worked out in a general relativistic context. This line of
development, coupled with a Newtonian treatment of modified gravity in the
Tree code, would be precious in the formulation of tests of gravity, because
relativistic effects may hide smoking-gun features of specific classes of
modified gravity theories.

##  Appendix A Code scaling

The code we presented in this work is the merging of two codes whose behaviour
in terms of run-time scaling is well-known and characterized; since we did not
modify the underlying algorithms, it is expected that the run-time scaling of
our code follows that of the parent codes.

However, the Libgevolution’s PM is obviously different from Gadget-4’s, and we
added the translation of particles data from the host code to the target
relativistic PM. Both this facts require that we establish the overall scaling
of GrGadget in its fully-relativistic configuration and the overhead
associated to both the relativistic PM and the interface between the two
codes.

In figure 12 we show the fraction of time spent in the PM in both the original
and relativistic configurations as a function of the grid cell size (see the
caption for details). The relativistic PM is an order of magnitude more
expensive than the original Gadget-4’s Newtonian PM, although in absolute
sense it is still either negligible or secondary in the simulation sets that
have been tested (it reaches a maximum value of $16\%$ at highest resolution,
i.e. in the $N=512$,$L=$250\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$).
However, it scales with both the resolution and the grid number as the
original Newtonian PM does.

Figures 13 and 14 report the scaling of run time in strong and weak scaling
tests respectively for the total run time, the tree time and the PM time
(left. middle and right panels in both figures; see the captions for details).
As inferred from 12, the run-time and hence its scaling, are dominated by the
Gadget-4’s Tree section.

Figure 12: The fraction of PM time to the total running time. Relativistic
runs are shown in blue while Newtonian runs are shown in red, whereas symbols
distinguish the value of grid points per dimension $N$ (squares and circles
for $N=256$ and $512$ respectively). We plot the time fraction on the $y$–axis
(log scale) vs the mesh resolution $L/N$ on the $x$–axis.

Figure 13: Strong-scaling test. We present the code scaling as the number $P$
of MPI tasks is increased while running the same simulation set-up. All the
results refer to GrGadget, i.e. to the configuration with fully-relativistic
PM. On the $x$–axis $P$ increases from 24 to 192, by $\times 2$ steps. On the
$y$–axis we report the speed-up (normalized so that the ideal speed-up for
$P=1$ is 1) for the total running time, the time spent in the PM and the time
spent in the Tree on the Left, Middle and Right panels respectively. Note that
the ideal behaviour (black dotted line) would result in a linear speed-up. The
PM data includes the translation of particles data from Gadget-4 to
Libgevolution. We show the results for $N=128$, $256$ and $512$ (solid, dashed
and dot–dashed lines respectively) for 4 different box sizes (i.e. mass
resolutions), $L=250$, $500$, $1000$ and
$$2000\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$ (circles, squares and stars
respectively). See the discussion in Appendix A for details.

Figure 14: Weak-scaling test. We present the code scaling as the number $P$ of
MPI tasks is increased for a proportionally increasing problem, then keeping
constant the particles–per–task occupancy. All the results refer to GrGadget,
i.e. to the configuration with fully-relativistic PM. On the $x$–axis $P$
increases from 24 to 192 with only 2 test cases. On the $y$–axis we report the
speed-up for the total running time, the time spent in the PM and the time
spent in the Tree on the Left, Middle and Right panels respectively. Note that
the ideal behaviour would result in a constant running time (horizontal dotted
black line). The PM data includes the translation of particles data from
Gadget-4 to Libgevolution. We show the results for two cases: from $N=128$, to
$N=256$ (solid lines with circles), and from $N=256$, to $N=512$ (dashed lines
with squares). Each of the two cases has been run for three different box
sizes (i.e. mass resolutions): $L=250\rightarrow
L=$500\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$, $L=500\rightarrow
L=$1000\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$ and $L=1000\rightarrow
L=$2000\text{\,}\mathrm{M}\mathrm{p}\mathrm{c}$/h$ (red, green and blue colors
respectively). See the discussion in Appendix A for details.

## Acknowledgements

We thank Julian Adamek for many fruitful discussions on gevolution, Volker
Springel for his comments on an early draft, Francesca Lepori, Marco Bruni,
Marco Baldi and Emilio Bellini for discussions. Simulations were performed
with the HOTCAT system of INAF (Taffoni et al., 2020; Bertocco et al., 2020).
PM acknowledges partial support by a Fondo di Ricerca di Ateneo grant of
University of Trieste.

## Data Availability

The simulation codes presented in this paper are publicly available on github
in the following path: https://github.com/GrGadget.

## References

  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 061102
  * Adamek et al. (2014) Adamek J., Durrer R., Kunz M., 2014, Classical and Quantum Gravity, 31, 234006
  * Adamek et al. (2016) Adamek J., Daverio D., Durrer R., Kunz M., 2016, Journal of Cosmology and Astroparticle Physics, 2016, 053
  * Adamek et al. (2017) Adamek J., Durrer R., Kunz M., 2017, J. Cosmology Astropart. Phys., 2017, 004
  * Adamek et al. (2022) Adamek J., et al., 2022, arXiv e-prints, p. arXiv:2211.12457
  * Alam et al. (2021) Alam S., et al., 2021, J. Cosmology Astropart. Phys., 2021, 050
  * Barnes & Hut (1986) Barnes J., Hut P., 1986, Nature, 324, 446
  * Barrera-Hinojosa & Li (2020) Barrera-Hinojosa C., Li B., 2020, Journal of Cosmology and Astroparticle Physics, 2020, 007
  * Barrera-Hinojosa et al. (2020) Barrera-Hinojosa C., Li B., Bruni M., hua He J., 2020, Vector modes in $\Lambda$CDM: the gravitomagnetic potential in dark matter haloes from relativistic $N$-body simulations (arXiv:2010.08257)
  * Bartelmann & Schneider (2001) Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291
  * Bertocco et al. (2020) Bertocco S., et al., 2020, in Pizzo R., Deul E. R., Mol J. D., de Plaa J., Verkouter H., eds, Astronomical Society of the Pacific Conference Series Vol. 527, Astronomical Data Analysis Software and Systems XXIX. p. 303 (arXiv:1912.05340)
  * Borzyszkowski et al. (2017) Borzyszkowski M., Bertacca D., Porciani C., 2017, MNRAS, 471, 3899
  * Bruni et al. (2014) Bruni M., Thomas D. B., Wands D., 2014, Phys. Rev. D, 89, 044010
  * Capozziello & De Laurentis (2012) Capozziello S., De Laurentis M., 2012, Annalen der Physik, 524, 545
  * Chisari & Zaldarriaga (2011) Chisari N. E., Zaldarriaga M., 2011, Physical Review D, 83
  * DESI Collaboration et al. (2016) DESI Collaboration et al., 2016, arXiv e-prints, p. arXiv:1611.00036
  * Daverio et al. (2015) Daverio D., Hindmarsh M., Bevis N., 2015, Latfield2: A c++ library for classical lattice field theory (arXiv:1508.05610)
  * Doré et al. (2014) Doré O., et al., 2014, arXiv e-prints, p. arXiv:1412.4872
  * Event Horizon Telescope Collaboration et al. (2019) Event Horizon Telescope Collaboration et al., 2019, ApJ, 875, L1
  * Ivezić et al. (2019) Ivezić Ž., et al., 2019, ApJ, 873, 111
  * Krause et al. (2017) Krause E., et al., 2017, arXiv e-prints, p. arXiv:1706.09359
  * Laureijs et al. (2011) Laureijs R., et al., 2011, arXiv e-prints, p. arXiv:1110.3193
  * Lepori et al. (2020) Lepori F., Adamek J., Durrer R., Clarkson C., Coates L., 2020, Monthly Notices of the Royal Astronomical Society, 497, 2078–2095
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Puchwein et al. (2013) Puchwein E., Baldi M., Springel V., 2013, MNRAS, 436, 348
  * Sefusatti et al. (2016) Sefusatti E., Crocce M., Scoccimarro R., Couchman H. M. P., 2016, mnras, 460, 3624
  * Silvestri & Trodden (2009) Silvestri A., Trodden M., 2009, Reports on Progress in Physics, 72, 096901
  * Spergel et al. (2015) Spergel D., et al., 2015, arXiv e-prints, p. arXiv:1503.03757
  * Springel (2005) Springel V., 2005, Monthly Notices of the Royal Astronomical Society, 364, 1105
  * Springel et al. (2021) Springel V., Pakmor R., Zier O., Reinecke M., 2021, MNRAS, 506, 2871
  * Taffoni et al. (2020) Taffoni G., Becciani U., Garilli B., Maggio G., Pasian F., Umana G., Smareglia R., Vitello F., 2020, in Pizzo R., Deul E. R., Mol J. D., de Plaa J., Verkouter H., eds, Astronomical Society of the Pacific Conference Series Vol. 527, Astronomical Data Analysis Software and Systems XXIX. p. 307 (arXiv:2002.01283)