Nbody simulations of γ gravity
^{1} Instituto de Física, Universidade Federal do Rio de Janeiro C P 68528, CEP 21941972, Rio de Janeiro, RJ, Brazil
email: vargas@if.ufrj.br; ioav@if.ufrj.br
^{2} Institute of Theoretical Astrophysics, University of Oslo, Postboks 1029, 0315 Oslo, Norway
email: d.f.mota@astro.uio.no
^{3} Astrophysics, University of Oxford, DWB, Keble Road, Oxford, OX1 3RH, UK
email: hans.a.winther@gmail.com
Received: 26 October 2015
Accepted: 11 January 2016
We have investigated structure formation in the γ gravity f(R) model with Nbody simulations. The γ gravity model is a proposal which, unlike other viable f(R) models, not only changes the gravitational dynamics, but can in principle also have signatures at the background level that are different from those obtained in ΛCDM (Cosmological constant, Cold Dark Matter). The aim of this paper is to study the nonlinear regime of the model in the case where, at late times, the background differs from ΛCDM. We quantify the signatures produced on the power spectrum, the halo mass function, and the density and velocity profiles. To appreciate the features of the model, we have compared it to ΛCDM and the HuSawicki f(R) models. For the considered set of parameters we find that the screening mechanism is ineffective, which gives rise to deviations in the halo mass function that disagree with observations. This does not rule out the model per se, but requires choices of parameters such that  f_{R0}  is much smaller, which would imply that its cosmic expansion history cannot be distinguished from ΛCDM at the background level.
Key words: gravitation / dark energy / galaxies: clusters: general / largescale structure of Universe / galaxies: halos
© ESO, 2016
1. Introduction
Since the discovery in 1998 (Riess et al. 1998; Perlmutter et al. 1999) that the Universe is speeding up instead of slowing down (as would be expected if gravity is always attractive), considerable effort has been devoted to understanding the physical mechanism behind this cosmic acceleration. The two main theoretical approaches considered in the literature to explain this phenomena are (1) to assume the existence of a new component with a sufficiently negative pressure (p< −ρ/ 3), generically denoted dark energy; and (2) to consider that general relativity has to be modified at large scales, or, more accurately, at low curvature (modified gravity). The simplest dark energy candidate is Einstein’s cosmological constant (Λ) with an equation of state w_{DE} ≡ p_{DE}/ρ_{DE} = −1. However, in spite of its very good accordance with current observations, Λ has some theoretical difficulties such as its tiny value as compared with theoretical predictions of the vacuum energy density, the cosmic coincidence problem, and related finetuning. This situation has motivated the search for alternatives like modifiedgravity theories. The simplest modifiedgravity candidates are the socalled f(R)theories, in which the Lagrangian density ℒ = R + f(R) is a nonlinear function of the Ricci scalar R.
As is well known, metric f(R)theories can be thought of as a special case of a scalartensor theory; a BransDicke model with a coupling constant ω_{BD} = 0. An accelerated expansion appears naturally in these theories. The very first inflationary model, proposed by Starobinsky more than three decades ago (Starobinsky 1980), is driven by a term of the type f(R) = αR^{2} (α> 0) and is still in excellent accordance with observations (Planck Collaboration XXII 2014). More recently, the idea of an acceleration driven by latetime curvature has also been explored in Capozziello et al. (2003) and Carroll et al. (2004). These authors considered a theory in which f(R) = −αR^{− n} (n> 0 and α> 0). However, these models do not have a regular matterdominated era and are incompatible with structure formation (Planck Collaboration XXII 2014).
To build a cosmologically viable f(R) theory, some stability conditions have to be satisfied (Pogosian & Silvestri 2008): (a) f_{RR} ≡ d^{2}f/ dR^{2}> 0 (no tachyons); (b) 1 + f_{R} ≡ 1 + df/ dR> 0 (the effective gravitational constant, (G_{eff} = G_{N}/ (1 + f_{R})) does not change sign (no ghosts)); (c) after inflation, lim_{R → ∞}f(R) /R = 0 and lim_{R → ∞}f_{R} = 0 (General Relativity is recovered at early times) ; (d)  f_{R}  is small at recent times, to satisfy solar system and galactic scale constraints. In addition to these conditions, there are some desirable characteristics that a viable cosmological model has to satisfy (Amendola et al. 2007). It should have a radiationdominated era at early times and a saddlepoint matterdominated phase followed by an accelerated expansion as a final attractor. By using the parameters and r≡−R(1 + f_{R}) / (R + f) , it can be shown that an early matterdominated epoch of the Universe can be achieved if and . Furthermore, a necessary condition for a latetime accelerated attractor is .
There are viable f(R) gravity theories that satisfy all the criteria mentioned above (Hu & Sawicki 2007; Starobinsky 2007; Appleby & Battye 2007; Cognola et al. 2008; Linder 2009; O’Dwyer et al. 2013). However, there is a generic difficulty from which all these “viable” f(R) theories (Thongkool et al. 2009) suffer: the curvature singularity in cosmic evolution at a finite redshift (Frolov 2008). It can be shown that this type of singularity problem can be cured, for instance, by adding a highcurvature term proportional to R^{2} (Appleby et al. 2010) to the density Lagrangian. Therefore, it is not possible to have cosmic acceleration with a totally consistent f(R) theory modifying gravity only at low curvatures. We remark that we do not address this problem here and only consider modifications at low curvatures.
We consider the specific case of a viable f(R) theory called γ gravity (O’Dwyer et al. 2013). Generically, in almost all viable f(R) theories, structure formation imposes such strong constraints on the parameters of the models that the effective equation of state parameter cannot be distinguished from that of a cosmological constant. In γ gravity the steep dependence on the Ricci scalar R facilitates the agreement with structure formation. O’Dwyer et al. (2013) showed that, in principle, the parameter that controls the steepness in γ gravity allows measurable deviations from ΛCDM (Cosmological constant, Cold Dark Matter) at both linear perturbation and background levels, while still compatible with both current observations. The main goal of this paper is to study the effects of γ gravity on the structure formation at nonlinear scales for choices of parameters where the model has observable signatures^{1} on the background expansion history of our Universe.
We go one step further and analyze the nonlinear evolution of structures computed from numerical simulations. The code that this paper is based on is a slight modification of ISIS (Llinares et al. 2014), which in turn is a modification of the RAMSES hydrodynamic Nbody code (Teyssier 2002). See also Zhao et al. (2011), Li et al. (2012a, 2011), Oyaizu et al. (2008), Mota et al. (2008), Boehmer et al. (2010), Zumalacarregui et al. (2013), Puchwein et al. (2013) for other codes that have implemented and performed simulations of f(R) gravity and Winther et al. (2015) for a recent code comparison between these codes.
We only consider the modifications made to implement the γ gravity field equations in the modified gravity part of the Nbody code. For more details on the implementation of the scalar fields and other technicalities we refer to Llinares et al. (2014). For this purpose, we focus on simple observables such as the matter power spectrum, halo mass function, density profiles, and velocity profiles to investigate modified gravity signatures that were previously studied in Li et al. (2012b, 2013), Hammami et al. (2015), Schmidt et al. (2009), Gronke et al. (2015), Lombriser et al. (2012, 2014), Winther et al. (2012), Terukina & Yamamoto (2012), Shi et al. (2015), Pujol et al. (2014), Hellwing et al. (2013), He et al. (2013, 2015), Schmidt (2010), Brax et al. (2013), Tessore et al. (2015), Gronke et al. (2015) that can also be observed (Berti et al. 2015; Zivick et al. 2015; Mak et al. 2012; Cai et al. 2015; Song et al. 2015; Wilcox et al. 2015; Jain & Khoury 2010; Schmidt 2010; Hellwing et al. 2014).
This paper is organized as follows: in Sect. 2 we revisit the γ gravity model and show the main properties of the background and linear perturbation evolution. Section 3 details the dynamics equations of scalaron field (f_{R}) and particle movement equations for γ gravity, which must be solved by our code during the simulations. The method for solving these equations is briefly explained in Sect. 4, and the code is tested in Sect. 5. Finally our results are shown in Sect. 7, and we conclude in Sect. 8.
2. γ gravity review
We investigate spatially flat cosmological models in the context of γ gravity (O’Dwyer et al. 2013), a viable f(R) theory defined by the following ansatz: (1)where dt is the incomplete Γfunction and α, n and R_{∗} are free positive constants. In reality, γ gravity can be thought of as a simple generalization of exponential gravity (Linder 2009) (2)obtained by fixing n = 1 in Eq. (1). We emphasize that γ gravity can satisfy all the stability and viability conditions. As discussed in O’Dwyer et al. (2013), for fixed n, there is a minimum value (α_{min}) of the parameter α such that for values α>α_{min} a latetime accelerated attractor is achieved. We consider this case throughout. From Eq. (1) we obtain the following derivatives: We note from Eq. (3) that with increasing n, the steepness of the f(R) function increases. Higher n means smaller  f_{R0}  , and the departures from GR will be smaller accordingly.
Although there is no cosmological constant, f(0) = 0, it follows from Eq. (1) that GR with Λ is recovered at high curvatures. Therefore, for R ≫ R_{∗} the models behave like ΛCDM. Since we are mainly interested in phenomena that occurred after the beginning of the matterdominated era, we neglect radiation and write the effective cosmological constant (the cosmological constant of the reference ΛCDM model) as (5)In the equation above, denotes the present value of the matter density parameter that a ΛCDM model would have if it had the same matter density today () as the modified gravity f(R) model. represents the Hubble constant in the reference ΛCDM model. Therefore, we have , where Ω_{m0} and H_{0} are the present value of the matter energy density parameter and Hubble parameter in the f(R) model, respectively. It is useful to rewrite R_{∗} as (6)where . To compute the background evolution, we start from the f(R) field equation for a FLRW metric (7)where ′ ≡ d / dy (y = lna), H ≡ ȧ/a is the Hubble parameter (a dot denotes the derivative with respect to cosmic time), which is related to R by (8)To solve these equations, we introduce the new variables With these definitions we obtain where R′ is given by Eq. (7). It is straightforward to verify that, as defined, x_{1} and x_{2} are always zero during the ΛCDM phase. We here focus on the three cases summarized in Table 1.
Furthermore, in terms of x_{1}(y) and x_{2}(y), the effective dark energy equation of state (w_{DE}) is given by,(13)For the considered models, the evolution of w_{DE} as a function of the redshift z is shown in Fig. 1.
Fig. 1 Effective equationofstate parameter w_{DE} as a function of redshift z for the parameters given in Table 1. The strongest deviation from − 1 is lower than 4%. 

