Excursion set peaks: a self-consistent model of dark halo abundances and clustering

Abstract

We describe how to extend the excursion set peaks framework so that its predictions of dark halo abundances and clustering can be compared directly with simulations. These extensions include: a halo mass definition which uses the TopHat filter in real space; the mean dependence of the critical density for collapse δ_c on halo mass m; and the scatter around this mean value. All three of these are motivated by the physics of triaxial rather than spherical collapse. A comparison of the resulting mass function with N-body results shows that, if one uses δ_c(m) and its scatter as determined from simulations, then all three are necessary ingredients for obtaining ∼10 per cent accuracy. For example, assuming a constant value of δ_c with no scatter, as motivated by the physics of spherical collapse, leads to many more massive haloes than seen in simulations. The same model is also in excellent agreement with N-body results for the linear halo bias, especially at the high mass end where the traditional peak-background split argument applied to the mass function fit is known to underpredict the measured bias by ∼10 per cent. In the excursion set language, our model is about walks centred on special positions (peaks) in the initial conditions – we discuss what it implies for the usual calculation in which all walks contribute to the statistics.

large-scale structure of Universe

1 INTRODUCTION

The evolution of the abundance of virialized clusters is a powerful probe of the combined effects of the nature of the initial conditions (Gaussian or not?), the nature of gravity (general relativity or not?) and the expansion history of the Universe (cosmological constant or more complex?). This has motivated studies to understand and parametrize this dependence, so that data sets can more precisely constrain the physical models.

There are three widespread approaches to modelling cluster abundances. The most accurate of these is the brute force numerical simulation method; one simply asserts that the dark matter haloes which form in these simulations represent clusters in the Universe. However, this method is computationally expensive owing to the non-linear nature of gravitational structure formation. Furthermore, while simulations can provide insights into structure formation, they cannot replace a deep understanding of the physical mechanisms. The other two are simpler analytic models which consider the statistics of the initial fluctuation field to model halo formation. One is based on peaks, and attempts to model the comoving number density dn/dm of objects of a given mass m (Bardeen et al. 1986; Appel & Jones 1990; Manrique et al. 1998; Hanami 2001). The other, known as the excursion set approach, seeks to count the mass fraction f(m) that is bound-up in objects of a given mass, and, from this, to infer the abundances of objects themselves (Press & Schechter 1974; Peacock & Heavens 1990; Bond et al. 1991; Sheth & Tormen 2002; Maggiore & Riotto 2010a; Musso & Sheth 2012). The two are related by the simple fact that mass weighting each halo must yield the mass fraction in haloes, so

\begin{equation} \int _M^\infty {{\rm d}m} \frac{{\rm d}n(m)}{{\rm d}m} m = \bar{\rho }\int _M^\infty f(m) {{\rm d}m}. \end{equation}

(1)

However, this subtle difference has led to rather different formulations of the problem. Only recently have the two approaches been written in the same formalism (Paranjape & Sheth 2012), allowing a direct comparison of how they differ. Whereas all previous excursion set calculations of f were based on the statistics of all positions in the initial density fluctuation field; Paranjape & Sheth (2012) showed that the excursion set peaks (ESP) calculation of f explicitly uses the statistics of a small subset of positions in the initial field (those associated with peaks).

Working with a special subset of positions is motivated by measurements extracted from simulations showing that the excursion set prediction for the mass of the object in which a randomly chosen particle ends up is almost always incorrect (White 1996), whereas the prediction for the subset of particles which are at the centre of mass of the protohalo is almost always quite accurate (Sheth, Mo & Tormen 2001). Moreover, for this special subset of (centre-of-mass) particles, the predicted mass depends on whether one assumes a spherical or a triaxial collapse – and the prediction is more accurate if one uses the more physically appropriate (though by no means perfect!) triaxial collapse model. In other words, the physics and the statistics are simpler, and tell a consistent story only if one works directly with the special subset of positions around which collapse occurs. More recent work has shown that protohalo positions do indeed coincide with peaks in the initial field (Ludlow & Porciani 2011), and their initial velocities are biased with respect to that of the dark matter (Elia, Ludlow & Porciani 2012), as is expected for peaks (Desjacques & Sheth 2010). This provides additional motivation for the ESP approach.

However, this approach raises the following technical problem. The peaks approach uses the statistics of the smoothed density field, as well as its first and second derivatives. As a result, the rms values of these quantities play an important role. In the theory, these rms values are related to integrals over the initial power spectrum. However, they depend on the form of the smoothing filter, and, for the TopHat smoothing filter which is expected to be most directly related to the physics of gravitational collapse, some of these integrals do not converge. For Gaussian filters there is no such problem, and Paranjape & Sheth (2012) showed that, when expressed in suitably scaled (and natural) units ν, the ESP prediction for halo abundances is in excellent agreement with simulations, when the measured halo abundances are also expressed in scaled units. Namely, if one sets

\begin{equation} \frac{m}{\bar{\rho }}\frac{{\rm d}n(m)}{{\rm d}\ln m} = \nu f(\nu ) \left|\frac{{\rm d}\ln \nu }{{\rm d}\ln m}\right|, \end{equation}

(2)

then the ESP νf(ν) is in good agreement with the distribution of ln ν measured in simulations. However, this agreement is only a partial success because the relation between halo mass and ν depends on the smoothing filter, and the ESP model employed a Gaussian filter whereas simulations use a TopHat. Therefore, when expressed as a function of mass, the Gaussian smoothed ESP prediction is not a good description of dn/dln m in simulations (because the Jacobian dln ν/dln m depends on the filtering kernel). Moreover, the Gaussian ESP prediction assumed that haloes form through a spherical collapse, even though this is an incorrect description of the physics of halo formation (Sheth et al. 2001). In fact, there is substantial scatter even in the triaxial collapse model (Sheth et al. 2001; Dalal et al. 2008; Robertson et al. 2009; Sheth 2009; Ludlow & Porciani 2011), so the latter also provides an incomplete description of collapse. There has been some discussion of how to include stochasticity into the predictions (Bond & Myers 1996; Chiueh & Lee 2001; Sheth & Tormen 2002; Desjacques 2008; Maggiore & Riotto 2010b; Corasaniti & Achitouv 2011; Paranjape, Lam & Sheth 2012; Sheth, Chan & Scoccimarro 2012), but none of these self-consistently combines the model for the stochasticity with the fact that estimates of this stochasticity, from simulations or from theory, are for the special subset of positions around which collapse occurs.

Therefore, the main goals of this work are to include TopHat smoothing in the ESP formalism and incorporate the effects of more complicated collapse models which are seen in simulations. All this is the subject of Section 2. Some testable predictions which elucidate the relation between the ESP approach and the usual excursion set analysis of all (rather than special) positions in space are discussed in Section 3. Section 4 presents implications for the spatial distribution of these objects – the ESP predictions for Lagrangian halo bias – which can be used to further test the model. A final section summarizes our results. We will show numerical results at redshift z = 0 for two Λ cold dark matter (ΛCDM) cosmologies: one with |$(\Omega _{\rm m},\Omega _\Lambda ,\sigma _8,n_s,h) = (0.3,0.7,0.9,1.0,0.7)$| which we denote Wilkinson Microwave Anisotropy Probe 1 (WMAP1) and the other with (0.25, 0.75, 0.8, 1.0, 0.7) which we denote WMAP3.

2 THE UNCONDITIONAL MASS FUNCTION

2.1 Notation

