Halos and Voids in f(R) Gravity

In this paper, we study the distribution of dark matter halos and voids using high resolution simulations in f(R) gravity models with the chameleon mechanism to screen the fifth force in dense environment. For dark matter halos, we show that the semi-analytic thin shell condition, with a suitably-defined environment, provides a good approximation to describe the mass and environmental dependence of the screening of the fifth force in halos. Due to stronger gravity, there are far more massive halos and large voids in f(R) models compared with the \Lambda CDM model. The numbers of voids with an effective radius of 15Mpc/h are twice and four times as many as those in \Lambda CDM for f(R) models with |f_{R0}|=1e-5 and 1e-4 respectively. This provides a new means to test the models using the upcoming observational data. We also find that halos inside voids are all unscreened in our simulations, which are ideal objects for the gravity test.


INTRODUCTION
The biggest problem in cosmology is to explain the observed accelerated expansion of the universe. Within the framework of General Relativity (GR), the acceleration originates from dark energy (Copeland et al. 2006). Alternatively, a largedistance modification to GR may account for the late-time acceleration of the universe.
It has been recognised that usually once we modify Einstein gravity on large scales, there can appear a new scalar degree of freedom in gravity which mediates a fifth force. Without a mechanism to suppress this additional force, most modified gravity models are excluded by stringent constraints on deviations from GR on the solar system scale. One way to evade these constraints is to exploit a chameleon mechanism (Khoury & Weltman 2004). The new scalar degree of freedom couples to the energy density of matter. By tuning the coupling and potential for the scalar mode, it is possible to realise a situation that in dense environments such as the solar system, the scalar field has a large mass and it essentially does not mediate the fifth force. On the other hand, on cosmological scales, this scalar mode is light and modifies gravity significantly.
In models with the chameleon mechanism, there ap-⋆ E-mail: baojiu.li@durham.ac.uk † E-mail: gong-bo.zhao@port.ac.uk ‡ E-mail: kazuya.koyama@port.ac.uk pears environmental dependence on the properties of dark matter distributions. In dense environment such as in clusters, the chameleon works efficiently and the modification of gravity is suppressed. On the other hand, in underdense regions such as in voids, the chameleon mechanism does not work and gravity is significantly modified. As is shown in our previous paper (Zhao, Li & Koyama 2011b), this environmental dependence is a smoking gun for the modification of gravity in models with the chameleon mechanism.
In this paper, we study the properties of dark matter halos and voids in models with the chameleon mechanism to reveal the environmental dependence of dark matter distributions. We use f (R) gravity as an example and exploit recent results from high resolution simulations described in Zhao, Li & Koyama (2011a). In the f (R) gravity, the fifth force can enhance gravity by a factor of 1/3 but this enhancement is suppressed by the chameleon mechanism in the overdense regions. Some other numerical simulations performed for models with chameleon mechanism are those in Oyaizu (2008); Oyaizu et al. (2008); Schmidt et al. (2008); Li & Zhao (2009);Zhao et al. (2010);  and . Although we will not study those simulations directly, we expect that the results found here will be true for them as well.
The paper is organised as follows. In Section II, we summarise f (R) gravity models that we shall study in this paper. In Section III, we study the probability distribution function of the smoothed density field and show how the properties of halos and voids in f (R) gravity are modified compared with the standard ΛCDM model. Section IV is devoted to the study of dark matter halos. We study how the difference between dynamical and lensing masses, which arises due to the fifth force, depends on the mass and environment of halos. We find that the semi-analytic thin shell condition that determines the efficiency of the chameleon can well describe those dependence found in simulations. In Section V, we study the underdense regions by identifying voids in our simulations. We study the properties of halos inside and near the voids. We show that the number density of large voids is significantly modified in f (R) gravity models.