Open with DEXTER 
Overview of the model parameters for γ gravity.
For a general f(R) model the differential equation for the matter density contrast (δ_{m}) in the linear regime for subhorizon scales is given by (Pogosian & Silvestri 2008; Zhang 2006; de la CruzDombriz et al. 2008) (14)where (15)In GR, f_{RR} = Q = 0, and there is no scale dependence for the density contrast in the linear regime. For ΛCDM the growing mode can be expressed in terms of hypergeometric function as (Silveira & Waga 1994) (16)We solved Eq. (14) numerically and obtained the growing mode for the γ gravity. By using (16), we then obtained the fractional change in the matter power spectrum P(k) relative to ΛCDM. Figure 2 shows ΔP_{k}/P_{Λ} at z = 0 for the three choices of parameters shown in Table 1.
Fig. 2 Fractional difference in the matter power spectrum with respect to ΛCDM for different values of n and α, as indicated in Table 1. These results are used to compare with the nonlinear power spectrum from our numerical simulations. 

Open with DEXTER 
3. Nbody equations
f(R)models are equivalent to a scalartensor theory (Brax et al. 2008), where the first derivative of the f(R) function, f_{R}. This field propagates according the equation (17)where κ^{2} = 8πG/c^{4} and T is the trace of energymomentum tensor, T = g_{μν}T^{μν}. In the quasistatic limit (see, e.g., Noller et al. 2014; Bose et al. 2015; Llinares & Mota 2014, 2013) this equation becomes (18)where (19)and Δ_{R}(a) = x_{2}(a) + 12x_{1}(a). The Ricci scalar R in function of f_{R} is given by inverting Eq. (3) (20)The geodesic equation, needed to update the particle positions, reads (21)where Φ is the Newtonian potential, which the dynamics is given by the Poisson equation (22)When implementing these equations in the Nbody code, we need to rewrite them in codeunits given by
(23)Here B_{0} is the size of the simulation box. In terms of the evolution equations becomes These are the only equations we need to implement and solve in the Nbody code.
For comparison we also need the linearized field equation. Simulations with this equation compared to the full f_{R} equation is a good measure of the amount of screening that takes place in the model. The linearized f_{R} equation is simply (27)where δf_{R} = f_{R}−f_{R}(a) and . In code units, taking , we obtain (28)and the geodesic equation becomes (29)We have (30)
4. Implementation in the ISIS code
Implementing scalartensor theories of gravity in Nbody code is rather straightforward because the scalartensor theories all contribute as a fifth force and because RAMSES, which ISIS is based on, has been widely used, thoroughly tested, and optimized. In this section we describe how the equations we need to solve are implemented in ISIS. For more details see Llinares et al. (2014).
To solve for f_{R} directly is not numerically stable since the solution can potentially vary over several orders of magnitude when going from deep voids to massive clusters in our simulation. We therefore introduce a field redefinition  f_{R}  = A(u) → ∇  f_{R}  = b(u)∇u where b(u) = dA(u) / du. The general field equation for f_{R} discretized on a grid with the fieldredefinition , where u is the field we solve for, can be written as (31)where (32)where h is the grid spacing and and where we have defined .
The equations are solved using NewtonGaussSeidel relaxation (with multigrid acceleration). The method consists of going through the grid and updating the solution using (33)For this we also need ∂ℒ(u_{i,j,k}) /∂u_{i,j,k}, which is given by (34)where (35)and c_{i,j,k} = c(u_{i,j,k}), where . Some problems arise related to how boundary conditions are handled on the refined grids in the code. This is discussed in Llinares et al. (2014) and Li et al. (2012b).
Fig. 3 Scalar field from tests. Different colors depict the different refinement levels. The left panel shows the field from a spherical density distribution, the right side shows a 1D sine field obtained in the second test. 

