Nonlinear evolution of $f(R)$ cosmologies. I. Methodology

We introduce the method and the implementation of a cosmological simulation of a class of metric-variation $f(R)$ models that accelerate the cosmological expansion without a cosmological constant and evade solar-system bounds of small-field deviations to general relativity. Such simulations are shown to reduce to solving a nonlinear Poisson equation for the scalar degree of freedom introduced by the $f(R)$ modifications. We detail the method to efficiently solve the nonlinear Poisson equation by using a Newton-Gauss-Seidel relaxation scheme coupled with the multigrid method to accelerate the convergence. The simulations are shown to satisfy tests comparing the simulated outcome to analytical solutions for simple situations, and the dynamics of the simulations are tested with orbital and Zeldovich collapse tests. Finally, we present several static and dynamical simulations using realistic cosmological parameters to highlight the differences between standard physics and $f(R)$ physics. In general, we find that the $f(R)$ modifications result in stronger gravitational attraction that enhances the dark matter power spectrum by $\sim 20\%$ for large but observationally allowed $f(R)$ modifications. A more detailed study of the nonlinear $f(R)$ effects on the power spectrum are presented in a companion paper.

Figure 1

The fractional discretization error (left), shown as $|f_{R,\mathrm{code}}/f_{R,\mathrm{exact}}-1|$, and the fractional residual error (right), $\mathcal{L}(f_R)$, for the analytically tractable case of Sec. 4a. The discretization error for the highest resolution box is $\sim {10}^{-6}$, and it allows us to define a sensible stopping criteria for the multigrid iterations. The fractional residual error is floored by limits on machine accuracy, which is $\sim {10}^{-15}$.

Figure 2

The computed and expected solutions for the point mass test. The analytic solution is arbitrarily normalized such that it and the numerical solutions match at scales where the simulation is expected to best match the analytic solution. This scale corresponds to $x=10h^{-1}\mathrm{Mpc}$. The background effective field mass is $0.13h{\mathrm{Mpc}}^{-1}$.

Figure 3

Orbit test. A point mass is placed at the center of a ${128}^3$ grid, and the trajectory of test particles are simulated without $f(R)$ modification. The inner particle is placed at two grid cells away from the center, while the outer particle is placed seven grid cells away. Both particles are given initial velocities such that they would be in circular orbit without $f(R)$ modifications.

Figure 4

Zeldovich 1D plane wave collapse with no $f(R)$ modifications. Particle momenta ($p$) is plotted against particle positions ($x$). The points represent the simulation output and the line represents the exact solution given by the Zeldovich approximation. The snapshot is taken at $a=a_{\mathrm{cross}}=1$.

Figure 5

Power spectra computed from cosmological simulations without $f(R)$ modifications. The Smith et al. fit is plotted in a solid line, showing good agreement with our simulations. The vertical lines, from left to right, are the particle half-Nyquist scales, defined as $k_{\mathrm{hN}}=\pi N_p/(4L_{\mathrm{box}})$, for $L_{\mathrm{box}}=256$, 128, and $64h^{-1}\mathrm{Mpc}$ boxes.

Figure 6

The ratio of the time derivative of the $f_R$ field to the spatial gradient. The 2-norm is defined simply as $\sqrt{⟨x^2⟩}$, where $x$ is either ${\partial }_t^2{f}_R$ or ${\nabla }^2{f}_R$. The factor of $c^2$ comes from the definition of the FRW metric.

Figure 7

The residual error, defined as the 2-norm of the left-hand side minus the right-hand side of Eq. (43) divided by the 2-norm of the left-hand side, is plotted as a function of CPU time spent in the iterative partial differential equations (PDE) solvers. The NGS scheme is only efficient at reducing low wavelength error modes and thus the convergence rate stalls as the large wavelength mode begins to dominate the error budget. The multigrid scheme can reduce error modes of all sizes equally well, and thus retains a good convergence rate throughout the iterations.

Figure 8

Two-body acceleration as a function of separation for several $f_R$ field strengths. The accelerations are shown as relative differences from the exact Newtonian two-body acceleration. The large scatter at $x\sim 2$ is because of truncation errors resulting from differencing the potential at scales comparable to the grid size. At large radii from the center, the periodic boundary condition invalidates the comparison of forces in the simulation to those due to isolated two-body interactions. At such radii, the test particles are attracted to the neighboring periodic images of the central point mass.

Figure 9

Two-body acceleration as a function of separation for several $f_R$ field strengths. The spherical blob has a constant density of $\delta \sim 500$ and a radius of 20 grid cells ($15.6h^{-1}\mathrm{Mpc}$). The increased scatter in the accelerations at the surface of the sphere are due to the fact that a smoothness of the spherical surface is limited by the finite grid size of the simulation. As in Fig. 8, the acceleration field is affected by the periodic boundary condition at large radius.

Figure 10

The relative differences in particle momenta in Zeldovich pancake simulations with varying $f_R$ field strengths. The simulation snapshot is taken at $a=a_{\mathrm{cross}}=1$.

Figure 11

Relative differences between the power spectra of $|f_{R0}|={10}^{-5}$ simulations and unmodified $\Lambda \mathrm{CDM}$ simulations. The linear prediction is taken from Hu and Sawicki [12] and assumes no chameleon behavior. The power spectra for each box size are only plotted up to their respective half-Nyquist wave numbers, which are $0.79h{\mathrm{Mpc}}^{-1}$, $1.57h{\mathrm{Mpc}}^{-1}$, and $3.14h{\mathrm{Mpc}}^{-1}$ for L256, L128, and L64 simulations, respectively.