# Cosmological N-body simulations including radiation perturbations

Jacob Brandbyge, Cornelius Rampf, Thomas Tram, Florent Leclercq, Christian Fidler, Steen Hannestad

► To cite this version:

Jacob Brandbyge, Cornelius Rampf, Thomas Tram, Florent Leclercq, Christian Fidler, et al.. Cosmological N-body simulations including radiation perturbations. Monthly Notices of the Royal Astronomical Society: Letters, 2017, 466 (1), pp.L68-L72. 10.1093/mnrasl/slw235. hal-05347503

HAL Id: hal-05347503

https://hal.science/hal-05347503v1

Submitted on 5 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Cosmological N-body simulations including radiation perturbations

Jacob Brandbyge,<fnref n="1" />,<fnref n="2" /> Cornelius Rampf,<fnref n="3" />,<fnref n="4" />★ Thomas Tram,<fnref n="5" /> Florent Leclercq,<fnref n="5" /> Christian Fidler<fnref n="6" /> and Steen Hannestad<fnref n="1" />

Accepted 2016 November 17. Received 2016 November 17; in original form 2016 October 19

## ABSTRACT

Cosmological N-body simulations are the standard tools to study the emergence of the observed large-scale structure of the Universe. Such simulations usually solve for the gravitational dynamics of matter within the Newtonian approximation, thus discarding general relativistic effects such as the coupling between matter and radiation (≡ photons and neutrinos). In this Letter, we investigate novel hybrid simulations that incorporate interactions between radiation and matter to the leading order in General Relativity, whilst evolving the matter dynamics in full non-linearity according to Newtonian theory. Our hybrid simulations come with a relativistic space–time and make it possible to investigate structure formation in a unified framework. In this work, we focus on simulations initialized at z = 99 and show that the extracted matter power spectrum receives up to 3 per cent corrections on very large scales through radiation. Our numerical findings compare favourably with linear analytical results from Fidler et al., from which we deduce that there cannot be any significant non-linear mode-coupling induced through linear radiation corrections.

**Key words:** cosmology: theory – dark matter – large-scale structure of Universe.

## 1 INTRODUCTION

According to the Λ cold dark matter (CDM) model, we live in a Universe that is nowadays dominated by a cosmological constant (Λ) and CDM, whereas the impact of baryons and radiation (i.e. photons and neutrinos) is secondary. None the less, at sufficiently early times, baryons and radiation played yet a major role in the early gravito-electroweak dynamics and were, for example, responsible for the observed acoustic oscillations in the cosmic microwave background (CMB; Smoot et al. 1992; Netterfield et al. 2002). Cosmological structure formation is mainly the result of gravitational instability, with initial conditions set in the period of recombination that is around 380 000 yr after the big bang. Electro-weak interactions freeze out at recombination and baryons are released from their tight coupling to radiation. Subsequently, this freeze-out in interactions enables the matter to cluster significantly for the first time.

At early times such as recombination, a Newtonian approximation of structure formation is not appropriate. Instead, the evolution of the multifluid components (dark matter, baryons, neutrinos and photons) is governed by the coupled set of Einstein–Boltzmann equations, i.e. General Relativity (GR). At early times, the linearized Einstein–Boltzmann equations are an excellent description, but cosmological structure formation for matter becomes fairly quickly a non-linear problem and unfortunately, solving the coupled set of Einstein–Boltzmann equations in full non-linearity is not yet feasible. Instead, it is common to solve for the process of cosmological structure formation using cosmological N-body simulations (Stadel 2001; Teyssier 2002; Springel 2005). Such simulations usually demand the Newtonian approximation and evolve only the matter component. Therefore, these simulations neglect the evolution of radiation perturbations.