Open with DEXTER 
Our solver needs a starting guess, and for this we use the cosmological background solution, that is, the solution we would expect when there are no matter sources, (36)This is only needed when the simulation is stared as otherwise we can use the old solution as our guess. For γ gravity our choice for A and the related expressions for b are and the background value of u is (40)in terms of u we have
When the fifth force is implemented, we simply replace it with an effective force F_{eff} that includes the effects of modified gravity wherever the code normally works with the gravitational force, F_{N}, (43)The expression for the fifth force in codeunits is given in Eq. (24), and we use the same fivepoint stencil as RAMSES uses to compute the gravitational force to compute the fifth force .
5. Tests of the Nbody solver
To verify that the solver is implemented correctly, we tested it several times. We present some of these tests in this section. The density is provided to the code through a distribution of particles. The density estimation (CIC) and refinement criteria are the same as those used for the cosmological simulations. We also have the option to set the density field analytically in the code. To test that the treatment of the boundary of the refinement is correct, we included two levels of refinement. The model parameters used for the tests are the same as mentioned previously in Table 1.
The first test case corresponds to a sphere of radius R of constant density located in the center of the box and embedded in a uniform background: (44)where characterizes the density contrast between the inside and outside of the sphere, is the mean density, R is the radius of the sphere, and B_{0} the size of the box. The value of δ chosen for the test is 5000. For the f(R) test we used R = 25 Mpc/h and B_{0} = 250 Mpc/h. A spherical symmetric configuration is effectively onedimensional, therefore the field equation reduces to an ODE, which we solved by using Mathematica and used this to compare with.
For the second test we analytically computed a density field ρ(x,y,z) ≡ ρ(x) (i.e., a 1D configuration), using the field equation, so that the solution is given by a sine: u ∝ 2 + sin(2πx).
Figure 3 shows the result of both these tests and the different colors depict the different refinement levels. We see that the curves are smooth when going from one level to another, which demonstrates that boundary conditions are handled properly. The tests were performed using the serial version of the code, and both tests give the expected results, which demonstrates that the code works properly.
6. Simulations
The γ gravity simulations were run using 512^{3} dark matter particles with a box size of B_{0} = 250 Mpc/h, and we refined cells whenever it had more than eight particles in it. The highest refinement level reached in our simulation was eight. The background cosmology used for the simulation was computed using Eq. (7) with h = 0.71, Ω_{Λ} = 0.733 and Ω_{m0} = 0.267.
The model parameters were presented in Sect. 2, with some plots of various background quantities. We also compared these simulations with simulations of the HuSawicki f(R) model from Llinares et al. (2014).
7. Results
7.1. Power spectrum
The nonlinear matter power spectrum is an important observable and could be used to distinguish among different models of structure formation. As we showed above, γ gravity can have a strong effect on the growth rate of the linear perturbations. We expect these signatures to be detectable in the nonlinear matter power spectrum.
To compute the power spectrum we used a public code, POWMES (Colombi & Novikov 2011), which uses folding methods to compute the power spectrum.
Figure 4 displays the difference of the matter power spectrum with respect to ΛCDM, defined as ΔP/P_{ΛCDM} ≡ P(k) /P_{ΛCDM}(k)−1 for γ gravity together with the corresponding predictions from linear perturbations theory and the HuSawicki f(R) model for comparison. We focus on the presentday epoch, which corresponds to z = 0.
Figure 4 shows that the differences from ΛCDM are lower than 5−10% for all of our runs in the range 0.05 h Mpc^{1} ≲ k ≲ 0.1 h Mpc^{1} and in agreement with the linear perturbation result seen in Fig. 2. The deviation from ΛCDM approaches zero for larger scales. This is because the range of fifth force is smaller than the horizon, therefore the modifications of gravity are not felt on the largest scales. On smaller, nonlinear scales, the full effect of the fifth force acts, and we see larger deviations. However, ΔP/P_{ΛCDM} continues to grow with k in γ gravity, in contrast to the HuSawicky model, where the chameleon screening mechanism is in play on these scales, which reduces the power enhancement. This is an indication that the screening mechanism does not work very well in our simulation. We discuss this in more detail in Sect. 7.4.
Fig. 4 Fractional deviation of the dark matter power spectrum with respect to the ΛCDM model. The linear predictions of γ gravity are represented by dashed lines. For comparison we also show the results from simulations with the HuSawicky f(R) model with  f_{R0}  = { 10^{4},10^{5},10^{6} } from Llinares et al. (2014). 

