A Stochastic Theory of the Hierarchical Clustering. II. Halo Progenitor Mass Function and Large-scale Bias

In the standard picture of structure formation, dark matter (DM) halos grow hierarchically, progressively evolving toward larger and more massive units via mergers and accretion (see, e.g., textbooks by Cimatti et al. 2020; Mo et al. 2010). The leading-order halo statistics is constituted by the halo mass function N(M, t), which provides the abundance of halo masses at any given cosmic epoch. However, to better characterize the growth history of halos, a higher-order statistics is needed, namely the progenitor mass function ${ \mathcal N }({M}_{1},{z}_{1}| {M}_{2},{z}_{2});$ the latter specifies, for a descendant halo mass M₂ at final redshift z₂, the conditional distribution of progenitor halo masses M₁ < M₂ at any redshift z₁ > z₂.

Such a quantity can be exploited for a variety of purposes, such as constructing numerical realization of the merging process known as merger trees (see Lacey & Cole 1993; Somerville & Kolatt 1999; Cole et al. 2000; Parkinson et al. 2008; Zhang et al. 2008; Jiang & van den Bosch 2014), computing halo formation time distributions and formation rates (see Lacey & Cole 1993; Kitayama & Suto 1996; Percival & Miller 1999; Percival et al. 2000; Miller et al. 2006; Neistein & Dekel 2008; Zhang et al. 2008; Moreno et al. 2009; Lapi et al. 2013), deriving the average growth history of halos (see Nusser & Sheth 1999; van den Bosch 2002; Miller et al. 2006; Neistein et al. 2006), and investigating the relation of halos with their large-scale environment (see Mo & White 1996; Sheth & Tormen 1999; Musso et al. 2012; Paranjape et al. 2012, 2018).

Clearly, the halo progenitor distribution can be directly extracted from N-body simulations (see Kauffmann et al. 1999; Benson et al. 2000; Wechsler et al. 2002; Helly et al. 2003; Kang et al. 2005; Springel et al. 2005; McBride et al. 2009; Behroozi et al. 2013, 2019). However, this procedure is demanding in terms of time and storing capacity since it requires running a full simulation and analyzing its outputs; moreover, it comes together with some tricky and systematic issues, concerning not only the limited mass and force resolution, but also the algorithms used to identify and relate halos among different simulation snapshots (see, e.g., Harker et al. 2006; Fakhouri & Ma 2008; Fakhouri et al. 2010; Jiang & van den Bosch 2014). Thus a first principles theoretical understanding of the halo progenitor mass function is also very important.

The standard theoretical framework to address such an issue is constituted by the excursion set formalism and can be traced back to the seminal work by Bond et al. (1991) and to the ensuing developments by Lacey & Cole (1993) and Mo & White (1996). In this approach, the issue of counting halos is remapped in finding the first crossing distribution of a random walk that hits a suitable barrier. Specifically, the random walk is executed by the over-density field δ around a given spatial location as a function of the smoothed mass scale σ²(M), which encapsulates the power spectrum of perturbations. When the smoothing is performed with a k-sharp filter in Fourier space, the problem greatly simplifies since the random walk is Markovian. In addition, the barrier is provided by the linear collapse threshold δ_c (t, σ), which depends on cosmic time and possibly on the mass scale (in the latter case, one deals with a so-called 'moving barrier'). When the mass-independent threshold δ_c (t) from the spherical collapse theory is used, the resulting framework is also known as extended Press and Schechter theory (EPS). ⁵ The excursion set framework can be exploited not only to derive the halo mass function, but also other higher-order statistics like the progenitor distribution and bias. These quantities are related to the first crossing distributions for walks that do not start from the origin, and cross the barrier farther away from it on a larger mass scale.

Remarkably, the ensuing progenitor mass distribution is found to be in overall good agreement with the N-body outcomes, especially when a moving barrier with shape inspired by the ellipsoidal collapse is adopted (see Sheth & Tormen 1999, 2002). A further bonus is that the progenitor distribution can be computed via closed-form analytic expressions (see Lacey & Cole 1993; Sheth & Tormen 2002; Mahmood & Rajesh 2005) or via fast semi-analytic algorithms (see Zhang et al. 2008). This is useful both theoretically to transparently understand the physics underlying the hierarchical clustering of halos, and practically to investigate, much faster than in simulations, the halo merger histories under diverse initial conditions and cosmological backgrounds.

Despite these successes and advantages, the excursion set framework is known to hide some pitfalls and drawbacks: the mere idea that one can relate the merging history of halos to the first upcrossing in the abstract δ_c − σ² space constitutes a leap of faith (see Mo et al. 2010); the merging kernel associated to the excursion set theory is not symmetric, and this causes a mathematical inconsistency or at least an ambiguity in the definition of the merger rates (see Benson et al. 2005; Neistein & Dekel 2008; Zhang et al. 2008). The relation between the probability of first upcrossing and the progenitor mass function cannot be strictly correct for individual mass elements (see discussion by Mo et al. 2010, their Section 7.2.2b).

In Paper I (Lapi & Danese 2020) we submitted a new theory of the hierarchical clustering based on stochastic differential equations in real space, which constitutes a change of perspective with respect to the excursion set formalism and can alleviate the above shortcomings. Our previous work focused on the (unconditional) halo mass function. In the present work, we show how to generalize the treatment to the (conditional) progenitor mass distribution and large-scale halo bias.

The plan of the paper is as follows. First, we present a stochastic differential equation that describes fluctuations in the mass growth of progenitor halos of a given descendant mass and redshift, as driven by a multiplicative Gaussian white noise involving the spherical collapse thresholds and the mass variances of the progenitor and descendant halos. In Section 2 we demonstrate that, as cosmic time passes, the noise yields an average drift of the progenitors toward larger masses, which quantitatively renders the expectation from the EPS theory. Then, in Section 2.1 we solve the Fokker–Planck equation associated with the stochastic dynamics, and obtain as a solution the EPS progenitor mass function. We also point out that the solution is stationary when the original equation is written in terms of a convenient scaled variable. In Section 2.2 we highlight how the formalism based on the Fokker–Planck equation allows us to define unambiguously the total derivative of the (conditional and unconditional) mass function in terms of halo collapse versus merging rates. In Section 3 we introduce a minimal modification of the stochastic equation in terms of a mass-dependent collapse threshold modulating the noise, and solve analytically the associated Fokker–Planck equation for the progenitor mass function. The latter is found to be in excellent agreement with the outcomes of N-body simulations. Even more remarkably, all of that is achieved with the same mass-dependent threshold shape used in Paper I to reproduce the unconditional halo mass function. In Section 4 we exploit the progenitor mass function to compute the large-scale halo bias, and show that it reproduces very well the simulation results. Finally, in Section 5 we summarize and briefly discuss our findings. In the Appendix, for the reader's convenience, we recap and discuss the main assumptions and results of Paper I, and we present an additional test concerning the unconditional mass function in warm DM cosmologies.

Unless otherwise explicitly specified, throughout this work, we adopt the standard flat cold dark matter (ΛCDM) cosmology (Planck Collaboration 2020) with rounded parameter values: matter density Ω_M ≈ 0.3, dark energy density Ω_Λ ≈ 0.7, baryon density Ω_b ≈ 0.05, Hubble constant H₀ =100 h km s⁻¹ Mpc⁻¹ with h ≈ 0.7, and mass variance σ₈ ≈ 0.8 on a scale of 8 h⁻¹ Mpc. The most relevant equations are highlighted with a box.

Analogously to the proposal put forward in Paper I (see the Appendix for a summary), we now write down a nonlinear stochastic differential equation aimed at creating the conditional probability of the EPS theory. This reads (overdot means differentiation with respect to t₁)

$$\begin{eqnarray}&&\displaystyle \frac{d}{{{dt}}_{1}}\,\mathrm{ln}{\left({\rm{\Delta }}{\sigma }^{2}\right)}^{-1/2}=\displaystyle \frac{{\left({\rm{\Delta }}{\sigma }^{2}\right)}^{1/2}}{{\rm{\Delta }}{\delta }_{c}}{\left|\displaystyle \frac{{\rm{\Delta }}\dot{{\delta }_{c}}}{{\rm{\Delta }}{\delta }_{c}}\right|}^{1/2}\,\eta ({t}_{1}),\end{eqnarray} \tag{ 1 }$$

or equivalently, in terms of mass ⁶