F (R) GRAVITY AND SIMULATIONS
The f (R) gravity, in which the Ricci scalar R in the Einstein-Hilbert action is generalised to a function of R, was designed to explain the observed cosmic acceleration without introducing dark energy. In such theories, the structure formation is determined by the following equations, where Φ denotes the gravitational potential, fR ≡ df (R) dR is the scalaron, the extra scalar degree of freedom, δR = R −R, δρM = ρM −ρM, and the quantities with overbar take the background values. In GR, the gravitational potential is solely determined by the distribution of matter, say, ∇ 2 Φ = 4πGa 2 δρM. This is a linear Poisson equation, which is much easier to solve. In f (R), however, the scalar field complicates the Poisson equation, making the effective Newton's constant vary with the local density: in underdense regions, the δR(fR) term in Eq. (1) vanishes thus two equations decouple, making the effective Newton's constant enhanced by 1/3. On the other hand, in the dense region, δfR in Eq.
(2) is negligible, so δR(fR) = −8πGδρM, which means that GR is locally restored. This is the chameleon mechanism making f (R) evade the stringent solar system tests, thus is important for the cosmological viability of the f (R) gravity. One could rewrite Eq. (1) as, where the effect of the scalar field is absorbed into the definition of the effective energy density δρ eff . Then the dynamical mass MD(r) of a halo is defined as the mass contained within a radius r, inferred from the gravitational potential felt by a test particle at r. It is given by MD ≡ a 2 δρ eff dV , in which the integral is over the extension of the body. On the other hand, the lensing mass is the true mass of the halo, i.e., ML ≡ a 2 δρMdV . Comparing the lensing mass with the dynamical mass of the same halo can in principal be a easy way to test GR. This is because the lensing mass and the dynamical mass are identical in GR, but quite different in MG scenarios. To quantify the difference, we calculate the relative difference ∆M between ML and MD for each halo, ∆M ≡ MD/ML −1.
Note that in f (R), ∆M 1/3 (Zhao, Li & Koyama 2011b). The probability distribution of matter density contrast field δ ≡ ρ(x; R)/ρ − 1. The δ field was filtered by a top-hat window with radius R = 2h −1 Mpc. In the plot we offset δ by 1 for the ease of using the logarithmic scale and the points plotted homogeneously in lg(1 + δ). The black squares, red circles, green triangles and blue diamonds are from the ΛCDM simulation and f (R) simulations with |f R0 | = 10 −6 , 10 −5 , 10 −4 respectively. Each curve represents the averaged result over ten realisations and is normalised so that the integration of P (1 + δ) is 1.
The presence of the chameleon effect indicates that Eqs (1) and (2) are highly nonlinear, so that the system cannot be solved without using N-body simulations. In this work, we we shall use the high-resolution N -body simulation catalogue (Zhao, Li & Koyama 2011a) for a f (R) gravity model, f (R) = αR/(βR + γ) (Hu & Sawicki 2007) where α = −m 2 c1, β = c2, γ = −m 2 , m 2 = H 2 0 ΩM and c1, c2 are free parameters. The expansion rate of the universe in this f (R) model is determined by c1/c2, and the structure formation depends on |fR0|, which is the value of |df /dR| at z = 0, and is proportional to c1/c 2 2 . We tune c1/c2 to obtain the same expansion history as that in a ΛCDM model, and choose values for |fR0| so that those models cannot be ruled out by current solar system tests. To satisfy these requirements, we set c1/c2 = 6ΩΛ/ΩM and simulate three models with |fR0| = 10 −4 , 10 −5 , 10 −6 . In this paper, we study the distribution of dark matter halos and voids in these simulations at z = 0. The method to identify halos is described in Li & Barrow (2011);Zhao, Li & Koyama (2011a).

