Part V: Structure Formation | Chapter 1

Large-Scale Structure Formation

From tiny quantum fluctuations to the cosmic web: how gravitational instability amplifies primordial density perturbations into galaxies, clusters, filaments, and voids

1. Introduction

The cosmic microwave background reveals that the universe at recombination ($z \approx 1100$) was astonishingly smooth, with density fluctuations of order $\delta\rho/\bar{\rho} \sim 10^{-5}$. Yet today we observe a rich tapestry of structure: galaxies, galaxy clusters, filaments spanning hundreds of megaparsecs, and vast empty voids. The central question of structure formation is: how did gravity amplify these tiny seeds into the cosmic web we observe?

The Standard Picture

During inflation, quantum vacuum fluctuations in the inflaton field were stretched to macroscopic scales, imprinting a nearly scale-invariant spectrum of density perturbations with spectral index $n_s \approx 0.965$. After inflation, these perturbations evolved under gravity: dark matter perturbations grew continuously from matter-radiation equality ($z_{\rm eq} \approx 3400$), while baryon perturbations were frozen by radiation pressure until recombination ($z \approx 1100$), after which baryons fell into the pre-existing dark matter potential wells.

This hierarchical picture — small structures forming first, then merging into larger ones — is the foundation of the $\Lambda$CDM model and has been spectacularly confirmed by galaxy surveys, weak lensing measurements, and CMB observations.

The mathematical framework proceeds in stages of increasing complexity: linear perturbation theory (valid for $\delta \ll 1$), the Jeans analysis for gravitational instability, the transfer function encoding scale-dependent evolution, the matter power spectrum as the central statistical descriptor, nonlinear collapse models, and finally the Press-Schechter formalism for predicting the abundance of collapsed structures.

2. Linear Perturbation Theory

2.1 The Density Contrast

We define the density contrast (or overdensity) as the fractional deviation of the local density from the cosmic mean:

$$\delta(\mathbf{x}, t) \equiv \frac{\rho(\mathbf{x}, t) - \bar{\rho}(t)}{\bar{\rho}(t)}$$

By construction, $\langle \delta \rangle = 0$. The density field is$\rho = \bar{\rho}(1 + \delta)$. Linear theory applies when $|\delta| \ll 1$.

Similarly, we decompose the velocity field into the background Hubble flow plus a peculiar velocity:$\mathbf{v} = H\mathbf{r} + \mathbf{v}_{\rm pec}$, where $\mathbf{v}_{\rm pec}$ is the deviation from uniform expansion.

2.2 Perturbed Fluid Equations in an Expanding Universe

We treat the matter as a pressureless (or pressure-bearing) fluid in an expanding background. Working in comoving coordinates $\mathbf{x} = \mathbf{r}/a(t)$ and using the peculiar velocity $\mathbf{u} = \mathbf{v}_{\rm pec}/a$, the three fundamental equations of fluid dynamics become:

Continuity Equation (mass conservation)

$$\frac{\partial \delta}{\partial t} + \frac{1}{a}\nabla \cdot \left[(1+\delta)\mathbf{u}\right] = 0$$

Linearizing (dropping terms of order $\delta \cdot \mathbf{u}$):

$$\frac{\partial \delta}{\partial t} + \frac{1}{a}\nabla \cdot \mathbf{u} = 0$$

Euler Equation (momentum conservation)