Recently, the first general relativistic cosmological simulations have been performed, either by using an N-body approach in the weak-field approximation (Adamek, Durrer & Kunz 2014; Adamek et al. 2016) or by assuming the fluid approximation for a pure CDM component in full numerical relativity (Bentivegna & Bruni 2016; Giblin, Mertens & Starkman 2016; Mertens, Giblin &

<footnote n="1">¹Department of Physics and Astronomy, University of Aarhus, Ny Munkegade 120, DK-8000 Aarhus C, Denmark</footnote>

<footnote n="2">²Centre for Star and Planet Formation, Niels Bohr Institute & Natural History Museum of Denmark, University of Copenhagen, Øster Voldgade 5-7, DK-1350 Copenhagen, Denmark</footnote>

<footnote n="3">³Department of Physics, Israel Institute of Technology – Technion, Haifa 32000, Israel</footnote>

<footnote n="4">⁴Institute for Theoretical Physics (ITP), Philosophenweg 16, Heidelberg University, D-69120 Heidelberg, Germany</footnote>

<footnote n="5">⁵Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK</footnote>

<footnote n="6">⁶Catholic University of Louvain – Center for Cosmology, Particle Physics and Phenomenology (CP3), Chemin du Cyclotron 2, B-1348 Louvain-la-Neuve, Belgium</footnote>

<footnote>* E-mail: rampf@physics.technion.ac.il</footnote>
Starkman 2016). While the former also includes (massive) neutrinos, generally these simulations do not take into account the evolution of all multifluid species and furthermore, also rely on the validity regime of the weak-field and fluid approximation, respectively. Here we follow a different strategy that comes, of course, with different approximations and assumptions. We present a modification to the Newtonian N-body code GADGET-2 such that it solves for the full multispecies to first order in cosmological perturbation theory, whilst evolving the matter component in full non-linearity. Essentially, to achieve this modification, we use a first-order Einstein-Boltzmann solver to compute the force exerted on matter through radiation and update the momentum conservation of the matter in the simulation accordingly. By pairing the Einstein-Boltzmann solver with the N-body code, we obtain hybrid simulations with mutual benefits whilst minimizing the disadvantages from both Newtonian and relativistic worlds. First and foremost, cosmological N-body codes in the Newtonian approximation have been tested and improved over a long time period (e.g. Schneider et al. 2016), which gives us confidence that the non-linear matter is evolved to sufficient accuracy. Secondly, the N-body output of our relativistic simulations is in accordance with the first-order cosmological perturbation theory and thus comes with a solid approximation of the underlying relativistic space-time. This should serve as a solid basis for e.g. investigating ray tracing.

This Letter is organized as follows. In the following section, we report the relativistic fluid equations that we aim to solve numerically by modifying the existing N-body code GADGET-2. We pair the simulation with the Einstein-Boltzmann solver CLASS and call the resulting code COsmological Simulations with RAdiation (COSIRA). Details on the specific implementation can be found in Section 3. We discuss our numerical results and confront them with theoretical predictions in Section 4. Finally, we summarize and discuss our results in Section 5.

## 2 EQUATIONS OF MOTION IN GR

In a multifluid Universe, the relativistic equations of motion for CDM in an arbitrary gauge contain general relativistic corrections at multiple instances. Here we work in the N-body gauge that minimizes the appearance of such GR corrections. As was shown in Fidler et al. (2015), a fully relativistic analysis of the multifluid dynamics reveals that momentum and mass conservation for CDM are, to the leading order, given by

$$ \partial_{\eta} \mathbf{v}_{\text{cdm}} + \mathcal{H} \mathbf{v}_{\text{cdm}} = \nabla \Phi + \nabla \gamma, \quad (1a) $$

$$ \partial_{\eta} \delta_{\text{cdm}} + \nabla \cdot \mathbf{v}_{\text{cdm}} = 0, \quad (1b) $$

respectively. Here, $\eta$ is the conformal time, $\mathcal{H} = \dot{a}/a$ the conformal Hubble parameter with $a$ being the cosmic scalefactor, $\mathbf{v}_{\text{cdm}}$ and $\delta_{\text{cdm}}$ are, respectively, the (peculiar) velocity and the density contrast of the dark matter component; $\Phi$ is the cosmological potential satisfying the Poisson equation

$$ \nabla^2 \Phi = -4\pi G a^2 \sum_{\alpha} \bar{\rho}_{\alpha} \delta_{\alpha}, \quad (1c) $$

where $\bar{\rho}_{\alpha}$ is the background density of the species labelled with $\alpha$ and the summation on the right-hand side runs over all the relevant species in the Universe, i.e. dark matter, baryons, neutrinos and photons.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Ratio of $|\gamma|$ compared to the cosmological potential $\Phi$ in Fourier space, illustrating the impact of radiation contaminants on the momentum conservation equation (1a) of the dark matter component. On scales $k \ge 10^{-3}h\,\text{Mpc}^{-1}$, the relevance of residual radiation is continuously decreasing in time. Figure taken from Fidler et al. (2015).</figcaption></figure>

