A Nonlinear equations

# Spherical collapse and halo mass function in \mathboldf(R) theories

## Abstract

We compute the critical density of collapse for spherically symmetric overdensities in a class of modified gravity models. For the first time we evolve the Einstein, scalar field and non-linear fluid equations, making the minimal simplifying assumptions that the metric potentials and scalar field remain quasi-static throughout the collapse. Initially evolving a top hat profile, we find that the density threshold for collapse depends significantly on the initial conditions imposed, specifically the choice of size and shape. By imposing ‘natural’ initial conditions, we obtain a fitting function for the spherical collapse as a function of collapse redshift, mass of the overdensity and , the background scalar field value at . By extending into drifting and diffusing barrier within the context of excursion set theory, we obtain a realistic mass function that might be used to confront this class of scalar-tensor models with observations of dark matter halos. The proposed analytic formula for the halo mass function was tested against Monte Carlo random walks for a wide class of moving barriers and can therefore be applied to other modified gravity theories.

Modified gravity, spherical collapse, halo mass function
###### pacs:
04.70.Bw, 04.25.dc, 98.80.Cq

## I Introduction

Einstein’s theory of General Relativity (GR) Einstein (1915) has withstood nearly one century of experimental testing. Many of its predictions have been confirmed through high precision laboratory and solar system experiments, and more recently with astrophyiscal and cosmological data (see Will (2006) for a review). Testing GR on the largest scales is a field still in its infancy, as probing cosmological distances is technically challenging. However current and future surveys mapping the large scale structure of the universe already provide a powerful tool to constrain gravity models Eisenstein et al. (2011); DES (); Pan-STARRS (); Laureijs et al. (2011).

In spite of its success, there are theoretical reasons to believe that GR is not a fundamental theory of gravity. In particular one expects corrections to the Einstein-Hilbert action through 1-loop corrections induced by matter Davies et al. (1977). In addition, there are puzzling cosmological observations such as the recently discovered accelerated expansion Riess et al. (1998); Perlmutter et al. (1999); Albrecht et al. (2006) which, although still consistent with GR and a cosmological constant term, might require a new gravity theory or the existence of an exotic form of matter known as dark energy.

In Starobinsky (1980) it was noticed that the local 1-loop corrections to the Einstein-Hilbert action are quadratic in curvature, and lead to a period of inflation in the very early universe. This motivates the introduction of additional functions of the curvature invariants, that allow for a dynamical late time acceleration Capozziello (2002); Starobinsky (2007) without the need for additional scalar fields or cosmological constant . The observed value of the cosmological constant is so extraordinarily small that one needs a finely tuned bare cosmological constant in the Einstein Hilbert action in order to cancel the large quantum vacuum energy and classical contributions related to Standard Model phase transitions . This severe “old” cosmological constant problem Weinberg (1989) is not resolved in theories. In this paper we assume that effectively due to some other physical process acting at the scale such that the scale enters the function naturally. In addition, the unknown physics leading to the removal of the vacuum energy might be accompanied by an effective scalar degree of freedom. There exist some ideas of screening Dvali et al. (2007); de Rham et al. (2008); Berkhahn et al. (2011) and relaxing Dolgov (1985); Charmousis et al. (2011); Emelyanov and Klinkhamer (2011); Bauer et al. (2010) the vacuum energy.

In this work we focus on a simple class of modified gravity models, so called gravity, where an arbitrary function of the Ricci scalar is introduced into the standard Einstein Hilbert action. It is well known that certain functional forms can give rise to an expansion history that exactly mimics a Universe governed by a cosmological constant and cold dark matter (CDM). However even for such models, the modified gravity contribution will still affect the growth of structure, both in the linear and nonlinear regime via a fifth force mediated by a scalar degree of freedom, the “scalaron” Starobinsky (1980) .

To be consistent with local and astrophysical gravity tests, functions must be used that essentially approach a constant value (usually ) in regions of high curvature. This ensures that in these high curvature regions the local scalaron mass becomes large enough to shut down the fifth force on much smaller scales. This effect is called the chameleon mechanism Khoury and Weltman (2004); Brax et al. (2004). At high curvature, such as in regions with deep potential wells in the late universe, e.g. galaxies and the solar system or in the early universe, at recombination, these theories are pushed towards GR. However at small curvature, typically well after matter-radiation-equality and on large scales, detectable deviations from GR are possible and still allowed observationally. This is exactly the regime where linear and non-linear large scale structure formation takes place. While theories should not be viewed as fundamental theories of gravity, they offer a self-consistent tool to scrutinize GR and look for new physics in the currently mapped large scale structure.

