Non-linear evolution of f(R) cosmologies I: methodology

TL;DR

This work tackles non-linear structure formation in viable $f(R)$ gravity, focusing on the Hu–Sawicki model with chameleon screening. It develops a non-linear N-body framework that solves the coupled $f_R$ and gravitational potential $\phi$ equations using a Newton–Gauss–Seidel relaxation with multigrid acceleration, stabilizing the solver by rewriting $f_R$ as $f_R=\bar{f_R} e^u$. Validation via analytic tests, point-mass and pancake dynamics, and cosmological runs shows that gravity is enhanced relative to $\Lambda$CDM and that non-linear screening reduces this enhancement on small scales, with typical power-spectrum boosts around $20\%$ for strong HS fields. The results demonstrate a robust, scalable methodology for simulating non-linear $f(R)$ cosmologies and offer insights into observational signatures, with a companion paper providing a more detailed cosmological analysis.

Abstract

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 non-linear Poisson equation for the scalar degree of freedom introduced by the f(R) modifications. We detail the method to efficiently solve the non-linear Poisson equation by using a Newton-Gauss-Seidel relaxation scheme coupled with 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 ~20% for large but observationally allowed f(R) modifications. More detailed study of the non-linear f(R) effects on the power spectrum are presented in a companion paper.

Figures (11)

• Figure 1: The fractional discretization error ( left), shown as $f_{R, {\rm code}} / f_{R, {\rm exact}} - 1$, and the fractional residual error ( right), $\mathcal{L}(f_R)$, for the analytically tractable case of § \ref{['subsection:atc']}. The discretization error for the highest resolution box is $\sim 10^{-6}$, and allows us to define a sensible stopping criteria for the multigrid iterations. The fractional residual error is floored by limits on numerical precision, 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 = 10 h^{-1}$ Mpc. The background effective field mass is $0.13 h$ 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 in 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_{\rm cross} = 1$.
• Figure 5: Power spectra computed from cosmological simulations without $f(R)$ modifications. The Smith et. al. fit is plotted solid line, showing good agreement with our simulations. The vertical lines, from left to right, are the particle half-Nyquist scales, defined as $k_{\rm hN} = \pi N_{p} / (4 L_{\rm box})$, for $L_{\rm box}=256$, $128$, and $64$$h^{-1}$ Mpc boxes.