The above equations have been derived from GR in the N-body gauge with corresponding line element

$$ ds^2 = a^2 \left[ -(1 + 2A)d\eta^2 - 2v_i dx^i d\eta + (\delta_{ij} - 2D_{ij}H_T) dx^i dx^j \right], \quad (2) $$

where $D_{ij} \equiv \partial_i \partial_j - \delta_{ij} \nabla^2 / 3$, A is a perturbation in the time–time component of the metric (A is sourced by radiation pressure and anisotropic stress), $v_i = \partial_i v$ is the scalar part of the total velocity of all fluid components, where $v$ is the potential of the said scalar part, and $H_T$ is a perturbation in the trace-free part of the spatial metric. These linear perturbations are determined by GR and can be obtained from conventional Einstein-Boltzmann codes. Also, the relativistic correction

$$ \gamma \equiv \ddot{H}_T + \frac{\dot{a}}{a} \dot{H}_T - 8\pi G a^2 \sum_{\alpha} p_{\alpha} \Pi_{\alpha} \quad (3) $$

can be obtained from such codes. Here, $p_{\alpha}$ and $\Pi_{\alpha}$ are, respectively, the pressure and anisotropic stress of photons and neutrinos.

Equations (1a)–(1c) are, apart from the additional function $\gamma$ on the right-hand side in (1a) and the non-CDM source terms in (1c), identical with the Newtonian fluid equations. The function $\gamma$ is non-zero when there is a significant amount of radiation, whereas it is vanishing when radiation becomes negligible. Chronologically, this is precisely the case during cosmological structure formation. As can be seen from Fig. 1, at early times ($z \ge 50$), the $\gamma$ term yields corrections of up to 1 per cent–20 per cent to the momentum conservation of the CDM component, whereas at late times ($z < 50$), $\gamma$ becomes negligible on all scales. Furthermore, at earlier times when there is still significant radiation, the Poisson equation (1c) is sourced by non-vanishing density perturbations from all species, whereas at late times, it is effectively sourced only by matter perturbations. Summing up, equations (1a)–(1c) coincide with the Newtonian fluid equations to the leading order, but only at late times when $\gamma$ and the non-matter density perturbations are vanishing. Thus, only in that limit, conventional Newtonian N-body simulations solve for the dynamics of dark matter in accordance with GR. In the following section, we will introduce a modification to N-body simulations that makes it possible to run conventional simulations at early times, and in accordance with GR. Having implemented the relativistic corrections in the Newtonian simulation
the respective N-body output satisfies automatically first-order GR for all fluid species and thus comes with an approximation of the underlying space-time, here to be interpreted in the N-body gauge.

### 3 NUMERICAL IMPLEMENTATION

While conventional N-body simulations solve for the gravitational collapse of CDM in full non-linearity, they do not incorporate relativistic terms that stem from the presence of other fluid species. Here we introduce the code COSIRA. COSIRA makes it possible to perform hybrid simulations that are compatible with GR.

To achieve this, we separate the force term in the momentum conservation for CDM (see the right-hand side of equation 1a) into parts that are determined either in the N-body simulation or in the Einstein-Boltzmann code,

$$ \Phi + \gamma = \Phi_{\text{sim}} + \Phi_{\text{GR}}, \quad (4a) $$

where

$$ \nabla^2 \Phi_{\text{sim}} = -4\pi G a^2 \bar{\rho}_{\text{cdm}} \delta_{\text{sim}}, \quad (4b) $$

$$ \nabla^2 \Phi_{\text{GR}} = -4\pi G a^2 \bar{\rho}_{\text{cdm}} \delta_{\text{GR}}, \quad (4c) $$

where $\bar{\rho}_{\text{cdm}} \delta_{\text{GR}} = \sum_{\alpha \neq \text{cdm}} \bar{\rho}_{\alpha} \delta_{\alpha} - (4\pi G a^2)^{-1} \nabla^2 \gamma$ and $\delta_{\text{sim}}$ is the non-linear CDM density contrast obtained from the N-body simulation. The effective GR perturbations, $\delta_{\text{GR}}$, are determined by using a modified version of the Einstein-Boltzmann code CLASS (Blas, Lesgourgues & Tram 2011) and this information is fed to the Poisson solver in the N-body simulation at each time step. We note that in the present approach, for reasons of consistency, the information is only passed from the Einstein-Boltzmann system to the N-body code, but not vice versa. This means that we include the linear impact of the relativistic species on the non-linear evolution of CDM, but not the backreaction of the non-linear clustering of matter on the metric potentials or the relativistic species.