Open with DEXTER 
7.2. Halo mass function
The halo mass function, the number density of halos with a given mass, is a useful tool for investigating the efficiency of a model in forming halos of different masses. To locate halos in the simulation outputs, we used the AHF (Amiga Halo Finder, Knollmann & Knebe 2009). We determined the halo mass function by binning halos in logarithmic mass intervals.
In Fig. 5 we show the fractional difference with respect to ΛCDM of halo mass function computed from our simulations. Our measurement of the halo mass function itself is limited by statistics and to a lesser extent, by the resolution in the high and low mass end. However, we can reduce the impact of these two effects by considering the relative difference between the halo mass functions measured in modified gravity and ΛCDM simulations with the same initial condition.
Figure 5 shows that we find an excess of halos in the range 10^{13} ≲ M/M_{⊙} ≲ 10^{15} probed by our simulation in γ gravity, and the signatures are very similar to the  f_{R0}  ≳ 10^{5} HuSawicki model. For the most massive halos (M/M_{⊙} ≳ 10^{15}) the effect of screening is much more pronounced in HuSawicki models, and the massfunction approaches ΛCDM. This does not occur for γ gravity and is again an indication that the chameleon screening mechanism does not work very efficiently in our simulations.
For γ gravity, the number of halos increases significantly, especially at the high mass end, by up to 40−100% for clustersized halos. For HuSawicki models, when the value of the f_{R} field becomes comparable to the cosmological potential wells, the chameleon effect starts to operate. This can be seen in the mass function, where deviations from ΛCDM approach zero in the high mass end for models with  f_{R0}  = 10^{5} and  f_{R0}  = 10^{6}.
The large deviations we find Δn/n_{ΛCDM} ~ 0.5−1 in the highmass end are probably already ruled out by present cluster counts (Cataneo et al. 2015; Mak et al. 2012). This does not rule out the model per se, but requires choices of parameters where  f_{R0}  is much smaller today, leading to no signatures in the background evolution (i.e., in the Hubble factor) of the Universe.
Fig. 5 Fractional difference of the halo mass function for γ gravity (data points) with respect to the ΛCDM model. For comparison we also show the results from simulations with the HuSawicky f(R) model with  f_{R0}  = { 10^{4},10^{5},10^{6} } (dashed lines) from Llinares et al. (2014). 

Open with DEXTER 
7.3. Halo profiles
Fig. 6 Fractional difference in the halo density profiles with respect to ΛCDM for four different mass bins. For comparison we also show the results from simulations with the HuSawicky f(R) model with  f_{R0}  = { 10^{4},10^{5},10^{6} } (dashed lines) from Llinares et al. (2014). 