In Oyaizu et al. (2008); Martinelli et al. (2011) a constraint on theories using the galaxy power spectrum was calculated , where it was found that the chameleon mechanism supresses deviations from GR on small scales where perturbations become nonlinear. However the formation of galaxy clusters involves both linear and nonlinear growth and is thus better suited to probe both the fifth force and the chameleon effects present in scalar-tensor theories Li and Hu (2011).

The aim of this paper is to predict the number density of dark matter halos with a given mass and observed redshift as a function of the model parameter of the Hu-Sawicki-Starobinsky model Hu and Sawicki (2007); Starobinsky (2007). Building on the work of Corasaniti and Achitouv (2011); Achitouv et al. (2012), we employ the spherical collapse model to obtain a critical density contrast , from which we construct a realistic mass function. Similar semi-analytical formalisms were already applied to theories in previous works Schmidt et al. (2009); Brax et al. (2010); Borisov et al. (2011); Li and Efstathiou (2011); Li and Lam (2012); Lam and Li (2012); Lombriser et al. (2013), see Clampitt et al. (2013) for voids. In Schmidt et al. (2009) only two limiting cases were considered: either the collapsing overdensity was considered to be fully chameleon screened or fully unscreened. As these two cases correspond to spherical collapse in CDM with a rescaled Newton constant in the unscreened case, an initial top-hat density profile retains its shape during collapse. This allowed the authors to study the spherical collapse of a top-hat profile by simply comparing the evolution of a closed patch of Friedmann-Robertson-Walker (FRW) spacetime. However comparison to -body simulations showed that this model was too simple and missed the interesting regime between the two limiting cases Schmidt et al. (2009).

The above restrictions were dropped in Borisov et al. (2011), where the collapse of the top-hat was studied numerically by solving the modified gravity field equations. An important result of this work was the discovery that a top-hat profile develops shell crossing during its evolution. In order to alleviate this problem, a smooth transition region between the top-hat and the ‘background’ FRW spacetime was introduced. However the resulting values of showed dependence on the shape of the transition region, and did not lie in the expected range obtained in Schmidt et al. (2009).

In this paper we improve the spherical collapse calculation for models by using as initial condition the average density profile around a density peak. This is completely determined by the input cosmology, thus removing any ambiguity in the choice of initial profile and making the spherically symmetric setup as physically accurate as possible. The profile is calculated using peaks theory Bardeen et al. (1986) and the linear matter transfer function Lewis et al. (2000).

The paper is organized as follows. In Section II and III we review and obtain the relevant equations for spherical collapse in theories and explain the method of their numerical solution. Here special care is taken of the initial conditions and applicability of the quasistatic approximation, which is crucial due to the breakdown of Birkhoff’s theorem. We exhibit one of the main results of this work; a fitting function for the critical overdensity which the linearly extrapolated matter density has to reach in order to form a halo of mass at redshift within the model parameter range . In Section IV we review the excursion set formalism and extending into a drifting and diffusing barrier, modeling this way aspherical collapse, in order to obtain a realistic mass function in terms of the collapse redshift , mass of the object and modified gravity parameter . The accuracy of this mass function function is checked against Monte Carlo random walks in Section IV and we conclude in Section V.

## Ii Review of the Hu-Sawicki-Starobinsky \mathboldf(R) model

The action is given by

 S=Sm+12κ2∫d4x√−g(R+f(R)), (1)

where and a minimally coupled matter action. Variation with respect to the metric gives the Einstein field equations

 Gμν=e−φ(κ2Tμν−12gμνV(φ)+(eφ);μν−gμν□eφ), (2)

with energy momentum tensor and Einstein tensor . We have introduced the notation

 eφ≡1+f,R,V(φ)≡Reφ−R−f, (3)

where writing requires a form of such that can be inverted to give . The condition is required to ensure that the model remains ghost free Starobinsky (2007). Using (3), the trace of the field equations

 3□eφ=−2V+V,φ+κ2T (4)

allows us to interpret the fourth order differential equations for , Eqs. (2) and (3), as second order field equation (2) plus a second order equation (4) for the scalar . One can view this procedure as a first step in the direction of a Hamiltonian formulation Olmo and Sanchis-Alepuz (2011); Saltas and Hindmarsh (2011); Deruelle et al. (2009). It is also possible to make the replacement (3) already in the action (1)

 S=Sm+12κ2∫d4x√−g(eφR−V(φ)), (5)

where variation with respect to and the scalar leads to (2) and (3). The action (5) is kown as O’Hanlon theory O’Hanlon (1972). Since matter is minimally coupled, the energy momentum tensor of matter is conserved

 Tμ ν;μ=0, (6)

implying that the Euler and continuity equation for a perfect fluid take the same form as in GR. The same is true for the geodesic equation and thus for the general relativistic virial theorem Jackson (1970).1 In the following we consider only cold dark matter in the single stream approximation, such that the energy momentum tensor takes the form of pressureless perfect fluid

 Tμν=ϱuμuν,uμuμ=−1,T=−ϱ, (7)