$$\begin{eqnarray}&&\displaystyle \frac{d}{{{dt}}_{1}}\,{M}_{1}={\left|\displaystyle \frac{{\rm{\Delta }}{\sigma }^{2}}{{\rm{\Delta }}{\delta }_{c}}\right|}^{3/2}\,\displaystyle \frac{2\,| {\rm{\Delta }}{\dot{\delta }}_{c}{| }^{1/2}}{| d{\sigma }_{1}^{2}/{{dM}}_{1}| }\ \eta ({t}_{1}).\end{eqnarray} \tag{ 2 }$$

Here t is the cosmic time, and η(t) is a Gaussian white noise (physical dimension $1/\sqrt{\mathrm{time}}$) with ensemble-average properties 〈η(t)〉 = 0 and $\langle \eta (t)\eta (t^{\prime} )\rangle =2\,{\delta }_{{\rm{D}}}(t-t^{\prime} )$, where δ_D is the Dirac δ-function. In the above, the subscript "2" refers to the quantities at the final cosmic time t₂ of the descendant halo, and the subscript "1" refers to those at the cosmic time t₁ of the progenitor halos; keeping this in mind, we have defined Δδ_c ≡ δ_c (t₁) − δ_c (t₂) and ${\rm{\Delta }}{\sigma }^{2}\equiv {\sigma }^{2}({M}_{1})-{\sigma }^{2}({M}_{2})\,={\sigma }_{1}^{2}-{\sigma }_{2}^{2}$. The quantity δ_c (z) = δ_c0 D(0)/D(z) is the critical threshold for collapse extrapolated from linear perturbation theory; in a flat universe (see Weinberg 2008; Mo et al. 2010), one can use the approximations ${\delta }_{c0}\,\simeq \tfrac{3}{20}\,{(12\pi )}^{2/3}$ $\left[1+0.0123{\mathrm{log}}_{10}{{\rm{\Omega }}}_{M}(z)\right]\approx 1.68$ and $D(z)\,\approx \displaystyle \frac{5}{2}\,\displaystyle \frac{{{\rm{\Omega }}}_{M}(z)}{1\,+\,z}{\left[\displaystyle \frac{1}{70}+\tfrac{209}{140}{{\rm{\Omega }}}_{M}(z)-\tfrac{1}{140}{{\rm{\Omega }}}_{M}^{2}(z)+{{\rm{\Omega }}}_{M}^{4/7}(z)\right]}^{-1}$ with ${{\rm{\Omega }}}_{M}{(z)\equiv {{\rm{\Omega }}}_{M}(1+z)}^{3}/[{{\rm{\Omega }}}_{{\rm{\Lambda }}}+{{\rm{\Omega }}}_{M}{(1+z)}^{3}]$. The quantity σ is the mass variance filtered on the mass scale M:

$$\begin{eqnarray}&&{\sigma }^{2}(M)=\displaystyle \frac{1}{{\left(2\pi \right)}^{3}}\,\int {d}^{3}k\,P(k){\tilde{W}}_{M}^{2}(k),\end{eqnarray} \tag{ 3 }$$

in terms of the power spectrum P(k) and of a Fourier transformed window function ${\tilde{W}}_{M}^{2}(k)$ whose volume in real space encloses the mass M. σ(M) is a deterministic function of M, which, for standard CDM power spectra (e.g., Bardeen et al. 1986), features an inverse, convex, slowly varying behavior. In the present framework, M(t) and σ(t) = σ(M(t)) constitute stochastic variables that fluctuate over cosmic time t under the influence of the noise. A discussion of the key assumptions underlying Equation (1) and a summary of the ensuing stochastic framework are presented in Paper I and summarized in the Appendix.

Assigned the descendant halo mass M₂ at cosmic time t₂ and an initial condition at early times, the above equation can be integrated via the Euler–Maruyama or Euler–Heun methods (e.g., Kloeden & Platen 1992) on a discretized time grid. In Figure 1 we show the resulting evolution of $\mathrm{ln}\sigma ({M}_{1}(t))$ and of $\mathrm{log}{M}_{1}(t)$ as a function of cosmic time, with initial condition M₁ ∼ 10⁴ M_⊙ or $\mathrm{ln}{\sigma }_{1}\sim 2.2$ at z ∼ 10 (reasonable initial conditions do not significantly change the outcome). The colored solid lines surrounded by shaded areas illustrate the median and quartiles over 3000 realizations of the noise, for different descendant masses M₂ ≈ 10⁹–10¹¹–10¹³–10¹⁵ M_⊙ at redshift z₂ ≈ 0. The corresponding colored dashed lines show the expected evolution when the quantity ${\nu }_{12}\equiv {\rm{\Delta }}{\delta }_{c}/\sqrt{{\rm{\Delta }}{\sigma }^{2}}$ stays constant, as expected in the EPS framework. It is found that the noise induces a drift of M₁ (σ₁), making it increase (decrease) and, after a burn-in period needed to erase memory of the initial condition, converge toward the evolution predicted by the EPS theory.

Zoom In Zoom Out Reset image size

2.1. Fokker–Planck Equation and the EPS Mass Function

The probability density ${ \mathcal P }({M}_{1},{t}_{1}| {M}_{2},{t}_{2})$ of finding a progenitor halo with mass between M₁ and M₁ + dM₁ at cosmic time t₁ given a descendant halo with mass M₂ at cosmic time t₂ can be derived by solving the Fokker–Planck equation associated with Equation (2). This reads (see the Appendix in Paper I for details; also textbooks by Risken 1996 and Paul & Baschnagel 2013)

$$\begin{eqnarray}&&{\partial }_{{t}_{1}}{ \mathcal P }=-{{ \mathcal T }}^{2}\,{\partial }_{{M}_{1}}\left[{ \mathcal D }\,({\partial }_{{M}_{1}}{ \mathcal D }){ \mathcal P }\right]+{{ \mathcal T }}^{2}\,{\partial }_{{M}_{1}}^{2}\left[{{ \mathcal D }}^{2}\,{ \mathcal P }\right],\end{eqnarray} \tag{ 4 }$$

where we have defined the two quantities ${ \mathcal D }({M}_{1},{M}_{2})\,\equiv 2{({\rm{\Delta }}{\sigma }^{2})}^{3/2}/| d{\sigma }^{2}/{{dM}}_{1}|$ and ${ \mathcal T }{({t}_{1},{t}_{2})\equiv | \dot{{\rm{\Delta }}}{\delta }_{c}{| }^{1/2}/({\rm{\Delta }}{\delta }_{c})}^{3/2}$ such that in Equation (2), the mass and time dependencies are factorized as ${\dot{M}}_{1}={ \mathcal D }({M}_{1},{M}_{2}){ \mathcal T }({t}_{1},{t}_{2})\eta ({t}_{1})$. The Fokker–Planck equation may also be written as a pure continuity equation ${\partial }_{{t}_{1}}{ \mathcal P }+{\partial }_{{M}_{1}}\,{ \mathcal J }=0$ in terms of a probability current ${ \mathcal J }=-{{ \mathcal T }}^{2}\,{ \mathcal D }\,{\partial }_{{M}_{1}}\,[{ \mathcal D }\,{ \mathcal P }]$. The natural boundary conditions ${\mathrm{lim}}_{{M}_{1}\to \infty }{ \mathcal P }=0$, ${ \mathcal P }({M}_{1},0| {M}_{2},{t}_{2})={\delta }_{{\rm{D}}}({M}_{1})$ and the constraint ${ \mathcal P }({M}_{1},{t}_{1}| {M}_{2},{t}_{2})=0$ whenever M₁ < 0 or M₁ > M₂ must apply; the latter corresponds to reflective barrier conditions ${ \mathcal J }{| }_{{M}_{1}=0,{M}_{1}={M}_{2}}=[{ \mathcal D }\,{\partial }_{{M}_{1}}({ \mathcal D }\,{ \mathcal P })]{| }_{{M}_{1}=0,{M}_{1}={M}_{2}}=0$ at the M₁ = 0 and M₁ = M₂ points (no net probability currents through M₁ = 0 or M₁ = M₂). To solve the Fokker–Planck equation, we employ the transformations:

$$\begin{eqnarray}&&\left\{\begin{array}{l}X\equiv \displaystyle \int \displaystyle \frac{{{dM}}_{1}}{{ \mathcal D }({M}_{1},{M}_{2})}=\displaystyle \frac{1}{{\left({\rm{\Delta }}{\sigma }^{2}\right)}^{1/2}}\\ Y\equiv \displaystyle \int {{dt}}_{1}\,{{ \mathcal T }}^{2}({t}_{1},{t}_{2})=\displaystyle \frac{1}{2\,{\left({\rm{\Delta }}{\delta }_{c}\right)}^{2}}\\ { \mathcal W }\equiv { \mathcal D }\,{ \mathcal P }.\end{array}\right.\end{eqnarray} \tag{ 5 }$$

Then Equation (4) turns into

$$\begin{eqnarray}&&{\partial }_{Y}\,{ \mathcal W }={\partial }_{X}^{2}\,{ \mathcal W },\end{eqnarray} \tag{ 6 }$$

which is a standard diffusion equation with boundary conditions ${\mathrm{lim}}_{X\to \infty }{ \mathcal W }=0$, ${ \mathcal W }(X,0)={\delta }_{{\rm{D}}}(X)$ and $({\partial }_{X}\,{ \mathcal W }){| }_{X=0,X=\infty }\,=0$, whose solution just writes ${ \mathcal W }={(\pi Y)}^{-1/2}\,{e}^{-{X}^{2}/4Y}$. Coming back to the original variables, we get the EPS conditional probability

$$\begin{eqnarray}&&{ \mathcal P }({M}_{1},{t}_{1}| {M}_{2},{t}_{2})=\displaystyle \frac{1}{\sqrt{2\pi }}\,\displaystyle \frac{{\rm{\Delta }}{\delta }_{c}}{{\left({\rm{\Delta }}{\sigma }^{2}\right)}^{3/2}}\left|\displaystyle \frac{d{\sigma }^{2}}{{{dM}}_{1}}\right|\,{e}^{-{\left({\rm{\Delta }}{\delta }_{c}\right)}^{2}/2{\rm{\Delta }}{\sigma }^{2}}.\end{eqnarray} \tag{ 7 }$$

Note that in the limit δ_c (t₂) → 0 and σ(M₂) → 0, the above probability converges to

$$\begin{eqnarray}&&{ \mathcal P }(M,t)=\sqrt{\displaystyle \frac{2}{\pi }}\,\displaystyle \frac{{\delta }_{c}}{{\sigma }^{2}}\left|\displaystyle \frac{d\sigma }{{dM}}\right|\,{e}^{-{\left({\delta }_{c}\right)}^{2}/2{\sigma }^{2}},\end{eqnarray} \tag{ 8 }$$

which is related to the unconditional mass function by $N(M,t)=(\bar{\rho }/M){ \mathcal P }(M,t)$ in terms of the co-moving matter density $\bar{\rho };$ as pointed out in Paper I, this is exactly the Press and Schechter formula.

We can derive the same results in another way, just by recasting Equation (1) in terms of the scaled (or self-similar) variable ${\nu }_{12}\equiv {\rm{\Delta }}{\delta }_{c}/\sqrt{{\rm{\Delta }}{\sigma }^{2}}$ as

$$\begin{eqnarray}&&\displaystyle \frac{d}{{{dt}}_{1}}\,{\nu }_{12}=-\left|\displaystyle \frac{{\rm{\Delta }}{\dot{\delta }}_{c}}{{\rm{\Delta }}{\delta }_{c}}\right|\,{\nu }_{12}+{\left|\displaystyle \frac{{\rm{\Delta }}{\dot{\delta }}_{c}}{{\rm{\Delta }}{\delta }_{c}}\right|}^{1/2}\,\eta ({t}_{1}).\end{eqnarray} \tag{ 9 }$$

The corresponding Fokker–Planck equation for the probability ${ \mathcal P }({\nu }_{12},t)$ reads

$$\begin{eqnarray}&&{\partial }_{{t}_{1}}\,{ \mathcal P }=\left|\displaystyle \frac{{\rm{\Delta }}{\dot{\delta }}_{c}}{{\rm{\Delta }}{\delta }_{c}}\right|\,{\partial }_{{\nu }_{12}}\left[{\nu }_{12}\,{ \mathcal P }+{\partial }_{{\nu }_{12}}\,{ \mathcal P }\right],\end{eqnarray} \tag{ 10 }$$

with boundary conditions ${\mathrm{lim}}_{{\nu }_{12}\to \infty }{ \mathcal P }=0$ and ${ \mathcal J }{| }_{{\nu }_{12}=0}\,=-| \dot{{\rm{\Delta }}}{\delta }_{c}/{\rm{\Delta }}{\delta }_{c}| \,[{\nu }_{12}\,{ \mathcal P }+{\partial }_{{\nu }_{12}}\,{ \mathcal P }]{| }_{{\nu }_{12}=0}=0$, implying ${\int }_{0}^{\infty }d{\nu }_{12}{ \mathcal P }=1.$

Under the ergodic hypothesis, we can focus on the stationary solution ${ \mathcal P }({\nu }_{12},{t}_{1})=\bar{{ \mathcal P }}({\nu }_{12});$ this is determined by ${\partial }_{{t}_{1}}\,{ \mathcal P }({\nu }_{12},{t}_{1})=0$ or

$$\begin{eqnarray}&&{\nu }_{12}\,\bar{{ \mathcal P }}+\displaystyle \frac{d}{d{\nu }_{12}}\,\bar{{ \mathcal P }}=0,\end{eqnarray} \tag{ 11 }$$

where the constant on the right-hand side (rhs) must be zero to satisfy the no-current boundary condition ${ \mathcal J }{| }_{{\nu }_{12}=0}=0$. The ensuing solution is just $\bar{{ \mathcal P }}({\nu }_{12})=\sqrt{2/\pi }\,{e}^{-{\nu }_{12}^{2}/2}$. The conditional probability function reads

$$\begin{eqnarray}&&{ \mathcal P }({M}_{1},{t}_{1}| {M}_{2},{t}_{2})=\left|\displaystyle \frac{d{\nu }_{12}}{{{dM}}_{1}}\right|\,\bar{{ \mathcal P }}({\nu }_{12})\,;\end{eqnarray} \tag{ 12 }$$

since $| d{\nu }_{12}/{{dM}}_{1}| ={\rm{\Delta }}{\delta }_{c}/2{({\rm{\Delta }}{\sigma }^{2})}^{3/2}\times | d{\sigma }_{1}^{2}/{{dM}}_{1}|$, this again boils down to the EPS probability of Equation (7).

Note that often in the literature, the quantity ${\nu }_{12}\,\bar{{ \mathcal P }}({\nu }_{12})$ is referred to as the conditional multiplicity function, while the quantity $f(\mathrm{log}{M}_{1}| {M}_{2})\equiv { \mathcal P }({M}_{1},{t}_{1}| {M}_{2},{t}_{2}){M}_{1}\,\mathrm{ln}(10)$ is the mass fraction in progenitor halos per logarithmic mass bin. Finally, the (conditional) progenitor mass function ${ \mathcal N }({M}_{1},{t}_{1}| {M}_{2},{t}_{2})$, which represents the ensemble-averaged number of progenitors with mass in the range M₁ and M₁ + dM₁ at time t₁ for a halo of mass M₂ at current time t₂ is given by

$$\begin{eqnarray}&&{ \mathcal N }({M}_{1},{t}_{1}| {M}_{2},{t}_{2})=\displaystyle \frac{{M}_{2}}{{M}_{1}}\,{ \mathcal P }({M}_{1},{t}_{1}| {M}_{2},{t}_{2})=\displaystyle \frac{{M}_{2}}{{M}_{1}}\left|\displaystyle \frac{d{\nu }_{12}}{{{dM}}_{1}}\right|\,\bar{{ \mathcal P }}({\nu }_{12}).\end{eqnarray} \tag{ 13 }$$

In Figures 2 and 3 we illustrate, as dotted lines, the conditional multiplicity function and the mass fraction in progenitors derived above. As is well known from EPS theory, such conditional mass functions tend to exceed that the outcomes of N-body simulations at the low-mass end, and to fall short of them at the high-mass end.

Zoom In Zoom Out Reset image size

2.2. Drift versus Diffusion

Before proceeding further toward a better agreement with simulations, it is interesting to highlight the relative role played by the two terms appearing on the rhs of the Fokker–Planck equation in determining the total derivative of the (conditional and unconditional) mass function. To this purpose, after Equations (4) and (13), one can write

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\partial }_{{t}_{1}}\,{ \mathcal N }={\left({\partial }_{{t}_{1}}{ \mathcal N }\right)}_{\mathrm{DRIFT}}+{\left({\partial }_{{t}_{1}}{ \mathcal N }\right)}_{\mathrm{DIFF}},\\ \quad {\left({\partial }_{{t}_{1}}{ \mathcal N }\right)}_{\mathrm{DRIFT}}=-\displaystyle \frac{{{ \mathcal T }}^{2}}{{M}_{1}}\,{\partial }_{{M}_{1}}\left[{ \mathcal D }\,({\partial }_{{M}_{1}}{ \mathcal D }){M}_{1}\,{ \mathcal N }\right],\\ \quad {\left({\partial }_{{t}_{1}}{ \mathcal N }\right)}_{\mathrm{DIFF}}=+\displaystyle \frac{{{ \mathcal T }}^{2}}{{M}_{1}}\,{\partial }_{{M}_{1}}^{2}\left[{{ \mathcal D }}^{2}\,{M}_{1}\,{ \mathcal N }\right].\end{array}\right.\end{eqnarray} \tag{ 14 }$$