We will designate Gaussian filtered quantities with the subscript G. By contrast, TopHat filtered objects will have no subscript. The Fourier transforms of the TopHat and Gaussian filters are W(kR) = (3/kR)j₁(kR) and |$W_{\rm G}(kR_{\rm G}) = {\rm e}^{-k^2R_{\rm G}^2/2}$|⁠, respectively, with j₁(y) being a spherical Bessel function.

We will frequently need the following integrals over the linearly extrapolated dimensionless matter power spectrum Δ²(k) ≡ k³P(k)/2π².

TopHat-filtered variance:

\begin{eqnarray} s \equiv \sigma _0^2 & \equiv \displaystyle\left\langle\delta ^2 \right\rangle = \int {\rm d}\ln k \Delta ^2(k)W(kR)^2. \end{eqnarray}

(3)

We will exclusively use the notation ν ≡ δ_c/σ₀, where δ_c = 1.686 is the usual spherical collapse threshold.

Gaussian-filtered spatial moments:

\begin{eqnarray} \sigma _{j{\rm G}}^2 & \equiv \displaystyle\int {\rm d}\ln k \Delta ^2(k) k^{2j}W_{\rm G}(kR_{\rm G})^2 \ \ ,\ j\ge 1, \end{eqnarray}

(4)

so that σ²_1G = 〈(∇δ_G)² 〉 and σ²_2G = 〈(∇²δ_G)²〉.

Mixed moment:

\begin{eqnarray} \sigma _{1{\rm m}}^2 & \equiv & \displaystyle\left\langle -\delta \nabla ^2\delta _{\rm G}\right\rangle \nonumber \\ &=& \displaystyle\int {\rm d}\ln k \Delta ^2(k) k^2 W_{\rm G}(kR_{\rm G})W(kR). \end{eqnarray}

(5)

Characteristic volume:

\begin{equation} R_{\ast } \equiv \sqrt{3}\sigma _{1{\rm G}}/\sigma _{2{\rm G}}; \qquad \qquad V_{\ast } = (2\pi )^{3/2} R_{\ast }^3. \end{equation}

(6)

These are the same as defined by Bardeen et al. (1986).

Spectral parameter:

\begin{eqnarray} \gamma \equiv \frac{\sigma _{1{\rm m}}^2}{\sigma _0\sigma _{2{\rm G}}} = \frac{\left\langle -\delta \nabla ^2\delta _{\rm G} \right\rangle }{\sqrt{\left\langle \delta ^2 \right\rangle \left\langle (\nabla ^2\delta _{\rm G})^2 \right\rangle }} = \left\langle x\mu \right\rangle , \end{eqnarray}

(7)

where we defined the standardized variables

\begin{equation} \mu \equiv \delta /\sigma _0\qquad {\rm and}\qquad x\equiv -\nabla ^2\delta _{\rm G}/\sigma _{2{\rm G}}. \end{equation}

(8)

Note that this definition of γ is similar but not identical to the corresponding one in Bardeen et al. (1986), since our peak heights are defined using TopHat smoothing.

2.2 Matching filter scales

The technical problem which we address in this subsection is how to ensure that the peaks in the Gaussian filtered density field δ_G will have an overdensity δ_TH = δ_c when smoothed with a TopHat, since essentially all measurements in simulations use TopHats only. In practice, we need a mapping between the scales R_TH and R_G of the two filters. All previous work accomplishes this either by matching the volumes: (4π/3) R³_TH = (2π)^3/2 R_G³, or by matching the variances: 〈δ²_TH〉 = 〈δ²_G〉. In this second case, the relation between R_TH and R_G depends on the shape of the power spectrum. But neither of these conditions guarantee that a peak identified on scale R_G satisfies δ_TH = δ_c.

For this reason, we construct the mapping between R_TH and R_G by finding that R_G (at a given R_TH) for which 〈 δ_G|δ_TH 〉 = δ_TH. This is equivalent to 〈δ_Gδ_TH〉 = 〈δ²_TH〉. This definition circumvents the technical complication arising from the variance of the Laplacian of δ_TH (which involves a divergent integral). For the ΛCDM P(k) we consider in this paper, R_G ≈ 0.46R_TH with a mild mass dependence which we account for, in the range we explore, 3 × 10¹⁰ < m/(h^{− 1} M_⊙) < 3 × 10¹⁵.

Our choice of matching between R_TH and R_G, namely 〈 δ_G|δ 〉 = δ, means that p(δ_G|δ) is a Gaussian with mean δ and variance 〈δ²_G〉 − s. However, 〈δ²_G〉 ≈ s to better than 5 per cent and, thus, p(δ_G|δ) ≈ δ_D(δ_G − δ) (since the cross-correlation between δ_G and δ is very close to unity). This significantly simplifies our ESP analysis for the following reason. In principle, the ESP framework also requires us to keep track of the variable |$\tilde{x}\equiv \delta ^{\prime }/\sqrt{\langle\delta ^{\prime 2}\rangle }$|⁠, where δ^′ = dδ/ds. This scalar (isotropic) variable correlates with μ and x, but not with the other (anisotropic) variables used by Bardeen et al. (1986). Therefore, one should in principle work with the three-dimensional Gaussian vector |$(\mu ,x,\tilde{x})$|⁠. However, because p(δ_G|δ) ≈ δ_D(δ_G − δ), it turns out that setting |$\tilde{x}=x$| yields an excellent approximation (note that this relation is exact for the Gaussian filter). We have indeed checked that our final answers for dn/dln m only change by a few per cent if we switch between the approximate and exact treatments. In the following, we will therefore always set |$\tilde{x}=x$| and only work with the two-dimensional Gaussian vector (μ, x).

Even though we smooth the density field with a TopHat filter, there is no compelling reason to use such a filter except to make a connection with the spherical collapse approximation. Numerical simulations indeed suggest that, while the Lagrangian volume occupied by the protohaloes is rather compact, it is more diffuse than a TopHat window (Porciani, Dekel & Hoffman 2002; Dalal et al. 2008; Despali, Tormen & Sheth 2013). However, since a deep understanding of halo collapse is still lacking, we will stick to the TopHat (and Gaussian) filters for simplicity.

2.3 Excursion set peaks with a constant barrier

Apart from the fact that the peak height is defined using TopHat filtering, there is no formal change in the derivation by Paranjape & Sheth (2012) of the number density of ESP for a constant barrier (and this derivation is formally the same as that in Appel & Jones 1990). In this case, we get

\begin{eqnarray} f_{\rm ESP}(\nu ) &=& \displaystyle V \mathcal {N}_{\rm ESP}(\nu )\nonumber \\ &=&(V/V_{\ast })({\rm e}^{-\nu ^2/2}/\sqrt{2\pi })\nonumber \\ && \displaystyle\times \int _0^\infty {\rm d}x \frac{x}{\gamma \nu } F(x)p_{\rm G}(x-\gamma \nu ;1-\gamma ^2), \end{eqnarray}

(9)

where |$V=m/\bar{\rho }=4\pi R^3/3$| is the Lagrangian volume associated with the TopHat smoothing filter, F(x) is the peak curvature function from equation A15 of Bardeen et al. (1986) and |$p_{\rm G}(y-\bar{y};\Sigma ^2)$| is a Gaussian distribution in y with mean |$\bar{y}$| and variance Σ². (As noted in Paranjape & Sheth 2012, F(x) ≠ 1 reflects the fact that the ESP calculation averages over a special subset of positions.) Since we have been careful to define ν and hence the mass using TopHat filtering, the mass function follows in the usual way from equation (2).