$$\frac{\partial \mathbf{u}}{\partial t} + H\mathbf{u} + \frac{1}{a}(\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{a}\nabla\Phi - \frac{1}{a\bar{\rho}}\nabla p$$

Linearizing:

$$\frac{\partial \mathbf{u}}{\partial t} + H\mathbf{u} = -\frac{1}{a}\nabla\Phi - \frac{c_s^2}{a}\nabla\delta$$

where $c_s^2 = dp/d\rho$ is the sound speed squared, and we used $\nabla p = c_s^2 \bar{\rho}\nabla\delta$.

Poisson Equation (gravity)

$$\nabla^2\Phi = 4\pi G \bar{\rho} a^2 \delta$$

Here $\Phi$ is the peculiar gravitational potential (the perturbation to the Newtonian potential beyond the background expansion), and $\nabla^2$ is the comoving Laplacian.

2.3 The Growth Equation

We combine the three linearized equations into a single second-order ODE for $\delta$. Take $\partial/\partial t$ of the continuity equation and $\nabla \cdot$ of the Euler equation, then eliminate $\nabla \cdot \mathbf{u}$ and use Poisson's equation to replace $\nabla^2\Phi$:

The Growth Equation

$$\boxed{\ddot{\delta} + 2H\dot{\delta} - 4\pi G\bar{\rho}\,\delta = \frac{c_s^2}{a^2}\nabla^2\delta}$$

This is a damped oscillator equation with a gravitational source term. The $2H\dot{\delta}$ term is Hubble friction from the expanding background.

For a single Fourier mode $\delta_k \propto e^{i\mathbf{k}\cdot\mathbf{x}}$ with comoving wavenumber $k$:

$$\ddot{\delta}_k + 2H\dot{\delta}_k + \left(\frac{c_s^2 k^2}{a^2} - 4\pi G\bar{\rho}\right)\delta_k = 0$$

The competition between the pressure term ($c_s^2 k^2/a^2$) and gravity ($4\pi G\bar{\rho}$) determines whether a mode grows or oscillates.

2.4 Growing and Decaying Modes

For pressureless matter ($c_s = 0$) in a matter-dominated universe ($\bar{\rho} = \bar{\rho}_0 a^{-3}$, $a \propto t^{2/3}$,$H = 2/(3t)$), the growth equation becomes:

$$\ddot{\delta} + \frac{4}{3t}\dot{\delta} - \frac{2}{3t^2}\delta = 0$$

where we used $4\pi G\bar{\rho} = 2/(3t^2)$ in the Einstein–de Sitter background.

Trying a power-law ansatz $\delta \propto t^n$, we get $n(n-1) + \frac{4}{3}n - \frac{2}{3} = 0$, yielding $n = 2/3$ or $n = -1$. The general solution is:

$$\delta(t) = A\,\delta_+(t) + B\,\delta_-(t) = A\,t^{2/3} + B\,t^{-1}$$

Since $a \propto t^{2/3}$ in the matter era: the growing mode grows as $\delta_+ \propto a$, while the decaying mode$\delta_- \propto t^{-1} \propto a^{-3/2}$ rapidly becomes negligible.

2.5 The Growth Factor D(a)

In a general $\Lambda$CDM cosmology, the growing mode is no longer simply proportional to $a$. We define the linear growth factor $D(a)$as the growing solution normalized so that $D(a) \propto a$ during matter domination:

$$D(a) = \frac{5\Omega_{m,0}}{2}\,H(a)\int_0^a \frac{da'}{[a' H(a')]^3}$$

This is the Heath (1977) integral formula. In the matter era, $D \propto a$. In the $\Lambda$-dominated era, the growth factor saturates (growth freezes) because dark energy accelerates expansion faster than gravity can amplify perturbations.

Key result: The growth suppression factor $g(a) = D(a)/a$decreases from $g = 1$ during matter domination to $g \approx 0.78$ today for $\Omega_{m,0} = 0.315$. This 22% suppression is one of the primary observational signatures used to constrain dark energy through weak lensing and galaxy clustering.

3. The Jeans Instability

3.1 The Jeans Length and Jeans Mass

Returning to the Fourier-space growth equation, a mode grows when the gravitational term overcomes the pressure term, i.e., when $4\pi G\bar{\rho} > c_s^2 k^2/a^2$. The critical scale separating growing from oscillating modes is the Jeans length:

Jeans Length

$$\boxed{\lambda_J = c_s\sqrt{\frac{\pi}{G\bar{\rho}}}}$$

Perturbations with $\lambda > \lambda_J$ (i.e., $k < k_J = 2\pi/\lambda_J$) undergo gravitational collapse. Perturbations with $\lambda < \lambda_J$ oscillate as pressure-supported sound waves.

The corresponding Jeans mass is the mass enclosed in a sphere of diameter $\lambda_J$:

$$M_J = \frac{4\pi}{3}\bar{\rho}\left(\frac{\lambda_J}{2}\right)^3 = \frac{\pi}{6}\bar{\rho}\,\lambda_J^3$$

Only perturbations containing mass $M > M_J$ can collapse under self-gravity.

3.2 Baryons Before Recombination

Before recombination ($z > 1100$), baryons are tightly coupled to photons via Thomson scattering, forming a single baryon-photon fluid. The sound speed is:

$$c_s = \frac{c}{\sqrt{3(1 + R)}}, \qquad R \equiv \frac{3\bar{\rho}_b}{4\bar{\rho}_\gamma} = \frac{3\Omega_b}{4\Omega_\gamma}\,a$$

At high redshift ($R \ll 1$), the sound speed approaches the relativistic limit $c_s \approx c/\sqrt{3} \approx 1.7 \times 10^5$ km/s.

With this large sound speed, the Jeans length is comparable to the Hubble radius — essentially no baryonic perturbation can grow before recombination. Instead, perturbations oscillate as acoustic waves, producing the famous acoustic peaks in the CMB angular power spectrum.

3.3 Baryons After Recombination

At recombination, photons decouple from baryons. The baryon sound speed drops precipitously:

$$c_s = \sqrt{\frac{k_B T}{m_p}} \approx 5\;\text{km/s} \quad \text{at } T \approx 3000\;\text{K}$$

A drop by a factor of $\sim 3 \times 10^4$! The Jeans mass plummets from $\sim 10^{16}\;M_\odot$ (cluster-scale) to $\sim 10^5\;M_\odot$(globular cluster scale).

After decoupling, baryons quickly fall into the potential wells already established by dark matter, rapidly catching up to the dark matter density contrast. This is why dark matter is essential: it provides a gravitational scaffolding into which baryons can collapse.

3.4 Dark Matter: Collisionless Growth

Cold dark matter is collisionless — it has no electromagnetic interactions and essentially zero pressure ($c_s \approx 0$). Therefore, the Jeans length for dark matter is effectively zero, and dark matter perturbations can grow on all sub-horizon scales. Growth begins at matter-radiation equality ($z_{\rm eq} \approx 3400$), giving dark matter a significant head start over baryons. By recombination, dark matter overdensities have grown by a factor of $\sim a_{\rm rec}/a_{\rm eq} \approx 3$.

4. The Transfer Function

Not all Fourier modes evolve identically. Modes that enter the Hubble horizon during different cosmic epochs experience different growth histories. The transfer function $T(k)$ encodes this scale-dependent processing:

$$T(k) \equiv \frac{\delta(k, z=0)}{\delta_{\rm primordial}(k)} \cdot \frac{\delta_{\rm primordial}(k_{\rm ref})}{\delta(k_{\rm ref}, z=0)}$$

Normalized so that $T(k) \to 1$ for $k \to 0$ (super-horizon modes that entered the horizon after matter domination).

4.1 The Meszaros Effect

The key physical effect is the Meszaros effect(Meszaros 1974): during the radiation-dominated era, the expansion rate is set by radiation ($H^2 \propto \rho_r \propto a^{-4}$), but the gravitational source for matter perturbations is only the matter density ($4\pi G\bar{\rho}_m$). The ratio$\bar{\rho}_m/\bar{\rho}_r \propto a$ is small, so the gravitational term is subdominant compared to the Hubble friction. As a result, sub-horizon matter perturbations stagnate during radiation domination — they grow only logarithmically rather than as a power law.

4.2 The Characteristic Scale

The critical scale is the comoving horizon at matter-radiation equality:

$$k_{\rm eq} = a_{\rm eq}H(a_{\rm eq}) = \frac{H_0\sqrt{2\Omega_{m,0}}}{a_{\rm eq}^{1/2}} \approx 0.01\;\text{Mpc}^{-1}$$

Modes with $k \ll k_{\rm eq}$ entered the horizon during matter domination (unaffected). Modes with $k \gg k_{\rm eq}$ entered during radiation domination (suppressed).

4.3 Asymptotic Forms

Large scales (entered during matter era):

$$T(k) \to 1 \qquad \text{for } k \ll k_{\rm eq}$$

Small scales (entered during radiation era):

$$T(k) \approx \frac{A}{(k/k_{\rm eq})^2}\ln\!\left(\frac{B\,k}{k_{\rm eq}}\right) \qquad \text{for } k \gg k_{\rm eq}$$

The logarithmic factor arises from the slow (logarithmic) growth during radiation domination. The constants A and B depend on the baryon fraction.

A widely used analytic fitting formula for CDM (Bardeen, Bond, Kaiser & Szalay 1986, “BBKS”) is:

$$T(q) = \frac{\ln(1 + 2.34\,q)}{2.34\,q}\left[1 + 3.89\,q + (16.1\,q)^2 + (5.46\,q)^3 + (6.71\,q)^4\right]^{-1/4}$$

where $q = k/(\Omega_{m,0}\,h^2\;\text{Mpc}^{-1})$ and $h = H_0/(100\;\text{km/s/Mpc})$. More accurate fits (Eisenstein & Hu 1998) include baryon oscillation wiggles.

5. The Matter Power Spectrum

5.1 Definition and Physical Content

The matter power spectrum is the Fourier transform of the two-point correlation function and encodes the amplitude of density fluctuations as a function of scale:

The Matter Power Spectrum

$$\boxed{P(k, z) = A\,k^{n_s}\,T^2(k)\,D^2(z)}$$

Three ingredients: the primordial spectrum ($k^{n_s}$ from inflation), the transfer function ($T^2(k)$ from post-inflationary processing), and the growth factor ($D^2(z)$ encoding time evolution).

The power spectrum is formally defined via:

$$\langle \delta(\mathbf{k})\,\delta^*(\mathbf{k}')\rangle = (2\pi)^3\,P(k)\,\delta_D^{(3)}(\mathbf{k} - \mathbf{k}')$$

Statistical homogeneity ensures $P$ depends only on $|\mathbf{k}|$(isotropy), and the Dirac delta enforces translational invariance.

5.2 Shape of P(k)

Combining the primordial spectrum with the transfer function gives the characteristic shape:

Large scales ($k \ll k_{\rm eq}$)

$$P(k) \propto k^{n_s}$$

The primordial spectrum is preserved on scales that entered the horizon after matter-radiation equality.

Small scales ($k \gg k_{\rm eq}$)

$$P(k) \propto k^{n_s - 4}\ln^2\!\left(\frac{k}{k_{\rm eq}}\right)$$

The $k^{-4}$ suppression from $T^2 \propto k^{-4}$ (Meszaros effect) gives $P(k) \propto k^{-3}$ for $n_s \approx 1$.

The turnover occurs at $k_{\rm eq} \approx 0.01\;\text{Mpc}^{-1}$, corresponding to a comoving scale of $\sim 100\;h^{-1}$ Mpc.

5.3 The Dimensionless Power Spectrum

It is often useful to work with the dimensionless power spectrum, which gives the variance per logarithmic interval in $k$:

$$\Delta^2(k) = \frac{k^3\,P(k)}{2\pi^2}$$

The variance of the density field is $\langle\delta^2\rangle = \int \Delta^2(k)\,d\ln k$. Nonlinearity sets in when $\Delta^2(k) \sim 1$, which occurs at$k \sim 0.1 - 0.3\;h$ Mpc$^{-1}$ today.

5.4 Normalization: $\sigma_8$

The power spectrum is conventionally normalized by the parameter $\sigma_8$, the rms density fluctuation in spheres of radius $R = 8\;h^{-1}$ Mpc:

$\sigma_8$ Definition

$$\boxed{\sigma_8^2 = \int_0^\infty \frac{dk}{k}\,\frac{k^3 P(k)}{2\pi^2}\,\left|\tilde{W}(kR_8)\right|^2}$$

where $\tilde{W}(x) = 3(\sin x - x\cos x)/x^3$ is the Fourier transform of a top-hat window function, and $R_8 = 8\;h^{-1}$ Mpc.

Planck 2018 value: $\sigma_8 = 0.811 \pm 0.006$

The $S_8$ tension: Galaxy weak lensing surveys (KiDS, DES, HSC) consistently find $S_8 \equiv \sigma_8\sqrt{\Omega_m/0.3} \approx 0.76$, which is$2$–$3\sigma$ lower than the Planck CMB prediction of $S_8 \approx 0.83$. This $S_8$ tension is one of the most active areas of research in modern cosmology.

6. Nonlinear Gravitational Collapse

6.1 Spherical Top-Hat Collapse Model

When $\delta$ approaches unity, linear theory breaks down and we must solve the full nonlinear problem. The simplest analytic model is a uniform spherical overdensity (a “top-hat”) embedded in a flat Einstein–de Sitter background.

Consider a sphere of radius $r$ and uniform overdensity $\delta_i \ll 1$at some early time. By Birkhoff's theorem, the sphere evolves as a closed Friedmann universe with $k > 0$:

$$r(\theta) = A(1 - \cos\theta), \qquad t(\theta) = B(\theta - \sin\theta)$$

where $A^3 = GMB^2$, and $\theta$ is the development angle.

The sphere reaches maximum expansion at $\theta = \pi$ (turnaround), then collapses to a singularity at $\theta = 2\pi$. The key result is the linear overdensity at the moment of collapse:

Critical Overdensity for Collapse

$$\boxed{\delta_c = \frac{3}{20}(12\pi)^{2/3} \approx 1.686}$$

This is the linear-theory overdensity extrapolated to the time of actual collapse. Despite $\delta_c \approx 1.7$, the true nonlinear overdensity at collapse is infinite (singularity). In practice, the object virializes at a finite overdensity.

6.2 Virialization

In reality, the collapsing sphere does not reach a singularity — it virializes through violent relaxation. By the virial theorem ($2K + U = 0$) and energy conservation:

$$r_{\rm vir} = \frac{1}{2}\,r_{\rm ta}$$

The virial radius is half the turnaround radius.

$$\Delta_{\rm vir} = 18\pi^2 \approx 178$$

The mean density within the virial radius is $\sim 178$ times the background density. In $\Lambda$CDM, $\Delta_{\rm vir}(z=0) \approx 337$ with respect to the mean matter density. A common convention uses $\Delta = 200$ times the critical density.

6.3 The Zel'dovich Approximation

For mildly nonlinear evolution, the Zel'dovich approximation(1970) provides an elegant Lagrangian description. Particles at initial (Lagrangian) position$\mathbf{q}$ move to Eulerian position:

$$\mathbf{x}(\mathbf{q}, t) = \mathbf{q} + D(t)\,\mathbf{\Psi}(\mathbf{q})$$

where $\mathbf{\Psi}(\mathbf{q}) = -\nabla\phi(\mathbf{q})/4\pi G\bar{\rho}a^2$ is the displacement field determined by the initial potential. This is exact in 1D and correctly predicts the formation of pancakes (sheets), filaments, and nodes.

The density in the Zel'dovich approximation is given by the Jacobian:

$$\rho(\mathbf{x}, t) = \frac{\bar{\rho}}{|(1 - D\lambda_1)(1 - D\lambda_2)(1 - D\lambda_3)|}$$

where $\lambda_1 \geq \lambda_2 \geq \lambda_3$ are the eigenvalues of the deformation tensor $\partial\Psi_i/\partial q_j$. Collapse occurs first along the direction of largest $\lambda$, forming a “pancake” (Zel'dovich pancake).

6.4 N-Body Simulations

For fully nonlinear structure formation, numerical N-body simulations are indispensable. Modern simulations (Millennium, Bolshoi, IllustrisTNG, AbacusSummit) track$10^{10}$–$10^{12}$ particles in boxes of$\sim$1–3 Gpc, using tree-PM or adaptive mesh refinement algorithms.

These simulations confirm the hierarchical structure formation paradigm: small halos form first, then merge into progressively larger ones. They reveal a universal dark matter halo density profile (the NFW profile: $\rho(r) \propto r^{-1}(1 + r/r_s)^{-2}$) and provide the calibration data for the halo mass function used in semi-analytic models of galaxy formation.

7. The Cosmic Web

7.1 Morphology of Large-Scale Structure

Galaxy redshift surveys reveal that matter on scales of 1–100 Mpc is organized into a striking network called the cosmic web (Bond, Kofman & Pogosyan 1996). It consists of four distinct morphological elements:

Nodes (Clusters)

The densest regions, where filaments intersect. Galaxy clusters with masses$10^{14}$–$10^{15}\;M_\odot$, containing hundreds to thousands of galaxies and hot intracluster gas at $T \sim 10^7$–$10^8$ K.

Filaments

Elongated structures connecting clusters, typically 10–100 Mpc long. They contain roughly 40% of all matter. Galaxies within filaments tend to have their spin axes aligned with the filament direction.

Walls (Sheets)

Two-dimensional structures bounding voids, containing ~5% of all galaxies. The Sloan Great Wall, discovered in SDSS data, extends over 400 Mpc.

Voids

Vast underdense regions with $\delta \lesssim -0.8$, occupying ~80% of the volume. Typical diameters: 20–50 Mpc. Nearly empty of galaxies. The Bootes Void spans ~100 Mpc.

7.2 Baryon Acoustic Oscillations (BAO)

Before recombination, the baryon-photon fluid supports acoustic oscillations. When photons decouple, the acoustic wavefronts freeze in place, leaving a characteristic imprint at the sound horizon scale:

BAO Standard Ruler

$$\boxed{r_s = \int_0^{t_{\rm rec}} \frac{c_s(t)}{a(t)}\,dt = \int_{z_{\rm rec}}^\infty \frac{c_s(z)}{H(z)}\,dz \approx 147\;\text{Mpc}}$$

This scale appears as a bump in the galaxy correlation function at $\sim 150$ Mpc (comoving) and as oscillatory wiggles in $P(k)$. Since $r_s$ is precisely calibrated by CMB physics, BAO serve as a “standard ruler” for measuring the expansion history $H(z)$ and the angular diameter distance $d_A(z)$.

7.3 Galaxy Surveys

Survey	Era	Galaxies	Key Science
SDSS	2000–2014	$\sim 10^6$	First BAO detection (2005)
BOSS/eBOSS	2009–2019	$\sim 2 \times 10^6$	BAO at z = 0.1–2.2
DESI	2021–2026	$\sim 4 \times 10^7$	Dark energy equation of state w(z)
Euclid	2023–2029	$\sim 10^9$	Weak lensing + BAO synergy

8. The Press-Schechter Formalism

8.1 The Halo Mass Function

Press & Schechter (1974) developed a remarkably successful analytic framework for predicting the number density of collapsed dark matter halos as a function of mass and redshift. The central assumption is that a region collapses into a virialized halo when its linearly extrapolated overdensity exceeds $\delta_c = 1.686$.

The smoothed density field on mass scale M (corresponding to radius $R = (3M/4\pi\bar{\rho})^{1/3}$) has variance:

$$\sigma^2(M, z) = D^2(z)\int_0^\infty \frac{dk}{k}\,\frac{k^3 P(k)}{2\pi^2}\,|\tilde{W}(kR)|^2$$

The fraction of the volume (or mass) in regions with $\delta > \delta_c$is $F(M) = \text{erfc}(\nu/\sqrt{2})/2$ where $\nu \equiv \delta_c/\sigma(M)$is the peak height.

The differential halo mass function — the comoving number density of halos per unit mass — is obtained by differentiating:

Press-Schechter Mass Function

$$\boxed{\frac{dn}{dM} = \frac{\bar{\rho}}{M}\left|\frac{d\ln\sigma}{dM}\right|\sqrt{\frac{2}{\pi}}\;\nu\,e^{-\nu^2/2}}$$

where $\nu = \delta_c/\sigma(M, z)$. The “fudge factor” of 2 (the famous factor-of-2 problem) accounts for mass in underdense regions that eventually accretes onto overdense ones.

8.2 Extended Press-Schechter (Excursion Set Theory)

Bond et al. (1991) provided a rigorous derivation of the Press-Schechter result using excursion set theory. The smoothed density$\delta(R)$ at a given point is treated as a random walk in the variable$S = \sigma^2(R)$ as the smoothing scale R decreases (S increases). A halo of mass M forms when the walk first crosses the barrier $\delta_c$. For a sharp-$k$ filter, the walks are Markovian and the first-crossing distribution exactly reproduces the Press-Schechter formula, including the factor of 2.

8.3 Sheth-Tormen Correction

N-body simulations show that the Press-Schechter formula underpredicts the abundance of massive halos and overpredicts low-mass halos. Sheth & Tormen (1999) introduced an improved fitting function based on ellipsoidal (rather than spherical) collapse:

$$\frac{dn}{dM} = A\,\frac{\bar{\rho}}{M}\left|\frac{d\ln\sigma}{dM}\right|\sqrt{\frac{2a}{\pi}}\left[1 + \left(\frac{\sigma^2}{a\,\delta_c^2}\right)^p\right]\frac{\delta_c}{\sigma}\exp\!\left(-\frac{a\,\delta_c^2}{2\sigma^2}\right)$$

with $a = 0.707$, $p = 0.3$, and $A \approx 0.322$(fixed by requiring all mass is in halos). This gives 10–20% agreement with N-body simulations across a wide range of masses and redshifts.

9. Computational Exercise: Matter Power Spectrum and Growth Factor

The following Python code computes the linear matter power spectrum $P(k)$ using the BBKS transfer function, evaluates the growth factor $D(a)$ for$\Lambda$CDM, and computes $\sigma_8$. It produces publication-quality plots showing how dark energy suppresses structure growth.

Python

structure_formation.py154 lines

#!/usr/bin/env python3
"""
structure_formation.py
Compute the linear matter power spectrum P(k), growth factor D(a),
and sigma_8 for LCDM cosmology.
"""
import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt

# ============================
# Cosmological Parameters (Planck 2018)
# ============================
H0 = 67.4              # km/s/Mpc
h = H0 / 100.0         # Dimensionless Hubble parameter
Omega_m = 0.315         # Total matter
Omega_b = 0.0493        # Baryons
Omega_L = 0.685         # Dark energy
Omega_r = 9.1e-5        # Radiation
n_s = 0.965             # Spectral index
sigma8_target = 0.811   # Planck normalization

# Derived quantities
Omega_c = Omega_m - Omega_b  # CDM density
a_eq = Omega_r / Omega_m     # Scale factor at equality
k_eq = 0.0732 * Omega_m * h**2  # Mpc^{-1}, horizon at equality

# ============================
# 1. Transfer Function (BBKS fitting formula)
# ============================
def transfer_bbks(k):
    """
    BBKS transfer function for CDM.
    k in units of h/Mpc.
    """
    Gamma = Omega_m * h  # Shape parameter (no baryon correction)
    q = k / Gamma        # Dimensionless wavenumber
    T = (np.log(1.0 + 2.34*q) / (2.34*q) *
         (1.0 + 3.89*q + (16.1*q)**2 +
          (5.46*q)**3 + (6.71*q)**4)**(-0.25))
    return T

# ============================
# 2. Dimensionless Hubble Parameter
# ============================
def E(a):
    """E(a) = H(a)/H0"""
    return np.sqrt(Omega_r/a**4 + Omega_m/a**3 + Omega_L)

# ============================
# 3. Linear Growth Factor D(a) via Heath integral
# ============================
def growth_factor(a):
    """
    D(a) = (5/2) Omega_m H(a) int_0^a da'/(a'H(a'))^3
    Normalized so D(1) gives present-day value.
    """
    def integrand(ap):
        return 1.0 / (ap * E(ap))**3
    integral, _ = quad(integrand, 1e-8, a)
    return 2.5 * Omega_m * E(a) * integral

# Normalize D so that D(a) -> a in matter era
D0 = growth_factor(1.0)

def D_normalized(a):
    """Growth factor normalized to D(1) = 1 for convenience."""
    return growth_factor(a) / D0

# ============================
# 4. Matter Power Spectrum P(k)
# ============================
def power_spectrum_unnorm(k):
    """Unnormalized P(k) = k^{n_s} T^2(k), k in h/Mpc."""
    return k**n_s * transfer_bbks(k)**2

# ============================
# 5. Compute sigma_8 for normalization
# ============================
def window_tophat(x):
    """Fourier transform of top-hat window."""
    return np.where(x < 1e-4, 1.0,
                    3.0*(np.sin(x) - x*np.cos(x)) / x**3)

R8 = 8.0  # h^{-1} Mpc

def sigma_sq_integrand(lnk):
    """Integrand for sigma^2(R8) in d(ln k)."""
    k = np.exp(lnk)
    W = window_tophat(k * R8)
    return k**3 * power_spectrum_unnorm(k) * W**2 / (2.0*np.pi**2)

sigma8_sq_unnorm, _ = quad(sigma_sq_integrand, np.log(1e-4),
                           np.log(1e2), limit=200)
A_s = sigma8_target**2 / sigma8_sq_unnorm  # Normalization constant

def P_k(k):
    """Normalized matter power spectrum, k in h/Mpc."""
    return A_s * power_spectrum_unnorm(k)

# Verify sigma_8
def sigma_R(R):
    """RMS density fluctuation in top-hat of radius R (h^{-1} Mpc)."""
    def integrand(lnk):
        k = np.exp(lnk)
        W = window_tophat(k * R)
        return k**3 * P_k(k) * W**2 / (2.0*np.pi**2)
    result, _ = quad(integrand, np.log(1e-4), np.log(1e2), limit=200)
    return np.sqrt(result)

print(f"sigma_8 = {sigma_R(8.0):.4f}  (target: {sigma8_target})")

# ============================
# 6. Plot Results
# ============================
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# Panel 1: Matter Power Spectrum P(k)
k_arr = np.logspace(-4, 1.5, 500)  # h/Mpc
P_arr = np.array([P_k(k) for k in k_arr])
axes[0].loglog(k_arr, P_arr, 'purple', lw=2)
axes[0].axvline(k_eq, color='gray', ls='--', alpha=0.5,
                label=r'$k_{\rm eq}$')
axes[0].set_xlabel(r'$k\;[h\,{\rm Mpc}^{-1}]$', fontsize=13)
axes[0].set_ylabel(r'$P(k)\;[h^{-3}\,{\rm Mpc}^3]$', fontsize=13)
axes[0].set_title('Matter Power Spectrum', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].set_xlim(1e-4, 30)

# Panel 2: Growth Factor D(a)
a_arr = np.linspace(0.01, 1.0, 200)
D_arr = np.array([D_normalized(a) for a in a_arr])
axes[1].plot(a_arr, D_arr, 'purple', lw=2, label=r'$D(a)/D(1)$')
axes[1].plot(a_arr, a_arr/a_arr[-1], 'gray', ls='--', lw=1.5,
             label=r'$\propto a$ (EdS)')
axes[1].set_xlabel(r'$a$', fontsize=13)
axes[1].set_ylabel(r'$D(a)/D(1)$', fontsize=13)
axes[1].set_title(r'Growth Factor ($\Lambda$CDM vs EdS)',
                  fontsize=14)
axes[1].legend(fontsize=11)

# Panel 3: Growth suppression g(a) = D(a)/a
g_arr = D_arr / (a_arr / a_arr[-1])
axes[2].plot(a_arr, g_arr, 'magenta', lw=2)
axes[2].set_xlabel(r'$a$', fontsize=13)
axes[2].set_ylabel(r'$g(a) = D(a)/a$', fontsize=13)
axes[2].set_title('Growth Suppression by Dark Energy', fontsize=14)
axes[2].set_ylim(0.7, 1.05)
axes[2].axhline(1.0, color='gray', ls='--', alpha=0.5)

plt.tight_layout()
plt.savefig('structure_formation.png', dpi=150, bbox_inches='tight')
plt.show()
print("Plot saved to structure_formation.png")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Expected Output

The script prints $\sigma_8 = 0.8110$ (matching the target), and produces three panels:

Left: $P(k)$ showing the $k^{n_s}$ rise, turnover at $k_{\rm eq}$, and$k^{n_s-4}$ decline
Center: $D(a)$ compared to the Einstein–de Sitter prediction $D \propto a$, showing suppression at late times
Right: The suppression factor $g(a) = D(a)/a$, falling from 1.0 to ~0.78 as $\Lambda$ dominates

10. Summary

Key Equations of Structure Formation

Growth Equation

$$\ddot{\delta} + 2H\dot{\delta} - 4\pi G\bar{\rho}\,\delta = \frac{c_s^2}{a^2}\nabla^2\delta$$

Jeans Length and Jeans Mass

$$\lambda_J = c_s\sqrt{\frac{\pi}{G\bar{\rho}}}, \qquad M_J = \frac{\pi}{6}\bar{\rho}\,\lambda_J^3$$

Matter Power Spectrum

$$P(k, z) = A\,k^{n_s}\,T^2(k)\,D^2(z), \qquad \Delta^2(k) = \frac{k^3 P(k)}{2\pi^2}$$

Normalization

$$\sigma_8^2 = \int_0^\infty \frac{dk}{k}\,\Delta^2(k)\,|\tilde{W}(kR_8)|^2, \qquad \sigma_8 = 0.811$$

Spherical Collapse Threshold

$$\delta_c = \frac{3}{20}(12\pi)^{2/3} \approx 1.686, \qquad \Delta_{\rm vir} = 18\pi^2 \approx 178$$

Press-Schechter Mass Function

$$\frac{dn}{dM} = \frac{\bar{\rho}}{M}\left|\frac{d\ln\sigma}{dM}\right|\sqrt{\frac{2}{\pi}}\;\nu\,e^{-\nu^2/2}, \qquad \nu = \frac{\delta_c}{\sigma(M)}$$

BAO Standard Ruler

$$r_s = \int_{z_{\rm rec}}^\infty \frac{c_s(z)}{H(z)}\,dz \approx 147\;\text{Mpc}$$

The Physical Picture

1. Inflation generates nearly scale-invariant perturbations ($n_s \approx 0.965$), seeding $\delta\rho/\bar{\rho} \sim 10^{-5}$ at all scales.

2. During radiation domination, sub-horizon matter perturbations stagnate (Meszaros effect), while super-horizon modes are frozen.

3. After matter-radiation equality ($z_{\rm eq} \approx 3400$), dark matter perturbations grow as $\delta \propto a$ on all scales.

4. Baryons, trapped in acoustic oscillations before recombination, decouple at $z \approx 1100$ and fall into dark matter potential wells. The frozen acoustic scale becomes the BAO feature at 147 Mpc.

5. The transfer function imprints a characteristic turnover in $P(k)$ at $k_{\rm eq} \sim 0.01$ Mpc$^{-1}$, encoding the horizon at equality.

6. When $\delta \gtrsim 1$, perturbations go nonlinear. Spherical collapse gives $\delta_c = 1.686$; the Zel'dovich approximation explains the anisotropic formation of the cosmic web.

7. The Press-Schechter formalism predicts the halo mass function, connecting the linear power spectrum to the abundance of galaxies and clusters.

8. Dark energy ($\Omega_\Lambda \approx 0.685$) suppresses structure growth at late times, reducing $D(a)$ by ~22% compared to an Einstein–de Sitter universe.

Bibliography

Textbooks & Monographs

Dodelson, S. & Schmidt, F. (2020). Modern Cosmology, 2nd ed. Academic Press. — Standard reference for linear perturbation theory, the matter power spectrum, and BAO physics.
Peebles, P.J.E. (1980). The Large-Scale Structure of the Universe. Princeton University Press. — Foundational treatment of galaxy clustering, correlation functions, and structure formation.
Padmanabhan, T. (1993). Structure Formation in the Universe. Cambridge University Press. — Detailed coverage of gravitational instability, Jeans analysis, and nonlinear collapse.
Mo, H., van den Bosch, F. & White, S. (2010). Galaxy Formation and Evolution. Cambridge University Press. — Comprehensive treatment linking cosmological perturbations to galaxy formation.
Weinberg, S. (2008). Cosmology. Oxford University Press. — Rigorous derivation of perturbation growth, transfer functions, and the power spectrum.
Bernardeau, F., Colombi, S., Gaztañaga, E. & Scoccimarro, R. (2002). “Large-scale structure of the Universe and cosmological perturbation theory,” Physics Reports 367, 1–248. arXiv:astro-ph/0112551. — Comprehensive review of perturbation theory beyond linear order.

Key Papers

Jeans, J.H. (1902). “The Stability of a Spherical Nebula,” Philosophical Transactions of the Royal Society A 199, 1–53. — The original gravitational instability analysis defining the Jeans length.
Press, W.H. & Schechter, P. (1974). “Formation of Galaxies and Clusters of Galaxies by Self-Similar Gravitational Condensation,” Astrophysical Journal 187, 425–438. — The Press-Schechter formalism for the halo mass function.
Eisenstein, D.J. & Hu, W. (1998). “Baryonic Features in the Matter Transfer Function,” Astrophysical Journal 496, 605–614. arXiv:astro-ph/9709112. — Fitting formulae for the matter transfer function including baryon acoustic oscillations.
Eisenstein, D.J. et al. (2005). “Detection of the Baryon Acoustic Peak in the Large-Scale Correlation Function of SDSS Luminous Red Galaxies,” Astrophysical Journal 633, 560–574. arXiv:astro-ph/0501171. — First detection of BAO in galaxy surveys.
Navarro, J.F., Frenk, C.S. & White, S.D.M. (1997). “A Universal Density Profile from Hierarchical Clustering,” Astrophysical Journal 490, 493–508. arXiv:astro-ph/9611107. — The NFW profile for dark matter halos.
Springel, V. et al. (2005). “Simulations of the formation, evolution and clustering of galaxies and quasars,” Nature 435, 629–636. arXiv:astro-ph/0504097. — The Millennium Simulation of cosmic structure formation.
Planck Collaboration (2020). “Planck 2018 results. VI. Cosmological parameters,” Astronomy & Astrophysics 641, A6. arXiv:1807.06209. — Constraints on the matter power spectrum parameters.

← CMB Physics Dark Matter →

Survey	Era	Galaxies	Key Science
SDSS	2000–2014	\(\sim 10^6\)	First BAO detection (2005)
BOSS/eBOSS	2009–2019	\(\sim 2 \times 10^6\)	BAO at z = 0.1–2.2
DESI	2021–2026	\(\sim 4 \times 10^7\)	Dark energy equation of state w(z)
Euclid	2023–2029	\(\sim 10^9\)	Weak lensing + BAO synergy

Cosmology

Large-Scale Structure Formation

1. Introduction

The Standard Picture

2. Linear Perturbation Theory

2.1 The Density Contrast

2.2 Perturbed Fluid Equations in an Expanding Universe

Continuity Equation (mass conservation)

Euler Equation (momentum conservation)

Poisson Equation (gravity)

2.3 The Growth Equation

The Growth Equation

2.4 Growing and Decaying Modes

2.5 The Growth Factor D(a)

3. The Jeans Instability

3.1 The Jeans Length and Jeans Mass

Jeans Length

3.2 Baryons Before Recombination

3.3 Baryons After Recombination

3.4 Dark Matter: Collisionless Growth

4. The Transfer Function

4.1 The Meszaros Effect

4.2 The Characteristic Scale

4.3 Asymptotic Forms

5. The Matter Power Spectrum

5.1 Definition and Physical Content

The Matter Power Spectrum

5.2 Shape of P(k)

Large scales (\(k \ll k_{\rm eq}\))

Small scales (\(k \gg k_{\rm eq}\))

5.3 The Dimensionless Power Spectrum

5.4 Normalization: \(\sigma_8\)

\(\sigma_8\) Definition

6. Nonlinear Gravitational Collapse

6.1 Spherical Top-Hat Collapse Model

Critical Overdensity for Collapse

6.2 Virialization

6.3 The Zel'dovich Approximation

6.4 N-Body Simulations

7. The Cosmic Web

7.1 Morphology of Large-Scale Structure

Nodes (Clusters)

Filaments

Walls (Sheets)

Voids

7.2 Baryon Acoustic Oscillations (BAO)

BAO Standard Ruler

7.3 Galaxy Surveys

8. The Press-Schechter Formalism

8.1 The Halo Mass Function

Press-Schechter Mass Function

8.2 Extended Press-Schechter (Excursion Set Theory)

8.3 Sheth-Tormen Correction

9. Computational Exercise: Matter Power Spectrum and Growth Factor

Expected Output

10. Summary

Key Equations of Structure Formation

Growth Equation

Jeans Length and Jeans Mass

Matter Power Spectrum

Normalization

Spherical Collapse Threshold

Press-Schechter Mass Function

BAO Standard Ruler

The Physical Picture

Bibliography