$\texttt{fRevolution}$ $-$ Relativistic Cosmological Simulations in $f(R)$ Gravity I: Methodology

TL;DR

The paper introduces $ exttt{fRevolution}$, a new relativistic simulation code for non-linear structure formation in $f(R)$ gravity, extending previous Newtonian approaches with a comprehensive relativistic framework.

Contribution

It presents the development of a relativistic cosmological simulation code for $f(R)$ gravity, incorporating full gravitational perturbations, which is a novel advancement over existing Newtonian codes.

Findings

Results for point mass fields in $f(R)$ gravity are demonstrated.

Cosmological simulations in the Hu-Sawicki model are performed and compared to Newtonian codes.

The code's capabilities are validated through these comparisons.

Abstract

We present the new relativistic cosmological particle-mesh code $\texttt{fRevolution}$, based on $\texttt{gevolution}$, aimed at simulating non-linear structure formation in $f(R)$ gravity. We introduce the general framework and approximation scheme, and the set of equations used to solve for the full set of gravitational perturbations. We show results for a point mass field and for cosmological simulations in the Hu-Sawicki model, and compare them to those of existing Newtonian codes. A more detailed analysis and discussion of our solutions will be carried out in a following paper.

Tables

Table 1. Table 1 : Orders of approximation used. The various fields are defined in § 1.1 .

Quantity	Order
$\Phi ,\dot{\Phi },\ddot{\Phi },\delta f_R,\delta {\dot{f}}_R,\chi ,\dot{\chi },\ddot{\chi }$	$ϵ$
${\Phi }_{,i},{\dot{\Phi }}_{,i},\delta f_{R,i},\delta {\dot{f}}_{R,i},{\chi }_{,i},{\dot{\chi }}_{,i}$	$\sqrt{ϵ}$
${\Phi }_{,ij},\delta f_{R,ij},{\chi }_{,ij}$	1
$B_i$, ${\dot{B}}_i$, ${\ddot{B}}_i$, $B_{i,j}$, ${\dot{B}}_{i,j}$, ${\ddot{B}}_{i,j}$	$ϵ$
$h_{ij}$, ${\dot{h}}_{ij}$, ${\ddot{h}}_{ij}$, $h_{ij,k}$, ${\dot{h}}_{ij,k}$, $h_{ij,k\ell }$	$ϵ$
$\delta T_0^0/{\overline{T}}_0^0$, $\delta R/\overline{R}$	1
$T_i^0/{\overline{T}}_0^0$	$\sqrt{ϵ}$
${\Pi }_{ij}/{\overline{T}}_0^0$	$ϵ$
$v^i$, $q_i$	1

Equations

$$\mathrm{d}{s}^2=a^2(\tau )\left[-(1+2\Psi )\mathrm{d}{\tau }^2-2B_i\mathrm{d}{x}^i\mathrm{d}\tau +(1-2\Phi ){\delta }_{ij}\mathrm{d}{x}^i\mathrm{d}{x}^j+h_{ij}\mathrm{d}{x}^i\mathrm{d}{x}^j\right].$$

$${\delta }^{ij}{\partial }_i{B}_j=0,{\delta }^{ij}{h}_{ij}={\delta }^{ij}{\partial }_i{h}_{jk}=0.$$

$$\chi \equiv \Phi -\Psi ,$$

$${\mathcal{B}}_i\equiv a^{-2}{B}_i.$$

$$S=\frac{1}{16\pi G}\int d^{4}x\sqrt{-g}F(R)\equiv \frac{1}{16\pi G}\int d^{4}x\sqrt{-g}[R+f(R)]+S_m[\psi ;g_{\mu \nu }].$$

$$(1+f_R)R_{\mu \nu }-\frac{R+f}{2}{g}_{\mu \nu }+{\square }_{\mu \nu }{f}_R=8\pi G{T}_{\mu \nu },$$

$$f(R)=-m^2\frac{c_1(R/m^2{)}^n}{1+c_2(R/m^2{)}^n},$$

$$f$$

$$R\approx 8\pi G{\rho }_m-2f\approx 8\pi G{\rho }_m+2m^2\frac{c_1}{c_2},$$

$$\frac{c_1}{c_2}=\frac{16\pi G{\rho }_0{\Omega }_{\Lambda }}{m^2}$$

$$f_{R_0}\equiv f_R(R_0)\approx -\frac{n{c}_1}{c_2^2}{\left(\frac{m^2}{R_0}\right)}^{n+1}$$

$$\bar{\rho }={\bar{T}}_0^0=\sum _i{\Omega }_i{a}^{-3(1+w_i)}\bar{P}{\delta }_k^j={\bar{T}}_k^j={\delta }_k^j\sum _i{\Omega }_i{w}_i{a}^{-3(1+w_i)}$$