with energy density and four velocity . This is justified as we only consider times well after matter radiation equality and choose initial conditions where shell crossing does not occur.

The functional form of is strongly restricted due to theoretical and phenomenological constraints Hu and Sawicki (2007); Antonio De Felice (2010), meaning that the resulting theory is free of classical and quantum instabilities and consistent with all known gravity experiments and observations. We will consider the class of functions proposed by Hu and Sawicki Hu and Sawicki (2007) and also Starobinsky Starobinsky (2007)

 f(R)= −2Λ+ϵn(4Λ)n+1Rn (8) ≃ −2Λ1+2ϵ/n (4Λ/R)n (9) ≃ −2Λ+2Λ⎛⎝1+[(n/2)1/nϵR4Λ]2⎞⎠−n/2, (10)

which can fulfill the above-mentioned constraints and are thus viable. In the cosmologically important region the three models above exhibit essentially the same behaviour. The Hu-Sawicki model (9) and the Starobinsky model (10) interpolate between and , where steepness is controlled by and the position of the transition is at Here is a constant energy scale whose value coincides with the measured value and is a small positive deformation parameter which is related to the more commonly used , via

 fR0≡|f,R(R0)|=ϵ(1+14(Ω−1Λ−1))−(n+1). (11)

Galactic rotation curves require and cepheids , such that viable theories and CDM have virtually indistinguishable expansion histories Hu and Sawicki (2007); Jain et al. (2012). Stability of the de Sitter vacuum together with solar system tests demands Capozziello and Tsujikawa (2008). Using the expansion history and baryon acoustic oscillations (BAO) gives complementary constraints on Martinelli et al. (2011). Note also that quantum corrections to the scalaron mass place strong constraints on Upadhye et al. (2012). Enforcing effectively introduces two additional energy scales; and . As we will see, is the range over which the effective potential of varies and is the squared mass of the scalaron, corresponding to small fluctuations () around the background field during the cosmic late time acceleration.

Now we derive the explicit form of for model (8), which will be assumed in the rest of this paper. For this we note first that

 f,R=−ϵ(4ΛR)n+1, (12)

such that in and after matter domination

 |f,R|≤ϵ≪1⇒f,R≃φ. (13)

From Eqs. (12) and (13) it follows

 R=4Λ(|φ|ϵ)−1n+1,f=−2Λ+ϵn4Λ(|φ|ϵ)nn+1. (14)

Finally becomes

 Missing or unrecognized delimiter for \right (15)

We can write the scalar equation (4) as

 Veff,φ≡−e−φ13(2V+κ2ϱ−V,φ)=□φ, (16)

where is the effective potential the scalar tries to minimize. Using it is given by

 Veff=ϵ4Λ3⎛⎜⎝(1+κ2ϱ4Λ)|φ|ϵ−n+1n(|φ|ϵ)nn+1⎞⎟⎠. (17)

The minimum of is at the field value ,

 φmin=−ϵ(1+κ2ϱ4Λ)−(n+1), (18)

which approaches zero for , see Fig 1. If occupies the minimum in this limit then GR is restored: the effective Newtonian constant returns to its GR value and the potential becomes . For small fluctuations around , the scalaron has a mass

 m2≡Veff,φφ=43(n+1)Λϵ(|φ|ϵ)−n+2n+1, (19)

which diverges at the General Relativistic limit .2 In this limit the fifth force is turned off since the effective interaction range of a Yukawa type interaction is given by the Compton wavelength .

Whether actually sits at the potential minimum in high density regions depends on the magnitude of spatial gradients Khoury and Weltman (2004). For galaxy clusters, if the forming cluster is too small or to large, -gradients cost too much energy and will stay near its background value even within the cluster, preventing the chameleon mechanism from operating. Large clusters and small on the other hand leave enough space for the scalar to change from to and hence the chameleon mechanism becomes active above a certain cluster mass scale. A necessary condition for the chameleon mechanism to activate is that the Compton wavelength must be much smaller than the size of the overdensity Hu and Sawicki (2007). In this case there is only a “thin shell” at the edge of the object that can mediate a large distance fifth force, from which the interior is unaffected and thus behaves like a typical overdensity in GR.

For overdensities of different magnitude and shape, and for different values of the model parameters () we encounter a time dependent mixture of all the above cases.

## Iii Spherical collapse

The spherical collapse model Peebles (1967); Gunn and Gott (1972) is a deterministic criterion which allows one to map an initially small, spherically symmetric overdensity to the formation of a virialized dark matter halo. More concretely this spherical collapse of an initial density profile allows one to estimate the formation time of a halo as function of the initial density amplitude , which can be inverted to give the threshold

 δi(zc)≡δ(zi,r=0|zc) (20)

