Halo Formation and Mass Functions

A halo mass function turns a statistical field into an abundance: it predicts the comoving number of collapsed dark-matter halos in each mass interval. The answer is conditional on the smoothing window, collapse barrier, halo finder, mass definition, cosmology, and calibration range. This page derives the Press–Schechter baseline, explains how excursion sets address its counting problem, and verifies the normalization with a deterministic quadrature. It does not model baryonic selection or convert halo mass into an observed galaxy or cluster proxy.

Helpful background. The Einstein Equation supplies Einstein Equations and their weak-field limit; Relativistic Matter and Compact Objects supplies Relativistic Matter and Compact Objects kinetic and fluid closures; Relativistic Cosmology supplies the Relativistic Cosmology background; Computational Relativity supplies convergence and benchmark practice; and Probability and Statistical Inference supplies Gaussian fields, probability densities, and likelihood reasoning.

From a smoothed field to a mass scale

Let $\delta_L(\mathbf x,z)$ be the linearly extrapolated matter contrast at the epoch of interest. A spherical real-space top-hat of comoving Lagrangian radius $R$ has Fourier transform

$$W_{\rm TH}(x)=3\frac{\sin x-x\cos x}{x^3}, \qquad M=\frac{4\pi}{3}\bar\rho_{m0}R^3,$$

where $\bar\rho_{m0}$ is the present comoving mean matter density. Its variance is

$$\sigma^2(M,z)= \int_0^\infty d\ln k\, \Delta_m^2(k,z)\lvert W_{\rm TH}(kR)\rvert^2.$$

The apparent singularity at $x=0$ is removable: $W_{\rm TH}(x)=1-x^2/10+O(x^4)$. A numerical implementation should use this limit for small $x$ rather than subtract nearly equal sine and cosine terms.

The mass assignment above belongs to the real-space top-hat. A sharp cutoff in Fourier space produces a simpler Markov random walk in the variance $S=\sigma^2$, but it has no unique enclosed real-space volume. Mixing the sharp-$k$ walk with the top-hat mass relation without declaring a calibration is therefore a convention change, not an innocent algebraic shortcut.

For a constant spherical-collapse barrier $\delta_c(z)$, define the peak height

$$\nu(M,z)=\frac{\delta_c(z)}{\sigma(M,z)}.$$

The numerator and denominator must be extrapolated to the same epoch. Here $\sigma(M,z)$ is evaluated at collapse redshift $z$, so $\delta_c(z)$ is the threshold extrapolated to that redshift and is nearly constant in the usual matter-dominated calculation. An alternative convention uses the present-day $\sigma(M,0)$ and the threshold $\delta_c(z)/D(z)$ with $D(0)=1$. Applying the growth factor to both conventions at once would count it twice.

Large $\nu$ labels rare massive peaks. In Einstein–de Sitter spherical collapse $\delta_c\simeq1.686$, but the preceding page shows why this is an idealized reference rather than a universal halo-finder threshold. Desjacques, Jeong, and Schmidt 2018, §5.2.1, printed pp. 108–110, eqs. (5.1)–(5.9) (Open PDF) derive the collapse mapping and state its dependence on the halo definition.

The Press–Schechter abundance

If the smoothed linear field is Gaussian, the one-sided fraction above the barrier is

$$P(\delta_R\gt\delta_c) =\frac12\operatorname{erfc}\!\left(\frac{\nu}{\sqrt2}\right).$$

In the idealized limit $\sigma(M)\to\infty$ as $M\to0$, this tail accounts for only half the mass. Press and Schechter doubled it, interpreting initially underdense material as matter that is eventually incorporated into larger collapsed regions. Their original derivation and the explicit doubling appear in Press and Schechter 1974, printed pp. 429–430, eqs. (11)–(24) (NASA ADS PDF).

It is safest to declare the multiplicity convention before writing a mass function. Here $f(\nu)d\ln\nu$ is the fraction of matter in the interval $d\ln\nu$, with