As evident from Fig. 1, $\gamma$ is rapidly oscillating in $k$ (over most scales) and therefore also oscillates in time. Since the time-scale of these oscillations is very short, of the order of $1/k$, it is not feasible to track them in an N-body simulation. However, as was shown in Weinberg (2002) and more recently in Voruz, Lesgourgues & Tram (2014), CDM does not couple to these fast modes and thus we consider only the slowly varying mode of $\gamma$ in our implementation. To extract this averaged mode of $\gamma$, we first apply a cut-off of $k \le 3 \cdot 10^{-3} h \text{ Mpc}^{-1}$ below which an average procedure is not required. For $k$-scales larger than the cut-off, we resample the transfer function of $\gamma$ such that it has 40 points per period $T = 2\pi/k$. Then we compute the rolling mean with a window size of 40 points (i.e. exactly one period). This rather complicated method is needed as, at early times, the oscillator is driven, implying that the frequency is not constant.

Our hybrid simulations have a box size of $L = 16384 \text{ Mpc} h^{-1}$ and a resolution of 1024 particles per dimension. The initial conditions for the CDM component are given in terms of the Zel'dovich approximation (Zel'dovich 1970; Buchert 1992), and are initialized with a weighted CDM and baryon power spectrum at $z_{\text{ini}} = 99$, which can be obtained from e.g. CAMB (Lewis, Challinor & Lasenby 2000). For the cosmological parameters, we use $\Omega_m = 0.3133$, $\Omega_\Lambda = 0.6867$, $\Omega_b = 0.0490$, $h = 0.6731$, $n_s = 0.9655$, $A_s = 2.215 \cdot 10^{-9}$, $\sigma_8 = 0.845$, and we consider a cosmology with three massless neutrino species.

In this work, the effective GR perturbations are not represented by discretized N-body particles but in terms of a fluid description. These GR perturbations are updated whenever a long-range force is

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Relative difference of the matter power spectrum $P(k)$ showing the impact of the relativistic perturbations $\delta_{\text{GR}}$ when switched on and off, i.e. $P(\delta_{\text{GR}} = 0) / P(\delta_{\text{GR}} \neq 0)$. All results are initialized at $z_{\text{ini}} = 99$. Shown are results from our numerical simulations ('sim'; solid lines) and linear theory predictions ('linear'; dotted lines) for the evolved differences of the power spectra at different redshifts. Linear predictions are obtained by applying the methods of Fidler et al. (2016).</figcaption></figure>

calculated in GADGET-2, and then added to the CDM long-range perturbations in Fourier space. Since there is no small-scale tree part for the radiation perturbations, they are, unlike the CDM perturbations, not smoothed with a Gaussian factor. For more details, we refer to Brandbyge & Hannestad (2009) where similar techniques have been applied to include linear perturbations in N-body simulations, however, in a Newtonian set-up with massive neutrinos.

Note also that for reasons of better comparison with theoretical predictions, we use the same initial matter power spectrum and set of phases, for both runs with or without the GR perturbations. Furthermore, to account for the radiation in the background evolution of the simulations, we added the background component of the radiation density to the Friedmann equation in GADGET-2 in both simulations with and without the GR perturbations.

CDM density fields are obtained by assigning N-body particles to grid points by applying the cloud-in-cell (CIC) algorithm (Hockney & Eastwood 1981). When plotting pure estimates of the CDM density (right-hand panel in Fig. 4), we utilize a $512^3$ grid so as to reduce shot noise. For estimating the matter power spectrum, we use a $1024^3$ grid and the CIC kernel is deconvolved in Fourier space.

### 4 RESULTS

Observe that the linear fluid equations (1a)-(1c) can be combined into a single differential equation for the matter density,