for the initial density profile to collapse at redshift . Spherical collapse in an Einstein-de Sitter or CDM universe can be modeled analytically by assuming that the density is homogeneous within the perturbation. Due to Birkhoff’s theorem, the inner part is not influenced by the transition region and simply behaves as a closed FRW universe. The redshift of collapse of this patch measured in the flat exterior FRW then approximately equals the formation time of a bound virialized object Jenkins et al. (2001).

In the case of theories we actually need to solve the full field equations since Birkoff’s theorem does not apply. There are further complications which hamper the calculation in models. It was noticed in Li and Hu (2011) that halos are actually composed of subhalos, and this will increase the chameleon effect as screened subhalos attract each other less strongly than particles in a homogeneous dust cloud. Another environmental effect arises due to the fact that the forming cluster is itself a subcluster of a larger sized over/underdensity, enhancing/diminishing the chameleon effect Li and Efstathiou (2011); Lombriser et al. (2013) by diminishing/enhancing field gradients. In this work we do not consider these two effects, although part of the environmental dependence of the collapse threshold is taken into account by using the average density profile around a peak, which only depends on linear power spectrum . We comment further on the environment in III.3 and V.

As will be reviewed in Section IV the halo mass function depends on the spherical collapse derived quantity via

 νc≡δi(zc)σ(zi,R), (21)

which is -independent for an Einstein de Sitter universe and slightly -dependent in an CDM universe. The standard deviation is given by the linear matter power spectrum and filter function

 σ(zi,R)2=∫d3k(2π)3W(kR)2P(zi,k). (22)

For convenience one usually considers

 δc(zc)≡D(zc,zi)δi(zc), (23)

which defines the collapse threshold at redshift . This quantity is the linearly extrapolated density field, where the linear growth function was used to evolve from to . The introduced time evolution has no physical meaning but is convenient as is constant in an Einstein-de Sitter universe, and only weakly dependent on in CDM. The approximate - and -independence of leads to the universality of halo mass function if written as a function of . Despite the artificial time evolution introduced in the definition of (23), it is the collapse criterion (20) defined at the initial time that one should have in mind both for GR and gravity: the halo mass function is determined within the initial conditions and the information about formation time only enters via (20). Due to the practically identical expansion histories in CDM and the assumed model (8), the initial conditions can be assumed to be equivalent for both models. In particular is identical in both models if is chosen such that all relevant scales are still linear. Therefore in Eq. (21) only should be adjusted when we consider models.

Since the collapse criterion (20) is the quantity we wish to calculate, we can trivially rewrite the definition (21) of using the CDM growth function

 νc≡δi(zc,R)σ(zi,R)=δi(zc,R)Dσ(zi,R)D=δc(zc,R)σ(zc,R). (24)

This expression holds for both and CDM, however as we will see in Section IV, theories predict a that is a function of both and . This is already the case for ellipsoidal collapse in GR Sheth and Tormen (2002).

As a final point, let us emphasize that using the scale dependent linear growth factor instead of the GR equivalent in Eq. (24) would be both incorrect and inconvenient. Incorrect since linear growth in is scale dependent and therefore not multiplicative, and inconvenient as would no longer encode all of the -dependent deviations; instead one would need to provide . as defined in (20) and (23) is a convenient way to provide . Another convenient way would be to fix and to fold all the -dependence into a modified Li and Hu (2011).

### iii.1 Quasistatic equations

To obtain , we must calculate the collapse of a spherically symmetric pressureless matter distribution in an asymptotic FRW spacetime, such that the 3+1 dimensional problem simplifies to a 1+1 dimensional one. In GR the calculation is much simpler. Due to Birkhoff’s theorem an initially homogeneous (“top-hat”) overdensity retains its shape during collapse. This allows us to treat the size of the homogeneous overdense region as the scale factor of a closed FRW universe Weinberg (2008).

In theories the additional scalar degree of freedom allows for monopole radiation Sexl (1966), thus Birkoff’s theorem no longer applies. Another, more severe, problem is that in the linear regime of collapse the gravitational force is scale dependent due to mass of fluctuations. Finally, since the energy density becomes sufficiently large during collapse for the chameleon mechanism to take effect, the gravitational force will depend on the local density. As a result of these effects, an initial top-hat overdensity will not retain its shape during collapse and we cannot use a closed FRW to describe its collapse. Rather, we must solve the spherically symmetric field equations.

A spacetime which a admits a spherically symmetric spatial slicing has a metric which can written in the form Misner et al. (1973)

 ds2=−e2Φdt2+a2e−2Ψ(dr2+r2dΩ2), (25)