Open with DEXTER 
We also studied how modifications of gravity change the density and velocity profiles of dark matter halos. We focused on halos in the mass ranges 0.5−1 × 10^{13} Mpc/h, 1−5 × 10^{13} Mpc/h, 0.5−1 × 10^{14} Mpc/h, and 1−5 × 10^{14} Mpc/h. The density profiles were calculated by binning dark matter particles in annular bins for each halo. Our calculated density profiles are averages of all density profiles of the proper size, ranging from 10% of the virialization radius, r = 0.1R_{vir}, to ten times the virialization radius, r = 10R_{vir}. This range was chosen to properly include all behaviors of the fifth force on the dark matter halos while also avoiding the inner regions of the halos, where the resolution of our simulations is not sufficient.
Figure 6 shows the fractional difference with respect to ΛCDM in the density profiles. We first note that the inner regions (R<R_{vir}) of halos for γ gravity are significantly denser than in ΛCDM. This difference is compensated for in outer regions (R>R_{vir}). Moreover, the profiles between the different model parameters do not differ appreciably from each other, and this pattern repeats for all ranges. However, the density profiles for HuSawicki models in general show stronger clustering in the lowdensity regions in the outskirts of halos than in the inner regions.
In the velocities profiles, shown in Fig. 7, we expect the fifth force to increase the velocity dispersion. This effect is very similar for both models, except for most massive halos, for which the HuSawicki models are more screened, causing a substantial decrease in comparison to γ gravity. For the HuSawicki model the velocities are boosted by ~20% in the  f_{R0}  = 10^{4} case and by 5% in the  f_{R0}  = 10^{6} case for all the three halo mass ranges analyzed. Only for the  f_{R0}  = 10^{5} parameters a mixed behavior is seen, that is, for the lower two halo mass ranges the boost is ~20%, but for the heaviest halos there is no deviation from ΛCDM in the inner parts and ~15% higher velocities are found in the outer parts. For γ gravity, the difference between the models we simulated is more expressive for less massive halos (5 × 10^{12}<M/M_{⊙}< 10^{13} and 10^{13}<M/M_{⊙}< 5 × 10^{13}), for most massive halos (10^{14}<M/M_{⊙}< 5 × 10^{14}) only the case α = 1.5 is distinguishable from the others.
Fig. 7 Fractional difference in the velocity dispersion profiles with respect to ΛCDM for four different mass bins. For comparison we also show the results from simulations with the HuSawicky f(R) model with  f_{R0}  = { 10^{4},10^{5},10^{6} } (dashed lines) from Llinares et al. (2014). 

Open with DEXTER 
7.4. Fifth force and screening
General relativity is very well tested on very small scales, especially inside the solar system. To ensure that F_{φ} is not relevant at these scales, we need a screening mechanism to suppress the fifth force at smallscale, very high curvature regimes.
When γ gravity was proposed in O’Dwyer et al. (2013), the authors explored the compatibility between the model and the solar system experiments using the chameleon mechanism for screening, as any other f(R) (Brax et al. 2008).
As we showed in the results above, the screening mechanism does not seem to be working very efficiently for the models simulated here. To check how much the fifth force is screened in our simulations, we compared the magnitude of the fifth and Newtonian force on the particles in our simulation box, as in Davis et al. (2012).
Figure 8 shows a scatter plot for this comparison at redshift z = 0. The dispersion for small F_{N}< 0.01 is expected (numerical scatter), here the forces are tiny, so the scatter here has a very weak effect on the simulation. The important part in this figure is the behavior for large F_{N}. If we have screening, then we should see F_{φ} ≪ F_{N}/ 3 in highdensity regions (which corresponds to high values of F_{N}). The result we find shows that the force ratio is roughly onethird – which is the linear prediction – everywhere, meaning that there is very little screening present in our simulations.
Fig. 8 Comparison between the magnitudes of Newtonian force (F_{N}) and fifth force (F_{φ}), the dashed line indicates F_{N}/F_{φ} = 1 / 3. The forces are in units of H_{0}/c^{2}. 

