EDP Sciences
Free Access
Issue
A&A
Volume 607, November 2017
Article Number A47
Number of page(s) 6
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/201731467
Published online 10 November 2017

© ESO, 2017

1. Introduction

Over the course of the past several decades, cosmologists using a large number of observations came up with a model describing the structure and evolution of the universe, dubbed ΛCDM model. In this model the Universe is constituted by cold dark matter (CDM), and vacuum energy (represented by the cosmological constant Λ). This model fits a large number of data (Del Popolo 2007; Komatsu et al. 2011; Del Popolo 2013, 2014; Planck Collaboration XIII 2016), but suffers from drawbacks on small scales (see Del Popolo & Le Delliou 2017), the fine tuning problem (Weinberg 1989; Astashenok & del Popolo 2012) and the cosmic coincidence problem. Another fundamental test that the ΛCDM model has to pass is to accurately predict the dark matter (DM) haloes distribution (i.e., the halo mass function (MF; see Del Popolo & Yesilyurt 2007; Hiotelis & Del Popolo 2006, 2013)). The high mass end of the MF at small redshift (z ≤ 2) is very sensitive to cosmological parameters like the Universe matter and dark energy (DE) content (Ωm and ΩΛ), the equation of state of the Universe, w, and its evolution (Malekjani et al. 2015; Pace et al. 2014). At redshifts higher than the previously quoted ones, the MF is of fundamental importance in the study of the reionization history of the universe (e.g., Furlanetto et al. 2006), quasar abundance (e.g., Haiman & Loeb 2001), and to study the distribution of DM.

Press & Schechter (1974, hereafter PS) proposed a very simple model based on the assumption of Gaussian distribution of the initial density perturbation, and the spherical collapse model. The quoted approach has the drawback of overpredicting the number of objects at small masses, and underpredicting those at high mass (e.g., Jenkins et al. 2001; White 2002). The extended-PS formalism, or excursion set approach (Bond et al. 1991; Bower 1991; Lacey & Cole 1993; Gardner 2001), introduced to overcome the quoted problems, was unable to solve them.

Extension of the quoted formalism (Del Popolo & Gambera 1998, 1999, 2000; Sheth et al. 2001), moving from the spherical collapse to non-spherical collapse gave much better agreement with N-body simulations (Sheth & Tormen 1999, hereafter ST). However, a deeper analysis of ST, and Sheth et al. (2001) showed that the ST MF overpredicts the halo number at large masses (Warren et al. 2006; Lukić et al. 2007; Reed et al. 2007; Crocce et al. 2010; Bhattacharya et al. 2011; Angulo et al. 2012; Watson et al. 2013), and when the redshift evolution is studied the situation worsen (Reed et al. 2007; Lukić et al. 2007; Courtin et al. 2011).

Another important issue is that of the universality of the MF, namely its independence on cosmology and redshift. Several studies (e.g., Tinker et al. 2008; Crocce et al. 2010; Bhattacharya et al. 2011; Courtin et al. 2011; Watson et al. 2013) showed that the MF is not universal nor in its z dependence or for different cosmologies.

In the present paper, we want to show how the excursion set approach can be improved to the extent that it can produce a MF in good agreement with N-body simulations like that of (Tinker et al. 2008) who showed clear evidences of the MF deviations from universality, calibrated the MF at z = 0 in the 1011<M< 1015h-1M mass range within 5%, and found the redshift evolution of the same.

The paper is organized as follows. In Sect. 2, we discuss the stochastic process, and Sect. 3 is devoted to results and discussion.

2. The stochastic process

As already reported, the excursion set model is based on the ideas of Press & Schechter (1974) and on their extensions which are presented in the pioneered works of Bond et al. (1991) and Lacey & Cole (1993). We will improve on these ideas in this paper but before that we will write some useful relations about the density fields and the smoothing filters which will be used in what follows.

The smoothed density perturbation at the center of a spherical region is (1)where is the density at distance r from the center of the spherical region and Wf is a smoothing filter. Reducing R, the variable δ executes a random walk that depends on the form of the density field and on the smoothing filter Wf. If the density field is Gaussian, then δ is a central Gaussian variable and its probability density is given by (2)For a spherically symmetric filter, the variance at radius R is given by (3)where is the Fourier transform of the filter and P is the power spectrum.