$$f_{\rm PS}(\nu) =\sqrt{\frac{2}{\pi}}\,\nu e^{-\nu^2/2}, \qquad \int_0^\infty f_{\rm PS}(\nu)d\ln\nu=1.$$

This mathematical normalization ranges over all $0\lt\nu\lt\infty$. A spectrum with a physical small-scale cutoff can have bounded variance and thus a positive minimum $\nu$; integrating its actual mass range need not place all matter in halos. That physical limitation is distinct from omitting the factor of two in the idealized construction.

The comoving number density per logarithmic mass is then

$$\frac{dn}{d\ln M} =\frac{\bar\rho_{m0}}{M} f(\nu)\left\lvert\frac{d\ln\nu}{d\ln M}\right\rvert.$$

This formula exposes every conversion: $\sigma(M)$ supplies the spectrum and window, $\delta_c$ supplies the barrier, and $\bar\rho_{m0}/M$ converts mass fraction into number density. Different papers sometimes attach the factor $\nu$ to a differently named multiplicity; compare the integration measure, not just the symbol $f$.

Why first crossings matter

Thresholding each smoothing scale independently double-counts a small region that lies inside a larger collapsed region—the cloud-in-cloud problem. In the excursion-set construction, $\delta(S)$ is followed as the smoothing variance increases, and a mass element is assigned to the first barrier crossing. For a constant barrier and a sharp-$k$ filter, independent increments make the walk Markovian and the first-crossing distribution recovers the normalized Press–Schechter form without an ad hoc doubling.

That clean result has sharp boundaries. A real-space top-hat correlates steps; ellipsoidal collapse motivates a moving, stochastic barrier; primordial non-Gaussianity changes the walk statistics; and environmental tides complicate the mapping from first crossing to a particular halo finder. Desjacques, Jeong, and Schmidt develop the original tail, factor-of-two failure, first-crossing resolution, and correlated-walk extensions in Desjacques, Jeong, and Schmidt 2018, §§5.2–5.7, printed pp. 108–129, eqs. (5.10)–(5.61) (Open PDF). Excursion sets organize these choices; they do not remove them.

A normalization and boundary benchmark

The benchmark first integrates

$$\int_0^\infty \sqrt{\frac{2}{\pi}}e^{-\nu^2/2}d\nu=1$$

after mapping $\nu=t/(1-t)$ to $0\le t\le1$. Composite Simpson refinement gives:

Grid intervals	Mass fraction with factor two	Absolute error	Without factor two	Observed order
16	$1.000534397997$	$5.343980\times10^{-4}$	$0.500267198998$	—
32	$0.999997584824$	$2.415176\times10^{-6}$	$0.499998792412$	$7.7896$
64	$0.999999997236$	$2.764010\times10^{-9}$	$0.499999998618$	$9.7711$
128	$0.999999999802$	$1.982032\times10^{-10}$	$0.499999999901$	$3.8017$
256	$0.999999999988$	$1.238598\times10^{-11}$	$0.499999999994$	$4.0002$
512	$0.999999999999$	$7.726042\times10^{-13}$	$0.500000000000$	$4.0028$

The pre-asymptotic orders are not a claimed superconvergence rate: endpoint mapping and cancellation make them nonuniform. The last refinements recover the fourth-order Simpson regime.

A separate mass-boundary test declares a scale-free $n=0$ variance with a real-space top-hat,

$$\sigma(M)=\delta_c(M/M_*)^{-1/2}, \qquad \delta_c=1.68647, \qquad M_*=10^{12}M_\odot.$$

The Simpson step in $\ln M$ is held fixed while both logarithmic boundaries are extended:

$\log_{10}(M/M_*)$ range	Intervals	With factor two	Absolute error	Without factor two
$[-4,4]$	512	$0.992021287366$	$7.978713\times10^{-3}$	$0.496010643683$
$[-8,8]$	1024	$0.999920211544$	$7.978846\times10^{-5}$	$0.499960105772$
$[-12,12]$	1536	$0.999999202115$	$7.978846\times10^{-7}$	$0.499999601058$
$[-16,16]$	2048	$0.999999992021$	$7.978845\times10^{-9}$	$0.499999996011$