where the quantities ${ \mathcal D }$ and ${ \mathcal T }$ have been defined below Equation (4). Here, the term ${({\partial }_{{t}_{1}}{ \mathcal N })}_{\mathrm{DRIFT}}$ stands for drift, and can be interpreted as the direct hierarchical collapse of larger and more massive perturbations; the term ${({\partial }_{{t}_{1}}{ \mathcal N })}_{\mathrm{DIFF}}$ stands for diffusion, and can be interpreted as describing the effect of stochastic merging events.

In Figure 4 we illustrate the contributions (in absolute value) of the drift (solid lines) and diffusion (dashed line) terms to the time derivative of the conditional mass function ${\partial }_{{t}_{1}}{ \mathcal N };$ these are plotted as a function of the progenitor mass M₁ at various progenitor redshifts z₁ (color-coded), for different descendant masses M₂ ≈ 10¹¹ (top left panel), 10¹³ (top right panel), and 10¹⁵ M_⊙ (bottom left panel) at z₂ ≈ 0. The bottom right panel shows the same for the unconditional mass function.

Zoom In Zoom Out Reset image size

The drift term is always positive (despite the minus sign appearing in Equation (14)), while the diffusion term is positive at the high-mass end, and becomes negative toward lower masses. The overall effect of the former is to cause an upward drift of the mass function toward larger masses, while that of the latter is to reshape the mass function by subtracting halos from the low-mass end and incorporating them into more massive objects. This is indeed what is to be expected from the combined action of direct collapse and stochastic mergers.

All in all, the total time derivative of the (conditional and unconditional) mass function is positive and dominated by diffusion at the high-mass end (new halos originate from merging of smaller ones). It is positive and dominated by drift at intermediate masses (new halos collapse from the perturbation field), while it becomes negative and is again dominated by diffusion toward the low-mass end (small masses are destroyed by merging into larger ones). Finally, the impact of diffusion is more relevant toward relatively larger masses and later cosmic times, while the drift term becomes increasingly important over a progressively extended mass range at higher redshifts.

We stress that the separation of the total time derivative of the mass function in a drift versus a diffusion term, and the related physical interpretation in terms of direct collapse versus merging, is exclusively allowed by our stochastic theory based on the Fokker–Planck equation. Similar arguments in the standard EPS framework prove to be difficult (e.g., Benson et al. 2005; Neistein & Dekel 2008; Zhang et al. 2008; Lapi et al. 2013), and actually spawn tricky ambiguities in the definition of the halo collapse and merging rates. An interesting future development, which goes beyond the scope of the present paper, would be to construct refined merger trees based on the above Fokker–Planck separation between the direct collapse and merging rates. Speculating somewhat, this could help alleviate some tension between the canonic merger-tree-based approaches to galaxy formation and high-redshift observations (especially related to massive, strongly star-forming systems).

In Paper I we showed that the N-body unconditional mass function can be reproduced in our stochastic approach by including in the multiplicative noise term a mass-dependent collapse threshold ⁷ with shape

$$\begin{eqnarray}&&{\delta }_{c}(\sigma ,t)=\sigma \,\sqrt{q}\,\nu \left[1+\displaystyle \frac{\beta }{{\left(\sqrt{q}\nu \right)}^{2\gamma }}\right],\end{eqnarray} \tag{ 15 }$$

in terms of three parameters q, β, and γ. Comparison with the N-body outcomes (see Sheth & Tormen 1999; Bhattacharya et al. 2011; Watson et al. 2013) indicates the approximate values q ≈ 0.65, β ≈ 0.15, and γ ≈ 0.35. Now we generalize the approach to derive the conditional probability function and show that, with the same shape of the collapse threshold and parameter values given above, it performs very well in reproducing the progenitor distribution extracted from N-body simulations.

The basic idea is just to replace Δδ_c ≡ δ_c (t₁) − δ_c (t₂) with ΔB(ν₁₂) ≡ δ_c (σ₁, t₁) − δ_c (σ₂, t₂) in the multiplicative noise term of Equation (1); when formulated in terms of the scaled variable ${\nu }_{12}\equiv {\rm{\Delta }}{\delta }_{c}/\sqrt{{\rm{\Delta }}{\sigma }^{2}}$, one gets

$$\begin{eqnarray}&&\displaystyle \frac{d}{{{dt}}_{1}}\,{\nu }_{12}=-{\nu }_{12}\left|\displaystyle \frac{{\rm{\Delta }}{\dot{\delta }}_{c}}{{\rm{\Delta }}{\delta }_{c}}\right|+\displaystyle \frac{{\rm{\Delta }}{\delta }_{c}}{{\rm{\Delta }}B({\nu }_{12})}{\left|\displaystyle \frac{{\rm{\Delta }}{\dot{\delta }}_{c}}{{\rm{\Delta }}{\delta }_{c}}\right|}^{1/2}\,\eta ({t}_{1}),\end{eqnarray} \tag{ 16 }$$

where ${\rm{\Delta }}B{({\nu }_{12})\equiv {\rm{\Delta }}{\delta }_{c}[{B}_{0}+{B}_{1}(1+{B}_{2}/{\nu }_{12}^{2})}^{\gamma }]$ can be expressed in terms of the coefficients

$$\begin{eqnarray}\left\{\begin{array}{rcl}{B}_{0} & \equiv & \sqrt{q}-\beta \,\displaystyle \frac{\sqrt{q}\,{\delta }_{2}}{{\rm{\Delta }}{\delta }_{c}}{\left(\displaystyle \frac{{\sigma }_{2}}{\sqrt{q}{\delta }_{2}}\right)}^{2\gamma }\\ {B}_{1} & \equiv & \beta \,\displaystyle \frac{\sqrt{q}\,{\delta }_{1}}{{\rm{\Delta }}{\delta }_{c}}{\left(\displaystyle \frac{{\sigma }_{2}}{\sqrt{q}{\delta }_{1}}\right)}^{2\gamma }\\ {B}_{2} & \equiv & {\left(\displaystyle \frac{{\rm{\Delta }}{\delta }_{c}}{{\sigma }_{2}}\right)}^{2}.\end{array}\right.\end{eqnarray} \tag{ 17 }$$

Note that the spherical collapse threshold independent of mass is recovered in the limit q → 1 and β → 0, which readily imply B₀ ≃ 1 and B₁ ≃ 0; then ΔB(ν₁₂) ≃ Δδ_c holds, and Equation (16) collapses into Equation (9).

The corresponding Fokker–Planck equation reads:

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{{t}_{1}}\,{ \mathcal P } & = & \left|\displaystyle \frac{{\rm{\Delta }}{\dot{\delta }}_{c}}{{\rm{\Delta }}{\delta }_{c}}\right|\,{\partial }_{{\nu }_{12}}\left\{{\nu }_{12}\,{ \mathcal P }\right.\\ & & \left.+\displaystyle \frac{{\rm{\Delta }}{\delta }_{c}}{{\rm{\Delta }}B({\nu }_{12})}\,{\partial }_{{\nu }_{12}}\left[\displaystyle \frac{{\rm{\Delta }}{\delta }_{c}}{{\rm{\Delta }}B({\nu }_{12})}\,{ \mathcal P }\right]\right\},\end{array}\end{eqnarray} \tag{ 18 }$$

and the relevant stationary solution ${ \mathcal P }=\bar{{ \mathcal P }}({\nu }_{12})$ is determined by

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}{\delta }_{c}}{{\rm{\Delta }}B({\nu }_{12})}\,\displaystyle \frac{d}{d{\nu }_{12}}\left[\displaystyle \frac{{\rm{\Delta }}{\delta }_{c}}{{\rm{\Delta }}B({\nu }_{12})}\,\bar{{ \mathcal P }}\right]=-{\nu }_{12}\,\bar{{ \mathcal P }},\end{eqnarray} \tag{ 19 }$$

