Calculation of the critical overdensity in the spherical-collapse approximation

TL;DR

This paper compares methods for calculating the critical overdensity in spherical collapse models, emphasizing the importance of precise definitions, especially in modified gravity theories like $f(R)$, and provides analytical and numerical results for these calculations.

Contribution

It clarifies the accuracy of the constant-infinity method and introduces an analytical expression for density contrast in $f(R)$ gravity models.

Findings

The constant-infinity method requires careful definition of numerical infinities.

Relative errors in $\, ext{delta}_c$, $n_{ ext{ln} M}$, and $N_{ ext{bin}}$ depend on redshift.

An analytical formula for density contrast as a function of collapse redshift and matter density parameter is derived.

Abstract

Critical overdensity $\delta_c$ is a key concept in estimating the number count of halos for different redshift and halo-mass bins, and therefore, it is a powerful tool to compare cosmological models to observations. There are currently two different prescriptions in the literature for its calculation, namely, the differential-radius and the constant-infinity methods. In this work we show that the latter yields precise results {\it only} if we are careful in the definition of the so-called numerical infinities. Although the subtleties we point out are crucial ingredients for an accurate determination of $\delta_c$ both in general relativity and in any other gravity theory, we focus on $f(R)$ modified-gravity models in the metric approach; in particular, we use the so-called large ($F=1/3$) and small-field ($F=0$) limits. For both of them, we calculate the relative errors (between our method and the others) in the critical density $\delta_c$, in the comoving number density of halos per logarithmic mass interval $n_{\ln M}$ and in the number of clusters at a given redshift in a given mass bin $N_{\rm bin}$, as functions of the redshift. We have also derived an analytical expression for the density contrast in the linear regime as a function of the collapse redshift $z_c$ and $\Omega_{m0}$ for any $F$.

Figures

Equations

$$S=\int d^{4}x\sqrt{-g}\left[\frac{R+f(R)}{2\kappa }+{\mathcal{L}}_m\right],$$

$$G_{\mu \nu }+f_R{R}_{\mu \nu }-\frac{1}{2}f{g}_{\mu \nu }-[{\nabla }_{\mu }{\nabla }_{\nu }-g_{\mu \nu }\square ]f_R=\kappa T_{\mu \nu },$$

$$f_{R}R-2f+3\square f_R-R=\kappa T,$$

$$d{s}^2=-(1+2\varphi )d{t}^2+a^2(t)(1+2\psi ){\delta }_{ij}d{x}^{i}d{x}^j.$$

$$\dot{\delta }+\frac{1}{a}(1+\delta )\vec{\nabla }\cdot \vec{v}=0\mathrm{and}$$

$$\dot{\vec{v}}+\frac{1}{a}(\vec{v}\cdot \vec{\nabla })\vec{v}+H\vec{v}=-\frac{1}{a}\vec{\nabla }\varphi ,$$

$$\ddot{\delta }+2H\dot{\delta }-\frac{4}{3}\frac{{\dot{\delta }}^2}{(1+\delta )}=\frac{1}{a^2}{\nabla }^2\varphi (1+\delta ).$$

$${\nabla }^2\varphi =\frac{16\pi G}{3}{a}^2\delta {\rho }_m-\frac{a^2}{6}\delta R(f_R),$$

$${\nabla }^2\delta f_R=\frac{a^2}{3}[\delta R(f_R)-8\pi G\delta {\rho }_m].$$

$$\ddot{\delta }+2H\dot{\delta }-\frac{4}{3}\frac{{\dot{\delta }}^2}{(1+\delta )}=\frac{3}{2}(1+\delta )H^2{\Omega }_m(t)(1+F)\delta$$

$$H^2(y)\frac{d^2\delta }{d{y}^2}+\left[2H^2(y)+H(y)\frac{dH(y)}{dy}\right]\frac{d\delta }{dy}-$$

$$-\frac{4H^2(y)}{3(1+\delta (y))}{\left(\frac{d\delta }{dy}\right)}^2=\frac{3}{2}(1+\delta )H^2(y){\Omega }_m(y)(1+F)\delta .$$

$$q(y\equiv \ln a)\equiv \frac{r}{r_i}-\frac{a}{a_i},$$