$$ \ddot{\delta}_{\text{cdm}} + \mathcal{H}\dot{\delta}_{\text{cdm}} - 4\pi G a^2 \bar{\rho}_{\text{cdm}} \delta_{\text{cdm}} = 4\pi G a^2 \bar{\rho}_{\text{cdm}} \delta_{\text{GR}}. \quad (5) $$

We use the novel semi-analytical methods of Fidler et al. (2016) to find the linear solution for $\delta_{\text{cdm}}$ and its power spectrum and confront it with the numerical solutions of the present approach. We note that for $\delta_{\text{GR}} = 0$, equation (5) is the usual differential equation for linear density fluctuations. The homogeneous solution of (5) is well known in the literature, so we will not report it here (see e.g. Rampf, Villone & Frisch 2015 and references therein).

In Fig. 2, we show the relative matter power spectrum when excluding/including the GR perturbations, for various redshifts. As evident, there is a striking agreement between the linear theory predictions (dotted lines) and our numerical results (solid lines). On
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Shown are the slices of the effective density perturbations $\delta_{GR}$ (see equation 4c) at various redshifts, compared against $\delta_{cdm}$ at $z = 99$. The CDM density does not change visibly at later times. The thickness of the slices is $512$ Mpc $h^{-1}$ and the width is $16384$ Mpc $h^{-1}$. For better visibility, we have increased the amplitudes of $\delta_{GR}(z)$ with respect to $\delta_{cdm}(z = 99)$ by factors of 1099, 15265 and 53936, respectively, for the redshifts $z = 99, 9$ and 0.</figcaption></figure>

the scales we simulate, the impact of GR perturbations is negligible for $k > 10^{-2} h \text{ Mpc}^{-1}$, but yields up to about 3 per cent corrections on larger scales. We note that at very large scales, close to the fundamental mode of our simulations, our numerical results depart from the linear predictions; this is due to sample variance because of the lack of N-body particles and thus could easily be rectified if needed.

In Fig. 3, we show the slices of the GR perturbations $\delta_{GR}$ at various redshifts, compared to $\delta_{cdm}$ at $z_{\text{ini}} = 99$. The relative importance of the radiation and CDM perturbations scales as $\Omega_{\text{rad}}/\Omega_{\text{cdm}}$ and at the largest scales simulated at $z = 99$, the radiation transfer function is roughly 4/3 times the CDM counterpart. At smaller scales, the effect of radiation free-streaming is clearly visible and becomes more pronounced at low redshift. Note that the largest structures in the box are not directly visible, since the density slices are dominated by scales around the peak of the power spectrum, $k \sim 0.01 h \text{ Mpc}^{-1}$. Since the radiation component has free-streamed out of these scales at $z = 99$, the radiation and CDM slices look quite different at this redshift. The CDM density, by contrast, follows a nearly scale-invariant growth on the scales shown and thus there is no visible evolution, which is the reason why we only show the slice at $z = 99$.

In Fig. 4, we show two panels, the left one depicts the difference of the Lagrangian potentials $\Delta\phi_{\text{cdm}} \equiv \phi_{\text{cdm}}(\delta_{\text{GR}} \neq 0) - \phi_{\text{cdm}}(\delta_{\text{GR}} = 0)$ and the right-hand panel the difference of Eulerian density contrasts, $\Delta\delta_{\text{cdm}} \equiv \delta_{\text{cdm}}(\delta_{\text{GR}} \neq 0) - \delta_{\text{cdm}}(\delta_{\text{GR}} = 0)$. Both panels are evaluated at final time $z = 0$. From these panels, one can recognize that large gradients in $\Delta\phi_{\text{cdm}}$ cause clusters to move and the change of the final positions of these clusters manifest as 'dipoles' in the Eulerian panel.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** The left-hand panel shows the difference $\Delta\phi_{\text{cdm}} \equiv \phi_{\text{cdm}}(\delta_{\text{GR}} \neq 0) - \phi_{\text{cdm}}(\delta_{\text{GR}} = 0)$ of snapshots obtained from our hybrid simulations, where $\phi_{\text{cdm}}$ is the potential of the scalar part of the Lagrangian displacement field ($\Psi = \nabla\phi_{\text{cdm}}$). The right-hand panel shows the difference of Eulerian density contrasts, $\Delta\delta_{\text{cdm}} \equiv \delta_{\text{cdm}}(\delta_{\text{GR}} \neq 0) - \delta_{\text{cdm}}(\delta_{\text{GR}} = 0)$. Both panels are evaluated at final time $z = 0$ and show a region of $2048$ Mpc $h^{-1}$ side length. Clusters move (in Eulerian coordinates) whenever $\Delta\phi_{\text{cdm}}$ exhibits high gradients. The two simulations predict the same structure in regions of constant $\Delta\phi_{\text{cdm}}$.</figcaption></figure>