where the integration constant must be null to satisfy the no-current boundary condition ${ \mathcal J }{| }_{{\nu }_{12}=0}=0$. The above equation can be easily solved by multiplying both sides by Δδ_c /ΔB(ν₁₂) and recognizing that it becomes separable for the function $[{\rm{\Delta }}{\delta }_{c}/{\rm{\Delta }}B({\nu }_{12})]\,\bar{{ \mathcal P }};$ we find

$$\begin{eqnarray}&&\begin{array}{l}\bar{{ \mathcal P }}({\nu }_{12})={{ \mathcal A }}_{12}\,\displaystyle \frac{{\rm{\Delta }}B({\nu }_{12})}{{\rm{\Delta }}{\delta }_{c}}\\ \quad \times \exp \left\{-\displaystyle \int d{\nu }_{12}\,{\nu }_{12}{\left[\displaystyle \frac{{\rm{\Delta }}B({\nu }_{12})}{{\rm{\Delta }}{\delta }_{c}}\right]}^{2}\right\},\end{array}\end{eqnarray} \tag{ 20 }$$

where the normalization constant ${{ \mathcal A }}_{12}$ is determined by the condition ${\int }_{0}^{\infty }d{\nu }_{12}\,\bar{{ \mathcal P }}({\nu }_{12})=1$. Inserting the full shape of ΔB(ν₁₂) and performing explicitly the integration yields the closed form expression

$$\begin{eqnarray}&&\begin{array}{l}{\bar{{ \mathcal P }}}_{\mathrm{CMF}}({\nu }_{12})={{ \mathcal A }}_{12}\,{B}_{0}\left[1+\displaystyle \frac{{B}_{1}}{{B}_{0}}{\left(1+\displaystyle \frac{{B}_{2}}{{\nu }_{12}^{2}}\right)}^{\gamma }\right]\,\exp \left\{-\displaystyle \frac{{B}_{0}^{2}}{2}\,{\nu }_{12}^{2}\left[1+\right.\right.\\ \quad \left.+\,{\left(\displaystyle \frac{{B}_{1}}{{B}_{0}}\right)}^{2}{\left(\displaystyle \frac{1+{B}_{2}/{\nu }_{12}^{2}}{1+{\nu }_{12}^{2}/{B}_{2}}\right)}^{2\gamma }\,\displaystyle \frac{{}_{2}{F}_{1}(1-2\gamma ,-2\gamma ,2-2\gamma ;-{\nu }_{12}^{2}/{B}_{2})}{1-2\gamma }+\\ \qquad \left.+\,2\,\displaystyle \frac{{B}_{1}}{{B}_{0}}{\left(\displaystyle \frac{1+{B}_{2}/{\nu }_{12}^{2}}{1+{\nu }_{12}^{2}/{B}_{2}}\right)}^{\gamma }\,\displaystyle \frac{{}_{2}{F}_{1}(1-\gamma ,-\gamma ,2-\gamma ;-{\nu }_{12}^{2}/{B}_{2})}{1-\gamma }\right]\right\}\,.\end{array}\end{eqnarray} \tag{ 21 }$$

In the above, ₂ F₁ is the ordinary hypergeometric function, which is defined by the power series ${}_{2}{F}_{1}{(a,b,c;x)\equiv {\sum }_{n=0}^{\infty }{(a)}_{n}({b}_{n}){x}^{n}/(c)}_{n}\,n!$ in terms of the Pochammer's symbols (q)₀ = 1 and (q)_n = q(q + 1)...(q + n − 1) for any integer n > 0. Note that in the limit x → 0, one has ₂ F₁(a, b, c; x) ≃ 1 + a b x/c + ....

To recover the unconditional probability (linked to the halo mass function) from the above expression, it is sufficient to take the limits σ₂ → 0 and δ₂ → 0. In such a case, ${B}_{0}\simeq \sqrt{q}$, B₁ ≃ 0, B₂ → ∞ but ${B}_{1}\,{B}_{2}^{\gamma }\simeq \beta {(\sqrt{q})}^{1-2\gamma }$, to imply ν₁₂ = δ₁/σ₁ ≡ ν and

$$\begin{eqnarray}&&\begin{array}{l}{\bar{{ \mathcal P }}}_{\mathrm{MF}}(\nu )\simeq { \mathcal A }\,\sqrt{q}\left[1+\displaystyle \frac{\beta }{{\left(\sqrt{q}\nu \right)}^{2\gamma }}\right]\\ \quad \times \,\exp \left\{-q\,\displaystyle \frac{{\nu }^{2}}{2}\left[1+\displaystyle \frac{2\,\beta }{1-\gamma }\,\displaystyle \frac{1}{{\left(\sqrt{q}\nu \right)}^{2\gamma }}\right.\right.\\ \qquad \left.\left.+\,\displaystyle \frac{{\beta }^{2}}{1-2\,\gamma }\,\displaystyle \frac{1}{{\left(\sqrt{q}\nu \right)}^{4\gamma }}\right]\right\},\end{array}\end{eqnarray} \tag{ 22 }$$

which is the result derived in Paper I (see also the Appendix).

In Figure 2 we illustrate the conditional multiplicity function ${\nu }_{12}\,\bar{{ \mathcal P }}({\nu }_{12})$ as a function of the scaled, self-similar variable ν₁₂. The expectation from our stochastic theory based on the Fokker–Planck equation is compared to the outcomes from the Millennium N-body simulations for FoF-identified halos by Cole et al. (2008). We find a reasonable agreement, and stress that this is obtained with the same mass-dependent threshold (i.e., the same parameters q, β, and γ) used in Paper I to reproduce the halo mass function. As a reference, in Figure 2 we also illustrate the results expected from the excursion set theory for a spherical collapse (EPS) and for a mass-dependent collapse threshold (moving barrier). As is well known, the former exceeds the N-body results at the low-mass end, and falls short at the high-mass end; the latter performs quite well in the low-/intermediate-mass range, but exceeds somewhat at the high-mass end.

In Figure 3 we show the mass fraction in progenitor halos $f(\mathrm{log}{M}_{1}| {M}_{2})\equiv { \mathcal P }({M}_{1},{t}_{1}| {M}_{2},{t}_{2})$ $\times {M}_{1}\,\mathrm{ln}(10)$ as a function of the progenitor to descendant mass ratio M₁/M₂, for different descendant masses M₂ ∼ 10¹², 3 × 10¹³, 10¹⁵ M_⊙ at z₂ ≈ 0 and different progenitor redshifts z₁ ∼ 0.5, 1, 2, 4. As above, we find an overall good agreement between the expectations from our stochastic theory based on the Fokker–Planck equation and the outcomes by the Millennium N-body simulations by Cole et al. (2008). Note that in the latter, the merger trees have been self-consistently constructed by linking FoF groups that contain one or more sub-halos between successive time steps (more details can be found in Springel et al. 2005).

Zoom In Zoom Out Reset image size

Two remarks are in order here. First, all of the analytic results have been computed with the same cosmological parameters of the Millennium simulations for fair comparison. Second, for self-consistency with the comparison pursued in Paper I for the unconditional mass function, we confront the stationary solutions of Equation (18) from our stochastic theory to the outcomes of simulations where halos are identified via FoF-based algorithms, since both tend to produce closely universal functions in terms of the scaled, self-similar variable ν₁₂. As mentioned in Paper I, more general (yet not analytic) solutions from our stochastic theory constitute transitional states that, though converging to the stationary one in the long time limit, could be possibly related to deviation of the mass function from the self-similar shape. A detailed comparison of these with the outcomes of numerical simulations based on spherical over-density algorithm for halos identification would be welcome but is beyond our scope here.

Given the conditional distribution derived in Section 3, it is straightforward to compute the large-scale halo bias, which is routinely exploited to evaluate the correlation between halo abundances and the surrounding environment. The (Eulerian) halo bias can be written as (see Mo & White 1996; Sheth & Tormen 1999)