PROBABILITY DISTRIBUTION OF DENSITY FIELD
In the standard cold dark matter paradigm, structures grow from the small inhomogeneities in the initial matter density field due to the pull of gravity. As a result, initial overdense (underdense) regions become more and more overdense (underdense). In f (R) gravity, gravity can be enhanced, so that the fifth force helps to pull more matter into overdense regions, and the underdense regions can be evacuated more efficiently.
In Fig. 1 we show the probability distribution of the matter density contrast field measured from our f (R) and ΛCDM simulations. This is calculated by filtering the density field by top-hat windows with as an example radius R = 2h −1 Mpc centred at each cell of the simulation grid, and counting how many such windows fall into a given density band.
As shown in this figure, the fifth force can tremendously increase the chance of creating extremely low-density regions in the Universe. For example, only ∼ 0.01% of the space in the ΛCDM paradigm has a density of 1 + δ = 0.05, while the f (R) models predict 5 (for |fR0| = 10 −6 ), 15 (for |fR0| = 10 −5 ) and 30 (for |fR0| = 10 −4 ) times as much. The effect becomes smaller for increasing δ, and for 1 + δ ∼ 0.2 − 0.3 the probability becomes roughly the same for all models. For 0.2 − 0.3 < 1 + δ <∼ 10 the fifth force actually decreases the probability and then for windows with 1 + δ >∼ 10 the fifth force makes it more likely to be found again by making matter cluster more strongly. Similar effects have been found for other models (Hellwing & Juszkiewicz 2009;Li 2011).
The peaks of the density distribution shifts towards low values as |fR0| increases, which shows that the Universe in f (R) gravity may look emptier overall. Meanwhile, Fig. 1 confirms that a f (R) universe will more likely host very big voids and very massive dark matter halos, both of which are rarer in a ΛCDM universe. We will come back to this point later. Fig. 1 also clearly shows the effect of the chameleon mechanism, which is known to work better for smaller |fR0| and for high density fields (Zhao, Li & Koyama 2011a). For |fR0| = 10 −6 , the deviation from LCDM is suppressed for high density fields while there is still a sizable deviation in the probability distribution for under-density fields. This shows that the modification of gravity is more prominent for voids for small |fR0|. This fact is important when we perform observational tests of modified gravity models with a realistic value of |fR0| compatible with the solar system constraints.