Fig. 1 compares the ESP result for a constant barrier B(s) = δ_c = 1.686 (solid red) with the fit to N-body simulations from Tinker et al. (2008) (dashed green). For the latter we used parameters from table 2 of Tinker et al. (2008) appropriate for haloes identified using the spherical overdensity (SO) definition at 200 times the mean density of the universe. We used cosmological parameters for the WMAP1 cosmology mentioned in Section 1, for which the Tinker et al. (2008) function gives a good fit (cf. their table 2). Later, we will also show results for the WMAP3 cosmology which was also explored by Tinker et al. (2008) – changing cosmology does not affect our results significantly. For comparison, we also show the functional form of Sheth & Tormen (1999) (short dashed blue) with their (q, p) = (0.7, 0.26) which gives a good fit to the MICE simulations (Crocce et al. 2010).

Figure 1.

ESP mass function for the mixed filtering approach discussed in the text, using a constant barrier (solid red). This is compared with the fit to N-body simulations given by Tinker et al. (2008), appropriate for SO haloes that enclose 200 times the mean density of the universe (dashed green). For comparison, we also show the Sheth & Tormen (1999) mass function (short dashed blue) with their (q, p) = (0.7, 0.26), which gives a good fit to the MICE simulations (Crocce et al. 2010). The ESP+constant barrier mass function overpredicts the abundance at all scales of interest. All curves used the WMAP1 cosmology – we have checked that changing the cosmology does not significantly alter the comparison.

Open in new tab Download slide

Clearly, equation (9) predicts far too many haloes at all interesting scales, including the high mass end (ν² = 10 corresponds to m = 1.3 × 10¹⁵h^{− 1} M_⊙ for this cosmology).

2.4 Beyond the spherical-collapse constant barrier

It is tempting to attribute the mismatch shown in Fig. 1 to the effects of non-spherical collapse. While these effects are certainly important at smaller mass scales, it is not obvious that they also matter at the largest scales, where the spherical collapse approximation should be quite accurate. To see that they are in fact relevant, note that the triaxial collapse model predicts that the critical density for collapse scales approximately as

\begin{equation} \delta _{\rm ec} \approx \delta _{\rm c}\left[1 + 0.4 (\sigma _0/\delta _{\rm c})^{1.2}\right] \end{equation}

(10)

(Sheth et al. 2001). Simulations indeed show that this expression for δ_ec provides a reasonable approximation to the mass dependence of the critical density required for collapse (Sheth et al. 2001). However, this implies that when ν = δ_c/σ₀ is as large as ∼ 4, then departures from the spherical δ_c are already of the order of 10 per cent. Therefore, we must include deviations from the spherical collapse barrier in order to have a fair comparison with the Tinker et al. (2008) simulations.

Before explaining how this can be done, we note that, although equation (10) is a reasonable approximation to the mass dependence of the critical density required for collapse (Sheth et al. 2001), there is substantial scatter around this value (Dalal et al. 2008; Sheth 2009). This scatter has been quantified by Robertson et al. (2009)¹ who found that, over the range log₁₀(ν²) ∈ (−0.3, 0.9), the value of B was Lognormally distributed around a mean value of ≃ 1.51 + 0.49σ₀ with variance (0.09σ²₀) (cf. their table 1 and fig. 3). The physical origin of this scatter is unclear. Triaxial collapse predicts some scatter (Sheth et al. 2001; Despali et al. 2013), which is related to the statistics of the tidal field, but the measured scatter is considerably larger. In any case, the physical origin of this scatter is unimportant for what follows – though we believe that this is an issue which deserves further investigation.

2.5 Excursion set peaks with a square-root barrier and scatter

To include the mean trend and scatter in barrier heights shown by Sheth et al. (2001) and quantified by Robertson et al. (2009), we will adopt the following model, first discussed in Paranjape et al. (2012). We consider a family of deterministic ‘square-root’ barriers B given by

\begin{equation} B = \delta _{\rm c}+ \beta \sqrt{s} , \end{equation}

(11)

with δ_c = 1.686 and β being a stochastic variable. Therefore, if f_ESP(ν|β) denotes the ESP solution for a given value of β, then our full solution will be

\begin{equation} f_{\rm ESP}(\nu ) = \int {\rm d}\beta p(\beta ) f_{\rm ESP}(\nu |\beta ). \end{equation}

(12)

Below, we use simulations to motivate our choice of p(β), and provide an explicit expression for f_ESP(ν|β). Still, it is easy to understand qualitatively what happens.

For ESP, the square-root barrier (11) with β > 0 dramatically lowers the halo counts at large masses (fig. 3 in Paranjape & Sheth 2012). While this goes in the right direction, introducing scatter in the value of β will counteract this decrease. This is because a scatter in β induces a scatter in the predicted mass of the halo. The steepness of the mass function implies that the scatter preferentially moves small mass objects to larger masses (the same effect that leads to Eddington bias), which in turn increases the abundance at large masses.

The importance of the scatter depends on the distribution of β, our choice for which is motivated by the Lognormal distribution of B seen in simulations for haloes of a fixed mass. Note that what we need as input is the ‘prior’ distribution p(β), knowing which our analysis then predicts f_ESP(ν) using equation (12). However, the distribution reported by Robertson et al. (2009) corresponds to the conditional distribution p(β|ν) which we can only construct a posteriori using the relation

\begin{equation} p(\beta |\nu ) = f_{\rm ESP}(\nu |\beta )p(\beta )/f_{\rm ESP}(\nu ), \end{equation}

(13)

which is a simple application of Bayes’ theorem.

In practice, we choose p(β) to be Lognormal with mean and variance chosen such that the mass function f_ESP(ν) as well as the conditional distribution p(β|ν) are in reasonable agreement with the results of simulations. Having done this, we will see later that the clustering of haloes is also correctly predicted by our model, which is a non-trivial test of the formalism. We set 〈 β 〉 = 0.5 and Var(β) = 0.25 for the Lognormal p(β); we have checked that varying these numbers within a ∼10 per cent range does not affect our results significantly. We avoided fitting these numbers more accurately because, as mentioned earlier, their physical origin is at present unclear and a more accurate fit would not lead to any physical insights. Moreover, recent work, in which haloes were identified using an ellipsoidal rather than spherical halo finder, suggests that the mean value is closer to 1.68 (1 + 0.2σ₀) with variance around this mean of 0.04σ²₀ (Despali et al. 2013).

For square-root barriers with fixed β (equation 19 of Paranjape & Sheth 2012, for the barrier in equation 11 above),

\begin{eqnarray} f_{\rm ESP}(\nu |\beta ) &=& (V/V_\ast ) ({\rm e}^{-(\nu +\beta )^2/2}/\sqrt{2\pi })\nonumber \\ && \displaystyle\times \int _{\beta \gamma }^\infty {\rm d}x \frac{x-\beta \gamma }{\gamma \nu } F(x)\nonumber \\ && \times p_{\rm G}(x-\beta \gamma - \gamma \nu ;1-\gamma ^2). \end{eqnarray}

(14)

We have assumed β > 0, which is true for the Lognormal model we will consider. If β is allowed to take negative values as well, then the lower limit of the integral over x must be replaced by x_min = max{0, βγ}, so as to enforce the peak constraint. While there is no closed form expression for f_ESP(ν), the integrals involved can be easily computed numerically.