$$\begin{eqnarray}&&b(M,t)=1+\displaystyle \frac{1}{{\delta }_{c0}}\left[\displaystyle \frac{N(M,t| {M}_{0},{t}_{0})}{N(M,t)V}-1\right],\end{eqnarray} \tag{ 23 }$$

in terms of the ratio between the number of halos at time t (density contrast δ_c ) that will end up at time t₀ (density contrast δ_c0) in an environment with volume V and mass ${M}_{0}\equiv \bar{\rho }V$.

Recalling from Paper I (see also the Appendix) the expression of the halo mass function $N(M,t)=(\bar{\rho }/M)$ ${\bar{{ \mathcal P }}}_{\mathrm{MF}}(\nu )| d\nu /{dM}|$ and using Equation (13), one can recast Equation (23) in the form

$$\begin{eqnarray}\begin{array}{rcl}b(M,z) & = & 1+\displaystyle \frac{1}{{\delta }_{c0}}\left[\displaystyle \frac{{\nu }_{\mathrm{CMF}}\,{\bar{{ \mathcal P }}}_{\mathrm{CMF}}({\nu }_{\mathrm{CMF}})}{\nu \,{\bar{{ \mathcal P }}}_{\mathrm{MF}}(\nu )}-1\right]\\ & \simeq & 1-\displaystyle \frac{1}{{\delta }_{c}}\left[1+\displaystyle \frac{\nu \,{\bar{{ \mathcal P }}}_{\mathrm{MF}}^{\prime} (\nu )}{{\bar{{ \mathcal P }}}_{\mathrm{MF}}(\nu )}\right],\end{array}\end{eqnarray} \tag{ 24 }$$

where ${\nu }_{\mathrm{CMF}}\equiv ({\delta }_{c}-{\delta }_{c0})/\sqrt{{\sigma }^{2}(M)-{\sigma }^{2}({M}_{0})}$ and ν ≡ δ_c (t)/σ(M). The last equality approximately holds on large scales where one can assume δ_c0 ≪ δ_c , σ(M₀) ≪ σ(M) and ν₀ ≡ δ_c0/σ(M₀) ≫ 1, so that at leading order, ν_CMF ≃ ν − (δ_c0/δ_c )ν and ${\bar{{ \mathcal P }}}_{\mathrm{CMF}}({\nu }_{\mathrm{CMF}})\simeq {\bar{{ \mathcal P }}}_{\mathrm{MF}}(\nu )-({\delta }_{c0}/{\delta }_{c})\nu \,{\bar{{ \mathcal P }}}_{\mathrm{MF}}^{\prime} (\nu )$, where the prime means differentiation with respect to ν. All in all, using the explicit expression of the unconditional probability Equation (22), we obtain

$$\begin{eqnarray}\begin{array}{rcl}b(M,z) & \simeq & 1+\displaystyle \frac{1}{{\delta }_{c}}\left\{-1+\displaystyle \frac{2\beta \gamma }{\beta +{\left(\sqrt{q}\nu \right)}^{2\gamma }}\right.\\ & & \left.+{\left(\sqrt{q}\nu \right)}^{2}{\left[1+\displaystyle \frac{\beta }{{\left(\sqrt{q}\nu \right)}^{2\gamma }}\right]}^{2}\right\}\,.\end{array}\end{eqnarray} \tag{ 25 }$$

In Figure 5 we illustrate the large-scale halo bias as a function of the variable ν ≡ δ_c /σ, and compare it with the outcomes from N-body simulations by Tinker et al. (2005) for halos identified by an FoF algorithm, and by Tinker et al. (2010) for halos identified by a spherical over-density algorithm with nonlinear threshold Δ = 200. We find a remarkable agreement with both determinations. In comparison, the EPS outcome for spherical collapse exceeds the simulation results for large ν (high masses and/or early times) and falls short for small ν. The excursion set with a moving barrier performs quite well for large ν but exceeds somewhat the N-body outcomes for ν ≲ 1.

Zoom In Zoom Out Reset image size

We have generalized the stochastic theory of hierarchical clustering presented in Paper I by Lapi & Danese (2020), so as to derive the (conditional) halo progenitor mass function and the related large-scale bias. In Section 2 we presented a stochastic differential equation that describes fluctuations in the mass growth of progenitor halos of given descendant mass and redshift, as driven by a multiplicative Gaussian white noise involving the spherical collapse thresholds and the mass variances of the progenitor and descendant halos. By numerically integrating the stochastic equation, we have demonstrated that, as cosmic time passes, the noise yields an average drift of the progenitors toward larger masses, which quantitatively renders the expectation from the standard EPS theory. In Section 2.1 we solved the Fokker–Planck equation associated with the stochastic dynamics, and obtained, as an exact solution, the EPS progenitor mass function. We have also pointed out that the solution is stationary when written in terms of a convenient, scaled variable. In Section 2.2 we highlighted how the formalism based on the Fokker–Planck equation allows us to define unambiguously the total derivative of the (conditional and unconditional) mass function in terms of halo direct collapse versus merging rates, naturally solving an otherwise tricky issue in the excursion set formalism.

In Section 3 we introduced a modification of the stochastic equation in terms of a mass-dependent collapse threshold modulating the noise. We solved analytically the associated Fokker–Planck equation for the progenitor mass function. The latter is found to be in excellent agreement with the outcome of N-body simulations. Even more remarkably, all of that was achieved with the same threshold as that used in Paper I to reproduce the halo mass function. In Section 4 we exploited the progenitor mass function and the halo mass function to compute the large-scale halo bias, finding it in pleasing agreement with the N-body outcomes.

All in all, in the present paper we have illustrated that the stochastic theory of hierarchical clustering introduced in Paper I can describe effectively not only the halos' abundance, but also their progenitor distribution and their correlation with the large-scale environment across cosmic times. Future developments will be focused on applying our stochastic theory to derive the statistics of cosmic voids, to build up refined merger-tree algorithms, and to investigate the distribution and clustering of sub-halos. It would be also worth exploring whether the same stochastic approach can be usefully applied to describe other phenomena in the modern astrophysical and cosmological landscape, such as galaxy formation or dark energy phenomenology.

We acknowledge the anonymous referee for a constructive report. We warmly thank C. Baccigalupi and T. Ronconi for helpful comments and critical reading. This work has been partially supported by PRIN MIUR 2017 prot. 20173ML3WW 002, "Opening the ALMA window on the cosmic evolution of gas, stars and supermassive black holes." A.L. acknowledges the MIUR grant "Finanziamento annuale individuale attivitá base di ricerca" and the EU H2020-MSCA-ITN-2019 Project 860744 "BiD4BEST: Big Data applications for black hole Evolution STudies."

In this appendix we recap the assumptions and results from Paper I (Lapi & Danese 2020), and provide a further test concerning the unconditional halo mass function in warm dark matter (WDM) cosmologies.

According to the standard cosmological framework, DM halos are thought to originate by the collapse of patches from an initial, closely Gaussian perturbation field. However, as demonstrated by many extensive N-body simulations (e.g., see textbook by Mo et al. 2010 and references therein), the detailed evolution of perturbations and proto-halo patches ultimately depends on a variety of effects. First, the role of initial conditions is crucial, in that a perturbation is more prone to collapse if resides within a sufficiently over-dense region of the initial density field. This was actually the foundational idea of pioneering estimates for the halo abundance (see Press and Schechter 1976), subsequently refined in terms of the excursion set approach (see Bond et al. 1991; Lacey & Cole 1993; Mo & White 1996) to avoid the double counting of over-dense regions overlapped with, or embedded within, larger collapsing ones (the so-called cloud-in-cloud issue). In addition, the shape of proto-halo patches may be also relevant, in that they tend to be ellipsoidal in shape, and this may influence the collapse efficiency and timescale at different mass scales (see Sheth & Tormen 1999, 2002). Moreover, collapse locations may be special points in the initial perturbation field (see Bardeen et al. 1986; Sheth et al. 2001; Dalal et al. 2008; Ludlow & Porciani 2011; Musso & Sheth 2012; Paranjape & Sheth 2012; Lapi & Danese 2014; Hahn & Paranjape 2014), such as peaks in density or in energy (see Musso & Sheth 2019). Finally, other nonlinear effects may influence the collapse of perturbations, such as the local environment, tidal forces, angular momentum, velocity fields, clumpiness, and dynamical friction.