where both and are functions of and . For convenience we factor out , which will be the scale factor of the asymptotic flat FRW spacetime, where we choose without loss of generality. Note that this metric is fully nonlinear. We present the nonlinear field equations in Appendix A and derive the conditions under which one can assume that , and remain small even when the density becomes non-linear . Under these conditions, the set of relevant field and fluid equations reduces to the

 Poisson equation a−2\upDeltaΦ=23κ2¯ϱδ−16δV,φ, (26a) the nonlinear scalar field equation a−2\upDeltaφ=13(δV,φ−κ2¯ϱδ), (26b) the energy conservation \large.δ+1ar2∂r(r2(1+δ)v)=0 (26c) and the Euler equation \large.v+vH+vav′=−1aΦ′, (26d)

with as the radial velocity and the perturbation in the Ricci curvature. Note that we cannot assume in (26b) since we would then miss the effect of the chameleon mechanism. It is important to treat (26b) as a nonlinear partial differential equation, even though . This non-linearity makes solving the system (26) a nontrivial task. We explain details of the numerical methods in Section III.4.

### iii.2 Initial conditions

In addition to the relevant dynamical equations, we must also specify the initial conditions of the problem. Well after matter-radiation equality and well before the late time accelerated expansion, the Universe is in a state that is well described by a linearly perturbed Einstein-de Sitter spacetime on all scales relevant for large scale structure formation. At such early times modified gravity effects due to the scalar field are completely negligible due to the temporal chameleon effect as seen in Fig. 1.

We choose as our initial time for the spherical collapse. At this redshift, radiation is already sub-dominant relative to matter by a factor of order . Also the subhorizon assumption implicit in eq. (28) at is acceptable; the largest masses considered in this work are of the horizon size.3 Note that at this initial redshift, the radiation component is subdominant to matter but non-negligible. To evade any potential problems with normalisation of the power spectrum, we use CAMB Lewis et al. (2000) to obtain with the choice of cosmological parameters

 σ8 =0.8, Ωm =0.27, h =0.7, ns =0.96

and then evolve the general relativistic growth equation to obtain to the collapse redshift, neglecting radiation. Throughout this paper, we stress that all quantities calculated using linear theory are obtained from the standard general relativistic equations.

We discuss our choice of the initial density profile in the following section. Given , we can use the constraint and Poisson equations

 2Φ′H=−κ2¯ϱ(1+δ)av (27)
 2a−2\upDeltaΦ=κ2¯ϱδ (28)

to obtain and at ;

 vi(r)=−aiHir2∫r0δi(r′)r′2dr′. (29)

We use the the following natural boundary conditions at all times; at , and at the outer boundary .

### iii.3 Density profile

It was observed in Borisov et al. (2011) that the shape of an initial tophat density profile will evolve in theories, in contrast to the shape preserving evolution obtained in GR. Specifically, they found that a tophat profile develops a large spike near the boundary between the overdensity and background FRW spacetime. We confirmed this behaviour when using the initial density profile

 δi(r)=δi,02(1−tanh(r/rb−1s)), (30)

where is the size of the tophat-like function and determines the steepness of the transition, with leading to . Decreasing has the effect of forming a steeper spike at an earlier redshift. Fig. 2 shows the normalized density profiles for two different steepness parameters but same at different instances during collapse.

The formation of a spike, which signals shell crossing, prohibits the use of tophat like functions for numerical studies. More importantly, it is clear that the shape of the density profile can dictate whether the chameleon mechanism becomes active; this is indicated in Fig .3. Depicted here are the potential and field at different times and positions. We see that the minimum of the effective potential only determines the position of far outside the overdensity , where gradients are always small. In the center it is prevented from settling into the minimum because the necessary field gradients are energetically too costly. However for the steep profile (blue dots), finally turns around and moves to the potential minimum. Hence the collapse time of an overdensity with fixed and will depend on the shape of the density profile. Fig. 4 shows that the growth rates of the two profiles start to deviate once becomes nonlinear. While the scalaron is nearly screened for , it remains unscreened for , enhancing the growth. It is thus clear that the collapse time depends on the shape of the initial profile.

Due to the above subtleties, in our numerical calculations we use a physically motivated mean density profile around a peak of height

 δi(r,R)=⟨δ(zi,x,R)|peak,ν⟩ (31)

which is completely determined by the gaussian statictics of the smoothed linear density field Bardeen et al. (1986). In Appendix C we derive and display the explicit shape function, and in Fig. 5 we exhibit the function for various initial values but fixed . The mass contained in a spherical tophat

 M=4π3¯ϱ0R3, (32)

is used to define the mass of the final halo, where is the dark matter density at the present time.