Fig. 2 shows the results for the WMAP1 (solid red) and WMAP3 (dotted red) cosmologies. Note that the value of σ₈ is ∼10 per cent smaller in the WMAP3 cosmology. We plot the ratio of the ESP prediction in equation (12) and the Tinker et al. (2008) fit. For comparison we also indicate the ratio of the MICE fit with the Tinker et al. fit (for the WMAP1 cosmology) as the short dashed blue curve. The horizontal dotted lines mark 5 per cent deviations from the Tinker et al. fit; the typical scatter in their measurements in this mass range is 5–10 per cent (cf. fig. 6 a of Tinker et al. 2008).

Figure 2.

Ratio of the ESP mass function for a square-root barrier with scatter in the barrier slope (equation 12) to the Tinker et al. (2008) fit, for WMAP1 (solid red) and WMAP3 (dotted red) cosmologies. In each case, we set δ_c = 1.686 and p(β) to be Lognormal with 〈 β 〉 = 0.5 and Var(β) = 0.25, for the reasons discussed in the text. For comparison, we also show the ratio of the MICE fit with the Tinker et al. fit for the WMAP1 cosmology (short dashed blue). Horizontal dotted lines mark 5 per cent deviations from the Tinker et al. fit; the typical scatter in their measurements in this mass range is 5–10 per cent.

Open in new tab Download slide

Fig. 3 shows the predicted mean 〈 β|ν 〉 and scatter Σ_β(ν) of the conditional distribution p(β|ν) (equation 13) for the WMAP1 (red) and WMAP3 (blue) cosmologies. The black curves show the corresponding quantities measured by Robertson et al. (2009), assuming the fit B_fit = 1.506 + 0.487σ from their table 1 and using 〈 β|ν 〉_Robertson + = (B_fit − 1.686)/σ in order to compare with our square-root barrier model in equation (11). The arrows indicate the approximate range over which Robertson et al. performed these fits for 〈 β|ν 〉 (dashed) and Σ_β(ν) (solid). The ESP predictions are within 30 per cent of the measured quantities. Together with Fig. 2, this shows that our model works rather well.

$The conditional mean 〈 β|ν 〉 and scatter $\Sigma _\beta (\nu )\equiv \sqrt{{\rm Var}(\beta |\nu )}$ of the distribution p(β|ν) (equation 13) for the WMAP1 (red) and WMAP3 (blue) cosmologies. The black curves show the corresponding quantities measured by Robertson et al. (2009). The ESP curves used δc = 1.686, fESP(ν|β) from equation (14) and assumed p(β) to be Lognormal with 〈 β 〉 = 0.5 and Var(β) = 0.25. The curve for 〈 β|ν 〉 from Robertson et al. (2009) assumed the ‘linear’ fit from their table 1, Bfit = 1.506 + 0.487σ0, and shows 〈β|ν〉 = (Bfit − 1.686)/σ0. The arrows indicate the approximate range over which Robertson et al. performed these fits for 〈 β|ν 〉 (dashed) and Σβ(ν) (solid). The ESP predictions are within 30 per cent of the measured quantities.$

Figure 3.

The conditional mean 〈 β|ν 〉 and scatter |$\Sigma _\beta (\nu )\equiv \sqrt{{\rm Var}(\beta |\nu )}$| of the distribution p(β|ν) (equation 13) for the WMAP1 (red) and WMAP3 (blue) cosmologies. The black curves show the corresponding quantities measured by Robertson et al. (2009). The ESP curves used δ_c = 1.686, f_ESP(ν|β) from equation (14) and assumed p(β) to be Lognormal with 〈 β 〉 = 0.5 and Var(β) = 0.25. The curve for 〈 β|ν 〉 from Robertson et al. (2009) assumed the ‘linear’ fit from their table 1, B_fit = 1.506 + 0.487σ₀, and shows 〈β|ν〉 = (B_fit − 1.686)/σ₀. The arrows indicate the approximate range over which Robertson et al. performed these fits for 〈 β|ν 〉 (dashed) and Σ_β(ν) (solid). The ESP predictions are within 30 per cent of the measured quantities.

Open in new tab Download slide

3 IMPLICATIONS FOR AVERAGES OVER ALL, RATHER THAN A SPECIAL SUBSET OF POSITIONS

The ESP model of the previous section is the first to self-consistently integrate the complex physics of non-spherical collapse with the fact that the physics is almost certainly simplest around the centre of mass of the protohalo. Although the agreement with the halo counts represents a significant and non-trivial success, and the next section shows that the predicted spatial distribution is also accurate, it makes a number of other testable predictions which are related to the logic of the excursion set approach.

For example, one can construct the centre-of-mass random walk associated with each halo and explicitly check, object by object, whether or not the corresponding barrier was crossed at a larger scale. This is, in fact, precisely the kind of test which was performed by Sheth et al. (2001). Revisiting this test is left to future work.

3.1 Distribution of δ for all particles in protohaloes of a given mass

We remarked in Section 1 that our ESP model explicitly averages only over special positions in the universe (see the text following equation 9 for exactly where this happens). If we treat these positions as the fundamental quantities – for which the physics of collapse is most easily applied – then by accounting for their initial density profiles, it is straightforward to estimate what would have happened had we averaged over all positions. Because our model has been calibrated for centre-of-mass positions solely, predicting the right distribution of δ values within each protohalo provides another test of our approach. We illustrate the idea below.

We begin by noting that because δ_pk ≡ δ_c + βσ₀, the distribution of β at fixed ν ≡ δ_c/σ₀ implies a distribution of δ_pk:

\begin{equation} p(\delta _{\rm pk}|S) {\rm d}\delta _{\rm pk} = p(\beta |\nu ) {\rm d}\beta . \end{equation}

(15)

Recall that fixed ν means fixed halo mass m, smoothing scale R or variance S and simulations suggest p(β|ν) is approximately Lognormal.

Now, suppose that there is an excursion set peak of height δ_pk identified on smoothing scale R_pk. Let p(δ|r, δ_pk, S_pk) denote the probability that the overdensity on smoothing scale R_pk at distance r from this peak is δ. This distribution is Gaussian, with mean and variance given by equations 7.8– 7.10 in Bardeen et al. (1986). Therefore, if we choose a random position within the set of all initial peak patches of scale R_pk, then it will have overdensity δ with probability

\begin{equation} p(\delta |S_{\rm pk}) = \int {\rm d}\delta _{\rm pk} p(\delta _{\rm pk}|S_{\rm pk}) \int _0^R 3\frac{{\rm d}r}{r} \frac{r^3}{R^3} p(\delta |r,\delta _{\rm pk},S_{\rm pk}). \end{equation}

(16)

This expression is particularly easy to evaluate using Monte Carlo methods: generate (r/R)³ and δ_pk from uniform and Lognormal distributions, respectively, and then generate a Gaussian variate whose mean and variance depend on r/R and |$\delta _{\rm pk}/\sqrt{S_{\rm pk}}$|⁠, following Bardeen et al. (1986). The histogram in Fig. 4 shows the result when S_pk = 0.9 (this corresponds to halo masses of 7.8 × 10¹³h^{− 1}M_⊙). The smooth solid curve which provides an excellent fit shows a Gaussian with the same mean and variance as δ. We have included it because the mean and variance of p(δ|r, δ_pk, S_pk) are known analytically (from Bardeen et al. 1986), so the mean and variance of p(δ|S_pk) are simply related to volume averages of correlation functions and the mean and variance of p(δ_pk|S_pk).

Figure 4.

Predicted distribution of δ for all particles in protohaloes of a given mass. Histogram shows our model, equation (16), which uses as input the smooth Lognormal distribution shown for p(δ_pk|S_pk). Our model is well approximated by a Gaussian (solid curve); it has a slightly different mean and variance than do the measurements in simulations (the dashed Gaussian).