The correlation of values of δ between scales is given by the autocorrelation function that is (4)Since S is a decreasing function of R and R an increasing function the mass M contained in the sphere of radius R, then, S can be considered as a function of mass.

The most interesting filter, because of its obvious physical meaning, is the top-hat in real space given by (5)where H is the Heaviside step function. The Fourier transform of the filter is given by, (6)In this paper we assume a Gaussian density field and the top-hat filter in real space. We used a flat model for the Universe with present day density parameters Ωm,0 = 0.3 and , where Λ is the cosmological constant and H0 is the present day value of Hubble’s constant. We have used the value H0 = 100 h KMs-1 Mpc-1 and a system of units with munit = 1012Mh-1, runit = 1 h-1 Mpc and a gravitational constant G = 1. In these units, H0/Hunit = 1.5276. Regarding the power spectrum, we employed the ΛCDM formula proposed by (Smith et al. 1998).

The stochastic process is defined as follows. We assume that(7)where K is a kernel and dW is the usual Wiener process. Thus, in the plane (S,δ) we have a random walk. We assume a kernel of the simple form (8)for uS and zero otherwise. Substituting in Eq. (7) and integrating by parts we have (9)Thus δ is a linear combination of a Wiener process, W(S), and an average integrated Wiener process,. For a = 0 the variable δ describes a Wiener process in which the value of δ(S + ΔS) depends only on the value of δ(S) and not on previous values, δ(S + ΔS) = δ(S) + c(1−aW(s). This is because, according to the definition of Wiener process, at every step the increment ΔW(s) is chosen from a central Gaussian with variance ΔS. Thus the steps of a walk on the (S,δ) plane are uncorrelated. For a ≠ 0 the second term in the right hand side of Eq. (9) includes information from all positions of the walk up to S and results to a correlation between steps. Obviously δ is a central Gaussian as a sum of central Gaussians.

The autocorrelation between scales is found multiplying δ(S) by δ(S′) and finding the expected value of the product taking into account the following property of Wiener integration, see for example Jacobs (2010) (10)For S′ ≤ S we have (11)Then, (12)From Eqs. (2) and (3) we have the condition c2 [a2/ 3−a + 1] = 1. Then, Eq. (11) can be written as (13)where .

It is reported in Eq. (90) of Maggiore & Riotto (2010), and is confirmed by our calculations that for the top-hat filter the predictions of Eq. (13) are in good agreement with those of Eq. (4) for values of λ close to 0.5. In Fig. 1 we give an example. The prediction of Eq. (4) and the prediction of Eq. (13) for λ = 0.45 that corresponds to a = 1.185 or a = 0.796 are plotted. This choice for the values of a, and consequently of c, gives the correct correlation between scales which is very important in describing the accurate evolution of random walks. This evolution corresponds completely to the power spectrum and the smoothing filer used and first distributions which will be presented below are fully consistent with the idea of the excursion set model. Since the excursion set model and the first crossing distribution are inseparably linked, an accurate evaluation of the first crossing of a barrier by the above random walks is essential.

An interesting quantity which correlates the past the present and the future of a time evolving stochastic process x(t) is defined by (14)In our problem this quantity gives a correlation between the steps of random walks for various values of S. Using Eq. (13) we have (15)Obviously when λ = 0 also C = 0 and the steps are uncorrelated. This corresponds to a Wiener process. For positive values of λ, steps are positive correlated for Δ < 2S (persisting walks) and negative correlated for Δ > 2S (anti-persisting walks). We note that for the fractional Brownian motion, which is a procedure with correlated steps, the persisting case corresponds to a Hurst exponent H> 1 and the anti-persisting to H< 1, (see for example Hiotelis & del Popolo 2013). However, roughly speaking, the procedure studied above looks like a fractional Brownian motion with varying H.

thumbnail Fig. 1

Predictions of Eqs. (4) and (13), solid line and dashed line respectively, for λ = 0.45.

Open with DEXTER

3. Results and discussion

We discretize Eq. (9) by dividing the mass interval [Mmin,Mmax] = [10-3,105.5] Munit into n intervals of equal length in logarithm spacing. N tracer particles are considered and ΔWi,i = 1,2...n values for each tracer particle are chosen from central Gaussians with respective variances S(Mi−1)−S(Mi). Then, δS(i) is calculated according to Eq. (9). The first crossing of the constant barrier δc = δ(z) is found for every tracer particle. We recall that δ(z) is the linear extrapolation up to present of the overdensity of a spherical region which collapses at redshift z (Peebles 1980). The number ni of particles which have their first upcrossing of the barrier between Si−1Si are grouped and the first crossing distribution is calculated by f(Si) = ni/ [N(Si−1Si)]. Finally, the mass function is calculated by 2Sif(Si).

thumbnail Fig. 2

Role of the kernel in amplifying the walk of a tracer particle. Smooth solid, noisy-dashed and noisy-solid lines correspond to a = 1, a = 0 and a = 3/2 respectively.

Open with DEXTER

In Fig. 2 we show the paths of the same tracer particle for a = 1 and for the cases a = 0 and a = 3/2. These paths are represented in the figure. It is clear that the path with a = 1 is much smoother, as expected. Consequently values of a control the degree of smoothness which is related to the distribution of values of δ(S + dS) for given δ(S). So the values of a define the correlation between various scales and, as we have shown above, a proper choice of these values results to autocorrelation functions which approximate very satisfactory the results of Eq. (4).

Before studying mass functions we have checked the reliability of our Monte Carlo approximation by testing our results with analytical solutions. For a = 0 or a = 3/2 the procedure is a Wiener process and the first crossing distribution is given by the inverse Gaussian (16)In Fig. 3, the prediction of Eq. (16) is plotted together withthe prediction of our Monte carlo approximation. The horizontal axis is ν = δc(z)/ for z = 0. We note that this is a test with only numerical interest. The results presented in this figure are derived for n = 1000 and N = 5 × 105. Our results are also compared with the interesting analytical predictions of Maggiore & Riotto (2010) and Musso & Sheth (2012).

In Maggiore & Riotto (2010) the authors use a path integral approach to estimate the first crossing distributions and their results have the form of infinite series which converge slowly. Their approximation is fully consistent with the idea of the excursion set model. The resulting mass function is approximated by the formula,

thumbnail Fig. 3

Comparison of the results of our Monte Carlo method with the exact solution for the cases a = 0 or a = 3/2. The agreement is satisfactory.

Open with DEXTER

(17)where , (see Eq. (120) in Maggiore & Riotto 2010) .

On the other hand, in the approximation Musso & Sheth (2012) the condition of the first up-crossing is replaced by a condition of any up-crossing, while these two conditions are obviously not equivalent. Additionally, a bivariate joint distribution between δ and v ≡ dδ/ dS is assumed, a choice which is unjustified. The mass functions is approximated by (18)where (19)where Γ depends on the power spectrum and the kernel used to smooth the density field.

In Fig. 4, we compare the predictions of our results, derived for a = 0.796, with those of analytical formulae of Eqs. (17) and (18) with the predictions of N-body simulations at z = 0. In both snapshots, squares are the predictions of N-body simulations of Tinker et al. (2008).

thumbnail Fig. 4

Comparison of the results of N-body simulations, squares, with those of analytical formulae of Eqs. (17) and (18), and with our predictions. Left snapshot: squares are the results of N-body simulations. The predictions of Eq. (17) are represented by the smooth solid line while the predictions of Eq. (18) for Γ = 1/3 and Γ = 1/2 are represented by the dashed lines (large dashes and small dashes respectively). Right snapshot: squares are the results of N-body simulations. The results of our Monte Carlo simulations are given by the solid line.

Open with DEXTER

The prediction of Eq. (17) and the predictions of Eq. (18) for Γ = 1/3 and for Γ = 1/2 are plotted in the figure. In the right snapshot the prediction of Monte Carlo approximation. Analytical formulae of Eqs. (17) and (18) result to smaller numbers for heavy haloes and larger numbers for smaller haloes compared to the results of N-body simulations, while our approximation gives the correct behavior of the mass function for small haloes, ν ≤ 1. This is an interesting result. It shows that the excursion set model works very satisfactory for small haloes with ν ≤ 1. This agreement has not been reported elsewhere. On the other hand, in agreement with the analytical formulae studied above, our approximation fails to produce the correct number of heavier haloes but for ν ≥ 1.1 our results coincide with those of Eq. (17) and those of Eq. (19) (for Γ = 1/3). Definitely, the problem of the approximation of the correct first crossing distribution by a simple analytical formula has not been solved yet but we believe that the simplicity of the approximation formula is a secondary issue. The important issue is that of the accurate evaluation of the first crossing distributions at various scales. Our results and those of the analytical formulae of Eqs. (17) and (18) indicate that the predictions of the excursion set model are, in any case, for heavy haloes far from the results of N-body simulations at least for the case of the top-hat filter and the constant barrier.

It seems possible that at large scales, additional parameters may be taken into account, as for example the ellipticity or the angular momentum of the structures. This could lead to think that the use of moving barriers as those proposed in the literature, (Bond & Myers 1996; Del Popolo & Gambera 1998; Sheth et al. 2001; Sheth & Tormen 2002) is necessary. Moving barrier models try to solve the problem of the overestimation of the number of small haloes and the underestimation of the number of small haloes using a critical threshold of collapse which varies with mass (S). The choice of a moving barrier is based on physical arguments since larger haloes appear with the larger ellipticity and larger angular momentum. A smaller critical threshold of collapse for large haloes rearranges first crossing distributions at various scales. Studying various cases of moving barriers we found that an agreement with the predictions of N-body simulations can be achieved using a constant, lower threshold for collapse. We used a barrier of the form δ(z) = pδc(z) where p is a constant. In Fig. 5 we present a comparison with the results of N-body simulations for two different values of p at z = 0. An agreement is shown. It is interesting to note the sensitivity of the distribution of heavy haloes to the values of p. We note that the use of a lower threshold results to an increase of the fraction ncross/N where N is the total number of walks studied and ncross is the number of walks which have passed the threshold in the range SminSmax, but this increase is larger for small values of S (larger haloes).

We used p = 0.866 as the best value. We calculated mass functions for redshifts z = 0,z = 1.25 and z = 2.5 and we presents them in Fig. 6. Our results are plotted inthe figure as are those derived from the fitting formula of Tinker et al. (2008) given in their Eq. (3) which is (20)where and the z-dependence is given by (see Eqs. (5)–(8) in Tinker et al. 2008). We used Δvir = 200.

In Fig. 7 we show the fractional error, defined by where MF is the mass function derived by our model, for z = 0 and z = 2.5.

We note that according to Tinker et al. (2008) their model described by Eq. (20)is valid for 0 ≤ z ≤ 2.5. However in order to find the cause of increasing difference between our results and those of Tinker et al. (2008), shown in Fig. 7, we present comparisons with some other analytical mass functions available in the literature, such that of

A. Watson et al. (2013) which is valid for 0 ≤ z ≤ 30. This is given by (23)where (24)

B. The formula of Warren et al. (2006), which is, (25)and C. that of Sheth et al. (2001)

(26)where A = 0.3222,as = 0.707 and ps = 0.3.

thumbnail Fig. 5

Comparison of the results N-body simulations, squares, with our results. The thin solid line solid corresponds to p = 0.894 and the thick one to p = 0.836.

Open with DEXTER

thumbnail Fig. 6

Comparison of our results, open squares with those of N-body simulations. The results of N-body simulations are represented by the fitting formula of Tinker et al. (2008) and are shown by the thick solid lines. Our results have been derived for a = 0.796 and p = 0.866.

Open with DEXTER

thumbnail Fig. 7

Fractional error [log (MF)−log (MFT)]/log (MFT) for two redshifts z = 0 (dotted-line) and z = 2.5 (solid-line). MFT is given by the fitting formula of Tinker et al. (2008) while MF represents our results.

Open with DEXTER

thumbnail Fig. 8

Upper snapshot: mass functions at redshift z = 5. The noisy solid line shows our results (a = 0.796,p = 0.866). The smooth solid line shows the model of Sheth et al. (2001). Dashed line represents the results of Warren et al. (2006) while dotted line shows the results of Watson et al. (2013). Lower snapshot: fractional errors between our results and the model of Sheth et al. (2001), solid line, between our results and the model of Warren et al. (2006), dashed line, and between our results and the model of Watson et al. (2013), dotted line.

Open with DEXTER

The comparisons for redshift z = 5 are shown in Fig. 8. We show that the agreement remains satisfactory.

We also note that for large redshifts resolution problems appear. This is because δc(z) is an increasing function of z and thus the percentage of walks which pass the barrier becomes smaller for large z. Large structures are more rare and the mass function appears noisy (see at the right side of Fig. 8). However, it becomes difficult to check if a disagreement is due to the physical process or to the poor resolution. This difficulty rises the challenge of finding an analytical solution for the first crossing distribution of the process of Eq. (9).

In Fig. 9 we present a comparison of our results with the formula of Watson et al. (2013) and for z = 10. These results are derived for n = 400 and N = 5 × 106. Resolution problems are obvious since only 4365 from N tracer particles cross the barrier, but the agreement remains satisfactory.

thumbnail Fig. 9

Mass functions at redshift z = 10. Noisy solid line shows our results (a = 0.796,p = 0.866). The dashed line shows the model of Watson et al. (2013).

Open with DEXTER

It is well known that the process of structure formation is a very complex one. It has been studied extensively in the literature and more than thirteen formulae for the mass function has been proposed for various cosmological models, various halo finding algorithms and various mass scales, see for example Watson et al. (2013) and references therein. Consequently, the probability of constructing an analytical approach that predicts the results of N-body simulations is extremely small. However, our results show that the stochastic process studied here, is not a N-body simulation, which is however able to shed more light to the physical process during the formation of structures. Since it is a process which describes accurately the correlation between scales for the realistic top-hat filter and is able to produce results close to these of N-body simulations for a constant barrier, deserves a more profound study. Any alternative approximation of first crossing distributions, resulting from the above described stochastic process, should be very interesting.

Acknowledgments

We acknowledge Dr. Andromachi Koufogiorgou for her kind help and J.Tinker for making available the results of their N-body simulations. A.D.P. was supported by the Chinese Academy of Sciences and by the President’s international fellowship initiative, grant No. 2017 VMA0044.

References

All Figures

thumbnail Fig. 1

Predictions of Eqs. (4) and (13), solid line and dashed line respectively, for λ = 0.45.

Open with DEXTER
In the text
thumbnail Fig. 2

Role of the kernel in amplifying the walk of a tracer particle. Smooth solid, noisy-dashed and noisy-solid lines correspond to a = 1, a = 0 and a = 3/2 respectively.

Open with DEXTER
In the text
thumbnail Fig. 3

Comparison of the results of our Monte Carlo method with the exact solution for the cases a = 0 or a = 3/2. The agreement is satisfactory.

Open with DEXTER
In the text
thumbnail Fig. 4

Comparison of the results of N-body simulations, squares, with those of analytical formulae of Eqs. (17) and (18), and with our predictions. Left snapshot: squares are the results of N-body simulations. The predictions of Eq. (17) are represented by the smooth solid line while the predictions of Eq. (18) for Γ = 1/3 and Γ = 1/2 are represented by the dashed lines (large dashes and small dashes respectively). Right snapshot: squares are the results of N-body simulations. The results of our Monte Carlo simulations are given by the solid line.

Open with DEXTER
In the text
thumbnail Fig. 5

Comparison of the results N-body simulations, squares, with our results. The thin solid line solid corresponds to p = 0.894 and the thick one to p = 0.836.

Open with DEXTER
In the text
thumbnail Fig. 6

Comparison of our results, open squares with those of N-body simulations. The results of N-body simulations are represented by the fitting formula of Tinker et al. (2008) and are shown by the thick solid lines. Our results have been derived for a = 0.796 and p = 0.866.

Open with DEXTER
In the text
thumbnail Fig. 7

Fractional error [log (MF)−log (MFT)]/log (MFT) for two redshifts z = 0 (dotted-line) and z = 2.5 (solid-line). MFT is given by the fitting formula of Tinker et al. (2008) while MF represents our results.

Open with DEXTER
In the text
thumbnail Fig. 8

Upper snapshot: mass functions at redshift z = 5. The noisy solid line shows our results (a = 0.796,p = 0.866). The smooth solid line shows the model of Sheth et al. (2001). Dashed line represents the results of Warren et al. (2006) while dotted line shows the results of Watson et al. (2013). Lower snapshot: fractional errors between our results and the model of Sheth et al. (2001), solid line, between our results and the model of Warren et al. (2006), dashed line, and between our results and the model of Watson et al. (2013), dotted line.

Open with DEXTER
In the text
thumbnail Fig. 9

Mass functions at redshift z = 10. Noisy solid line shows our results (a = 0.796,p = 0.866). The dashed line shows the model of Watson et al. (2013).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.