The physical reason for the above mentioned shape dependence of spherical collapse is the same as the environment dependence taken into account in Li and Efstathiou (2011); Lam and Li (2012); Li and Lam (2012); Lombriser et al. (2013). In both cases the effectiveness of the chameleon mechanism is influenced by size of gradients in , which in turn depends on the density profile and its environment. The mean density profile (31) is an approximate way to take into account both effects: the actual mean shape close to the “size” of the collapsing protohalo and its mean environment for .

### iii.4 Numerical method

To numerically evolve the system of equations we start from the initial time slice and evolve the energy conservation and Euler equations over a single timestep, where we use -foldings as our time variable with a staggered leapfrog method to decompose the temporal and spatial derivatives. For sufficiently small and coarse grained radial coordinate (we define as the radial coordinate, with ) the simple finite difference scheme remains stable. Between and , we set the timesteps to be relatively large; , but for we refine the timesteps to to ensure that we can accurately model the effect of the chameleon mechanism.

Once we have evolved the fluid equations to the next timestep, we solve the equation (26b) using a very similar relaxation algorithm as outlined in Borisov et al. (2011); decomposing the non-linear equation into discretized form and Taylor expanding around the previous timestep. By solving the resulting large, yet sparse matrix equation, we update the solution and repeat until convergence is achieved. Once the field is calculated on the new timestep, its contribution as a source in the Poisson equation equation is evaluated, and the metric potential is obtained by a simple numerical integration. The process is then repeated until collapse is reached. We cannot evolve the system of equations formally to collapse. Here we simply evolve our system to an arbitrary high value in the non-linear over density; . Once the collapse redshift has been determined, we use Eqs. (20) and (23) to obtain .

To test that the solution obtained with our code is accurate, we check (completely independently of the code) that the data output solves the redundant momentum constraint and Newtonian gauge equations as given in Eqs. (). We also test our code by attempting to reproduce the standard General Relativistic values of (by using the algebraic relation as opposed to solving the non-linear equation). We find an error of less than in for collapsing objects in the redshift range . Finally, we test that the code is unaffected by modifying the number of points used in the Poisson equation integration, changing the asymptotic boundary and decreasing the timesteps by a factor of ten. All tests produced a deviation of less than in the resulting .

We perform over runs of the code, varying over initial density , field value and scale of the perturbation , choosing such that the overdensity collapses between and such that the mass of the object lies in the range relevant to clusters, .

### iii.5 Spherical collapse threshold for \mathboldf(R) gravity

The main results of this work are the collapse threshold and a realistic halo mass function as a function of , and , which we present in the following sections. In Section III.5 we will obtain as a fit function to our numerical results, and in Section IV we use this adding a drifting and diffusing barrier in the excursion set theory to obtain a realistic mass function .

The threshold for collapse will be a non-linear function of the model parameters and , and also the initial density (or correspondingly the collapse redshift ), and the mass of the overdensity (or equivalently the size of the overdensity, fixed by as discussed in section III.3). We wish to construct a fitting function for using as input data the values of obtained from our numerical simulations. In what follows we fix for simplicity.

We exhibit the behaviour of as a function of , and in Fig. 6. These figures are instructive as a qualitative check of the influence of modified gravity on collapse. Each panel corresponds to runs with fixed and each set of data points of the same color/shape correspond to runs with the same average collapse redshift. We observe a clear linear dependence between and for small and a redshift and -dependent break in this behaviour. This break is determined by from eq. (III.5). The dashed vertical line shows the mass as defined in (35) for which at . We also observe an approach to the GR value for increasing . Similarly, for field values close to the General Relativistic limit , we find the correct limit . The full lines show the fit function (III.5). We clearly observe a non-trivial and dependence and a return to GR for large objects and those with a high collapse redshift. An approximately linear relationship between and is observed for large values of , in agreement with the results of Li and Efstathiou (2011); Lombriser et al. (2013).

To quantify the effect of modified gravity we provide an interpolation function to fit the data. From Fig. (6) we can impose the following ansatz for

 δc(z,M,fR0) =δΛc(z){1+b2(1+z)−a3(mb−√m2b+1)+ +b3(tanhmb−1)} (33) mb(z,M,fR0) =(1+z)a3(log10[M/(M⊙h−1)]−m1(1+z)−a4) m1(fR0) =1.99log10fR0+26.21 b2 =0.0166 b3(fR0) =0.0027⋅(2.41−log10fR0) a3(fR0) =1+0.99exp[−2.08(log10fR0+5.57)2] a4(fR0) =(tanh[0.69⋅(log10fR0+6.65)]+1)0.11

The fit function converges separately for and to its GR limit , which can be approximated by Nakamura and Suto (1997)

 δΛc(z)≃3(12π)2/320(1−0.0123log10[1+Ω−1m−1(1+z)3]). (34)