OVERDENSE REGIONS
In general, dark matter halos reside in high-density regions, which form their local environment. It is well known that the fifth force in f (R) gravity sensitively depends on the environment. Thus from a theoretical point of view, it is very important to understand how the environment changes the properties of the fifth force.
As mentioned above, f (R) gravity is a subclass of the chameleon scalar field theory, with fR = exp(γ √ κϕ) − 1, in which κ = 8πG = M −2 Pl and ϕ is the corresponding scalar field and γ = 2/3 is the constant coupling strength. The scalar field is governed by an effective potential (see, e.g., Li & Barrow (2007)) When the chameleon mechanism is at work, a spherical body will develop a thin shell, and the thickness of which is given by (Khoury & Weltman (2004); Li & Efstathiou (2011)), in which R is the radius of the body, ∆R is the thickness of the shell, ϕout, ϕin are the values of ϕ minimising V eff inside and outside the body, respectively. Similarly, ρin and ρout are the constant matter density inside and outside the body respectively. Note that only the matter inside the thin shell contributes to the fifth force exerting on a nearby test particle. From Eq. (5) it is evident that the fifth force could be suppressed if shell becomes thinner, which can be realised by the following two ways: (i) Increasing R and/or ρin, thereby making the body (in the case here the dark matter halo) more massive.
(ii) Decreasing ϕout, which involves increasing ρout or equivalently making the environment denser 1 .
As a result, massive halos in dense environments are strongly screened from the fifth force, while small halos in low-density environments are less screened and may experience the full fifth force.
In the f (R) gravity theory, the thin-shell expression Eq. (5) can be translated into the following equation, by using the relationship between √ κϕ and fR, and the fact that √ κϕ ∼ |fR| ≪ 1. The ratio between the magnitudes of the fifth force and gravity can be approximately estimated as (Li & Efstathiou 2011) ∆M has a maximum value of 1/3 as expected, and it can be analytically estimated by calculating ∆R/R from Eq. (6) as follows: (i) Given a halo's mass and virial radius we can compute the average ρin and therefore fR,in; (ii) ρout, the environmental density, can be estimated by computing the average density of a sphere with a radius Renv centring on the concerned halo (Li & Efstathiou 2011 at z = 0 whereρ out(in) ≡ ρ out(in) /ρ. So we get ∆M from Eq. (7) This could be measured for each halo from the N -body simulations, as is shown in Fig. 2. Note that the thin-shell condition for dark matter halos have been studied using N -body simulations by (Schmidt 2010), but here for the first time we have checked the environmental effects and compared with analytical results. The left panels of Fig. 2 show the dependence of ∆M on the halo mass (illustrated by the size and colour of the symbols) and the environment matter density (horizontal axis). The results are measured from our N -body simulations using two values of Renv, respectively 8h −1 Mpc (upper left) and 5h −1 Mpc (lower left). There are several interesting features: (i) Massive halos mostly reside in overdense regions, as could be seen from the correlation between the size of the symbols and the horizontal axis. This is as expected, because only in those regions there are enough particles to form large halos; (ii) With the same environmental density, small halos are less screened (have bigger ∆M ), which agrees with the analysis above; (iii) For halos with comparable mass, those in overdense regions are more strongly screened, because of the environmental effect; (iv) Only very few halos reside in underdense regions, and those are mostly small halos. This is because most particles in those regions have been pulled away. Note that the halos in those regions are essentially unscreened; (v) The dependence of the results on Renv is fairly weak, indicating that the exact value of Renv in the definition of the environment is not crucial.
In the right panels of Fig. 2 we have shown the analytical approximation obtained by using Eq. (10). We can find that the the analytical result agrees with the N -body simulation result quite well, and once again there is no sensitive dependence on Renv. This means that the analytical approximation Eq. (10) can well describe the nonlinear behaviour of the fifth force in f (R) gravity.
In Fig. 3 we show the same results for a different f (R) model, with |fR0| = 10 −5 . Again we could see that the analytical and numerical results agree well. Note that here all halos but the very massive ones are unscreened. Figs. 2, 3 lend supports to the simple excursion-set model of Li & Efstathiou (2011) for studying structure formation in the chameleon-type scalar field models.
Of course, the agreement between the analytical and simulation results is not perfect, although it is fairly good statistically. For |fR0| = 10 −5 , the analytic thin shell condition underestimates the fifth force especially with Renv = 5h −1 Mpc. This is probably because the Chameleon mechanism becomes weaker for larger |fR0| and the environmental effect becomes more non-local. Thus we need a larger Renv to capture the environmental effect correctly.
Another possible reason for the slight discrepancy between analytical and numerical results is that the dark matter halos are not rigorously spherical, which we have assumed when applying the thin-shell condition; we have tried to search for the possible correlation between the ellipticity of the halos and ∆M , but we didn't find any evidence, which means that using the spherical thin shell condition is indeed a good approximation.