In addition to the slight displacement of high-density regions in the presence of radiation/GR, we have investigated whether there are effects concerning the relative abundance of different cosmic web elements. One way to characterize such effects are the so-called mass filling fractions (MFFs) that measure how much of the overall mass content is contained in clusters, filaments, sheets and voids. (The MFF is defined such that the sum of the individual MFFs is unity.) Of course, the MFFs depend on the definition of such structures and in this work, we have chosen the Lagrangian classifier LICH (Leclercq et al. 2016). We find that radiation/GR has at most a $\mathcal{O}(10^{-5})$ relative effect on these numbers, leaving the MFFs effectively unchanged; in detail, they are MFF(CLUSTERS) $\simeq 0.097$, MFF(FILAMENTS) $\simeq 0.414$, MFF(SHEETS) $\simeq 0.389$ and MFF(VOIDS) $\simeq 0.100$. Thus, radiation/GR induces only a large-scale shift of cosmological structures (cf. Fig. 4), but does not change their relative abundance.

# 5 SUMMARY AND DISCUSSION

We have performed novel hybrid simulations to investigate cosmological structure formation for all fluid species within GR. By using our developed code COSIRA, we evolve photon and neutrino perturbations to first order within cosmological perturbation theory, whilst solving for the matter perturbations in full non-linearity in Newtonian gravity.

Relativistic computations within COSIRA are performed in the N-body gauge (Fidler et al. 2015), which is a gauge constructed such that at sufficiently late times, the equations of motion for CDM coincide exactly with the ones in Newtonian gravity. In that gauge, also the relativistic spatial volume element coincides at all times with the one in Newtonian gravity, which makes interpreting mass densities and velocites in the relativistic simulation a trivial exercise. We also note that at sufficiently late times, we have for the line element in the N-body gauge, equation (2), such that $A = 0$ and $H_T = 3\zeta$, where $\zeta$ is the gauge-invariant comoving curvature perturbation (see e.g. Fidler et al. 2016). This should enable straightforward post-analyses of the output of our hybrid simulations on a
fully relativistic space-time, e.g. to investigate weak lensing via ray
tracing techniques on the *N*-body output.

Our numerical code COSIRA delivers extremely stable results, es-
pecially on large scales. This is most clearly seen in Fig. 2 where
we compare our numerical findings against theoretical predictions.
The agreement between *N*-body simulation and linear theory re-
sults is excellent, except on very large scales, where our numerics
are mostly affected by the particle shot noise. It is because of this
agreement that the non-linear mode coupling induced through the
first-order GR corrections seems to be highly suppressed; however,
further studies with higher precision are required to assess this in
more detail.

In the present case study, we have focused on the regime of cosmological structure formation starting from z<sub>ini</sub> = 99 until today. COSIRA makes it however, possible to initialize relativistic simulations also at earlier times, possibly even as early as recombination. Furthermore, our findings enable us to perform comparisons to other relativistic simulations, where different approximations and assumptions have been applied, see e.g. Adamek et al. (2016), Bentivegna & Bruni (2016) and Hahn & Paranjape (2016). The limitation of the present approach is essentially by two effects that are of second order in cosmological perturbation theory. First, the radiation source is computed at the linear order using *CLASS*. Higher order corrections are expected to be very small, especially since the source is only relevant at early times. Secondly, our approach is based on the *N*-body gauge but, so far, this gauge is only defined at the linear order. Relativistic corrections to the CDM evolution are expected at second order in the metric potentials and are typically negligible. These corrections, however, could matter when studying observables that are suppressed at the linear order, such as the relativistic matter bispectrum (Tram et al. 2016).