We obtained our result (III.5) by considering as independent fit parameters for each value. For instance the various best fit parameters suggest a linear dependence on , see Fig. 7. In a similar fashion we obtain the other functional forms of . The parameter is of particular interest since it determines the position of the chameleon transition at , where changes its behavior from a linear growth in to a constant; see Fig. (6). Therefore roughly speaking the halo mass function at approaches CDM for masses larger than

 M1=1014.2(fR010−6)2M⊙h−1 (35)

due to the chameleon mechanism.

## Iv Halo Mass function: prediction for \mathboldf(R) gravity and deviation from GR

Dark matter halos result from the non-linear collapse of initial density perturbations. The abundance of these virialized structures depends on both the properties of the initial matter density field and the collapse threshold which leads to their formation. Following the seminal work of Press and Schechter (1974), the excursion set approach Bond et al. (1991) computes the abundance of dark matter halos as a function of their mass. The method involves smoothing the initial density field over different realisations and positing that once the overdensity encapsulated in a smoothing region is above a threshold criteria, the region will collapse. The key assumption is then to equate the fraction of collapsed comoving volume to the comoving density of halos . Thus the number density of haloes in the mass range , the halo mass function , is given by

 n(M)=f(σ)¯ϱ0M2dlnσ−1dlnM, (36)

where is the comoving background dark matter density and is related to fraction of collapsed volume. The fundamental quantity, which contains all information on the non-linear collapse dynamics, is . In what follows we first review the analytic derivation of in case of spherical GR collapse. We then extend this calculation to models with realistic collapse parameters. Having constructed the multiplicity function for these modified gravity models, we can provide an estimate of signatures in the cluster abundance. Note that our methodology is different to existing approaches in the literature Lam and Li (2012); Li and Lam (2012); Lombriser et al. (2013). In this work effects are taken into account by averaging the barrier over environments of the initial Lagrangian perturbations.

### iv.1 Halo mass function prediction for uncorrelated random walk and generic barrier

To estimate the fraction of collapsed volume, one has to compute the probability of having an overdensity smoothed on a scale . In the original Press-Schechter (PS) approach Press and Schechter (1974), assuming Gaussian initial conditions, the fraction of collapsed regions can be calculated analytically; it is given by

 F(R)=∫∞BΠ(δ,σ(R))dδ, (37)

where is the collapse threshold and the probability density function (PDF) is . However the PS approach suffers from the so called cloud in cloud problem: it requires an ad-hoc normalization of the mass function due to an incorrect counting of collapsed regions. To understand where the problem occurs let us review the standard excursion set procedure. We start by re-writing the smoothed overdensity on a scale at any random position as

 δ(R)=1(2π)3∫d3kW(k,R)~δ(z,k), (38)

where and are the Fourier transforms of the filter function and and the linearly extrapolated respectively. Since is a random quantity, it was shown in Bond et al. (1991) that its evolution follows a Langevin equation. Once we fix the filter, there is a one to one relationship between the smoothing scale , the mass of the halos and the variance defined as

 S≡σ2(z,R)=12π2∫dkk2P(z,k)W2(k,R). (39)

In the case of a sharp- filter and Gaussian initial conditions, the Langevin equation takes the form

 ∂δ∂S=ηδ(S), (40)

where is white Gaussian noise completely specified by its mean and variance . According to these equations, performs a random walk and its evolution between two scales and is determined by its previous step only. Since the system does not keep memory of previous steps, the dynamics corresponds to a Markovian random walk and the PDF follows a simple Fokker-Planck equation

 ∂Π∂S=12∂2Π∂δ2. (41)

The PDF is fully specified by two initial conditions. At , which corresponds to very large scales, the homogeneity of the universe implies that . If one relaxes the second condition then the solution of Eq. (41) is a Gaussian PDF corresponding to the original PS prediction.

In the Excursion set approach, when random walks cross the threshold collapse at scale a halo of mass is assumed to form. However random walks can cross more than once at different smoothing scales, and this can lead to double counting of halos. To evade this problem, one must remove walks when they cross for the first time. This can be encoded in an absorbing boundary condition; the PDF of uncollapsed objects is given by the solution of Eq. (41) with the second initial condition . An exact analytic solution for a barrier that is a generic function of the smoothing scale does not exist. However for a constant spherical collapse barrier the exact solution is given by Bond et al. (1991); Redner (2001):

 Π(δ,S)=1√2πS(e−δ2/(2S)−e−(2δc−δ)2/(2S)), (42)

where the first term on the right hand side is the previous Gaussian solution while the second term is known as the ‘anti-Gaussian’. The fraction of collapsed volume is then

 F(S)=1−∫δc−∞Π(δ,S). (43)

The first-crossing rate is given by . From the definition of the multiplicity function it follows

 f(σ)=√2πe−δ2c/(2σ2)δcσ, (44)

which is the original PS prediction with the correct normalisation.