$$w_i=\{\begin{array}{ll}0& \text{non relativistic matter: baryons, CDM}\\ 1/3& \text{relativistic species, radiation}\\ -1& \text{cosmological constant}\end{array}$$

$$(1+{\bar{f}}_R){\bar{R}}_{\mu \nu }+\frac{\bar{R}+\bar{f}}{2}{\bar{g}}_{\mu \nu }+{\bar{\square }}_{\mu \nu }{\bar{f}}_R=8\pi G{\bar{T}}_{\mu \nu },$$

$$3{\ddot{\bar{f}}}_R+6\mathcal{H}{\dot{\bar{f}}}_R+\left(2\bar{f}+\bar{R}-{\bar{f}}_R\bar{R}\right)a^2=-8\pi G{a}^2{\bar{T}}_{\mu }^{\mu },$$

$$3{\bar{f}}_{RR}\ddot{\bar{R}}+3{\bar{f}}_{RRR}{\dot{\bar{R}}}^2+6\mathcal{H}{\bar{f}}_{RR}\dot{\bar{R}}+(2\bar{f}+\bar{R}-{\bar{f}}_R\bar{R})a^2=-8\pi G{a}^2{\bar{T}}_{\mu }^{\mu }.$$

$$\bar{R}=-8\pi G\bar{T},$$

$${\mathcal{H}}_{\mathrm{in}}$$

$$R_{\mathrm{in}}$$

$${\dot{R}}_{\mathrm{in}}$$

$${\tau }_{\mathrm{osc}}^2\simeq a{\bar{f}}_{RR}^{-1}.$$

$$Q=\bar{Q}+\delta Q,$$

$$\delta f_R\equiv \delta f_R=f_R-{\bar{f}}_R.$$

$$(1+{\bar{f}}_R)\delta R_{\mu \nu }+{\bar{R}}_{\mu \nu }\delta f_R-\frac{\bar{R}+\bar{f}}{2}\delta g_{\mu \nu }-\frac{\delta R+\delta f}{2}({\bar{g}}_{\mu \nu }+\delta g_{\mu \nu })+$$

$$+(\delta {\square }_{\mu \nu })f_R+{\bar{\square }}_{\mu \nu }\delta f_R=8\pi G\delta T_{\mu \nu }$$

$$(1+2\Phi )\Delta \delta f_R-2\mathcal{H}\delta {\dot{f}}_R+\frac{a^2}{3}\left[f_R\delta R-2\delta f+\bar{R}\delta f_R\right]=\frac{a^2}{3}(\delta R+8\pi G\delta T),$$

$$\delta {\dot{f}}_R=\delta {\dot{f}}_R^{t-\frac{1}{2}}=\frac{\delta f_R^t-\delta f_R^{t-1}}{\mathrm{d}\tau }.$$

$$U[\delta f_R]=S(T_{\mu \nu },\delta f_R)\Rightarrow Y[\delta f_R,T_{\mu \nu }]\equiv U-S=0.$$

$$r^{(0)}\equiv Y[\delta f_R^{(0)}],$$

$$\delta f_R^{(n+1)}=\delta f_R^{(n)}+{\epsilon }^{(n)}\equiv \delta f_R^{(n)}-{\frac{Y}{\partial Y/\partial (\delta f_R)}|}_{\delta f_R=\delta f_R^{(n)}}.$$

$$||Y^{(n)}||<r_{\mathrm{c}},$$

$$||Y^{(n)}|{|}^2=\sum _{i,j,k}(Y_{i,j,k}^{(n)}{)}^2.$$

Full text

1 Introduction

The current and upcoming large scale structure surveys will be able to test cosmological structure formation with unprecedented precision. To properly understand and interpret this data, we therefore need to increase the precision of our theoretical prediction as well. The standard approach to investigate the non-linear process of structure formation is N-body simulations [3, 4, 5], which however ignore any relativistic effect. Nevertheless, recent works have shown that it is possible to overcome this intrinsic limitation by interpreting the predictions of a Newtonian simulation in a relativistic context [6, 7, 8, 9, 10, 11]. Such an approach is well defined in the weak field approximation [7, 12], and allows to rely on Newtonian N-body codes to predict the matter dynamics of a relativistic theory such as General Relativity (GR). It is worth to note that it is possible to include relativistic effect in the initial conditions [13, 14], include relativistic species [15, 16] and, finally, to reconstruct relativistic observables [11].

On the other hand, in the Newtonian framework, the scale factor is completely decoupled from the evolution of matter and therefore needs to be set by hand. Therefore when one wants to implement dynamical dark energy or modified gravity theories, even if such theory aims to modify the dynamic of the scale factor, the latter is de facto reconstructed and set by hand, or taken to be the one of the $\Lambda$–Cold Dark Matter concordance model ($\Lambda$CDM) [17, 18, 19, 20, 21, 22, 23, 24, 25]. While this is justified in first approximation for most viable alternative models, we cannot expect to be able to capture all the new dynamics following this approach, which thus might limit our ability to constrain these theories.

This motivates the creation of the code gevolution, based on GR and directly based on the weak field approximation [26, 27, 1]. In this paper, we propose a method to extend it to $f(R)$ gravity models and discuss its implementation in the new code fRevolution and the first results.

1.1 Metric

We consider the line element of a perturbed FLRW metric in the Poisson gauge

[TABLE]

where $a(\tau)$ is the cosmological scale factor, $\tau$ is the conformal time, and $x^{i}$ are the comoving Cartesian coordinates. As usual, Greek indices run through all spacetime dimensions, while Latin indices run only on space-like dimensions. An overdot will denote derivative with respect to $\tau$, a bar will denote background quantities and a tilde will denote Fourier transforms. The Hubble parameter is defined as $\mathcal{H}\equiv\dot{a}/a$.

The Poisson gauge corresponds to choosing $B_{i}$ divergenceless, and $h_{ij}$ traceless and divergenceless111We should point out that in the present version of the code we neglect the tensor perturbations $h_{ij}$ and their back-reaction on the particle evolution. Firstly, the effect of $h_{ij}$ is much smaller than that of $B_{i}$ which is itself much smaller than that of the standard scalar potentials $\Phi$ and $\Psi$. Secondly, resolving the full dynamics of massless tensor perturbations is computationally extremely expensive. One can still recover the approximate but pretty accurate configuration of $h_{ij}$ at each time by neglecting the time derivatives. See the original gevolution paper [1] for details on this point., that is

[TABLE]

We can maintain these even beyond the linear level provided that we remain in the weak field regime, where all the metric perturbations remain small $\ll 1$. Following the gevolution prescription, we also define the gravitational slip

[TABLE]

and where convenient we will use the auxiliary field

[TABLE]

1.2 $f(R)$ Gravity

Our goal is to study structure formation in $f(R)$ gravity. Among the possible modified gravity alternatives for cosmic acceleration [28, 29, 30], $f(R)$ gravity is one of the more popular and well-studied classes of theories. There is an extensive literature on $f(R)$ gravity and its cosmological implications, for a review see for instance [31, 32, 33] and references therein.

The action of the theory is222We choose to work in the Jordan frame, where matter is minimally coupled to gravity so that the geodesics are the same as in GR. Alternatively, one can perform a conformal transformation to the Einstein frame, where the gravitational action is Einstein-Hilbert $S\sim\int{\rm d}^{4}x\,\widetilde{R}(\widetilde{g}_{\mu\nu})$, but matter is coupled to a different metric than $\widetilde{g}_{\mu\nu}$, so that the additional complexity resides in computing non-standard geodesics instead of non-standard evolution for the metric perturbations. Both approaches are equally valid and must lead to the same observable predictions.

[TABLE]

where $f$ is a non-linear function of the Ricci scalar $R$. The field equations read

[TABLE]

where we denoted $\square_{\mu\nu}\equiv g_{\mu\nu}\square-\nabla_{\mu}\nabla_{\nu}$ for compactness.

Among the $f(R)$ models relevant for the cosmic acceleration, the Hu-Sawicki model [34] is likely the most studied and tested at the level of cosmological background, linear perturbations, and Solar System (post-Newtonian approximation) [35, 36, 37, 38, 39, 30], and at the non-linear level with the use of simulations [17, 18, 19, 21]. These studies include possible degeneracies with other effects produced, for instance, by massive neutrinos [40]. The model is given by

[TABLE]

where typically $m^{2}$ is of the order of the present curvature of the Universe. At large curvatures, this can be approximated by

[TABLE]

Moreover, in the limit $c_{1}/c_{2}^{2}\to 0$ at fixed $c_{1}/c_{2}$ the solutions are approximately [34]

[TABLE]

so choosing

[TABLE]

will produce the observed accelerated background expansion. In most situations, the single parameter

[TABLE]

where $R_{0}=8\pi G\rho_{0}(\Omega_{m}+4\Omega_{\Lambda})$ is the present cosmological curvature, is enough to essentially specify all the model dynamics, and it is generically found that for small enough values of $|f_{R_{0}}|\lesssim 10^{-5}$, the Hu-Sawicki model is in agreement with all existing observations333Possible issues related to the quasi-static approximation and to past singularities [41, 42, 43, 44, 45] suggest that we might be able to put even stronger constraints once these effects are taken into account; however, these issues are beyond the scope of this work and we will not consider them further.. Operationally, specifying $f_{R_{0}}$, $m^{2}$ and $n$ completely determines $c_{1}$ and $c_{2}$ once we impose the condition (1.10). In our simulations, we take $n=1$ and consider three values of $|f_{R_{0}}|=10^{-4},10^{-5},10^{-6}$ (see later for more details).

2 Cosmological Background

Because we work in the Jordan frame, the matter evolution of the different species remains the same as in GR, so that we only need the scale factor and the present (or initial) abundances in order to compute the background density and pressure:

[TABLE]

where as usual

[TABLE]

The background equations are given by

[TABLE]

and their trace is

[TABLE]

which we can rewrite as

[TABLE]

Although for $f(R)$ models of cosmic acceleration such as Hu-Sawicki the background is essentially that of $\Lambda$CDM at large curvatures, this is not necessarily true at very late times and surely not necessarily true for a generic $f(R)$ model. To facilitate extending our framework to other models and parameters, we decided to keep the background solver completely general. Moreover, because we are interested in looking for large scale relativistic effects, like back-reaction, we should make sure that the correct background is subtracted when calculating the perturbation equations lest we introduce unphysical homogeneous modes in the perturbations. Notice in particular that in general we will not have

[TABLE]

because of oscillatory solutions and/or of $\Lambda$-like components in the $f(R)$ solutions, which are not present explicitly in the matter energy-momentum tensor. In the code, we solve (2.5) using a Runge-Kutta-Fehlberg method [46] starting from the “GR” initial conditions

[TABLE]

deep in the matter-domination era, typically at redshift $z_{\rm in}=(1+a_{\rm in})^{-1}\sim 100$.

In order to resolve the curvature and Hubble oscillations, one must choose a time step smaller than the typical oscillation time, which is roughly

[TABLE]

3 Perturbations

In order to remove the background we expand each quantity $Q$ as

[TABLE]

keeping in mind that the background quantities of $f$ and its derivatives are computed at $R=\bar{R}$. For convenience, we define the scalaron field

[TABLE]

The perturbation equations then read

[TABLE]

So far these are fully general, and do not assume small perturbations. Naturally, we do expect some perturbations to be indeed small, which simplifies the problem greatly.

Matter and curvature perturbations, that is $\delta T^{0}_{0}/\bar{T}^{0}_{0}$ and $\delta R/\bar{R}$, and particle velocities are kept at all orders, which allows us to study non-linear structure formation and relativistic or quasi-relativistic particle motion, such as CDM that has undergone rare extreme accelerations, and also intrinsically relativistic species like massive neutrinos.

Metric perturbations are assumed very small and are normally kept at first order. However, we also keep terms quadratic in $\Phi$ provided that they contain two space derivatives (e.g. $\Phi\Delta\Phi$, etc.), which makes sense since the Poisson equation (not exact in GR nor $f(R)$, but qualitatively accurate nonetheless) dictates $\partial^{2}\Phi\propto\delta\rho$ and we are keeping $\delta\rho/\bar{\rho}$ at all orders.

Our approximation scheme is summarised in Table 1. We assume as usual that $\bar{f}_{R}\ll 1$ as well as $\delta f_{R}\ll 1$, although the latter can be of the same order of magnitude and even larger than $\bar{f}_{R}$, and could be a “large” perturbation in the same sense as $\Phi$. The scalaron appears in the equations in a way similar to that of $\Phi$, contributing to the fifth-force effects of $f(R)$, and also directly sources the gravitational slip $\chi$ (see §3.4.1), so we might risk losing potentially important features of the solutions by neglecting terms containing it.

Furthermore, we follow the standard approach and work in the quasi-static approximation, which consists in neglecting time derivatives. Therefore, we neglect $\delta\ddot{f}_{R}$, but unlike in other works which are based on the strictly Newtonian version of $f(R)$, we do keep the term proportional to $\mathcal{H}\,\delta\dot{f}_{R}$, as we do for $\Phi$ (see §3.2).

3.1 Scalaron Field: Trace Equation

We begin by considering the trace equation to update the scalaron field. Following the prescriptions in Tab.1, we obtain

[TABLE]

Keep in mind that $\delta R=\delta R(f_{R})$ and $\delta f=\delta f(f_{R})$. This equation is a convenient choice as the first equation to solve in the code (after updating the energy-momentum tensor) because metric perturbations (and in fact only $\Phi$) enter the equation only in a sub-leading term ($\Phi\Delta\delta f_{R}\ll\Delta\delta f_{R}$), so the error we make in using the old value of $\Phi$ will be negligible.

Note that the term containing $\delta\dot{f}_{R}$ is dealt with numerically by splitting it as

[TABLE]

The index $t$, though it has an obvious correspondence with the cosmological time, is simply a discrete index labelling the simulation steps.

3.1.1 Relaxation Solver

The difficulty in solving (3.4) comes from the fact that the relation between $\delta f_{R}$ and $\delta R$ is in general (highly) non-linear so we cannot use standard spectral methods (e.g. FFT), rather we have to rely on relaxation methods. Schematically, we start from an equation of the form

[TABLE]

In (3.6), $U$ is a non-linear differential operator acting on $\delta f_{R}$, and the source term $S$ can in principle contain terms depending on $\delta f_{R}$ as well. Starting from an initial guess $\delta f_{R}=\delta f_{R}^{(0)}$, associated with a residual $r^{(0)}$

[TABLE]

we implement a Newton-Raphson iterative method, defined by

[TABLE]

where the error $\varepsilon^{(n)}$ quantifies the correction between values of $\delta f_{R}$ at consecutive iterations. Convergence is reached comparing the residual of the equation with some (small) pre-determined constant $r_{\rm c}$:

[TABLE]

where typically $||\cdot||$ denotes $L_{2}$ norm taken over the whole grid:

[TABLE]

Notice that the index $n$ denotes a progression in the relaxation, not in cosmological time. To clarify, Eq. (3.8) defines a sequence which progressively approaches the solution of (3.6) at each time step. In this sense, an additional (fixed) index $t$ is implied in each quantity in (3.8).

There are many possible sweeping strategies, the most popular and one of the easier to parallelise is probably the red-black scheme, in which one solves the equation for cells of the same colour as in the colours of a chess board (straightforwardly generalised to 3 spatial dimensions), and then solves for the remaining half. The reason why this scheme is particularly useful for parallelisation is that when discretised on a lattice labelled by the indices $i,j,k$ and having cell size $\ell$, Eq. (3.4) only depends on the local value of $\delta f_{R}$ and of its Laplacian, which is computed from the nearest neighbours:

[TABLE]

Because the nearest neighbours of a black cell are red and vice versa, we can parallelise the update of all cells of one colour and afterwards update the remaining half.

Despite parallelising the relaxation, convergence often becomes increasingly slower as one approaches the exact solution. Formally, the issue is that the modes in the residual that have wavelengths longer than the grid size decrease more slowly than those with wavelengths comparable with the grid size. Therefore, especially for large grids, one often relies on multi-grid methods to speed up the convergence. Multi-grid algorithms speed up the convergence of these long wavelength modes by solving the equation on coarser grids (larger grid size), whose small-scale modes (in units of the grid size) correspond to larger scale modes in the finer grids.

For simplicity, we briefly illustrate the algorithm for two grids, but it can be easily generalised (see e.g. [47] for additional details). After a number of relaxation steps on the finer grid $\ell$, which produce a guess $\delta f_{R}^{\ell}$ with residual

[TABLE]

we move to the coarser grid having cell size $L$ (typically $L=2\ell$) using the restriction operator444The restriction/injection/fine-to-coarse and the prolongation/interpolation/coarse-to-fine operators define the mapping between fields on two grids with different cell sizes. We choose a tri-linear interpolation for the prolongation operator, and its adjoint or inverse (full-weighting) for the restriction. See e.g. [47] for details.

$\mathcal{R}_{\ell\to L}$ on the scalaron and on the residual:

[TABLE]

On the new grid, we perform additional relaxations steps solving the modified equation

[TABLE]

We then prolong the error $\varepsilon^{L}\equiv\delta f_{R,\rm new}^{L}-\delta f_{R,\rm old}^{L}$ from the coarser grid to the finer grid, and thus correct the guess for $\delta f_{R}$ on the latter:

[TABLE]

Additional relaxation steps are then performed on the finer grid, and if required the multigrid cycle can be repeated until the desired precision is achieved.

3.1.2 Change of variable

Depending on the specific $f(R)$ model, $\delta f_{R}$ might only have a finite range of “healthy” values555This is not a generic feature of $f(R)$ models, but it occurs in several models designed to produce cosmic acceleration, including the Hu-Sawicki model. Typically, this happens whenever $f(R)$ and/or the relation $R\leftrightarrow f_{R}$ are ill-defined for $R<0$, or when the unbounded interval $-\infty<R<\infty$ is mapped into a bounded interval for $f_{R}$. , and it may happen that the sequence (3.8) accidentally pushes $\delta f_{R}$ outside this range. For instance, in the Hu-Sawicki model [34], the relation $R(f_{R})$ is well-defined only for a definite sign of $f_{R}$, so clearly only a finite range of values of $\delta f_{R}\equiv f_{R}-\bar{f}_{R}$ is allowed. Where needed, as suggested in [17], we circumvent this problem by using the auxiliary variable $u$, defined as

[TABLE]

This is of course not the only possible choice, and in principle each model should be considered individually. Once we have a relation $R(u)$ that is well-defined on the whole real axis, we can convert $\delta f_{R}\to u$ on each lattice point, re-formulate the trace equation in terms of $u$, and apply a strategy analogous to (3.8), before converting back $u\to\delta f_{R}$. For the specific case (3.16), we have for example

[TABLE]

so expanding the right-hand side at first order in $\delta u$ (we drop the subscript $(n)$ for simplicity) we obtain

[TABLE]

In practice, we can skip converting between $\delta f_{R}$ and $u$ using the following sequence

[TABLE]

where for clarity we remind the reader that $\varepsilon^{(n)}$ is given by (3.8).

3.2 Gravitational Potential: $00$ Equation

Having solved the trace equation, we use the 00 equation to solve for $\Phi$, putting all terms containing the scalaron in the source term:

[TABLE]

where it is implied that $\Phi=\Phi_{t-1}$ and $\chi=\chi_{t-1}$ in the right-hand side. With this new source term, $\Phi_{t}$ is readily computed via

[TABLE]

which allows us to update both $\Phi_{t}$ and $\dot{\Phi}_{t-\frac{1}{2}}$ analogously to (3.5).

3.3 Vector Field (Elliptic Constraint): $0i$ Equation

The $0i$ equation

[TABLE]

can be used to evolve the vector mode $B_{i}$ via an elliptic constraint equation. Projecting (3.22) on the spin-1 component, using the operator

[TABLE]

we obtain

[TABLE]

3.4 Traceless $ij$ Equation

We finally consider the traceless part of the $ij$ equations, namely

[TABLE]

where

[TABLE]

and round brackets in indices denote symmetrisation:

[TABLE]

This equation will be used to evolve the gravitational slip $\chi$ (via its spin-0 projection), and possibly $B_{i}$ (via the spin-1 projection) through a parabolic equation.

As mentioned previously, we are neglecting the tensor perturbations $h_{ij}$. They would enter these equations through a term proportional to $\ddot{h}_{ij}+2\mathcal{H}\dot{h}_{ij}-\Delta h_{ij}$, unchanged from GR to $f(R)$, inside the square brackets.

We move all non-linear terms and terms containing $\Phi$ and $\delta f_{R}$ (already updated at this point of the cycle) to the right-hand side and project on the traceless part, obtaining

[TABLE]

We kept $\delta f_{R}$ separated from the rest of the source term in the right-hand side for reasons that will be clear shortly. Note also that

[TABLE]

so this is going to be the combination that is actually used to solve this equation. In Fourier space, we obtain

[TABLE]

3.4.1 Spin-0 Mode

We first calculate $\chi$ by projecting on the spin-0 part, using the projection operator

[TABLE]

which yields

[TABLE]

so we finally obtain

[TABLE]

Notably, we can avoid the computation of the scalaron-dependent terms in the source, and simply add $\delta f_{R}$ to the final result, which is why we chose to keep those terms explicit in (3.28).

Eq. (3.33) is essentially showing how the new scalar is a source of anisotropic stress in $f(R)$ gravity theories. In fact, if one assumes that $\Phi$ and $T_{\mu\nu}$ are essentially the same as in GR666This is obviously not a good quantitative approximation, but it helps in understanding the qualititative effect of modified gravity on $\chi$., then it is the difference $(\chi-\delta f_{R})$ that is roughly equal to $\chi_{\rm GR}$ (see also §4), or equivalently

[TABLE]

3.4.2 Spin-1 Mode

Projecting on the spin-1 component, using the projector (3.23) through the contraction

[TABLE]

where $P_{(1)}$ was defined in (3.23), finally yields, using the gauge condition $\widetilde{\mathcal{B}_{i}}k^{i}=0$,

[TABLE]

The field $\widetilde{\mathcal{B}}_{i}$ is then updated with a simple Euler criterion

[TABLE]

4 The Newtonian Limit

The Newtonian limit of (3.20) is essentially the $f(R)$ equivalent of the Poisson equation, which is

[TABLE]

where $\Phi_{N}$ is the Newonian potential. Furthermore, the geodesic motion of non-relativistic test particles can be approximated by

[TABLE]

Several fundamental assumptions are being made in the Newtonian limit, namely that $\Phi_{N}$ is small so that the first order terms suffice, and that we are in the deep sub-horizon, quasi-static regime so that

[TABLE]

We should also assume that $\delta P\ll\delta\rho$, as is the case if velocities are non-relativistic. Moreover, we are assuming that only the leading corrections to GR are relevant, which allows us to get rid of plenty of terms such as $f_{R}\ll 1$, $f_{RR}R\ll 1$ and so on. With these approximations, we find that the leading contribution to (3.20) is given by

[TABLE]

so that the correction to the field $\Phi$, assuming that $\delta\rho$ evolves practically as in GR, is roughly

[TABLE]

Moreover, as we have seen in §3.4.1, the gravitational slip $\chi$ (which vanishes identically in Newtonian gravity) is now sourced directly by the scalaron, so that

[TABLE]

These equations provide us with one the more intuitive ways to see how the additional scalar sources the gravitational potential and contributes as a fifth force, in fact the acceleration of a test particle will be

[TABLE]

Similarly, the Newtonian limit of the trace equation (3.4) reads

[TABLE]

which in combination with the previous results yields

[TABLE]

which are precisely the equations used in the pioneering [17].

In the Newtonian approximation, we replace (3.4) and (3.20) with (4.8) and (4.4), respectively. When computing the particle dynamics, we moreover neglect $B_{i}$ and $\chi$, as well as the standard relativistic corrections (typically of order $v^{2}/c^{2}\ll 1$). See the original gevolution paper [1] for details.

5 Results

In this section we will present some results from our simulations and their comparison with $\Lambda$CDM and existing modified gravity codes [21, 22]. A more detailed discussion will appear in a following paper [2].

5.1 Point Mass

As a first test, we consider the static field produced by a point mass located in the centre of a (256Mpc/$h$)3 cubic box (with periodic boundary conditions), solving on $128^{3}$ grid points (hence the spatial resolution is $\ell_{\rm cell}=2$Mpc/$h$), and compare these results with the analytical solutions obtained linearising (4.8):

[TABLE]

with a density field

[TABLE]

The formal solution for an actual point mass $\rho=m\,\delta^{(3)}(\mathbf{r})$ in an asymptotically flat Universe is trivially a Yukawa-like profile

[TABLE]

where

[TABLE]

The conversion between $m$ and (5.2) is given by

[TABLE]

where $V_{\rm cell}=\ell_{\rm cell}^{3}$ is the volume of a lattice cell. For $|f_{R_{0}}|=10^{-6}$, we replace the point mass with a Gaussian profile

[TABLE]

because $\delta f_{R}>|f_{R_{0}}|$ at small distances, and hence the linear approximation fails. In all cases, we expect solutions to deviate from the analytical result nearest to the overdensity and to the boundaries of the box, due to finite resolution and finite size effects, respectively. Results are shown in Fig. 1.

5.2 Cosmological Simulations

We performed three simulations of a (512 Mpc$/h)^{3}$ comoving box with $512^{3}$ grid points and the same number of CDM particles, for different values of $|f_{R_{0}}|=10^{-4},10^{-5},10^{-6}$. This allows us to study the limit in which the model reduces essentially to $\Lambda$CDM and when instead deviations become significant. The other cosmological parameters used in the simulations are [48]: $h=0.6736$, $\Omega_{b}h^{2}=0.02237$, $\Omega_{c}h^{2}=0.1200$, $n_{s}=0.9649$, $A_{s}=2.099\times 10^{-9}$ (at 0.05 Mpc*-1*). We should stress that these cosmological parameters are not the most up-to-date values available (see e.g. [48]) but were chosen for a direct comparison with the existing codes [21, 22]. We plan on using more recent values in an upcoming publication [2] in which we discuss our results in more detail.

5.2.1 Matter Power Spectrum

We present out results for the matter power spectrum in Fig. 2(a), where we overlay our power spectra to those of MG-gadget [22]. We observe excellent agreement with [22] at all scales except for $|f_{R_{0}}|=10^{-4}$ (where deviations from GR are most significant), where our solutions have a faster drop in the power excess around $k\gtrsim 1h/$Mpc. This is most likely due to the proximity to the Nyquist scale $k_{\rm Ny}$, and in general we do not expect our results to be accurate and competitive with those of MG-gadget so close to $k_{\rm Ny}$777Notice that the results of the original gevolution code also deviate from those of Gadget at small scales $k\simeq k_{\rm Ny}$, and that in particular less power is produced around those scales (see [1] for details)..

Altogether, we can still state that our results for $|f_{R_{0}}|=10^{-4}$ agree with those of existing codes very well, even in the non-linear regime, at scales larger than about a factor of 4-5 times the Nyquist scale; the agreement is even better and essentially perfect at any scale $k<k_{\rm Ny}$ for $|f_{R_{0}}|=10^{-5},10^{-6}$. Considering the intrinsic limitations of a fixed-grid approach compared to an adaptive mesh, we can consider this agreement very satisfactory.

5.2.2 Vector Modes

In figure 2(b) we present our results for the power spectrum of vector modes $B_{i}$. We also compare our results with those of [49], in which vector modes are computed in a Post-Friedmannian framework. For this reason, we used slightly different cosmological parameters (see [49]) and the values $|f_{R_{0}}|=10^{-4}$, $1.289\times 10^{-5}$ and $1.289\times 10^{-6}$.

We can see that the power excess compared to $\Lambda$CDM is larger by roughly 100%, 30% and a few percent for $|f_{R_{0}}|=10^{-4}$, $10^{-5}$ and $10^{-6}$ respectively, at $k\simeq 1h/$Mpc. We should stress again that this excess is not due to any additional term appearing directly in the evolution equation for $B_{i}$, but indirectly due to how the scalar potentials (and its gradient, which sources the vector modes) and the matter source are modified because of the $f(R)$ contributions. It is therefore remarkable that the power excess is even larger than for $\Phi$ which is sourced by $\delta f_{R}$ directly. The agreement with [49] is overall rather good but we do detect a slight excess of extra power. We plan to investigate this point further in a following publication.

5.2.3 Curvature

Next we present the solutions for the scalar curvature perturbations $\delta R$, in figure 3(a). The figure shows the ratio of power spectra instead of the relative power excess as in the previous cases, because the differences from the $\Lambda$CDM solution

[TABLE]

can be of several orders of magnitude and not at most of order unity as in the cases of matter and vector perturbations.

We see that the $\Lambda$CDM limit appears to be recovered in the appropriate limit as $|f_{R_{0}}|$ decreases, but we also notice that deviations are significant, especially at smaller scales, already at relatively small values of $|f_{R_{0}}|$ and rather strikingly for $|f_{R_{0}}|=10^{-4}$ at basically all scales. While deviations in $\delta R$ are not easily testable alone, as the main cosmological observable is the matter power spectrum, these results suggest that even in those cases in which deviations in $\delta\rho$ and $\Phi$ are relatively small, the linear expansion around some reference curvature888Note that $\widetilde{R}$ needs not be the cosmological background curvature $\bar{R}$, but can be any “sensible” choice, for instance the GR solution $\widetilde{R}=-8\pi GT$. $\widetilde{R}$ might give extremely inaccurate approximations to the real solution. For example, we can not assume

[TABLE]

if $|\delta R/\widetilde{R}|\sim 1$, and similarly for other derivatives. Moreover, this is typically exacerbated by the high non-linearity of the relation between $f_{R}$ and $R$ in $f(R)$ models relevant for cosmic acceleration

The conclusion is that one should be very careful when producing estimates for the effects of modified gravity in the approximation $R\simeq R^{\Lambda\text{CDM}}$, because this could be violated by many orders of magnitude even in those cases for which the gravitational potential and the matter power spectrum are not too different that in GR.

5.2.4 Gravitational Slip

We present our results for $\chi$ in figure 3(b). In this case we show the actual power spectra instead of the relative power enhancement. As we have seen in §3.4.1, the scalaron directly sources $\chi$ and so when $\delta f_{R}$ is much larger than $\chi$ would be in GR, it completely dominates and we have essentially $\chi\approx\delta f_{R}$. Because $\delta f_{R}$ can in principle be of the same order of magnitude as $|f_{R_{0}}|$ and even larger, we can easily see how $\chi$ can be many orders of magnitude larger than in $\Lambda$CDM.

Interestingly, this is the only quantity for which in the case $|f_{R_{0}}|=10^{-6}$ we do not recover $\Lambda$CDM plus very small corrections, but instead deviations remain large, of several orders of magnitude, even at large scales. The possibility of detecting signatures of modified gravity using the gravitational slip is an interesting topic (see e.g. the recent [50]), and we plan on discussing some of these possibilities in an upcoming work [2].

6 Conclusions

We have presented the framework and first results from the code fRevolution, based on the relativistic code gevolution [1]. We have discussed the approximation scheme which only relies on the weak field limit of GR with no further assumption on the smallness of density and scalar curvature perturbations, nor on the smallness of the scalaron $\delta f_{R}$ compared to $\bar{f}_{R}$. Moreover, we go beyond the Newtonian limit (see §4) and take into account the Hubble friction when solving for the scalaron dynamics.

Overall, our results agree very well with analytical predictions in the case of the field produced by a point mass, and with existing (Newtonian) modified gravity codes for the matter power spectrum. To our knowledge, we present for the first time direct results for the scalar curvature perturbations, which however can be computed even in a strictly Newtonian framework, and for the gravitational slip $\chi$ and frame dragging $B_{i}$ which instead are intrinsically relativistic effects and can be computed in Newtonian codes only a posteriori under the assumption that Newtonian and relativistic solutions for $\Phi$ and $\delta f_{R}$ are essentially the same.

For the chosen parameter values, we detect a power excess for vector modes of about a factor 2, and differences of several orders of magnitude for $\chi$ and $\delta R$ compared to the $\Lambda$CDM predictions. While the impact of modified structure formation on the gravitational slip might provide interesting new directions to test and constrain modified gravity, the observed deviations from GR in the solutions for $\delta R$ suggest that extra care should be taken, when formulating predictions in modified gravity based on the assumption that deviations from the GR curvature are small, namely that not only $|f_{R}|\ll 1$, but also that $|\delta R+8\pi G\delta T|\ll|8\pi G\delta T|$. We plan to carry out a more extensive discussion and analysis of our results in an upcoming publication [2].

Acknowledgments

The authors would like to thank J. Adamek, M. Kunz and I. Sawicki for useful comments and discussions during the development and testing of the code, and Á. de la Cruz-Dombriz and P. Dunsby for their support and contribution during the early phases of the project. LR would like to thank the DAMTP, University of Cambridge for its hospitality. Early tests of the code have been carried out on the clusters Zeus at the University of Cape Town and Baobab at the University of Geneva. The final simulations have been carried out on Koios at the Institute of Physics of the Czech Academy of Sciences in Prague. DD is supported by an Advanced Postdoc.Mobility grant of the Swiss National Science Foundation. LR is supported by ESIF and MEYS (Project CoGraDS-CZ.02.1.01/0.0/0.0/15_003/0000437).