Open in new tab Download slide

The dashed curve shows a Gaussian distribution with mean (1.67 + 0.075 S_pk)/1.35 and variance (0.35/1.35) S_pk, which describes the measurements in simulations of this quantity rather well (Achitouv et al. 2013). In effect, our model is trying to predict this dashed curve. Because we use as input the Lognormal distribution of peak heights p(δ_pk|ν) that is seen in simulations (shown as the other solid curve which peaks at larger δ: it has mean |$\langle\delta _{\rm pk}|S_{\rm pk}\rangle = 1.506 + 0.487\sqrt{S_{\rm pk}}$|⁠, and variance 0.09 S_pk), our model has no free parameters, so it is perhaps remarkable that the model does as well as it does. Although our prediction is narrower than the measured distribution (predicted variance is 0.18 S_pk whereas the measured value is 0.26 S_pk) its mean, 1.35, is within 5 per cent of the measured value of 1.29.

Fig. 5 shows a similar analysis for S_pk = 2, which corresponds to lower mass haloes (7.8 × 10¹²h^{− 1}M_⊙) for which the assumption that haloes form from peaks is less secure; clearly, the discrepancies between the model and the measurements are qualitatively similar even at these small masses. The predicted mean and variance are 1.43 and 0.136 S_pk, whereas the measured values are 1.36 and 0.26 S_pk.

Figure 5.

Same as previous figure, but now for S_pk = 2.

Open in new tab Download slide

3.2 Distribution of predicted mass for all particles in protohaloes of a given mass

The lower mean value implies that the mass which the excursion set approach predicts for the particles which are not at the centre of mass should be systematically smaller than the mass of the halo in which they end up – in qualitative agreement with one of the findings of Sheth et al. (2001, see their figs 2–4).

To see that this is indeed the case, we can use our ESP approach (walks centred on peaks in the initial field) to model what the usual excursion set estimate (walks centred on all positions in the initial field) really represents. Recall that, for the statistics of all positions, we know that if δ_c is independent of m then the excursion set prediction for the mass fraction in objects with mass m(s) is |$sf(s) = \nu \exp (-\nu ^2/2)/\sqrt{2\pi }$| with ν² = δ²_c/s. [Strictly speaking, to compare most directly with Sheth et al. (2001) we should account for the distribution in δ_c values at each s; but the constant deterministic δ_c discussion below captures the gist of the argument. Also, our expressions below are appropriate for a sharp k filter, but the logic, and so the qualitative trends, do not depend on this choice.]

In our ESP model, this fraction is given by taking each peak of mass M, and then accounting for the distribution of excursion set predictions for s(m) > S(M) which come from all the other mass elements which make up M (rather than just the one mass element at the centre of mass), and integrating the distribution of S over the range 0 < S < s, i.e.

\begin{equation} sf(s) = \int _0^s \frac{{\rm d}S}{S} Sf_{\rm ESP}(S) sf(s|S), \end{equation}

(17)

where

\begin{equation} Sf_{\rm ESP}(S) = \frac{M}{\bar{\rho }} \frac{{\rm d}n_{\rm ESP}(M)}{{\rm d}\ln M} \left|\frac{{\rm d}\ln M}{{\rm d}\ln S}\right| \end{equation}

(18)

(from equation 2) and f(s|S) is the excursion set prediction for the mass fraction of M which is incorrectly predicted to have mass m = M(s) instead of M(S). Because we know sf(s) as well as n_ESP(S), we can solve (by numerical back substitution) for f(s|S). However, we can get a simple, intuitive estimate of the result as follows.

For constant δ_c, figs 1 and 2 in Paranjape & Sheth (2012) show that Sf_ESP(S) is rather well approximated simply by |$\nu f(\nu ) = \nu \exp (-\nu ^2/2)/\sqrt{2\pi }$| with ν² = (Ξδ_c)²/S for some Ξ < 1. [Strictly speaking, they used the Sheth–Tormen functional form for νf(ν); however, at high masses, this form reduces to the simpler one given here with |$\Xi =\sqrt{0.7}$|⁠.] Since this is the same functional form for S that appears on the left-hand side of equation (17) for s, it is easy to check that sf(s|S) should also be the same function of ν, but with ν² = (δ_c − Ξδ_c)²/(s − S).

Thus, this calculation shows explicitly that the average over all walks will yield f which has a similar functional form to the one for the specially chosen subset, but the two will have different values of δ_c. Moreover, the distribution of predicted masses will be very different from a delta function centred on M. Rather, it will look just like the conditional excursion set prediction for the fraction of all walks which first cross δ_c on scale s given that they are known to pass through Ξδ_c < δ_c on scale S before first crossing δ_c. This latter point quantifies the explanation given by Sheth et al. (2001) when addressing the concerns raised by White (1996). It also expands on a point first made by Paranjape & Sheth (2012): the ESP framework naturally accounts for the fact that one must rescale |$\nu \rightarrow \sqrt{q}\nu$| with q ∼ 0.7 if one wants to use the usual excursion set expressions (based on the statistics of all positions) to model halo counts – but this rescaling is not necessary if one uses the statistics of the special positions for which the physics is simplest. Of course, because the physics of collapse refers to collapse around the halo centres of masses, and not around random positions in space, this discussion highlights the dangers of interpreting the barrier height in the random walk calculation with that for the physics of collapse around special positions in space.

4 HALO BIAS

In this section, we derive expressions for the predicted clustering of haloes in our model.

In the excursion set approach, the clustering of haloes relative to the initial dark matter field – or ‘Lagrangian halo bias’ – is captured by the conditional distribution f(s|δ₀, S₀) of walks that first cross the barrier at some small scale s having previously passed through the value δ₀ at some large scale S₀. As discussed by Musso, Paranjape & Sheth (2012) and Paranjape & Sheth (2012), a convenient way of organizing this bias, for both excursion sets and ESP, is through cross-correlations of the halo density 〈 ρ_h|δ₀ 〉 = f(ν|δ₀, S₀)/f(ν) with Hermite polynomials in the matter overdensity δ₀:

\begin{eqnarray} b_n & \equiv& S_{0}^{-n/2}\left\langle\rho _hH_n(\delta _{0}/\sqrt{S_{0}}) \right\rangle \nonumber \\ &=& \displaystyle S_{0}^{-n/2}\int _{-\infty }^{\infty }{\rm d}\delta _{0}p_{\rm G}(\delta _{0};S_{0}) \left\langle\rho _{\rm h}|\delta _{0} \right\rangle H_n(\delta _{0}/\sqrt{S_{0}}), \end{eqnarray}

(19)

where |$H_n(x)={\rm e}^{x^2/2}(-{\rm d}/{\rm d}x)^n {\rm e}^{-x^2/2}$| are the ‘probabilist’s’ Hermite polynomials, and the bias parameters b_n can be shown to have the structure

\begin{equation} b_n =\left(\frac{S_{\times }}{S_{0}}\right)^n \sum _{r=0}^n {n\choose r} b_{nr}\epsilon _{\times }^r, \end{equation}

(20)

where

\begin{eqnarray} S_{\times }&=& \displaystyle\int {\rm d}\ln k \Delta ^2(k)W(kR)W(kR_0),\nonumber \\ \epsilon _{\times }&=& 2 {\rm d}\ln S_{\times }/{\rm d}\ln s \end{eqnarray}

(21)

are cross-correlations between the small and the large scale.