The above calculation corresponds to spherically collapsing overdensities; the situation is considerably more complicated in the real Universe. The dynamics of collapse is aspherical and small over-dense regions require additional matter to collapse Sheth et al. (2001) since they are significantly affected by the surrounding shear field. Using ellipsoidal collapse in the excursion set approach introduces a stochastic barrier; this motivates the study of a generic barrier. In the CDM case a simple Gaussian distribution for the barrier with a mean value which drifts linearly as function of the variance is sufficient to reproduce the N-body halo mass function with high accuracy Corasaniti and Achitouv (2011, 2011); Achitouv and Corasaniti (2012, 2012). Furthermore this barrier is consistent with the overdensity required to collapse measured in the initial condition Achitouv et al. (2012) and has the advantage of admitting an exact solution for Markovian multiplicity function.

For gravity we have shown in Section III.5 that spherical collapse cannot be modeled using a linear barrier. To obtain an analytical prediction for using a generic barrier, we start by introducing the variable and assume that the barrier is described by a Gaussian PDF with mean value and variance , with constant . In such a scenario the Fokker-Planck equation for the variable is given by

 ∂Π(Y,S)∂S=1+DB2∂2Π(Y,S)∂Y2−d¯BdS∂Π(Y,S)∂Y (45)

In the special case where , the exact solution for is Corasaniti and Achitouv (2011, 2011)

 f(σ)=√2aπe−a¯B2/(2σ2)δcσ (46)

with . For generic , the solution of the Fokker-Planck equation without the absorbing boundary condition is simply given by a Gaussian with mean and variance . The crossing rate in this case would be given by

 F(S)=−ddS∫∞0√a2πSe−a(Y−¯B)2/(2S), (47)

leading to

 f(σ)=√2aπe−a¯B2/(2σ2)12σ(¯B−2Sd¯BdS). (48)

However, this expression does not have the correct normalisation since we did not solve the equations using an absorbing boundary condition. For constant barrier, one could correct this expression by multiplying by an ad-hoc factor two, however for a linear drift this would not be sufficient to recover the exact solution since there is no factor of two multiplying the first derivative of . In fact one can show that the factor of two in front of the first derivative of cancels once we add the anti-Gaussian term Corasaniti and Achitouv (2011). Thus the exact solution for a constant and linear barrier is given by

 f(σ)=√2aπe−a¯B2/(2σ2)1σ(¯B−Sd¯BdS). (49)

Similarly for a generic barrier one could approximate the exact solution by expanding in higher order derivatives of the barrier and dividing these terms by a factor of two. Thus we propose the following formula for a generic barrier

 f(σ)=√2aπe−a¯B2/(2σ2)1σ(¯B−σ2d¯Bdσ2++12∑n≥2(−σ2)nn!d¯Bnd(σ2)n). (50)

Note that this expression matches the first two terms of de Simone et al. (2011) which are equivalent to Sheth and Tormen (2002) for . In Sheth and Tormen (2002) it was proposed that the additional correction for a generic barrier is given by

 fST(σ)=√2πe−¯B2/(2σ2)1σ(5∑n=0(−σ2)nn!d¯Bnd(σ2)n), (51)

see also de Simone et al. (2011).

Note that (50) or (51) are approximations, and should only be applied to barriers well-described by algebraic functions for which all derivatives are well-defined. In order to test the robustness of our ansatz, we performed a series of Monte Carlo random walks for various barrier models following the procedure described in Bond et al. (1991). In the case where and , we exhibit some of the Monte Carlo random walks in Fig. 8. We find that for barriers which scale like and , the first order derivative term is sufficient to fit the Monte Carlo walks with high accuracy while for and we require the fourth/fifth term to obtain an accurate match. In such a case, is unable to reproduce the exact solution while Eq. (50) fits these Monte Carlo random walks with high accuracy. In what follows we will use Eq. (49) to model the spherical collapse barrier of but the method that we apply in section IV.2 to predict the halo mass function can be extended to non-standard GR where the barrier can be an arbitrary algebraic function of the form with and .

In the absence of N-body simulations, one way to evaluate the effect of gravity on the halo mass function is to use spherical collapse (ie: ) and measure the ratio between the GR and prediction for an uncorrelated walk (ie: sharp- filter). For that purpose we first need an accurate prediction for gravity. We run Monte Carlo walks for various parameters to test the accuracy of Eq. (49).

In Fig (9) we observe the halo mass function corresponding to the exact Monte Carlo solution (dot) and our prediction (full line) Eq. (49) with and given by Eq. (III.5). On the lower panel we show the relative difference between the Monte Carlo and theoretical prediction. We see that the difference is of order , confirming that Eq. (49) provides an excellent fit. The colours correspond to the different model parameters we test: blue, red, yellow and green correspond to ,