Finally, we note that in this Letter, for reasons of better compar-
ison with the theoretical model, we have chosen to initialize our
simulations with the actual matter power spectrum at zini, for both
runs when including and excluding radiation. While this procedure
is perfectly justified in the case when radiation is included, it is
rather not when neglecting radiation perturbations. Indeed, in the
latter case, it is a very common practice to initialize N-body sim-
ulations by using today's matter power spectrum and rescale its
amplitude by the growing mode of linear density fluctuations back
to when the simulation is initialized (e.g. Schneider et al. 2016;
Zennaro et al. 2016). This method, sometimes referred to as 'back
scaling', produces a matter power spectrum at zini of a designed
universe with vanishing radiation content. Of course, the back scal-
ing method neglects by construction the evolution of radiation (see
Fig. 3). This and many more details should be assessed in forthcom-
ing studies and confronted with the present approach to sufficient
precision.

ACKNOWLEDGEMENTS

We thank Julian Adamek for useful discussions. JB and SH ac-
knowledge support from the Villum Foundation. CR acknowl-
edges the support of the individual fellowship RA 2523/1-2 from
the Deutsche Forschungsgemeinschaft (DFG). FL acknowledges

funding from the European Research Council through grant 614030,
Darksurvey. CF is supported by the Wallonia-Brussels Federation
grant ARC11/15-040 and the Belgian Federal Office for Science,
Technical & Cultural Affairs through the Interuniversity Attraction
Pole P7/37.

**Statement of contribution:** JB implemented the general relativistic modification into the *N*-body code *GADGET-2*. CR was involved in the conception and design of this project and wrote the Letter. TT modified the Einstein–Boltzmann code *CLASS* such that it outputs the relativistic correction function *γ* and other related quantities. FL produced the Eulerian and Lagrangian plots (Fig. 4) and performed the cosmic web analysis using *LICH*. CF and SH contributed to the interpretation of the results and the overall clarity of this Letter. All authors read and approved the final manuscript.

REFERENCES

Adamek J., Durrer R., Kunz M., 2014, Class. Quantum Gravity, 31, 23, 234006

Adamek J., Daverio D., Durrer R., Kunz M., 2016, Nat. Phys., 12, 346

Bentivegna E., Bruni M., 2016, Phys. Rev. Lett., 116, no. 25, 251302

Blas D., Lesgourgues J., Tram T., 2011, J. Cosmol. Astropart. Phys., 1107, 034

Brandbyge J., Hannestad S., 2009, J. Cosmol. Astropart. Phys., 0905, 002

Buchert T., 1992, MNRAS, 254, 729

Fidler C., Rampf C., Tram T., Crittenden R., Koyama K., Wands D., 2015, Phys. Rev. D, 92, 12, 123517

Fidler C., Tram T., Rampf C., Crittenden R., Koyama K., Wands D., 2016, J. Cosmol. Astropart. Phys., 1609, 031

Giblin J. T., Mertens J. B., Starkman G. D., 2016, Phys. Rev. Lett., 116, 251301

Hahn O., Paranjape A., 2016, Phys. Rev. D, 94, 083511

Hockney R. W., Eastwood J. W., 1981, Computer Simulation Using Particles. McGraw-Hill, New York

Leclercq F., Jasche J., Lavaux G., Wandelt B., 2016, preprint (arXiv:1601.00093)

Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473

Mertens J. B., Giblin J. T., Starkman G. D., 2016, Phys. Rev. D, 93, 124059

Netterfield C. B. et al., 2002, ApJ, 571, 604

Rampf C., Villone B., Frisch U., 2015, MNRAS, 452, 1421

Schneider A. et al., 2016, J. Cosmol. Astropart. Phys., 1604, 047

Smoot G. F. et al., 1992, ApJ, 396, L1

Springel V., 2005, MNRAS, 364, 1105

Stadel J. G., 2001, PhD thesis, Univ. Washington

Teyssier R., 2002, A&A, 385, 337

Tram T., Fidler C., Crittenden R., Koyama K., Pettinari G. W., Wands D., 2016, J. Cosmol. Astropart. Phys., 1605, 058

Voruz L., Lesgourgues J., Tram T., 2014, J. Cosmol. Astropart. Phys., 1403, 004

Weinberg S., 2002, ApJ, 581, 810

Zel'dovich Y. B., 1970, A&A, 5, 84

Zennaro M., Bel J., Villaescusa-Navarro F., Carbone C., Sefusatti E., Guzzo L., 2016, preprint (arXiv:1605.05283)

This paper has been typeset from a *TEX/LATex* file prepared by the author.