VOIDS
Next let us turn to the underdense regions, or voids, in the f (R) gravity.
The voids are identified using VAMSUR (Voids As Merged Spherical Underdense Regions) code developed by Li (2011). The basic idea is to first find the low-density spherical regions (protovoids) and then merge them to form irregularlyshaped voids using a given algorithm. In this work we have chosen δ < −0.8 as the definition of voids.
In Fig. 4 we show the protovoids (left panel) and dark matter halos (right panel) identified in one of our simulation boxes for the model with |fR0| = 10 −6 . We can see that: (i) Most dark matter halos (in particular the more massive ones) distribute in regions where few protovoids can be identified, and vice versa, which is a trivial test of the code and numerical results.
(ii) Near the voids the dark matter halos are less screened (denoted in blue) while far from the voids they can be well screened (in red).
(iii) Halos inside voids are all unscreened as is also shown in Fig. 5.
These observations agree with our expectation very well.
As mentioned above, in the f (R) gravity the fifth force helps evacuate the low density regions, which results in more large voids than in ΛCDM. To see this more clearly, one could plot the number density of voids as a function of their effective volumes, and this is shown in Fig. 6 (see figure caption for details). Here we can see the clear trend: increasing |fR0|, which makes the fifth force less suppressed from ear- The horizontal axis of is ∆env, the matter overdensity of a sphere (the environment) centred at each halo, with comoving radius Renv as indicated in the panel. The vertical axis is ∆ M , the ratio between the magnitudes of the fifth force and gravity at the surface of a halo. Each circle represents a halo, and the halo's mass is illustrated by both the size (increasing size for increasing mass) and colour (from red to violet for increasing mass) of the circle. The left (right) panels are numerical (analytical) results (see text for a detailed description).  lier times, produces more large voids. For example, a f (R) universe with |fR0| = 10 −5 (10 −4 ) has twice (four times) as many voids with effective radius 15h −1 Mpc as a ΛCDM universe does, and the relative difference in the abundance for larger voids is even larger. Very large voids and very massive dark matter halos are rare objects in a ΛCDM Universe, and we have seen that both of them are more abundant in the f (R) universes. This is the reason why cluster abundance gives the strongest constraints on |fR0| with current observations (Schmidt et al. 2009;Ferraro et al. 2011;Lombriser et al. 2011). Compared with the halos, the modified gravity effect is more pronounced on the voids. This is because in f (R), gravity is maximumly enhanced in voids due to the presence of the fifth force. At z = 0, the model with |fR0| = 10 −5 predicts ∼30% more halos with mass ∼ 5 × 10 14 M⊙ than ΛCDM (Zhao, Li & Koyama 2011a), while it predicts twice as many voids of size ∼ 15000h −3 Mpc 3 .
Due to the limitation of our simulation box size, we do not have voids with radius larger than ∼ 15h −1 Mpc. However, there is an abundance of such large voids observed. For example, using the SDSS DR7, Pan et al. (2011) identified about 1000 voids in the northern galactic hemisphere with radii > 10h −1 Mpc; the largest and median radii in their void catalogue are 30h −1 Mpc and 17h −1 Mpc, respectively. Those voids have an edge density contrast of δ < −0.85. They find that their observations agree quite well with the ΛCDM simulations, which means that their data could place strong constraints on the f (R) gravity, making voids a promising tool to study the physics of the accelerated cosmic expansion and large-scale structure formation. Of course, their voids are identified by looking at galaxies in the survey, and to make direct comparison with their data we have to generate galaxy catalogues in the f (R) gravity. We will leave this to future work.

SUMMARY AND CONCLUSIONS
In this paper, we studied the over-and underdense regions, namely, the distribution of halos and voids, in f (R) gravity simulations. By comparing the probability distribution function of the density contrast δ for f (R) models with that for GR, we find that there are far more voids in f (R) gravity than that in GR. For example, the numbers of voids with an effective radius of 15h −1 Mpc are twice and four times as many as those in GR for f (R) models with |fR0| = 10 −5 and 10 −4 respectively. This in principle provides a new means to test GR observationally using the upcoming data. We also find that halos near the voids are less screened and experience stronger gravity. Especially, halos inside the voids are all unscreened in our simulations. This confirms the expectations that small galaxies inside voids provide us the best place for testing modification of gravity.
On the other hand, the overdense regions, ie, the distribution of dark matter halos, can provide important information for the GR test as well. In this work, we utilised the thin-shell condition developed in (Khoury & Weltman 2004) and (Li & Efstathiou 2011), and analytically predicted the fractional difference between the lensing mass and the dynamical mass of dark matter halos, ∆M , as a function of the environment. As we found, the analytic result agrees very well with the simulation result, which means that the thin-shell condition is a good approximation for f (R) gravity. This has important applications for the semi-analytic halo model building for f (R) gravity, which is crucial for realistic constraints of f (R) models using observations.