Press-Schechter Formalism

The trick: count noise, predict galaxies

Imagine you want to know how many galaxy clusters of mass 10¹⁵ solar masses exist in a cubic gigaparsec of the universe today, and how many of 10¹², and how many of every mass in between. The honest way is to take the measured initial conditions, evolve a billion particles forward under gravity in a supercomputer, and count the halos that condense out. The Press-Schechter formalism does it with a pencil. From two assumptions — that the initial density field is Gaussian noise, and that a patch collapses once its overdensity reaches δ_c = 1.686 — it writes down a closed-form expression for the abundance of halos at every mass. That expression is the halo mass function, and the leap that produces it is one of the most economical ideas in all of cosmology.

The intuition is almost embarrassingly simple. Take the density field of the early universe, δ(x) = δρ/ρ, and blur it with a filter that averages over a sphere containing a mass M. Now ask: what fraction of the volume sits above the collapse threshold? Those over-threshold patches are exactly the regions that will, by today, have pulled themselves together into a bound halo of at least mass M. Because the field is Gaussian, "the fraction above a threshold" is just the tail of a bell curve — an error function. Differentiate that collapsed fraction with respect to M and you have how many halos form at each mass. That is the whole machine.

How it works, ingredient by ingredient

Ingredient one — a Gaussian random field. Inflation predicts that the primordial density contrast is, to exquisite accuracy, a Gaussian random field: at any point, smoothed on any scale, δ is drawn from a normal distribution with zero mean and a variance σ²(M) you can compute. The variance comes from integrating the matter power spectrum P(k) against the filter:

σ²(M) = (1 / 2π²) ∫ P(k) Ŵ²(k R) k² dk ,   M = (4/3) π ρ̄ R³

Large regions average over more of the field, so σ(M) is a decreasing function of mass. This single curve, σ(M), carries all the cosmology — change σ₈ or the shape of P(k) and the whole prediction shifts.

Ingredient two — a collapse threshold. Spherical-collapse theory tells you when a region has gone nonlinear. Take a uniform overdense sphere; it expands, decelerates, turns around, and collapses. At the moment of collapse the real overdensity is formally infinite, but if you take the linear-theory growth of that same perturbation and extrapolate it to the collapse time, you get a finite number: δ_c = 1.686. Since the Gaussian field Press-Schechter counts is the linear field, this is the right yardstick. Any smoothed patch with linear δ > 1.686 has, in reality, already collapsed.

Putting them together. The fraction of mass that has collapsed into objects of mass ≥ M is the fraction of the field, smoothed on scale M, lying above δ_c:

F(>M) = 2 × ∫[δ_c → ∞] (1 / √(2π σ²)) exp(−δ² / 2σ²) dδ
      = 2 × (1/2) erfc( δ_c / (√2 σ(M)) )
      = erfc( δ_c / (√2 σ(M)) )

Note the factor of 2 out front — we will return to it, because that single number is the formalism's most famous wart. Differentiating −dF/dM and multiplying by ρ̄/M (to convert collapsed-mass fraction into a number density) gives the mass function:

dn/dM = √(2/π) (ρ̄ / M²) (δ_c / σ) |d ln σ / d ln M| exp(−δ_c² / 2σ²)

Everything interesting lives in the exponential. The dimensionless ratio ν = δ_c/σ(M) — "peak height" — controls the whole shape. Where σ ≫ δ_c (small M), ν is tiny, the exponential is ≈ 1, and the mass function is a near power law: low-mass halos are everywhere. Where σ drops below δ_c (cluster scale), ν exceeds 1 and the exponential slams the abundance down. That exponential cutoff is why clusters are exponentially rare and why their counts are such a sharp ruler.

Where 1.686 comes from

The threshold is not a fit — it is a derivation. In an Einstein-de Sitter universe (matter-dominated, flat), track a spherical top-hat overdensity with the exact parametric solution of the collapse:

collapse turnaround:  linear δ_ta = (3/20)(6π)^(2/3) ≈ 1.062
full collapse:        linear δ_c  = (3/20)(12π)^(2/3) ≈ 1.686
virialization:        nonlinear Δ_vir ≈ 178 (the "200×" of M_200)