The calculation of bias therefore requires the conditional mass function of ESP. Introducing scatter in the barrier parameter β is straightforward, as we discuss below. The ESP conditional multiplicity function for a square-root barrier at fixed β and a fixed density δ₀ at scale S₀ is given by

\begin{eqnarray} f_{\rm ESP}(\nu |\beta ,\delta _{0}) &=& (V/V_\ast ) p(B/\sqrt{s}|\delta _{0}) \nonumber \\ && \times \displaystyle\frac{1}{\gamma \nu } \int _{\beta \gamma }^\infty {\rm d}x (x-\beta \gamma )F(x)\nonumber \\ && \times p(x|B/\sqrt{s},\delta _{0}), \end{eqnarray}

(22)

where p(μ|δ₀) and p(x|μ, δ₀) are one-dimensional conditional Gaussian distributions that can be written using the covariance matrix of the trivariate Gaussian in (μ, x, δ₀).

The bias coefficients δⁿ_cb_n can be shown to be algebraically equivalent to the Taylor coefficients of the expansion of equation (22) in powers of |$\bar{\delta }_0\equiv \delta _{0}/\delta _{\rm c}$|⁠, in a specific limit (Musso et al. 2012; Paranjape & Sheth 2012). This leads to a straightforward prescription for calculating these quantities. In particular, one must evaluate equation (22) after formally setting |$\tilde{\sf{c}}=0$|⁠, where |$\tilde{\sf{c}}$| is the matrix that one must subtract from the unconditional covariance matrix of the variables (μ, x) to obtain the covariance matrix of the conditional Gaussian p(μ, x|δ₀). The resulting conditional multiplicity is

\begin{eqnarray} && \!\!\!\!f_{\rm ESP}(\nu |\beta ,\delta _{0},\tilde{\sf{c}}=0) \nonumber\\ &&= (V/V_\ast ) ({\rm e}^{-(\nu +\beta -\bar{\delta }_0(S_{\times }/S_{0})\nu )^2/2}/\sqrt{2\pi })\nonumber \\ && \displaystyle\times \frac{1}{\gamma \nu } \int _{\beta \gamma }^\infty {\rm d}x (x-\beta \gamma )F(x)\nonumber \\ && \times p_{\rm G}(x-\beta \gamma - \gamma \nu _1;1-\gamma ^2), \end{eqnarray}

(23)

where we have defined |$\nu _1\equiv \nu (1- \bar{\delta }_0(S_{\times }/S_{0})(1-\epsilon _{\times }))$|⁠, and S_× and ϵ_× were defined in equation (21).

Taylor expanding equation (23) in powers of |$\bar{\delta }_0$|⁠, the bias parameters at fixed β can be written as

\begin{eqnarray} \frac{\delta _{\rm c}b_1(\nu ,\beta )}{(S_{\times }/S_{0})} &= \mu _1(\nu ,\beta ) + (1-\epsilon _{\times })\lambda _1(\nu ,\beta ), \end{eqnarray}

(24)

\begin{eqnarray} \frac{\delta _{\rm c}^2 b_2(\nu ,\beta )}{(S_{\times }/S_{0})^2} &= \mu _2(\nu ,\beta ) + 2(1-\epsilon _{\times })\mu _1(\nu ,\beta )\lambda _1(\nu ,\beta )\nonumber \\ & \phantom{\mu _2} + (1-\epsilon _{\times })^2\lambda _2(\nu ,\beta ), \end{eqnarray}

(25)

where, using |$\Gamma \equiv \gamma /\sqrt{1-\gamma ^2}$|⁠,

\begin{eqnarray} \mu _n(\nu ,\beta ) &= \nu ^nH_n(\nu +\beta ), \end{eqnarray}

(26)

\begin{eqnarray} \lambda _n(\nu ,\beta ) &= (-\Gamma \nu )^n\left\langle H_n\left.\left(y-\beta \Gamma -\Gamma \nu \right)\right|\nu ,\beta \right\rangle _y, \end{eqnarray}

(27)

and for some function h(y, ν, β), we have

\begin{eqnarray} && \!\!\!\!\left\langle h(y,\nu ,\beta )|\nu ,\beta \right\rangle _y \nonumber \\ && \!= \displaystyle\frac{\int _{\beta \Gamma }^\infty {\rm d}y (y-\beta \Gamma )F(y\gamma /\Gamma )p_{\rm G}(y-\beta \Gamma -\Gamma \nu ;1)h(y,\nu ,\beta )}{\int _{\beta \Gamma }^\infty {\rm d}y (y-\beta \Gamma )F(y\gamma /\Gamma )p_{\rm G}(y-\beta \Gamma -\Gamma \nu ;1)}. \end{eqnarray}

(28)

Marginalizing over β gives the bias parameters as

\begin{eqnarray} \delta _{\rm c}b_1(\nu ) &=& \left(\frac{S_{\times }}{S_{0}}\right)\bigg [\left\langle \mu _1|\nu \right\rangle + (1-\epsilon _{\times })\left\langle \lambda _1|\nu \right\rangle \bigg ], \end{eqnarray}

(29)

\begin{eqnarray} \delta _{\rm c}^2 b_2(\nu ) &=& \displaystyle\left(\frac{S_{\times }}{S_{0}}\right)^2\bigg [\left\langle \mu _2|\nu \right\rangle + 2(1-\epsilon _{\times })\left\langle \mu _1\lambda _1|\nu \right\rangle \nonumber \\ &&+ (1-\epsilon _{\times })^2\left\langle\lambda _2|\nu \right\rangle \bigg ], \end{eqnarray}

(30)

where for some function g(ν, β), we have

\begin{eqnarray} && \!\!\!\left\langle g|\nu \right\rangle = \int {\rm d}\beta p(\beta |\nu )g(\nu ,\beta ), \end{eqnarray}

(31)

where p(β|ν) was given in equation (13).

The presence of ϵ_× in these expressions for real-space bias is an indication of k-dependent Fourier-space bias, as discussed at length by Musso et al. (2012).² At large scales R₀ → ∞, ϵ_× → 0. Consequently, the coefficients b_n0 in equation (20) correspond to the large-scale (k → 0) Fourier-space bias that is routinely measured in simulations. Further, Musso et al. (2012) pointed out that b_n0 correspond to the usual peak-background split expressions for bias, b_n0 = f(ν)^{− 1}( − ∂/∂δ_c)ⁿf(ν). In other words, the peak-background split bias parameters correspond to taking the large-scale limit R₀ → ∞, or equivalently, setting ϵ_× → 0 and S_×/S₀ → 1 in the coefficients b_n.

Fig. 6 shows our ESP prediction for the linear peak-background split bias δ_cb₁₀ (i.e. setting ϵ_× = 0, S_×/S₀ = 1 in equation 29) as the solid red curve. The dotted black curve is the spherical collapse prediction ν² − 1 (Mo & White 1996). The dashed green curve shows the fit presented by Tinker et al. (2010) for the same simulations analysed by Tinker et al. (2008). We have checked that in this case the results for the WMAP1 and WMAP3 cosmologies are nearly identical, so we only show the results for WMAP3. The ESP curve matches the N-body fit very well at large masses, which is where the ESP framework is expected to work well. Note that the curve for Tinker et al. (2010) is not simply the logarithmic derivative of the mass function calibrated by Tinker et al. (2008); that derivative would underestimate the large-scale bias by up to 10 per cent (Manera, Scoccimarro & Sheth 2010; Tinker et al. 2010). (For comparison, the short dashed blue curve shows the peak-background split prediction for the MICE mass function.)