This deterministic run used Node.js 24.18.0 on Darwin arm64, IEEE-754 binary64 arithmetic, and the JavaScript standard library only. Download the MIT-licensed benchmark, then use the standalone download-and-run command sequence:

curl -fL https://blackhole.org/code/cosmology/cosmology-benchmarks.mjs -o cosmology-benchmarks.mjs
node cosmology-benchmarks.mjs --case halo-mass-function

Changing $\delta_c$ or $M_*$ relabels the toy mass axis but cannot repair a missing normalization factor. Conversely, a normalized toy integral does not validate a physical mass function.

Simulation calibration is part of the definition

Analytic Press–Schechter theory captures a low-mass power-law region and a rare-peak cutoff, but precision abundances are calibrated on simulations. A fit is usable only together with its simulation cosmologies, volume and mass resolution, initial conditions, gravity solver, halo finder, overdensity definition, redshift outputs, fitting range, and covariance or residual model.

As a concrete historical example, Tinker et al. combined HOT, ART, and GADGET2 collisionless simulations near WMAP1 and WMAP3 cosmologies; identified spherical-overdensity halos relative to the mean density; and calibrated several overdensity definitions over roughly $10^{11}$–$10^{15}\,h^{-1}M_\odot$ and $0\le z\lesssim2$. Their results showed percent-level dependence on mass definition and non-universal redshift evolution. The simulation set, halo-finder and resolution tests, fitting functions, and evolution limits are documented in Tinker et al. 2008, §2, printed pp. 2–6; §3, pp. 7–12; §4, pp. 13–14 (Open PDF). Those fits are an example of a complete calibration contract, not an evergreen prescription for arbitrary cosmology, mass proxy, baryonic treatment, redshift, or gravity theory.

Observed cluster or galaxy counts add another layer: selection, projection, mass–observable scatter, calibration, and survey volume. They must be modeled in the likelihood rather than absorbed silently into $f(\nu)$. The next page develops the related bias and redshift-space mapping for tracers.

Exercises

For the toy model $\sigma(M)=\delta_c(M/M_*)^{-1/2}$, derive the mass-fraction integrand in $d\ln M$ and explain why omitting the Press–Schechter factor of two gives one half over infinite boundaries.

Solution

Here $\nu=(M/M_*)^{1/2}$, so

$$\frac{d\ln\nu}{d\ln M}=\frac12.$$

Using the declared multiplicity,

$$\frac{dF}{d\ln M} =f_{\rm PS}(\nu)\frac{d\ln\nu}{d\ln M} =\frac12\sqrt{\frac{2}{\pi}}\, \nu e^{-\nu^2/2}.$$

Because $d\nu=\nu\,d\ln\nu=(\nu/2)d\ln M$,

$$\int_{-\infty}^{\infty}\frac{dF}{d\ln M}d\ln M =\int_0^\infty\sqrt{\frac{2}{\pi}}e^{-\nu^2/2}d\nu =1.$$

Before the doubling, the integrand is smaller by exactly one half, so its integral is $1/2$. This is the normalization defect; changing the threshold or mass scale only changes variables and cannot remove it.

References

• Desjacques, V., D. Jeong, and F. Schmidt, “Large-Scale Galaxy Bias,” Physics Reports 733, 1–193 (2018), DOI, Open PDF.
• Press, W. H., and P. Schechter, “Formation of Galaxies and Clusters of Galaxies by Self-Similar Gravitational Condensation,” Astrophysical Journal 187, 425–438 (1974), DOI, NASA ADS PDF.
• Tinker, J. L., A. V. Kravtsov, A. Klypin, K. Abazajian, M. S. Warren, G. Yepes, S. Gottlöber, and D. E. Holz, “Toward a Halo Mass Function for Precision Cosmology: The Limits of Universality,” Astrophysical Journal 688, 709–728 (2008), DOI, Open PDF.