The same calculation hands you three iconic numbers at once: the linear threshold 1.686 used in Press-Schechter, and the virial overdensity ~178 (rounded to 200 in the M₂₀₀ halo-mass convention). In a ΛCDM universe the threshold drifts weakly with redshift because dark energy slows growth — at z = 0 it is closer to 1.674 — but the variation is at the percent level and δ_c = 1.686 is used almost universally as the canonical value.

The factor of 2 and the excursion set

Here is the wart. Press and Schechter's counting argument tallies only the mass sitting directly above δ_c after smoothing. But a Gaussian is symmetric about zero, and in the limit of a vanishing threshold it splits the mass exactly 50/50 — so their formula assigns only half of all the mass in the universe to collapsed objects, even as the threshold goes to zero where everything should have collapsed. To force the books to balance, Press and Schechter multiplied by 2 by hand. It worked, it matched the simulations of the day, and it bothered theorists for seventeen years.

The missing half is the "cloud-in-cloud" problem. A point can be underdense (δ < δ_c) when you smooth on a small scale, yet be embedded inside a much larger region that exceeds δ_c when smoothed on a big scale. That point is genuinely inside a halo; the naive counting throws it away.

The repair came from Bond, Cole, Efstathiou & Kaiser (1991) — the excursion-set picture. Fix a point in space and watch its smoothed overdensity δ(S) as you shrink the smoothing scale, where S = σ²(M) runs upward from zero. With a sharp filter in k-space, each shell of new k-modes is independent, so δ(S) executes a Brownian random walk against the variable S. A halo of mass M is identified with the first scale at which the walk crosses the barrier δ_c. The mathematics of first passage is classical: the first-crossing distribution of a Brownian walk against a constant barrier is exactly

f(S) dS = (δ_c / √(2π S^3)) exp(−δ_c² / 2S) dS

and integrating it reproduces the Press-Schechter mass function including the factor of 2 — now derived, not imposed. The factor of 2 is the contribution of trajectories that have already crossed the barrier on a larger scale and wandered back below it; the naive argument silently discarded them. After 1991, the formalism stood on a rigorous footing.

Worked example: putting numbers in

Take a ΛCDM universe with σ₈ = 0.81, Ω_m = 0.31, and evaluate the cluster-scale abundance today. At a mass scale M = 10¹⁵ h⁻¹ M_⊙, the smoothed variance is roughly σ(M) ≈ 0.55, so the peak height is

ν = δ_c / σ = 1.686 / 0.55 ≈ 3.07
exponential suppression = exp(−ν²/2) = exp(−4.7) ≈ 9 × 10⁻³

That factor of ~0.009 is what makes 10¹⁵-solar-mass clusters rare: a comoving number density of order 10⁻⁷ per (h⁻¹ Mpc)³. Now drop to M = 10¹² h⁻¹ M_⊙, a Milky-Way-scale halo. There σ(M) ≈ 1.8, so

ν = 1.686 / 1.8 ≈ 0.94
exponential suppression = exp(−0.44) ≈ 0.64

The suppression has all but vanished — galaxy-scale halos are common, with number densities a factor ~10⁵ higher than clusters. Three orders of magnitude in mass move you from "barely suppressed power law" to "exponentially crushed tail," entirely because σ(M) crossed δ_c. And notice the sensitivity: if you nudged σ₈ up by 5%, σ(M) at the cluster scale rises to ~0.58, ν falls to 2.91, and exp(−ν²/2) jumps to ~0.014 — a 50% increase in cluster abundance from a 5% change in σ₈. This exponential leverage is exactly why cluster counts are a premier σ₈ probe.

Variants and refinements

The bare Press-Schechter formula is the trunk of a whole family tree:

• Extended Press-Schechter (EPS). The conditional first-crossing distribution — the probability that a walk crossing δ_c on scale M₁ at redshift z₁ crossed an earlier barrier on scale M₂ > M₁ at an earlier redshift — gives halo merger rates and progenitor mass functions. EPS merger trees power semi-analytic models of galaxy formation.
• Sheth-Tormen (1999). Real perturbations collapse ellipsoidally, not spherically; the effective barrier becomes mass-dependent ("moving barrier"). Sheth, Mo & Tormen fit a form with this physics that matches N-body simulations far better — fixing PS's factor-~2 over-prediction at low mass and under-prediction at high mass.
• Tinker et al. (2008). A purely empirical fitting function calibrated directly against large N-body simulations, accurate to a few percent over the cluster-cosmology mass range, with explicit dependence on the chosen overdensity (Δ = 200, 500, etc.). This is what survey pipelines actually use.
• Non-Gaussian extensions. If primordial fluctuations have a skewness f_NL ≠ 0, the tail of the distribution above δ_c is enhanced or suppressed, shifting the cluster abundance. The high-mass tail is therefore a sensitive probe of primordial non-Gaussianity.