Scale-independent (peak-background split) linear Lagrangian halo bias δcb10 in the ESP framework for a square-root barrier with Lognormal scatter (solid red). This curve is in excellent agreement with the fit to N-body simulations from Tinker et al. (2010), the simulations being the same as those for which Fig. 2 showed a mass function comparison (see the text for a discussion). The dotted black curve is the spherical collapse prediction ν2 − 1 (Mo & White 1996), and the short dashed blue curve is the peak-background split prediction associated with the mass function fit to the MICE simulations.

Figure 6.

Scale-independent (peak-background split) linear Lagrangian halo bias δ_cb₁₀ in the ESP framework for a square-root barrier with Lognormal scatter (solid red). This curve is in excellent agreement with the fit to N-body simulations from Tinker et al. (2010), the simulations being the same as those for which Fig. 2 showed a mass function comparison (see the text for a discussion). The dotted black curve is the spherical collapse prediction ν² − 1 (Mo & White 1996), and the short dashed blue curve is the peak-background split prediction associated with the mass function fit to the MICE simulations.

Open in new tab Download slide

Fig. 7 shows the corresponding result for the quadratic peak-background split coefficient δ²_cb₂₀ for ESP (solid red) and the spherical collapse result ν²(ν² − 3) (Mo, Jing & White 1997, dotted black). The dashed blue curve shows the Mo et al. (1997) expression after rescaling |$\nu \rightarrow \sqrt{0.7}\nu$|⁠. In this case, there are no measurements in simulations with which to compare – this is the subject of ongoing work – but we have included the figure to show that b₂₀ behaves qualitatively like b₁₀, transitioning from the usual (sharp k excursion set) expression rescaled by |$\nu \rightarrow \sqrt{0.7}\nu$| at ν ∼ 1 to lying above this rescaled value at larger ν.

Figure 7.

Same as previous figure, but now for the second-order bias coefficient. In this case, there are no measurements from simulations with which to compare.

Open in new tab Download slide

These plots illustrate another point made by Paranjape & Sheth (2012): the ESP framework naturally accounts for the fact that the measured large-scale bias deviates from that predicted by a naive application of the peak-background split to a mass function that was fit using the rescaled variable |$\nu \rightarrow \sqrt{q}\nu$| with q ∼ 0.7.

4.1 Reconstructing the peak-background split parameters

In this subsection, we revisit the idea of reconstructing the peak-background split parameters b_n0 from suitable real-space measurements of the parameters b_n (Musso et al. 2012). In particular, we wish to address some subtleties that arise due to the stochasticity of the barrier height.

The theoretical quantities plotted in Figs 6 and 7 required us to set ϵ_× → 0 and S_×/S₀ → 1 by hand. In Fourier space, this would correspond in practice to taking the limit k ∼ 1/R₀ → 0, which is how Tinker et al. (2010), e.g., derived their fit for the linear bias. Musso et al. (2012), on the other hand, pointed out that the theoretical predictions for δⁿ_cb_n at finiteR₀ are equivalent to an average over the Hermite-transformed initial density field |$H_n(\delta _{0}/\sqrt{S_{0}})$| at the centres of haloes (equation 19). To be precise, if there are N haloes in a given mass bin, one estimates the bias parameters at this mass as

\begin{equation} b_n^{\rm (msd)} = S_{0}^{-n/2}\frac{1}{N}\sum _{i=1}^{N}H_n(\delta _{0i}/\sqrt{S_{0}}), \end{equation}

(32)

where δ_0i is the Lagrangian density contrast smoothed on scale R₀ and centred on the ith halo in the bin, and S₀ is the linearly extrapolated variance at R₀.

If we use a deterministic barrier in our theoretical model for the bias, then these measurements can be explicitly converted into estimates of the peak-background split parameters b_n0 (Musso et al. 2012; Paranjape & Sheth 2012). To see this, note that for the square-root barrier with fixed β, we can use equations (24) and (25) to write the peak-background split parameters in terms of b₁(ν, β) and b₂(ν, β) as

\begin{eqnarray} \delta _{\rm c}b_{10}(\nu ,\beta ) &=& \mu _1 + \lambda _1\nonumber \\ &=& \displaystyle\frac{1}{(1-\epsilon _{\times })}\bigg [\frac{\delta _{\rm c}b_1(\nu ,\beta )}{(S_{\times }/S_{0})} - \epsilon _{\times }\mu _1 \bigg ], \end{eqnarray}

(33)

\begin{eqnarray} \delta _{\rm c}^2 b_{20}(\nu ,\beta ) &=& \mu _2 + 2\mu _1\lambda _1 + \lambda _2\nonumber \\ &=& \displaystyle\frac{1}{(1-\epsilon _{\times })^2}\bigg [\frac{\delta _{\rm c}^2 b_2(\nu ,\beta )}{(S_{\times }/S_{0})^2} - 2\epsilon _{\times }\mu _1\left(\frac{\delta _{\rm c}b_1(\nu ,\beta )}{(S_{\times }/S_{0})} - \mu _1\right)\nonumber \\ &&- \epsilon _{\times }(2-\epsilon _{\times })\mu _2 \bigg ], \end{eqnarray}

(34)

which reduce to the expressions in equations (44) and (47) of Musso et al. (2012) when β = 0.

If β is stochastic (as we have seen it must be), then the measurements on the right-hand side of equation (32) correspond to the marginalized bias coefficients – b^(msd)_n↔〈b_n(ν, β)|ν 〉 – since we did not keep track of the value of β for each halo. For the linear coefficient b₁₀ this does not pose any problem, and we can estimate this coefficient using 〈 μ₁|ν 〉, which is calculable, and the measured linear bias,

\begin{eqnarray} \delta _{\rm c}b_{10}^{\rm (msd)}(\nu ) &= \displaystyle\frac{1}{(1-\epsilon _{\times })}\bigg [\frac{\delta _{\rm c}b_1^{\rm (msd)}}{(S_{\times }/S_{0})} - \epsilon _{\times }\left\langle\mu _1|\nu \right\rangle \bigg ] . \end{eqnarray}

(35)

For b₂₀, however, we must deal with a term involving 〈 b₁(ν, β)μ₁(ν, β)|ν 〉. In practice, measuring this term would involve knowing the value of β for each halo. While this is doable in an N-body simulation, it makes the reconstruction technique very cumbersome. In order to keep the algorithm simple, we will assume that this term can be split as 〈b₁(ν, β)μ₁(ν, β)|ν〉 = 〈b₁(ν, β)|ν〉〈μ₁(ν, β)|ν〉 = b^(msd)₁〈μ₁|ν〉. This assumption can be explicitly tested using N-body simulations. With this caveat, the reconstructed quadratic coefficient is

\begin{eqnarray} \delta _{\rm c}^2 b_{20}^{\rm (msd)}(\nu ) &=& \displaystyle\frac{1}{(1-\epsilon _{\times })^2}\bigg [\frac{\delta _{\rm c}^2 b_2^{\rm (msd)}}{(S_{\times }/S_{0})^2}\nonumber \\ &&- \displaystyle 2\epsilon _{\times }\left(\frac{\delta _{\rm c}b_1^{\rm (msd)}}{(S_{\times }/S_{0})}\left\langle\mu _1|\nu \right\rangle - \left\langle\mu _1^2|\nu \right\rangle \right)\nonumber \\ && \displaystyle - \epsilon _{\times }(2-\epsilon _{\times })\left\langle\mu _2|\nu \right\rangle \bigg ]. \end{eqnarray}