Thus it should appear evident that the collapse and evolution of DM (proto)halos constitute an inherently stochastic process, originating not only by a degree of randomness in the initial conditions, but also by the complexity of deterministic processes affecting the ensuing collapse. As a consequence, the fine details of the evolution for individual proto-halos at different spatial locations and cosmic times are, for all practical purposes, difficult to follow and/or model ab initio in (semi-)analytic terms. Nevertheless, in Paper I we submitted that, if one is interested mainly in the statistical properties of the halo population as a whole, an effective description could be conveniently obtained in terms of stochastic differential equations with appropriate noise.

The situation is somewhat analogous to the classic description of Brownian motion: a microscopic particle immersed in a fluid continuously undergoes collisions with the fluid molecules; the resulting motion x_t , despite being deterministic, appears to be random at the macroscopic level, especially to an external observer who has no access to the exact positions and velocities of the (innumerable) fluid molecules and to the initial conditions of the particle. In the way of a statistical description, the problem is effectively treated via a stochastic differential equation ${\dot{x}}_{t}\sim \kappa \,\eta (t)$, in terms of a diffusion constant κ and of a fluctuating white noise η(t), which allows us to implicitly account for the complex microscopic dynamics of the system. Note that often the system's state influences the intensity of the driving noise, like when the Brownian fluctuations of a microscopic particle near a wall are reduced by hydrodynamic interactions, so that the noise becomes multiplicative in terms of a nonuniform diffusion coefficient κ = κ(x_t ). Similar stochastic models with multiplicative noise have been employed to describe a wide range of physical phenomena, from Brownian motion in inhomogeneous media or in close approach to physical barriers, to thermal fluctuations in electronic circuits, to the evolution of stock prices, to computer science, to the heterogeneous response of biological systems and randomness in gene expression (e.g., Risken 1996; Reed & Jorgensen 2004; Mitzenmacher 2004; Paul & Baschnagel 2013).

Our proposal in Paper I was to adopt a similar description to capture the essence of the proto-halo collapse via the following nonlinear stochastic differential equation:

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\mathrm{ln}{\sigma }^{-1}({M}_{t})=\displaystyle \frac{\sigma ({M}_{t})}{{\delta }_{c}(t)}{\left|\displaystyle \frac{\dot{D}}{D}\right|}^{1/2}\,\eta (t),\end{eqnarray} \tag{ A1 }$$

where t is the cosmic time (corresponding to redshift z), η(t) is a Gaussian white noise with ensemble-averaged properties 〈η〉 = 0 and $\langle \eta (t)\eta (t^{\prime} )\rangle =2\,{\delta }_{D}(t-t^{\prime} )$ (the factor 2 is only a convention, and clearly it could be reabsorbed into the multiplicative term), δ_c (σ, t) is the critical threshold for collapse (possibly dependent on scale, see below), D(t) is the linear growth factor of perturbations, and σ(M) is the mass variance filtered on the mass scale M defined in Equation (3). The rationale naively followed to invent the above Equation (A1) is simple. On the left-hand side, it appears as the time derivative of an a-dimensional function of the mass that incorporates the power spectrum and the filtering scale; in choosing $\mathrm{ln}\sigma$, we have been inspired by a number of N-body simulations (e.g., Zhao et al. 2009), which suggest the mass growth of halos to be easily described in terms of such a quantity. On the rhs, it appears as a stochastic driving η(t), which, for dimensional consistency, must be multiplied by the (inverse) square root of a characteristic timescale. Since our aim here is to describe the growth of DM perturbations using quantities related to the linear regime, we find it natural to choose the timescale $| \dot{D}(t)/D(t)|$ that effectively describes the linear growth of perturbations under gravity in a given cosmological background.

We now discuss the meaning and relevance of the various terms in Equation (A1). The mass variance ${\sigma }^{2}(M)\,\propto \int {d}^{3}k\,P(k){\tilde{W}}_{M}^{2}(k)$ incorporates the initial conditions in terms of the power spectrum of perturbations P(k). We stress that in the present theory, the relation σ(M) between σ and M is purely deterministic; however, both M_t and σ(M_t ) are to be considered stochastic variables that fluctuate over cosmic time t under the influence of the noise. This is a different perspective with respect to the standard excursion set formalism, where the over-density field δ(σ) executes a random walk as a function of the mass variance σ, which in turn plays the role of a pseudo-time variable. We further note that in the excursion set approach, the choice of the filter function ${\tilde{W}}_{M}^{2}$ in the definition of σ(M) has a crucial impact, since the random trajectories δ(σ) are Markovian only when a k-sharp filter in Fourier space is adopted. In the present theory, assuming a different filter function (e.g., Gaussian or top-hat in real space) changes only the deterministic relation σ(M) but has otherwise no effect on the Markovianity (but see Paper I for a non-Markovian generalization) of the stochastic processes M_t and σ(M_t ). From a purely technical point of view, our approach has some degree of similarity to the formulation of the excursion sets in terms of stochastic differential equations and path integral representations by Maggiore & Riotto (2010a, 2010b). However, in our framework, the stochastic behavior is not associated with the trajectories in the abstract δ(σ) space but with the physical variables M_t or σ(M_t ); further key differences are also evident from the discussion below.

The fluctuating Gaussian white noise η(t) is meant to implicitly describe, at a macroscopic level, both the stochasticity in the initial conditions and all of the complex physical phenomena associated with nonlinearities in the gravitational evolution. In fact, the adopted noise is "multiplicative," in that its strength depends on the state of the system as expressed by the ratio σ/δ_c . Patches with σ ≳ δ_c tend to change their mass more abruptly, while the evolution is slower for σ ≲ δ_c ; this renders, in our stochastic theory, the naive expectation that over-dense regions in the initial perturbation field are more prone to collapse. Nevertheless, there is a crucial difference with the excursion set theory. In the latter, to avoid double counting of overlapping patches, only the first crossings of the walk δ(σ) with the barrier at δ_c must be considered; this implies that, for the relevant trajectories, the mass of proto-halos can only increase. In our approach, the fluctuating noise can induce both positive and negative variations of the mass M_t within the filtered region (with the obvious constraint M_t > 0). Positive variations of M_t can be reasonably related to mergers among collapsing proto-halos, or mass accretion from the field. Negative variations can be interpreted in terms of mass loss due to gravitational interactions, tidal forces, stripping, and fragmentation with overlapping or embedding proto-halo patches. On the other hand, the multiplicative nature of the noise induces, on statistical grounds, an average drift toward larger masses: as η(t) fluctuates, so does the random variable σ(M_t ), and hence the multiplicative factor σ/δ_c on the rhs of Equation (A1) varies; therefore, 〈σ η/δ_c 〉 is not null even if 〈η〉 is. It turns out, as shown in Paper I by the numerical integration of Equation (A1), that this noise-induced drift actually makes σ(M_t ) copy the decrease of δ_c (t) with time, so that given the inverse convex shape of the deterministic function σ(M), a net average increase in mass M_t is enforced. This is of the order of the correct amount to avoid the cloud-in-cloud issue and statistically recover, without the need to restrict to first crossing distributions, the excursion set result on halo abundance (at least when the mass-independent δ_c from spherical collapse is employed). In other words, the cloud-in-cloud issue is inherently treated and solved via the multiplicative noise structure of Equation (A1).

The collapse threshold δ_c (σ, t) is meant to represent the level of over-density required for a perturbation on given mass scale to undergo efficient collapse. In Paper I we adopted a dependence on the mass variance analogous to what is done in the excursion set

$$\begin{eqnarray}&&{\delta }_{c}(\sigma ,t)\approx \sqrt{q}{\delta }_{c}\left[1+\beta {\left(\displaystyle \frac{\sigma }{\sqrt{q}{\delta }_{c}}\right)}^{2\gamma }\right],\end{eqnarray} \tag{ A2 }$$

in terms of the spherical collapse threshold δ_c (t) ≈ 1.68 D(0)/D(t) and of three shape parameters (q, β, γ). In the excursion set approach, δ_c (σ, t) constitutes a deterministic barrier (but see Maggiore & Riotto 2010b for generalizations) that the random walk δ(σ) must hit to enforce collapse. Its mass dependence is ascribed, though a bit naively, to the fact that perturbations may be ellipsoidal in shape. In our stochastic framework, δ_c (σ, t) is just a quantity modulating the noise, and as such is meant to describe a bunch of effects not only limited to the shape of proto-halo patches, but also to other complex phenomena that can influence their collapse, like tidal torques and angular momentum, cosmological constant, dynamical friction, and possibly the specialty of collapse locations. In this vein, in Paper I we have set the values of the parameters (q, β, γ) in order to reproduce the unconditional halo mass function measured in N-body simulations.