PS versus its descendants

Formalism	Year	Collapse model	Factor of 2	N-body accuracy	Main use
Press-Schechter	1974	spherical, δ_c = 1.686	inserted by hand	tens of percent	analytic intuition
Bond et al. (excursion set)	1991	spherical, random walk	derived rigorously	same shape as PS	rigorous foundation
Extended PS (EPS)	1991–93	conditional crossings	derived	tens of percent	merger trees, SAMs
Sheth-Tormen	1999	ellipsoidal, moving barrier	built in	~10–20%	improved mass function
Jenkins et al.	2001	simulation fit	n/a	~10–20%	universal fit
Tinker et al.	2008	simulation fit, Δ-dependent	n/a	~few percent	cluster cosmology
Despali et al.	2016	simulation, virial Δ	n/a	~few percent	universal recalibration

Observational status and applications

Press-Schechter — or rather its calibrated descendants — is one of the workhorses of precision cosmology. The high-mass exponential tail maps directly onto the abundance of galaxy clusters, the most massive bound objects in the universe, and counting them as a function of mass and redshift is a clean test of structure growth:

• Sunyaev-Zel'dovich cluster counts. Planck, the South Pole Telescope (SPT), and the Atacama Cosmology Telescope (ACT) detect clusters by the SZ decrement they imprint on the CMB. Fitting their number counts to a Tinker-class mass function constrains σ₈ and Ω_m — and produced an early, mild tension with the CMB-inferred σ₈.
• X-ray surveys. eROSITA's all-sky X-ray survey detects ~10⁵ clusters; their mass function is among the strongest current handles on the growth of structure and dark energy.
• Weak-lensing-calibrated counts. DES, KiDS, and HSC weight cluster counts by lensing masses, sidestepping the messy baryonic mass-observable relations that limited earlier work.
• Reionization and the first halos. The same formalism, pushed to high redshift, predicts how many ~10⁸-solar-mass halos exist at z ~ 10–20 to host the first stars — the seeds JWST is now probing at the bright end.
• Galaxy formation. Through the extended (conditional) version, PS supplies the merger trees that semi-analytic models hang baryonic physics on, connecting the dark skeleton to observable galaxies.

Common pitfalls and misconceptions

• Thinking δ_c = 1.686 is a real, present-day overdensity. It is the linearly-extrapolated overdensity at collapse. The true nonlinear overdensity at collapse is infinite, and at virialization it is ~178–200. The 1.686 only makes sense as a threshold applied to the linear Gaussian field.
• Forgetting why the factor of 2 is there. It is not a normalization convenience — it is the cloud-in-cloud mass that the naive counting omits, recovered exactly by the random-walk first-crossing probability.
• Trusting the low-mass slope. Bare PS over-predicts small halos by roughly a factor of 2. For anything quantitative, use Sheth-Tormen or Tinker; reserve PS for intuition and scaling arguments.
• Confusing the sharp-k filter with a real-space top hat. The clean random-walk derivation needs a sharp k-space filter so increments are independent; but the mass-radius relation M ∝ R³ is usually quoted for a real-space top hat. Mixing the two introduces an order-unity ambiguity in how σ(M) maps to mass.
• Assuming universality is exact. The mass function expressed in terms of ν is nearly — but not perfectly — universal across cosmology and redshift. Percent-level work requires explicit recalibration (Tinker, Despali), and the choice of halo-finder and overdensity definition matters.
• Treating it as a derivation of structure formation. PS does not evolve anything dynamically; it is a statistical shortcut that maps linear initial conditions onto a final abundance via the collapse threshold. It says nothing about halo spatial clustering (that needs the bias / peak-background-split argument, also due to the excursion set).