Open with DEXTER 
To understand this better, we revisit the screening conditions. Considering a spherical symmetric body with constant density ρ_{c} embedded in the background of constant density ρ_{b}, the solutions to the field equation (see e.g. Brax et al. 2012) mean that the fifth force is given by where is the socalled screening factor (also called the thinshell) and f_{Rb} (f_{Rc}) is the value of f_{R} in the background (inside the body), where ρ = ρ_{b} (ρ_{c}). If the fifth force is screened and we recover General Relativity. If we instead find that and gravity is significantly modified.
When the field is located in the minimum of its effective potential in the background^{2}, we can calculate the screening factor analytically. Assuming this, we find that for the γ gravity models considered here, Earth and Sun are almost completely screened and pass local gravity constraints assuming the galaxy is screened.
However, for the parameter values considered in this paper, the screening condition gives that the galaxy is not screened which again implies that the Earth and the sun is not screened either^{3}. The only caveat to this is that screening might have survived from earlier times. At early redshift and almost all objects are screened. Then in the cosmological evolution very quickly evolves to high values . The field inside our galaxy might have been trapped, ensuring screening. This is not expected, but to rigorously rule out (or confirm) this possibility, we would need simulations beyond the quasistatic approximation that we assumed in our simulations. This is beyond the scope of this paper.
8. Summary and conclusions
We have investigated the nonlinear evolution in γ gravity, a f(R) theory of gravity that is a viable alternative to ΛCDM. The models we investigated use a screening mechanism to suppress the deviations from General Relativity at small (solar system) and large cosmological scales. Specifically, this is the chameleon screening mechanism. As a result of this screening mechanism, the strongest signatures in these models are expected to occur at the nonlinear regime of structure formation. Therefore, to unveil the imprints of such theories at astrophysical scales, we ran several cosmological Nbody simulations. We compared models with ΛCDM and the HuSawicki model and showed that several astrophysical observables (halo mass function, density profiles, and power spectra) show the signatures of the model significantly strongly.
For the matter power spectrum we found a small deviation, lower than 10%, on large scales (k ≲ 10^{1}h/Mpc), which is consistent with the predictions of linear perturbation theory. For small scales (k ≳ 10^{1}h/Mpc), on the other hand, we found a strong increase in power, the largest deviation (α = 1.18) reaches ~60%, and in the other cases (α = 1.05 and α = 1.5) it reaches ~50%. This is quite different from the results found in the HuSawicky model, where the screening of the fifth force is much more active in suppressing the enhancement of power on small scales, the largest deviation is ~30%, much lower than for γ gravity.
For the halo mass function we found that all of our runs gave very similar result in the mass range 11.5 < log (M/M_{⊙}) < 14.7, and these results are very similar to what is found in the HuSawicky model for  f_{R0}  ≥ 10^{5}; a 10−60% increase in the halo abundance. However, for γ gravity we found an enhancement of 40−100% in the abundance for log (M/M_{⊙}) > 14.5, which can be compared to the HuSawicki model, where the halo mass function approaches that of ΛCDM.
For halo density profiles we found that the different runs with γ gravity were not distinguished well in any mass range; in all situations we saw that the inner regions (R < R_{vir}) of halos for γ gravity are significantly thicker than the halos of ΛCDM, around 10% for less massive halos, 5 × 10^{12} < M/M_{⊙} < 10^{13}, increasing with the halo mass until they reached ~30% for the most massive halos, 10^{14} <M/M_{⊙} < 5 × 10^{14}. This difference is compensated for in the outer regions (R > R_{vir}). For the HuSawicki model the same effect acts for less massive halos 5 × 10^{12} <M/M_{⊙} < 10^{13}, but it changes for more massive halos, 5 × 10^{13} < M/M_{⊙} < 10^{14} and 10^{14} < M/M_{⊙} < 5 × 10^{14} when the screening mechanism suppresses the fifth force. This then leads to a suppression in the accretion of mass concentration in the inner regions (compared to the nonscreened case).
Finally, we computed the velocity profiles, and as we expected, the velocities are significantly enhanced in γ gravity compared to ΛCDM. The other cases are comparable in all mass ranges, but with some important signatures. The difference between the different γ gravity models is larger for lowmass halos, 5 × 10^{12} < M/M_{⊙} < 10^{13}, and we can distinguish them from each other; the differences reach ~30% close to the boundary region, R ~ R_{vir}. This difference decreases with the mass of the halos, and for the most massive halos 10^{14} < M/M_{⊙} < 5 × 10^{14}. Only the case α = 1.5 is lower than ~10% from the others, which are practically identical.
The chameleon mechanism – the screening mechanism that makes f(R) gravity viable – is not very effective for the parameter choices considered in this paper. This explains the large deviations from ΛCDM in the observables we considered. Especially the cluster count signatures of 40−100% in the highmass end M ≳ 10^{14.5}M_{⊙} disagree with current observations. This does not rule out the model per se, but require choices of parameters where  f_{R0}  is much smaller today, which implies that the model has no observable signatures in the background evolution (i.e., in the Hubble factor) of the Universe.
With observable signatures we mean that the equation of state differs enough from the ΛCDM value of w = −1 at low redshifts that such a deviation could be detected by nearfuture experiments like WFIRST (Spergel et al. 2015), for example.
Acknowledgments
M.V.S. thanks the Brazilian research agency CAPES and the University of Oslo for support and Max B. Grönke and Amir Hammami for useful discussions. H.A.W. is supported by BIPAC and the Oxford Martin School. D.F.M. would like to thank the Research Council of Norway for funding. The simulations were performed on the NOTUR cluster HEXAGON, which is the computing facility at the University of Bergen.
References
 Amendola, L., Gannouji, R., Polarski, D., & Tsujikawa, S. 2007, Phys. Rev. D, 75, 083504 [NASA ADS] [CrossRef] [Google Scholar]
 Appleby, S. A., & Battye, R. A. 2007, Phys. Lett. B, 654, 7 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Appleby, S. A., Battye, R. A., & Starobinsky, A. A. 2010, J. Cosmol. Astropart. Phys., 1006, 005 [NASA ADS] [CrossRef] [Google Scholar]
 Berti, E., Barausse, E., Cardoso, V., et al. 2015, Class. Quant. Grav. 32, 243001 [Google Scholar]
 Boehmer, C. G., Burnett, J., Mota, D. F., & Shaw, D. J. 2010, J. High Energy Phys., 07, 053 [NASA ADS] [CrossRef] [Google Scholar]
 Bose, S., Hellwing, W. A., & Li, B. 2015, J. Cosmol. Astropart. Phys., 2, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Brax, P., van de Bruck, C., Davis, A.C., & Shaw, D. J. 2008, Phys. Rev. D, 78, 104021 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Brax, P., Davis, A.C., Li, B., & Winther, H. A. 2012, Phys. Rev. D, 86, 044015 [NASA ADS] [CrossRef] [Google Scholar]
 Brax, P., Davis, A.C., Li, B., Winther, H. A., & Zhao, G.B. 2013, J. Cosmol. Astropart. Phys., 4, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Cai, Y.C., Padilla, N., & Li, B. 2015, MNRAS, 451, 1036 [NASA ADS] [CrossRef] [Google Scholar]
 Capozziello, S., Cardone, V. F., Carloni, S., & Troisi, A. 2003, Int. J. Mod. Phys. D, 12, 1969 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Carroll, S. M., Duvvuri, V., Trodden, M., & Turner, M. S. 2004, Phys. Rev. D, 70, 043528 [NASA ADS] [CrossRef] [Google Scholar]
 Cataneo, M., Rapetti, D., Schmidt, F., et al. 2015, Phys. Rev. D, 92, 044009 [NASA ADS] [CrossRef] [Google Scholar]
 Cognola, G., Elizalde, E., Nojiri, S., et al. 2008, Phys. Rev. D, 77, 046009 [NASA ADS] [CrossRef] [Google Scholar]
 Colombi, S., & Novikov, D. 2011, Astrophysics Source Code Library [record ascl:1110.017] [Google Scholar]
 Davis, A.C., Li, B., Mota, D. F., & Winther, H. A. 2012, ApJ, 748, 61 [NASA ADS] [CrossRef] [Google Scholar]
 de la CruzDombriz, A., Dobado, A., & Maroto, A. L. 2008, Phys. Rev. D, 77, 123515 [NASA ADS] [CrossRef] [Google Scholar]
 Frolov, A. V. 2008, Phys. Rev. Lett., 101, 061103 [NASA ADS] [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Gronke, M., Llinares, C., Mota, D. F., & Winther, H. A. 2015, MNRAS, 449, 2837 [NASA ADS] [CrossRef] [Google Scholar]
 Hammami, A., Llinares, C., Mota, D. F., & Winther, H. A. 2015, MNRAS, 449, 3635 [NASA ADS] [CrossRef] [Google Scholar]
 He, J.h., Li, B., & Jing, Y. P. 2013, Phys. Rev. D, 88, 103507 [NASA ADS] [CrossRef] [Google Scholar]
 He, J.h., Li, B., & Hawken, A. J. 2015, Phys. Rev. D, 92, 103508 [NASA ADS] [CrossRef] [Google Scholar]
 Hellwing, W. A., Cautun, M., Knebe, A., Juszkiewicz, R., & Knollmann, S. 2013, J. Cosmol. Astropart. Phys., 1310, 012 [NASA ADS] [CrossRef] [Google Scholar]
 Hellwing, W. A., Barreira, A., Frenk, C. S., Li, B., & Cole, S. 2014, Phys. Rev. Lett., 112, 221102 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, W., & Sawicki, I. 2007, Phys. Rev. D, 76, 064004 [NASA ADS] [CrossRef] [Google Scholar]
 Jain, B., & Khoury, J. 2010, Annals Phys., 325, 1479 [NASA ADS] [CrossRef] [Google Scholar]
 Knollmann, S. R., & Knebe, A. 2009, ApJS, 182, 608 [NASA ADS] [CrossRef] [Google Scholar]
 Li, B., Mota, D. F., & Barrow, J. D. 2011, ApJ, 728, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Li, B., Zhao, G.B., Teyssier, R., & Koyama, K. 2012a, J. Cosmol. Astropart. Phys., 1, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Li, B., Zhao, G.B., Teyssier, R., & Koyama, K. 2012b, J. Cosmol. Astropart. Phys., 1201, 051 [NASA ADS] [CrossRef] [Google Scholar]
 Li, B., Hellwing, W. A., Koyama, K., et al. 2013, MNRAS, 428, 743 [NASA ADS] [CrossRef] [Google Scholar]
 Linder, E. V. 2009, Phys. Rev. D, 80, 123528 [NASA ADS] [CrossRef] [Google Scholar]
 Llinares, C., & Mota, D. 2013, Phys. Rev. Lett., 110, 161101 [NASA ADS] [CrossRef] [Google Scholar]
 Llinares, C., & Mota, D. F. 2014, Phys. Rev. D, 89, 084023 [NASA ADS] [CrossRef] [Google Scholar]
 Llinares, C., Mota, D. F., & Winther, H. A. 2014, A&A, 562, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lombriser, L., Schmidt, F., Baldauf, T., et al. 2012, Phys. Rev. D, 85, 102001 [NASA ADS] [CrossRef] [Google Scholar]
 Lombriser, L., Koyama, K., & Li, B. 2014, J. Cosmol. Astropart. Phys., 1403, 021 [NASA ADS] [CrossRef] [Google Scholar]
 Mak, D. S. Y., Pierpaoli, E., Schmidt, F., & Macellari, N. 2012, Phys. Rev. D, 85, 123513 [NASA ADS] [CrossRef] [Google Scholar]
 Mota, D. F., Shaw, D. J., & Silk, J. 2008, ApJ, 675, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Noller, J., von BraunBates, F., & Ferreira, P. G. 2014, Phys. Rev. D, 89, 023521 [NASA ADS] [CrossRef] [Google Scholar]
 O’Dwyer, M., Joras, S. E., & Waga, I. 2013, Phys. Rev. D, 88, 063520 [NASA ADS] [CrossRef] [Google Scholar]
 Oyaizu, H., Lima, M., & Hu, W. 2008, Phys. Rev. D, 78, 123524 [NASA ADS] [CrossRef] [Google Scholar]
 Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565 [NASA ADS] [CrossRef] [Google Scholar]
 Pogosian, L., & Silvestri, A. 2008, Phys. Rev. D, 77, 023503; Erratum: 2010, 81, 049901 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration XXII. 2014, A&A, 571, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puchwein, E., Baldi, M., & Springel, V. 2013, MNRAS, 436, 348 [NASA ADS] [CrossRef] [Google Scholar]
 Pujol, A., Gaztañaga, E., Giocoli, C., et al. 2014, MNRAS, 438, 3205 [NASA ADS] [CrossRef] [Google Scholar]
 Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, F. 2010, Phys. Rev. D, 81, 103002 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, F., Lima, M. V., Oyaizu, H., & Hu, W. 2009, Phys. Rev. D, 79, 083518 [NASA ADS] [CrossRef] [Google Scholar]
 Shi, D., Li, B., Han, J., Gao, L., & Hellwing, W. A. 2015, MNRAS, 452, 3179 [NASA ADS] [CrossRef] [Google Scholar]
 Silveira, V., & Waga, I. 1994, Phys. Rev. D, 50, 4890 [NASA ADS] [CrossRef] [Google Scholar]
 Song, Y.S., Taruya, A., Linder, E., et al. 2015, Phys. Rev. D, 92, 043522 [NASA ADS] [CrossRef] [Google Scholar]
 Spergel, D., Gehrels, N., Baltay, C., et al. 2015, ArXiv eprints [arXiv:1503.03757] [Google Scholar]
 Starobinsky, A. A. 1980, Phys. Lett. B, 91, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Starobinsky, A. A. 2007, JETP Lett., 86, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Terukina, A., & Yamamoto, K. 2012, Phys. Rev. D, 86, 103503 [NASA ADS] [CrossRef] [Google Scholar]
 Tessore, N., Winther, H. A., Metcalf, R. B., Ferreira, P. G., & Giocoli, C. 2015, J. Cosmol. Astropart. Phys., 10, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Thongkool, I., Sami, M., Gannouji, R., & Jhingan, S. 2009, Phys. Rev. D, 80, 043523 [NASA ADS] [CrossRef] [Google Scholar]
 Wilcox, H., Bacon, D., Nichol, R. C., et al. 2015, MNRAS, 452, 1171 [NASA ADS] [CrossRef] [Google Scholar]
 Winther, H. A., Mota, D. F., & Li, B. 2012, ApJ, 756, 166 [NASA ADS] [CrossRef] [Google Scholar]
 Winther, H. A., Schmidt, F., Barreira, A., et al. 2015, MNRAS, 454, 4208 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, P. 2006, Phys. Rev. D, 73, 123504 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, G.B., Li, B., & Koyama, K. 2011, Phys. Rev. D, 83, 044007 [NASA ADS] [CrossRef] [Google Scholar]
 Zivick, P., Sutter, P. M., Wandelt, B. D., Li, B., & Lam, T. Y. 2015, MNRAS, 451, 4215 [NASA ADS] [CrossRef] [Google Scholar]
 Zumalacarregui, M., Koivisto, T. S., & Mota, D. F. 2013, Phys. Rev. D, 87, 083010 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Effective equationofstate parameter w_{DE} as a function of redshift z for the parameters given in Table 1. The strongest deviation from − 1 is lower than 4%. 

Open with DEXTER  
In the text 
Fig. 2 Fractional difference in the matter power spectrum with respect to ΛCDM for different values of n and α, as indicated in Table 1. These results are used to compare with the nonlinear power spectrum from our numerical simulations. 

Open with DEXTER  
In the text 
Fig. 3 Scalar field from tests. Different colors depict the different refinement levels. The left panel shows the field from a spherical density distribution, the right side shows a 1D sine field obtained in the second test. 

Open with DEXTER  
In the text 
Fig. 4 Fractional deviation of the dark matter power spectrum with respect to the ΛCDM model. The linear predictions of γ gravity are represented by dashed lines. For comparison we also show the results from simulations with the HuSawicky f(R) model with  f_{R0}  = { 10^{4},10^{5},10^{6} } from Llinares et al. (2014). 

Open with DEXTER  
In the text 
Fig. 5 Fractional difference of the halo mass function for γ gravity (data points) with respect to the ΛCDM model. For comparison we also show the results from simulations with the HuSawicky f(R) model with  f_{R0}  = { 10^{4},10^{5},10^{6} } (dashed lines) from Llinares et al. (2014). 

Open with DEXTER  
In the text 
Fig. 6 Fractional difference in the halo density profiles with respect to ΛCDM for four different mass bins. For comparison we also show the results from simulations with the HuSawicky f(R) model with  f_{R0}  = { 10^{4},10^{5},10^{6} } (dashed lines) from Llinares et al. (2014). 

Open with DEXTER  
In the text 
Fig. 7 Fractional difference in the velocity dispersion profiles with respect to ΛCDM for four different mass bins. For comparison we also show the results from simulations with the HuSawicky f(R) model with  f_{R0}  = { 10^{4},10^{5},10^{6} } (dashed lines) from Llinares et al. (2014). 

Open with DEXTER  
In the text 
Fig. 8 Comparison between the magnitudes of Newtonian force (F_{N}) and fifth force (F_{φ}), the dashed line indicates F_{N}/F_{φ} = 1 / 3. The forces are in units of H_{0}/c^{2}. 

Open with DEXTER  
In the text 