Specifically, the probability density function ${ \mathcal P }(M,t)$ of finding the system in a state with mass M at cosmic time t can be evaluated by writing down the Fokker–Planck equation associated with the stochastic Equation (A1), with the collapse threshold given by Equation (A2); in terms of the self-similar, scaled variable ν ≡ δ_c (t)/σ(M), this can be put in the compact form

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}\,{ \mathcal P }(\nu ,t) & = & \left|\displaystyle \frac{\dot{D}}{D}\right|\,{\partial }_{\nu }\left\{\nu \,{ \mathcal P }(\nu ,t)\right.\\ & & \left.+\displaystyle \frac{1}{2{ \mathcal P }(\nu ,t)}\,{\partial }_{\nu }{\left[\displaystyle \frac{{ \mathcal P }(\nu ,t)}{\sqrt{q}\left[1+\beta /{\left(\sqrt{q}\nu \right)}^{2\gamma }\right]}\right]}^{2}\right\}.\end{array}\end{eqnarray} \tag{ A3 }$$

As discussed in Paper I, the halo mass function N(M, t) can be related to the stationary solutions $\bar{{ \mathcal P }}(\nu )$ of the above equation; the result turns out to be

$$\begin{eqnarray}&&N(M,t)=\displaystyle \frac{\bar{\rho }}{M\,\sigma }\left|\displaystyle \frac{d\sigma }{{dM}}\right|\,\nu \,\bar{{ \mathcal P }}(\nu ),\end{eqnarray} \tag{ A4 }$$

where

$$\begin{eqnarray}&&\begin{array}{l}\bar{{ \mathcal P }}(\nu )={ \mathcal A }\,\sqrt{q}\left[1+\displaystyle \frac{\beta }{{\left(\sqrt{q}\nu \right)}^{2\gamma }}\right]\exp \left\{-q\,\displaystyle \frac{{\nu }^{2}}{2}\right.\\ \quad \left.\times \,\left[1+\displaystyle \frac{2\,\beta }{1-\gamma }\,\displaystyle \frac{1}{{\left(\sqrt{q}\nu \right)}^{2\gamma }}+\displaystyle \frac{{\beta }^{2}}{1-2\,\gamma }\,\displaystyle \frac{1}{{\left(\sqrt{q}\nu \right)}^{4\gamma }}\right]\right\},\end{array}\end{eqnarray} \tag{ A5 }$$

and ${ \mathcal A }$ is a constant determined by the normalization condition ${\int }_{0}^{\infty }d\nu \,\bar{{ \mathcal P }}(\nu )=1$. When posing q = 1 and β = 0 corresponding to the mass-independent, spherical collapse threshold δ_c (t), we exactly recover the Press & Schechter (1976) mass function $\bar{{ \mathcal P }}(\nu )=\sqrt{2/\pi }\,{e}^{-{\nu }^{2}/2}$. On the other hand, comparison with the outcomes of N-body simulations (see Sheth & Tormen 1999; Bhattacharya et al. 2011; Watson et al. 2013) yields the optimal parameter values q ≈ 0.65, β ≈ 0.15, and γ ≈ 0.35. Remarkably, the asymptotic behavior for ν >> 1, corresponding to large masses and/or early cosmic times, is seen to produce a shape akin to the empirical fit of N-body simulations adopted since Sheth & Tormen (1999); however, for finite ν, the terms in the exponential are important and must be taken into account.

A.1. Halo Mass Function in Warm Dark Matter Cosmologies

We now provide a further test of our stochastic theory, concerning the halo mass function in warm dark matter (WDM) cosmologies. In fact, if our framework has captured the macroscopic essence of the proto-halo evolution and collapse, it should be able to describe the halo mass function for cosmologies different from the standard CDM, and in particular for WDM models, without any adjustment of the parameters entering the expression of the halo mass function (e.g., Hahn & Paranjape 2014). To wit, the halo mass function should be still expressed by Equations (A4) and (A5) with the same fitted parameter values (q, β, γ) given above, and the only modification with respect to CDM would be in the deterministic relation σ(M) that depends explicitly on the power spectrum of perturbations.

Specifically, we employ a WDM power spectrum (Bode et al. 2001; Viel et al. 2005) given by

$$\begin{eqnarray}&&{P}_{\mathrm{WDM}}(k)={P}_{\mathrm{CDM}}(k){\left[1+{\left(\alpha k\right)}^{2\mu }\right]}^{-10/\mu }\end{eqnarray} \tag{ A6 }$$

with P_CDM(k) the standard CDM spectrum (e.g., Bardeen et al. 1986), μ ≈ 1.12 and

$$\begin{eqnarray}&&\alpha \approx 0.049{\left(\displaystyle \frac{{m}_{\mathrm{WDM}}}{\mathrm{keV}}\right)}^{-1.11}{\left(\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{WDM}}}{0.25}\right)}^{0.11}{\left(\displaystyle \frac{h}{0.7}\right)}^{1.22}\,{h}^{-1}\,\,\mathrm{Mpc},\end{eqnarray} \tag{ A7 }$$

in terms of the WDM mass m_WDM for a thermally produced particle; WDM free-streaming effects introduce a characteristic mass scale of suppression ${M}_{\mathrm{WDM}}\,=(4\pi /3)\bar{\rho }{\left[\pi \alpha {\left({2}^{\mu /5}-1\right)}^{-1/2\mu }\right]}^{3}$ in the power spectrum. We consider WDM masses m_WDM ≈ 0.25, 0.5, and 1 keV, corresponding to suppression scales M_WDM ≈ 10¹⁰, 10¹¹, and 10¹² M_⊙. We note that such particle masses are already ruled out by observations of the Lyα forest and high-z galaxy statistics (see Viel et al. 2013; Lapi & Danese 2015), and are considered here just for the sake of testing. From the expression of the power spectrum, we compute the mass variance ${\sigma }^{2}(M)\propto \int {d}^{3}k\,{P}_{\mathrm{WDM}}(k){\tilde{W}}_{M}^{2}(k)$ by choosing a k-sharp filter in Fourier space ${\tilde{W}}_{M}(k)={\theta }_{{\rm{H}}}({k}_{M}-k)$ with standard mass assignment ${k}_{M}={(6{\pi }^{2}\bar{\rho }/M)}^{1/3}={(9\pi /2)}^{1/3}/{R}_{M}$ and ${R}_{M}\,={(3M/4\pi \bar{\rho })}^{1/3};$ note that we do not find it necessary to change the multiplicative factor in the relation ${k}_{M}\propto {R}_{M}^{-1}$, as in Benson et al. (2013) and Schneider et al. (2013). As discussed by these authors, the choice of the k-sharp filter is motivated empirically, since it originates as a mass function that gets correctly suppressed below the free-streaming length of the WDM particles. Possibly a deeper reason lays in the real-space non-locality inherent in the filter, which may somehow capture the properties of the initial density environment near small-mass WDM peaks (see Hahn & Paranjape 2014). For completeness, we will also show results for a real-space top-hat window ${\tilde{W}}_{M}{(k)=3({{kR}}_{M})}^{-3}\,[\sin ({{kR}}_{M})-{{kR}}_{M}\,\cos ({{kR}}_{M})];$ a Gaussian filter produces indistinguishable outcomes.

The predictions on the WDM mass function from Equations (A4) and (A5) at z = 0 are illustrated in Figure 6 as solid (k-sharp filter) and dashed (real top-hat or Gaussian filter) lines, for three values of the WDM particle mass m_WDM = 0.25 (blue), 0.5 (green), and 1 (red) keV. The CDM mass function is also plotted for reference (black line). The WDM mass functions are found to deviate from the CDM one around the corresponding free-streaming mass scale M_WDM, with the k-sharp filter producing the expected downturn. Our results for the k-sharp filter agree well with the outcomes from the N-body simulations by Schneider et al. (2013). Two caveats are worth mentioning on the simulated WDM mass functions: (i) the steepening, which occurs far below the free-streaming mass scale, is well known to be spuriously caused by artificial fragmentation (see also discussion by Hahn & Paranjape 2014); and (ii) the effects of late-time thermal velocities are neglected or approximately treated, so that the shape of the halo mass function around the free-streaming mass scale is still debated (see Benson et al. 2013). Given these uncertainties, we consider the result of the test conducted here to be quite satisfying.

Zoom In Zoom Out Reset image size