(36)

The remaining averages over β simplify to some extent, and we have

\begin{eqnarray} \left\langle\mu _1|\nu \right\rangle &=& \nu \left(\nu +\left\langle\beta |\nu \right\rangle \right),\nonumber \\ \left\langle\mu _2|\nu \right\rangle &=& \left\langle\mu _1^2|\nu \right\rangle - \nu ^2\nonumber \\ &=& \nu ^2\left(\nu ^2-1-2\nu \left\langle\beta |\nu \right\rangle +\left\langle\beta ^2|\nu \right\rangle \right), \end{eqnarray}

(37)

where

\begin{equation} \left\langle\beta ^j|\nu \right\rangle = \int {\rm d}\beta p(\beta |\nu ) \beta ^j, \end{equation}

(38)

with p(β|ν) given in equation (13).

5 DISCUSSION

We have presented a derivation of equation (12), which represents the first analytic model of halo abundances which accounts self-consistently for the fact that the physics of collapse is best described around the halo centre of mass, is not spherical, and can vary substantially from one place in the universe to another. The comparison between equation (12) and halo counts in simulations (Fig. 2) is very encouraging. (Fig. 1 shows that accounting for the physics of non-spherical collapse is necessary.)

We also derived expressions for the associated large-scale bias factors b₁₀ and b₂₀ (these follow from setting ϵ_× → 0 and S_×/S₀ → 1 in equations 29 and 30); the former is in very good agreement with measurements from simulations. We are in the process of determining if b₂₀ is similarly accurate. Additionally, we showed how real-space measurements of b₁ and b₂ can be converted to practical estimates of b₁₀ and b₂₀ (equations 35 and 36), which extends the results of Musso et al. (2012) and Paranjape & Sheth (2012) to the case of stochastic barrier heights. Testing the accuracy of this reconstruction in N-body simulations is also work in progress.

Our ESP analysis highlighted the benefits of developing a model for the abundance of positions in space around which the physics of collapse is simplest. In Section 3, we showed that it can be used to understand what would have happened if one had instead studied all initial positions (Figs 4 and 5). This provides the basis for a number of other tests of the more formal aspects of our approach – as well as elucidating the relation between our ESP calculation and the usual excursion set approach.

Our model must take as input one physical quantity from simulations. This quantity describes how the typical density required for collapse depends on protohalo mass, and how much this density varies from one protohalo to the next. While the triaxial collapse model correctly predicts that this mean value and the scatter around it should increase as halo mass decreases, the predicted scatter is smaller than observed. We therefore hope that our work will motivate studies of the physical origin of p(β), as well as measurements in simulations of how this distribution evolves. If this distribution does not evolve, then our ESP model provides a particularly simple model for how the distribution of halo abundances, when expressed as a function of ν, should be nearly but not exactly universal (see discussion in Musso & Sheth 2012). Fig. 8 shows that the predicted departures from universality in this simple case are in reasonable agreement with measurements.

Figure 8.

Predicted departures from self-similarity in the ESP model, for the WMAP3 cosmology. Top panel shows the predicted halo abundances at z = 0 and 1, expressed as a function of y = [δ_cD(0)/D(z)σ₀(m)]². We have assumed that the distribution p(β) defined in the text does not evolve with redshift. Bottom panel shows the ratio of our predictions to the fitting formulae for these redshifts which Tinker et al. (2008) calibrated from halo abundances measured in their simulations.

Open in new tab Download slide

Returning to the issue of the distribution of β, an obvious direction of future work is that peaks have non-trivial deformation and inertia tensors. Non-spherical initial shapes mean the critical density for collapse will now depend not just on the tidal field but on the shape and the misalignment angle between the principal axes of the two tensors. In principle, our ESP model contains a prescription for the distribution of shapes (and misalignments) for which one should run the ellipsoidal collapse model.

This does raise one caveat: in peaks theory (on which our ESP model is based), the axes of the two tensors are expected to be well aligned, with the direction of maximum compression (due to the tidal field) aligned with the shortest axis (van de Weygaert & Bertschinger 1996). Measurements in simulations do show strong alignment, but the direction of maximum compression lies instead along the longest axis (Porciani et al. 2002; Despali et al. 2013). On the other hand, these measurements are dominated by masses for which δ_c(m)/σ₀(m) < 1, where our ESP model may not apply anyway: there is no physically compelling reason why a peak of (normalized) height 1 should dominate its surroundings. While this mass regime may be less interesting for studies which use the abundance of rich clusters to constrain cosmological models, understanding it is useful for understanding the formation and evolution of galaxies. So we hope our approach of identifying just what it is that is special about the positions around which gravitational collapse occurs, and then using only these positions to make inferences about halo assembly histories, will eventually provide insight into the low mass regime as well.

We thank Christophe Pichon for discussions, and for pointing out to us the work of Hanami (2001), and Ixandra Achitouv for discussions regarding barrier stochasticity. This work is supported in part by NSF-0908241 and NASA NNX11A125G (RKS and AP), and by the Swiss National Science Foundation (VD). AP thanks A. Réfrégier and A. Amara at ETH, Zürich for hospitality during a visit when this work was completed. VD also wishes to thank ICTP for hospitality when parts of this work were completed.

We use the results of Robertson et al. (2009) rather than those of Dalal et al. (2008) because the latter only show a plot whereas Robertson et al. provide a description of the full probability distribution. Both groups find a similar mean and rms for the barrier heights, however, which increase approximately linearly with σ₀.

The factor (S_×/S₀)ⁿ is a consequence of the fact that the excursion set calculation lives in real space rather than Fourier space, and would be present even if the Fourier-space bias were constant.

REFERENCES

Achitouv

Rasera

Sheth

R. K.

Corasaniti

P. S.

. ,

2013

preprint (arXiv:1212.1166)

Appel

Jones

B. J. T.

. ,

MNRAS

1990

, vol.

245

pg.

522

Bardeen

J. M.

Bond

J. R.

Kaiser

Szalay

A. S.

. ,

ApJ

1986

, vol.

304

pg.

Month:	Total Views:
January 2017	6
February 2017	2
April 2017	6
June 2017	2
July 2017	2
August 2017	1
September 2017	3
October 2017	4
November 2017	8
December 2017	6
January 2018	4
February 2018	5
March 2018	9
April 2018	10
May 2018	5
June 2018	11
July 2018	9
August 2018	8
September 2018	5
October 2018	10
November 2018	10
December 2018	5
January 2019	5
February 2019	2
March 2019	10
April 2019	14
May 2019	8
June 2019	9
July 2019	8
August 2019	10
September 2019	8
October 2019	7
November 2019	9
December 2019	2
January 2020	8
February 2020	2
March 2020	6
April 2020	6
May 2020	3
June 2020	5

Article Contents

Abstract

1 INTRODUCTION

2 THE UNCONDITIONAL MASS FUNCTION

2.1 Notation

2.2 Matching filter scales

2.3 Excursion set peaks with a constant barrier

2.4 Beyond the spherical-collapse constant barrier

2.5 Excursion set peaks with a square-root barrier and scatter

3 IMPLICATIONS FOR AVERAGES OVER ALL, RATHER THAN A SPECIAL SUBSET OF POSITIONS

3.1 Distribution of δ for all particles in protohaloes of a given mass

3.2 Distribution of predicted mass for all particles in protohaloes of a given mass

4 HALO BIAS

4.1 Reconstructing the peak-background split parameters

5 DISCUSSION

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only