$$q^{\prime \prime }+\frac{H^{\prime }}{H}{q}^{\prime }$$

$$\sigma \equiv {\left(\frac{1}{q{a}_i/a+1}\right)}^3(1+{\delta }_i)-1.$$

$$p\equiv p(F)\equiv \frac{5}{4}\left(-1+\sqrt{1+\frac{24}{25}F}\right).$$

$${\delta }_i=C{a}_i^{1+p},$$

$${\delta }_c=C{a}_c^{1+p}{}_2{F}_1\left(ϵ(p),b(p);c(p);-\frac{(1-{\Omega }_{m0})}{{\Omega }_{m0}}{a}_c^3\right)$$

$${\delta }_c(a)=C{a}_c^{1+p}.$$

$$\frac{d^2\delta }{d{y}^2}+\frac{1}{2}\frac{d\delta }{dy}=\frac{3}{2}(1+F){\delta }_m,$$

$${\delta }_l(y)=C{e}^{y(1+p)},$$

$${\delta }_i=C{e}^{y_i(1+p)},{\delta }_i^{\prime }=(1+p){\delta }_i.$$

$$\delta_{c}=\left\{\begin{array}[]{lcl}\frac{3}{5}\left(\frac{3\pi}{2}\right)^{2/3}&,&{\rm for\quad}F=0\\ {}&{}\hfil\\ 1.70605&,&{\rm for\quad}F=1/3,\end{array}\right.$$

$$C(y_{c})=\left\{\begin{array}[]{ccl}\frac{3}{5}\left(\frac{3\pi}{2}\right)^{2/3}e^{-y_{c}}&,&{\rm for\quad}F=0\\ {}\hfil&{}\hfil\\ 1.70605\,e^{-y_{c}(1+p)}&,&{\rm for\quad}F=1/3.\end{array}\right.$$

$$\mathrm{Inf}(y_c)\equiv {\delta }_m({\Omega }_{m0}=1,y_c),$$

$${\delta }_c({\Omega }_{m0},y_c)$$

$${\Delta }^i\equiv \frac{{\delta }_c^i}{{\delta }_c^{\mathrm{pw}}}$$

$$n_{\ln M}=\frac{dn}{d\ln M}(M,z)=\frac{{\rho }_{m0}}{M}f(\nu )\frac{d\nu }{d\ln M}$$

$$f(\nu )=A\sqrt{\frac{2}{\pi }a\nu }[1+(a{\nu }^2{)}^{-q}]\exp [-a{\nu }^2/2)]$$

$$\sigma (M,z)={\sigma }_{80}D(z){\left(\frac{M}{M_8}\right)}^{-\gamma (M)/3},$$

Full text

I Introduction

Since the discovery that the expansion of the Universe is speeding up in $1998$ ref4 ; ref5 , several attempts have been made to understand the physical mechanism responsible for this cosmic acceleration. One possibility is that an exotic new component with negative pressure (dubbed dark energy) would be responsible for it. The simplest dark energy candidate is the cosmological constant ($\Lambda$). As an alternative to dark energy, one considers modification of general relativity (GR). The simplest modified-gravity candidate is the so-called $f(R)$ gravity in which the Lagrangian density ${\cal L}=R+f(R)$ is a nonlinear function of the Ricci scalar $R$. Here, we focus on the metric approach.

There is a general belief (see, for instance, Ref. Hurtado ) that cosmological kinematical tests (those that depend only on the background evolution) alone cannot distinguish between the standard cosmological constant model ($\Lambda$CDM) and any viable $f(R)$ gravity theory. Indeed, one of the strongest constraints comes from the growth of perturbations in the nonlinear regime Schmidt ; Santos . In the linear regime, on the other hand, only weak constraints are obtained ODweyer .

In GR, Birkhoff’s theorem holds, and in the spherical collapse (SC) approximation, an initial top-hat density profile keeps being a top hat. That straightforwardness enables analytical results in the Einstein-de Sitter (EdS) background ($\Omega_{m}=1$) —and it is therefore the standard benchmark for all more realistic initial conditions. In modified theories, however, the fifth force mediated by the new scalar degree of freedom—the scalaron Starobinsky —and the so-called chameleon mechanism camaleon play a crucial role. Indeed, the chameleon mechanism is a key ingredient to hide the fifth-force effects in high-density environments such as the Solar System and at Galactic scales. Consequently, the validity of the SC approximation itself has been the subject of a large dispute in the current literature. Borisov et al. f1 numerically solved the full modified gravity equations for the model proposed by Hu and Sawicki Hu and found that an initial top-hat profile develops shell crossing during its evolution, and therefore, its shape changes. An improvement of the SC numerical calculation for again the same $f(R)$ model is found in Ref. koop , using as an initial condition the average density profile around a density peak. Using the results for the SC found in Ref. camaleon , Lombriser et al. Lombriser and subsequently Cataneo et al. cataneo have taken into account the chameleon suppression of modifications in high-density regions.

Precisely in order to circumvent such problems and to gain some insight on the role played by the chameleon mechanism, we work in the so-called large- and small-field limits (see Sec II). Indeed, such an approach has been proven effective before: In Refs. ref36 ; ferrano , the density profiles and the linear bias of the cluster halos were determined in such limits, showing good agreement with N-body simulations.

The key quantity in the SC is the critical overdensity $\delta_{c}(z_{c})$ at a given collapse redshift $z_{c}$. It is defined as the final value (i.e, at redshift $z_{c}$) of the linear evolution of a given spherical top-hat initial perturbation that actually collapses at $z_{c}$ according to the full nonlinear equations. Theoretically, the latter value should be infinite, but in practice one has to deal with numerical infinities, and here, is where lies a potential problem. As we see later on, different numerical infinities give rise to different results for $\delta_{c}(z)$. The main goal of this work is to discuss and carefully quantify these differences.

This paper is structured as follows. Section II reviews the basic equations that describe the SC model in $f(R)$ theories. In Sec. III we review the current techniques to calculate the critical density contrast and point out the delicate step in the constant-infinity method. In Sec. IV we compare and discuss the results. In particular, we show the relative errors between our method and the others, in the case of $f(R)$ theories, in both the so-called small-field limit ($F=0$) and the large-field limit ($F=1/3$), when calculating the critical density $\delta_{c}$, the comoving number density of halos per logarithmic mass interval $n_{\ln M}$, and the number of clusters at a given redshift in different mass bins $N_{\rm bin}$. Throughout this paper, we assume the background is given by the standard $\Lambda$CDM model, since we assume there is no observable difference between that and the actual viable $f(R)$ evolution. In the Appendix we present the evolution of the density contrast in the linear regime for arbitrary $\Omega_{m0}$.

II Spherical collapse in $f(R)$ theories

In $f(R)$ theories the Einstein-Hilbert action is modified to

[TABLE]

where ${\cal L}_{m}$ is the Lagrangian of the ordinary matter, $\kappa\equiv 4\pi G$, and throughout this paper, we use $c=\hbar=1$. Variation of Eq. (1) with respect to the metric $g^{\mu\nu}$ yields the modified Einstein equations:

[TABLE]

where $G_{\mu\nu}$ is the Einstein tensor and $f_{R}\equiv\frac{df}{dR}$. Taking the trace of Eq. (2), one gets

[TABLE]

where $T$ is the trace of the energy-momentum tensor $T_{\mu\nu}$. From the latter, equation one can see that $f_{R}$ represents an extra degree of freedom. Indeed, one can show m18 that upon a change to Einstein’s frame, one recovers standard GR plus an extra scalar field.

The SC model considers a homogeneous and isotropic region (a top-hat profile) with density $\rho(t)=\bar{\rho}(t)+\delta\rho(t)$, where $\bar{\rho}$ is the background fluid density. We suppose that this region contains only nonrelativistic matter (both the pressure $p_{m}$ and the effective sound speed $c_{\rm eff}^{2}$ are negligible). Such a region can be described as a perturbation in an otherwise homogeneous universe with density $\bar{\rho}(t)$, scale factor $a(t)$, and Hubble parameter $H\equiv\dot{a}/a$, whose metric is given by

[TABLE]

In GR, in the considered case (since the anisotropic stress vanishes), the gravitational potential $\phi$ would be equal to the negative of the second potential ($\phi+\psi=0$). In modified theories, however, the extra scalar field acts as a source of the deviation between them.

The nonlinear continuity and Euler equations Peebles for the density contrast $\delta(t)\equiv\frac{\delta\rho}{\bar{\rho}}$ and velocity-field perturbation $\vec{v}$ are, respectively,

[TABLE]

[TABLE]

in comoving spatial coordinates. Combining Eqs. (5) and (6), one obtains a second-order differential equation for $\delta$:

[TABLE]

When gravity is modeled by $f(R)$ theories, the potential $\phi$ is modified accordingly, as follows. A perturbation in the matter density produces a perturbation in the metric $g_{\mu\nu}$, which can be translated into a perturbation in the Ricci scalar $R$. Using Eqs. (2) and (4), one gets the equation for the modified potential $\phi$:

[TABLE]

where $\delta R\equiv R-\bar{R}$ and $\bar{R}$ is the background Ricci scalar. Accordingly, the function $f_{R}(\bar{R})\equiv\frac{df}{dR}(\bar{R})$ is perturbed by $\delta f_{R}\equiv f_{R}(R)-f_{R}(\bar{R})$. Using (3), we get

[TABLE]

To obtain Eqs. (8) and (9), we have considered $|f_{R}(\bar{R})|\ll 1$ and the quasi-static approximation $\dot{f}_{R}(\bar{R})\ll|\vec{\nabla}f_{R}(\bar{R})|$. The former condition indicates that the background is similar to $\Lambda$CDM, while the latter assumes that the time scale of the collapse is much smaller than the time scale of the expansion of the Universe. Therefore, any time variation of the (background) scalar field $f_{R}$ is negligible in a typical time scale of collapse. Such an approximation is equivalent to focus on the evolution of the perturbations inside the Hubble radius when the background evolution is close to $\Lambda$CDM ($f_{R}\ll 1$), as has been shown in Ref. ref35 . However, in Ref. PhysRevD.89.084023 , it was shown that the deviation in the global matter power spectrum between static and nonstatic simulations is only $~{}0.2\%$. Therefore, although in principle the static approximation is not supposed to be accurate, the corrections are actually small.

Two opposite regimes, inherent to any $f(R)$, can be represented by a single factor $F$ ref36 . When the curvature in the spherical region is large, the fluctuations of the field $\delta f_{R}$ are very small, so that the Laplacian in Eq. (9) can be neglected, yielding $\delta R\approx 8\pi G\delta\rho_{m}$ and the usual Poisson equation is recovered; this defines the so-called small-field limit. On the other hand, if the curvature $R$ in the spherical region is similar to the background curvature $\bar{R}$, its fluctuation $\delta R$ is small. Thus, it can be neglected in Eq. (8), which increases the gravitational potential by a global $4/3$ factor on the right-hand side of that equation. This regime is the so-called large-field limit, i.e., $|\delta f_{R}|\sim|\phi|$. In the former case, gravity is not modified due to the chameleon effect; the opposite situation occurs in the large-field limit where gravity is strengthened, becoming more attractive.

For both the limits above, Eq. (7) can be cast as

[TABLE]

where $F=0$ reproduces the small-field case and $F=1/3$ corresponds to the large-field limit. It is convenient to write the above equation in terms of $y\equiv\ln a$ (we take the present value of the scale factor $a_{0}=1$):

[TABLE]

III Calculating the critical density contrast $\delta_{c}$

The critical density can be calculated following two different procedures: one directly from the time evolution of $\delta(y)$ given by Eq. (11) and another from the difference in evolution between the bubble radius and the background scale factor, as we see below. In the former, one has to deal with numerical infinities. We start with the latter, where one can circumvent this problem.

III.1 The differential-radius method

We define (following Ref. ref36 ) the differential radius

[TABLE]

where $r_{i}$ is the initial bubble radius when its scale factor is $a_{i}$ . The full nonlinear evolution equation for $q(y)$ obtained from mass conservation and Eq. (11) is

[TABLE]

where

[TABLE]

The initial conditions—usually set at a high redshift, when matter dominates ($\Omega_{m}\sim 1$)—are $q(y_{i})\equiv q_{i}=0$ and $q^{\prime}_{i}=-\delta_{i}(1+p)/(3(1+\delta_{i}))$, with

[TABLE]

One then solves Eq. (III.1) requiring that the collapse, defined by $q(z_{c})=-a_{c}/a_{i}$, takes place at a given redshift $z_{c}\equiv 1/a_{c}-1$. Such requirement constraints $q_{i}^{\prime}$ and, consequently, the value of $\delta_{i}\equiv\delta_{m}(a_{i})$ which we parametrize as

[TABLE]

inspired by the linear growth of $\delta_{c}$ at such high redshifts [see Eq.(17)].

The critical density $\delta_{c}$ is the linear evolution—determined by Eq. (A)—of such initial density perturbation $\delta_{i}$, given by Eq. (46) at the collapse ($a=a_{c}$):

[TABLE]

for any $\Omega_{m0}$ and $F$—see Eqs. (A) for definitions of $\epsilon(p)$, $b(p)$ and $c(p)$.

The explicit dependence on $\Omega_{m0}$ is usually taken for granted because the full evolution (i.e, until the collapse) is supposed to happen while matter still dominates. In that case there is no dependence of the hypergeometric function on the collapse scale factor $a_{c}$ for any value of $F$, since ${{}_{2}}F_{1}(\epsilon,b;c;0)=1\;\forall\{\epsilon,b,c\}$ and one recovers the standard power-law dependence

[TABLE]

Of course, for $p(F=0)=0$ and $\Omega_{m0}=1$, we recover the standard value $\delta_{c}=1.68647$, as expected. For $F=1/3$ and $\Omega_{m0}=1$, our results yield a constant $\delta_{c}=1.70605$ for any $z_{c}$, which agrees with the result given in Ref. ref36 for $z_{c}=0$.

III.2 The constant-Infinity method and its mending

As mentioned before, the full solution of Eq. (11) should be infinite at the collapse. However, since this equation can be solved only numerically, it is necessary to establish the infinite as a very large number for a given collapse redshift $z_{c}$. As we see below, this numerical infinity depends both on the collapse redshift $z_{c}$ and on $\Omega_{m0}$.

Here, we consider an initial value for the redshift $z_{i}=1000$, where the Universe is completely dominated by nonrelativistic matter ($\Omega_{m}=1$). We then have to deal with the linear version of Eq. (11) when the background behaves like an EdS universe:

[TABLE]

whose linear growing solution is given by

[TABLE]

where $C$ is again a constant which clearly depends on the collapse redshift. The initial conditions are then given by

[TABLE]

The constant $C$ is obtained requiring the collapse occurs at $y=y_{c}$, i.e, $\delta(y_{c})\to\infty$. To fix the numerical value of infinity, we use the known values of $\delta_{c}$ for an EdS universe:

[TABLE]

where the latter value was obtained in the previous section. The constant $C$ is then given by

[TABLE]

Evolving Eq. (11) with the initial conditions (20) and (22), the numerical infinity is defined by

[TABLE]

which clearly depends on $z_{c}$. Some papers Ronaldo1 ; infinito do not take this dependence into account and assume a constant value $\widetilde{\rm Inf}\equiv{\rm Inf}(y_{c}=0)=10^{5}$ or $10^{8}$. While both $\widetilde{\rm Inf}$ and ${\rm Inf}$ have roughly the same order of magnitude ($\sim 10^{5}$) at $z_{c}=0$, ${\rm Inf}(z_{c})$ decreases monotonically with $z_{c}$ and can be as low as $10^{4}$ at $z\sim 3$ (for both $F=0$ and $F=1/3$). Later on (see Figs. 2, 3, and 6), we point out the numerical differences in the final outcome from this approximation.

Once the so-called infinity ${\rm Inf}(y_{c})$ is established, we then use Eq. (11) for different values of collapse redshifts $z_{c}$ and for different values of $\Omega_{m0}$ in order to calculate the constants $C_{j}\equiv C(y_{c},\Omega_{m0})$ that satisfy the condition for the collapse $\delta_{m}(\Omega_{m0},y_{c})={\rm Inf}(y_{c})$. The corresponding value of $\delta_{c}$ is given by Eq. (46). If $F=0$, we get

[TABLE]

If $F=1/3$, we get

[TABLE]

IV Comparing the results from different approaches

In this section we plot the relative errors between our method and the previously mentioned ones, for both values $F=0$ and $F=1/3$, when one calculates the critical density $\delta_{c}$ (as a function of the collapse redshift), the comoving number density of halos per logarithmic mass interval $n_{\ln M}$ [see Eq. (27 for the definition], and the number of clusters at a given redshift in mass bins $N_{\rm bin}$, as a function of either the redshift or the mass interval.

If the value of $z_{c}$ is fixed, our method describes the evolution of $\delta_{c}$ in terms of $\Omega_{m0}$ (Fig. 1, left panel). On the other hand, when we fix the value of $\Omega_{m0}$, the method describes the evolution of $\delta_{c}$ as function of redshift $z_{c}$ (Fig. 1, right panel).

In the following figures we compare the results of our method with the differential-radius equation and with the method where the numerical infinity assumes large constant values (we pick $\widetilde{\rm Inf}=10^{5}$ and $10^{8}$ for the sake of comparison with previous results in the literature Ronaldo1 ; infinito ) for both values $F=0$ and $F=1/3$. We define the ratio between other method and our as

[TABLE]

where the superscript $i$ stands for the different methods in the literature, specified in the figures, and ${\rm pw}$ stands for the approach introduced in the present work. For the figures we always assume $\Omega_{m0}=0.3$.

In Fig 2 we show $1-\Delta^{i}$ as a function of the collapse redshift $z_{c}$. Note that the differential-radius method and ours are equivalent. On the other hand, the relative differences between our method (or the differential-radius method) and the constant-infinity one increase (in magnitude) with the collapse redshift because the critical density contrast calculated with our method stabilizes at the EdS value for large $z_{c}$ (as it is supposed to), while it grows unbounded in the latter approaches.

Next, we calculate and compare the mass function using the results for $\delta_{c}$ obtained by our method and by the method where the infinite is fixed (we pick $10^{5}$ and $10^{8}$ again). We use the Sheth-Tormen mass function Sheth , which analytically determines the distribution of these objects as a function of their virial masses and redshifts. The virial mass $M$ is defined such that the average density inside the virial radius $r_{v}$ is $\Delta_{v}$ times the critical density. We should mention that in order to relate SC to virialized halos, the virial theorem has to be modified in $f(R)$ theories. See, for instance, Ref. Schmidt for more on this topic.” In this formalism, the comoving number density of halos per logarithmic interval (in the virial mass $M$) is

[TABLE]

where $\rho_{m0}$ is the present matter density of the Universe, the peak threshold $\nu=\delta_{c}(z)/\sigma(M,z)$,

[TABLE]

and $\sigma$ is the variance of the linear contrast density field in spheres of a comoving radius $r$ containing the mass $M$.Also, $A$ is a normalization constant such that $\int f(\nu)d\nu=1$. Here, we use the following approximation Viana for $\sigma$:

[TABLE]

where

[TABLE]

and shape parameter $\Gamma$ is given by Viana

[TABLE]

Above, $M_{8}=5.95$ x $10^{14}\Omega_{m0}h^{-1}M_{\odot}$ is the mass inside a sphere of radius $8h^{-1}\mbox{Mpc}$, $h$ is the current Hubble constant in units of $100{\rm~{}km~{}s}^{-1}{\rm Mpc}^{-1}$, $\Omega_{b}=0.02230/h^{2}$ is the baryonic density parameter, and $D(z)$ is the growing solution of Eq. (36) normalized in $z=0$. We assume that $\sigma_{80}=0.8159$ ref19 , $q=0.3$, and $a=0.707$ Sheth . We have checked that if, instead of the fit given by Eq. (29)—that is numerically simpler to deal with—we had used a more accurate approach, our results would not change significantly.

In Fig. 3, we show the relative error

[TABLE]

for $z=0$ between our method and the others when $\Omega_{m0}=0.3$ for both values $F=0$ and $F=1/3$.

In Fig. 4 we show relative errors as a function of the redshift when the virial mass takes the value $10^{13}h^{-1}M_{\odot}$ (upper panels) and $10^{15}h^{-1}M_{\odot}$ (lower panels). From these figures, it is clear that the relative errors are more significant for larger virial masses and higher redshift.

The observed quantity, however, is the number of clusters at a given redshift and in a given mass bin. It is defined as

[TABLE]

where $dV$ is the comoving volume at redshift $z$, $\frac{dV}{dz\,d\Omega}=\frac{r^{2}}{H(z)}$, and the comoving distance $r(z)$ is given by

[TABLE]

Figure 5, obtained with our method, shows $N_{\rm bin}$ for $F=0$ and virial masses in the ranges $10^{13}-10^{14}$ $h^{-1}M_{\odot}$ (left panel) and $10^{14}-10^{15}$ $h^{-1}M_{\odot}$ (right panel). Similar results are obtained by using $F=1/3$. Note that for higher-mass bins, the peak in the $N_{\rm bin}$ is lower and located at a smaller $z$. From this piece of information and from our results above—both $\Delta^{i}$ and $\Delta_{n}^{i}$ increase with $z$—we expect the bin $10^{13}-10^{14}$ $h^{-1}M_{\odot}$ to be the most sensitive one to the small differences in $\delta_{c}$ and $n_{\ln M}$ we have been pointing out.

We then define one last relative error, namely,

[TABLE]

and plot it, for virial masses between $10^{13}$ and $10^{14}$ $h^{-1}M_{\odot}$, in Fig. 6 for $F=0$ (left panel) and $F=1/3$ (right panel). Notice that, in the redshift range of interest, it is at most of the order of $0.01$.

V Conclusions

In summary, we have shown that our approach matches the results from the differential-radius method and pointed out that the so-called constant-infinity method does need corrections in the calculation of a key quantity, namely, the critical density $\delta_{c}$. We have also derived an analytical expression for the critical density $\delta_{c}$ as a function of $\Omega_{m0}$, $z_{c}$, and $F$.

In spite of being more rigorous and more accurate than the constant-infinity method used in the literature, we should mention that, for the the current stage of the observations, the procedure presented here does not yield observable differences in the cluster number count $N_{\rm bin}$. The small discrepancies pointed out (less than $1\%$) may be of use in the future, when more precise data become available.

We remark that our results are based on the validity of the SC approximation and of the Sheth and Tormen prescription. It would be interesting to compare the differences we found, due to distinct ways of calculating $\delta_{c}$, using a semianalytic approach, with those from N-body simulations. Besides, as shown in Ref. Burrage , departures from spherical symmetry affect chameleon screening and a detailed comparison of semianalytical methods and simulations are required to determine the correct functional form of the mass function. This is an important task to have in mind in the upcoming large-scale surveys.

ACKNOWLEDGEMENTS

D.H. acknowledges financial support from CAPES.

Appendix A SOLUTION OF THE DIFFERENTIAL

EQUATION FOR THE DENSITY CONTRAST IN THE LINEAR REGIME

Consider the linear equation for the perturbations

[TABLE]

where $\dot{~}{}\equiv d/dt$. Considering the change of variable $t\rightarrow a(t)$ and using the equation for the background in the standard $\Lambda$CDM scenarium (neglecting radiation), namely,

[TABLE]

then Eq. (36) is written as

[TABLE]

where ${}^{\prime}\equiv d/da$. In an EdS universe, the above expression becomes

[TABLE]

Suppose a solution to (39) of the form

[TABLE]

where $C$ is a constant. Replacing this solution in Eq. (39) yields

[TABLE]

whose nontrivial solution is given by

[TABLE]

Now consider a solution to Eq. (A) of the form

[TABLE]

Substituting this solution in Eq. (A), we get

[TABLE]

Making the change of variable

[TABLE]

one can recognize Eq. (A) as an hypergeometric differential equation. Thus, the growing solution of Eq. (A) is

[TABLE]

where ${{}_{2}}F_{1}(\epsilon,b;c;u)$ is the hypergeometric function, and

[TABLE]

For $F=0$, one obtains the solution ioav

[TABLE]

and for $F=1/3$, we get

[TABLE]