Galaxy Correlation Function and Clustering Statistics

Look up at a wide-field galaxy map and the eye immediately reads structure: galaxies knot into groups, thread into filaments, sheet into walls, and abandon vast rounded voids. That pattern is neither the smoothness of a uniform gas nor the featureless scatter of a dice roll — it is clustering, and the single number that captures it best is the two-point correlation function ξ(r)\xi(r), the excess probability, over a random Poisson distribution, of finding a second galaxy a distance rr from a first. Peebles carried this statistic from turbulence theory into cosmology in the 1970s and made it the workhorse of large-scale structure. This article is the flagship synthesis of the galaxy-clustering head: it takes the density contrast δ(x)\delta(\mathbf{x}) from the perturbation articles, the shape of P(k)P(k) from the transfer-function and primordial-spectrum articles, and growth from linear-growth theory, and asks what the resulting statistic is, how it is measured, and what it tells us about the Universe. Two facts organize everything that follows. First, ξ(r)\xi(r) is the autocorrelation of the density field, so it is the Fourier partner of the matter power spectrum P(k)P(k) — real-space clustering and kk-space power are the same information in two languages. Second, because the primordial field is nearly Gaussian, this one two-point statistic captures essentially all of the field's information on large scales. From these flow the small-scale power law, the linear tracing of dark matter, and the acoustic bump near 100 h−1 Mpc100\,h^{-1}\,\mathrm{Mpc} that turns a survey into a cosmic ruler.

Key takeaways
  • The two-point correlation function ξ(r)\xi(r) measures the excess probability, over a random (Poisson) distribution, of finding a pair of galaxies separated by a distance rr; equivalently it is the autocorrelation ⟨δ(x)δ(x+r)⟩\langle\delta(\mathbf{x})\delta(\mathbf{x}+\mathbf{r})\rangle of the density contrast.
  • ξ(r)\xi(r) and the matter power spectrum P(k)P(k) are a Fourier-transform pair, so clustering in real space and clustering in kk-space carry identical statistical content.
  • On scales 0.10.1–10 h−1 Mpc10\,h^{-1}\,\mathrm{Mpc} the measured function is a power law ξ(r)=(r/r0)−γ\xi(r)=(r/r_0)^{-\gamma} with r0≈5 h−1 Mpcr_0\approx 5\,h^{-1}\,\mathrm{Mpc} and γ≈1.8\gamma\approx 1.8; a scale-free P(k)∝knP(k)\propto k^{n} gives exactly γ=n+3\gamma=n+3.
  • Galaxies are biased tracers: on large linear scales ξgal(r)=b2 ξm(r)\xi_{\mathrm{gal}}(r)=b^{2}\,\xi_{\mathrm{m}}(r), where bb is the (tracer- and scale-dependent) galaxy bias.
  • A faint localized bump near r≈100 h−1 Mpcr\approx 100\,h^{-1}\,\mathrm{Mpc} is the baryon acoustic oscillation feature, a standard ruler fixed by the sound horizon at the drag epoch.
  • Estimators such as Landy–Szalay, ξ^=(DD−2DR+RR)/RR\hat\xi=(DD-2DR+RR)/RR, convert galaxy and random catalogs into ξ(r)\xi(r); redshift-space distortions and survey geometry must be modeled with care.
  • Key idea: ξ(r)\xi(r) is the excess pair probability over random and equals the Fourier transform of P(k)P(k); galaxies trace the mass with bias bb, so on linear scales ξgal=b2ξm\xi_{\mathrm{gal}}=b^{2}\xi_{\mathrm{m}}, while nonlinear collapse produces the observed power law with a BAO ruler at ∼100 h−1 \sim100\,h^{-1}\,Mpc.

Where this lives today

The correlation function is the daily currency of the great galaxy surveys. BOSS and its successor eBOSS mapped millions of galaxies and quasars to pin the BAO bump in ξ(r)\xi(r) across redshift, and the Dark Energy Spectroscopic Instrument (DESI), now returning tens of millions of redshifts, has measured the ruler to sub-percent precision — one of the tightest constraints on the expansion history and hence on dark energy. The anisotropy of the redshift-space ξ(s,μ)\xi(s,\mu), decomposed into Legendre multipoles, delivers the growth rate fσ8f\sigma_8 and tests general relativity against modified gravity. For faint populations where spectroscopy is infeasible, the projected angular correlation w(θ)w(\theta) — linked to ξ(r)\xi(r) through Limber's equation — extends clustering to imaging surveys such as the Dark Energy Survey and, soon, Rubin's LSST. Open problems remain sharp: the scale dependence and stochasticity of galaxy bias, nonlinear and redshift-space modeling at percent precision, and the tension between the clustering amplitude σ8\sigma_8 from low-redshift structure and from the cosmic microwave background — a discrepancy that ξ(r)\xi(r) measurements are central to resolving.

The physical picture

Stand at any galaxy and count neighbours in a thin shell of radius rr and thickness dr\mathrm{d}r. Under an unbiased Poisson process the expected count would be the mean density times the shell volume, nˉ 4πr2 dr\bar n\,4\pi r^2\,\mathrm{d}r, with no memory of the central galaxy. Real galaxies cluster, so the count runs high at small rr; the fractional excess, averaged over all galaxies as centres, is exactly ξ(r)\xi(r). A value ξ=0.5\xi=0.5 means fifty percent more pairs than random; ξ=0\xi=0 carries no clustering; and ξ<0\xi<0 — perfectly allowed — means anticorrelation, the signature of the underdense voids that balance the overdense knots. The first figure makes this concrete by contrasting a random field with a clustered one.

:::svg[A random Poisson field versus a clustered field: a shell at separation r around a galaxy holds more neighbours than the random expectation n̄ dV, and that fractional excess is ξ(r), with joint pair probability dP₁₂ = n̄²[1 + ξ(r)] dV₁ dV₂] ξ(r) is the excess of pairs over a random Poisson field: dP₁₂ = n̄²[1 + ξ(r)] dV₁ dV₂ random N = n̄ dV, ξ = 0 clustered N = n̄ dV (1 + ξ), ξ positive shell r Number conservation forces ∫ ξ(r) 4π r² dr = 0: close-pair excess is repaid by large-scale anticorrelation.

::: Think of ξ(r)\xi(r) and P(k)P(k) as two dialects for one story. The correlation function answers "how much does a galaxy at separation rr know about its neighbour?"; the power spectrum answers "how much variance lives in fluctuations of size ∼1/k\sim 1/k?". They are exact Fourier transforms, so no information is lost between them — but features look utterly different: the acoustic physics of the early Universe imprints a sinusoidal ripple on P(k)P(k) and a single localized bump on ξ(r)\xi(r) near the sound-horizon scale. The second figure is the punchline the subject builds toward: the measured ξ(r)\xi(r), a steep nonlinear power law giving way to a shallow linear regime carrying the acoustic ruler.

Figure
The measured correlation function on a log–log plot: a power law ξ = (r/r₀)^−γ with r₀ ≈ 5 h⁻¹ Mpc and γ ≈ 1.8 dominates below ~20 h⁻¹ Mpc and crosses ξ = 1 at r₀, while on large scales ξ falls steeply and carries a faint bump near 100 h⁻¹ Mpc — the baryon acoustic oscillation ruler set by the sound horizon
The punchline: ξ(r) — a power law (r/r₀)^−1.8 crossing ξ = 1 at r₀ = 5 h⁻¹ Mpc, plus the BAO bump near 100 0.001 0.01 0.1 1 10 100 0.1 1 10 100 separation r (h⁻¹ Mpc) correlation function ξ(r) r₀ = 5 h⁻¹ Mpc (ξ = 1) slope −γ = −1.8 nonlinear power law BAO bump ≈ 100 h⁻¹ Mpc

Two more facts round out the picture. Galaxies do not trace the dark matter faithfully: they light up only where the density field crests high enough to collapse a halo, so they clump more strongly than the mass — a bias b>1b>1 for bright galaxies, entering ξ\xi as b2b^2. And because the total galaxy number is fixed, the positive clustering piled up at small rr must be repaid by faint anticorrelation at large rr: the volume integral of ξ\xi vanishes, overdense knots being matter stolen from the surrounding voids.

The core equations

The cleanest definition is operational, phrased through pair counts. Consider two infinitesimal volumes separated by rr; the joint probability of finding one galaxy in each is the random expectation modulated by ξ\xi.

dP12=nˉ2 [1+ξ(r)] dV1 dV2(1)\mathrm{d}P_{12} = \bar{n}^{2}\,[1+\xi(r)]\,\mathrm{d}V_1\,\mathrm{d}V_2 \qquad (1)

where: dP12\mathrm{d}P_{12} — joint probability of a galaxy in each of two volume elements; nˉ\bar{n} — mean galaxy number density; ξ(r)\xi(r) — two-point correlation function; dV1,dV2\mathrm{d}V_1,\mathrm{d}V_2 — infinitesimal volume elements. Reading it: ξ=0\xi=0 makes the two events independent (Poisson); ξ>0\xi>0 enhances close pairs, ξ<0\xi<0 suppresses them. The same statistic is the autocorrelation of the density contrast δ\delta, the form connecting clustering to the field theory of perturbations; statistical isotropy collapses the dependence to the scalar separation r=∣r∣r=|\mathbf{r}|.

ξ(r)=⟨ δ(x) δ(x+r) ⟩(2)\xi(r) = \langle\, \delta(\mathbf{x})\,\delta(\mathbf{x}+\mathbf{r}) \,\rangle \qquad (2)

where: δ(x)=[ρ(x)−ρˉ]/ρˉ\delta(\mathbf{x})=[\rho(\mathbf{x})-\bar\rho]/\bar\rho — density contrast; ⟨⋯ ⟩\langle\cdots\rangle — ensemble (or spatial) average; r\mathbf{r} — separation vector. Reading it: ξ(0)=⟨δ2⟩\xi(0)=\langle\delta^2\rangle is the field variance, and the correlation decays as points decouple at large rr. Because ξ\xi is an autocorrelation, the Wiener–Khinchin theorem makes it the Fourier transform of P(k)P(k); the angular integral over the direction of k\mathbf{k} reduces the 3D transform to a 1D integral with a sinc kernel.

ξ(r)=1(2π)3∫P(k) e ik⋅r d3k=12π2∫0∞P(k) sin⁡(kr)kr k2 dk(3)\xi(r) = \frac{1}{(2\pi)^{3}}\int P(k)\,e^{\,i\mathbf{k}\cdot\mathbf{r}}\,\mathrm{d}^{3}k = \frac{1}{2\pi^{2}}\int_0^\infty P(k)\,\frac{\sin(kr)}{kr}\,k^{2}\,\mathrm{d}k \qquad (3)

where: P(k)P(k) — matter power spectrum, ⟨δ(k)δ∗(k′)⟩=(2π)3P(k) δD(k−k′)\langle\delta(\mathbf{k})\delta^{*}(\mathbf{k}')\rangle=(2\pi)^{3}P(k)\,\delta_D(\mathbf{k}-\mathbf{k}'); sin⁡(kr)/(kr)\sin(kr)/(kr) — isotropic sinc kernel from the angular average; kk — comoving wavenumber. Reading it: a sharp scale in P(k)P(k) becomes an oscillation in ξ(r)\xi(r) and vice versa. Observationally, on nonlinear scales the measured function collapses onto a strikingly simple power law over more than two decades of separation.

ξ(r)=(rr0)−γ,r0≈5 h−1 Mpc,  γ≈1.8(4)\xi(r) = \left(\frac{r}{r_0}\right)^{-\gamma},\qquad r_0\approx 5\,h^{-1}\,\mathrm{Mpc},\ \ \gamma\approx 1.8 \qquad (4)

where: r0r_0 — correlation length, the separation at which ξ=1\xi=1; γ\gamma — logarithmic slope of the clustering; hh — Hubble constant in units of 100 km s−1 Mpc−1100\,\mathrm{km\,s^{-1}\,Mpc^{-1}}. Reading it: r0r_0 sets the amplitude and γ\gamma the steepness of clustering in the halo regime. Galaxies are biased tracers, so the measured galaxy correlation is the matter correlation scaled up by a single constant on large linear scales.

ξgal(r)=b2 ξm(r)(5)\xi_{\mathrm{gal}}(r) = b^{2}\,\xi_{\mathrm{m}}(r) \qquad (5)

where: bb — linear bias factor, δgal=b δm\delta_{\mathrm{gal}}=b\,\delta_{\mathrm{m}}; ξm\xi_{\mathrm{m}} — matter correlation function; ξgal\xi_{\mathrm{gal}} — galaxy correlation function. Reading it: since ξ\xi is quadratic in δ\delta, the bias enters squared; redder, more luminous galaxies carry larger bb. To estimate ξ\xi from a real catalog, one compares observed pair counts against a synthetic random catalog sharing the survey's geometry and selection; the Landy–Szalay combination is the minimum-variance, edge-insensitive estimator.

ξ^(r)=DD−2 DR+RRRR(6)\hat{\xi}(r)=\frac{DD-2\,DR+RR}{RR} \qquad (6)

where: DDDD — normalized data–data pair count; DRDR — data–random; RRRR — random–random. Reading it: the −2DR+RR-2DR+RR structure cancels survey-boundary effects to first order, which is why Landy–Szalay beats the naive DD/RR−1DD/RR-1. Finally, the field's amplitude is summarized by its variance in spheres, the quantity that defines σ8\sigma_8 and normalizes the spectrum.

σ2(R)=12π2∫0∞P(k) W2(kR) k2 dk,W(x)=3[sin⁡x−xcos⁡x]x3(7)\sigma^{2}(R)=\frac{1}{2\pi^{2}}\int_0^\infty P(k)\,W^{2}(kR)\,k^{2}\,\mathrm{d}k,\qquad W(x)=\frac{3[\sin x - x\cos x]}{x^{3}} \qquad (7)

where: σ2(R)\sigma^{2}(R) — variance of mass fluctuations in top-hat spheres of comoving radius RR; σ8\sigma_8 uses R=8 h−1 MpcR=8\,h^{-1}\,\mathrm{Mpc}; W(kR)W(kR) — top-hat window. Reading it: σ8≈0.81\sigma_8\approx0.81 is the standard ΛCDM normalization; ξ(r)\xi(r), P(k)P(k), and σ(R)\sigma(R) are three windows on one field.

Derivation AFrom pair counts to the autocorrelation of the density contrast

The goal is to prove that the operational pair-count definition of ξ\xi in equation (1) is identical to the field-theoretic autocorrelation of the density contrast in equation (2), so that "excess pairs over random" and "⟨δδ⟩\langle\delta\delta\rangle" are literally the same object — the bridge on which the entire link between observation and theory rests. We start from the local number density, form the joint two-point density, and reduce it using the vanishing mean of δ\delta.

Step 1. Write the local number density in terms of the contrast, n(x)=nˉ [1+δ(x)]n(\mathbf{x})=\bar n\,[1+\delta(\mathbf{x})], the definition of δ\delta rearranged, with ⟨n⟩=nˉ\langle n\rangle=\bar n since ⟨δ⟩=0\langle\delta\rangle=0. Why: this hinge turns a statement about counting discrete galaxies into one about the continuous field.

Step 2. The joint probability of a galaxy in each of two volume elements is the expectation of the product of local densities, dP12=⟨n(x1) n(x2)⟩ dV1 dV2\mathrm{d}P_{12}=\langle n(\mathbf{x}_1)\,n(\mathbf{x}_2)\rangle\,\mathrm{d}V_1\,\mathrm{d}V_2. Why: a coincidence probability is the correlation of the two local rates, as in Poisson-process theory.

Step 3. Substitute and expand, ⟨n1n2⟩=nˉ2⟨(1+δ1)(1+δ2)⟩=nˉ2[1+⟨δ1⟩+⟨δ2⟩+⟨δ1δ2⟩]\langle n_1 n_2\rangle=\bar n^{2}\langle(1+\delta_1)(1+\delta_2)\rangle=\bar n^{2}[1+\langle\delta_1\rangle+\langle\delta_2\rangle+\langle\delta_1\delta_2\rangle]; the linear terms vanish because ⟨δ⟩=0\langle\delta\rangle=0. Why: the mean-zero contrast kills the cross terms, leaving only the genuine two-point correlation.

Step 4. What survives is dP12=nˉ2[1+⟨δ1δ2⟩] dV1 dV2\mathrm{d}P_{12}=\bar n^{2}[1+\langle\delta_1\delta_2\rangle]\,\mathrm{d}V_1\,\mathrm{d}V_2; comparing with equation (1) forces ξ(r)=⟨δ(x) δ(x+r)⟩\xi(r)=\langle\delta(\mathbf{x})\,\delta(\mathbf{x}+\mathbf{r})\rangle. Why: equal for all volume elements means equal coefficients — the pair-count excess is exactly the field autocorrelation.

Go deeper · Derivation B · 1 minThe Fourier transform relating the correlation function and the power spectrum
Derivation BThe Fourier transform relating the correlation function and the power spectrum

The goal is to derive the isotropic sinc-kernel integral of equation (3), showing that ξ(r)\xi(r) is the Fourier transform of P(k)P(k) and, crucially, that a localized feature in one domain becomes an oscillation in the other — the mathematical fact behind the BAO bump appearing as acoustic wiggles in P(k)P(k). We expand the contrast in Fourier modes, insert the power-spectrum definition to collapse the double integral, then do the angular integral to reach the radial form.

Step 1. Expand the contrast in Fourier modes, δ(x)=(2π)−3∫δ(k) e ik⋅x d3k\delta(\mathbf{x})=(2\pi)^{-3}\int\delta(\mathbf{k})\,e^{\,i\mathbf{k}\cdot\mathbf{x}}\,\mathrm{d}^3k. Why: clustering statistics are diagonal in Fourier space, so switching to modes makes the calculation tractable.

Step 2. Insert two copies into the autocorrelation: ξ(r)=(2π)−6∫ ⁣ ⁣∫⟨δ(k)δ∗(k′)⟩ e ik⋅xe−ik′⋅(x+r) d3k d3k′\xi(r)=(2\pi)^{-6}\int\!\!\int\langle\delta(\mathbf{k})\delta^{*}(\mathbf{k}')\rangle\,e^{\,i\mathbf{k}\cdot\mathbf{x}}e^{-i\mathbf{k}'\cdot(\mathbf{x}+\mathbf{r})}\,\mathrm{d}^3k\,\mathrm{d}^3k'. Why: computing ⟨δδ⟩\langle\delta\delta\rangle mode by mode exposes exactly the correlator P(k)P(k) encodes.

Step 3. Homogeneity locks the wavevectors together, ⟨δ(k)δ∗(k′)⟩=(2π)3P(k) δD(k−k′)\langle\delta(\mathbf{k})\delta^{*}(\mathbf{k}')\rangle=(2\pi)^{3}P(k)\,\delta_D(\mathbf{k}-\mathbf{k}'), and the Dirac delta collapses one integral, leaving ξ(r)=(2π)−3∫P(k) e−ik⋅r d3k\xi(r)=(2\pi)^{-3}\int P(k)\,e^{-i\mathbf{k}\cdot\mathbf{r}}\,\mathrm{d}^3k. Why: translation invariance means different Fourier modes do not correlate — this diagonalizes the problem.

Step 4. Do the angular integral over the sphere of fixed kk, ∫dΩk e−ik⋅r=4π sin⁡(kr)/(kr)\int\mathrm{d}\Omega_k\,e^{-i\mathbf{k}\cdot\mathbf{r}}=4\pi\,\sin(kr)/(kr), leaving ξ(r)=12π2∫0∞P(k) sin⁡(kr)kr k2 dk\xi(r)=\tfrac{1}{2\pi^{2}}\int_0^\infty P(k)\,\frac{\sin(kr)}{kr}\,k^{2}\,\mathrm{d}k. Why: the sinc kernel is the fingerprint of isotropy, and its oscillatory tail smears a knife-edge scale in P(k)P(k) into a single bump in ξ(r)\xi(r).

Go deeper · Derivation C · 1 minPower-law clustering from a scale-free power spectrum
Derivation CPower-law clustering from a scale-free power spectrum

The goal is to show that a scale-free spectrum P(k)∝knP(k)\propto k^{n} produces a scale-free correlation function ξ(r)∝r−(n+3)\xi(r)\propto r^{-(n+3)}, so that the measured slope γ≈1.8\gamma\approx1.8 directly reports an effective spectral index n≈−1.2n\approx-1.2 in the nonlinear regime — linking the observed power law of equation (4) to the spectrum's shape with no free normalization. We substitute the power law into the transform and extract the rr-dependence by a change of variables.

Step 1. Substitute P(k)=A knP(k)=A\,k^{n} into the radial transform, ξ(r)∝∫0∞kn sin⁡(kr)kr k2 dk=∫0∞kn+1 sin⁡(kr)r dk\xi(r)\propto\int_0^\infty k^{n}\,\frac{\sin(kr)}{kr}\,k^{2}\,\mathrm{d}k=\int_0^\infty k^{n+1}\,\frac{\sin(kr)}{r}\,\mathrm{d}k. Why: a single power law isolates the pure scaling behaviour, stripping away the transfer-function detail.

Step 2. Change variables to x=krx=kr, pulling every factor of rr outside: ξ(r)∝r−(n+3)∫0∞x n+1sin⁡(x) dx\xi(r)\propto r^{-(n+3)}\int_0^\infty x^{\,n+1}\sin(x)\,\mathrm{d}x, leaving a pure number (regularized). Why: scale-free inputs give scale-free outputs — the integral cannot reintroduce a scale it does not contain.

Step 3. Read off the slope γ=n+3\gamma=n+3; the measured γ≈1.8\gamma\approx1.8 implies an effective n≈−1.2n\approx-1.2, steeper than the primordial ns≈0.96n_s\approx0.96 because nonlinear collapse transfers power to small scales. Why: this converts a measured real-space slope into a statement about the spectrum, and its success is strong evidence for gravitational instability from Gaussian initial conditions.

Limiting cases

The correlation function wears a different face in each regime of separation, and the transitions between them carry the physics:

CaseResultMeaning
Small scales (nonlinear)ξ(r)≃(r/r0)−γ, r0≈5 h−1Mpc, γ≈1.8\xi(r)\simeq (r/r_0)^{-\gamma},\ r_0\approx 5\,h^{-1}\mathrm{Mpc},\ \gamma\approx 1.8Power-law clustering from nonlinear collapse and halo occupation; r0r_0 is the correlation length where ξ=1\xi=1.
Large scales (linear)ξgal(r)=b2 ξm(r)≪1\xi_{\mathrm{gal}}(r)=b^{2}\,\xi_{\mathrm{m}}(r)\ll 1Linear biasing; ξ\xi is small and traces the primordial spectrum shaped by the transfer function and linear growth.
BAO scalerBAO≈100 h−1Mpcr_{\mathrm{BAO}}\approx 100\,h^{-1}\mathrm{Mpc}A localized bump — the standard ruler at the comoving sound horizon; amplitude only ∼10−3\sim10^{-3}.
Very large scales∫0∞ξ(r) 4πr2 dr=0\int_0^\infty \xi(r)\,4\pi r^{2}\,\mathrm{d}r = 0Integral constraint: positive small-scale clustering is compensated by weak large-scale anticorrelation.
Zero lagξ(0)=⟨δ2⟩=σ2\xi(0)=\langle\delta^{2}\rangle=\sigma^{2}The correlation at zero separation is the variance of the field, set by σ8\sigma_8 when smoothed on 8 h−1Mpc8\,h^{-1}\mathrm{Mpc}.

Worked examples

1 · Power-law clustering: amplitudes, correlation length, and spectral slope (fully worked)

Assume the canonical fit ξ(r)=(r/r0)−γ\xi(r)=(r/r_0)^{-\gamma} with r0=5.0 h−1Mpcr_0=5.0\,h^{-1}\mathrm{Mpc} and γ=1.8\gamma=1.8. Evaluate ξ\xi at 11 and 10 h−1Mpc10\,h^{-1}\mathrm{Mpc}, interpret r0r_0, translate the excess into pair counts, and infer the effective spectral index.

Step 1. At r=1 h−1Mpcr=1\,h^{-1}\mathrm{Mpc}: ξ=(1/5)−1.8=51.8=18.1\xi=(1/5)^{-1.8}=5^{1.8}=\mathbf{18.1}. At r=10 h−1Mpcr=10\,h^{-1}\mathrm{Mpc}: ξ=(10/5)−1.8=2−1.8=0.287\xi=(10/5)^{-1.8}=2^{-1.8}=\mathbf{0.287}. Clustering is overwhelmingly strong at 1 h−1Mpc1\,h^{-1}\mathrm{Mpc} — eighteen times more pairs than random — and already modest by 10 h−1Mpc10\,h^{-1}\mathrm{Mpc}.

Step 2. By definition r0r_0 is where ξ=1\xi=1: at 5 h−1Mpc5\,h^{-1}\mathrm{Mpc} a galaxy is twice as likely to have a neighbour as a random point (1+ξ=21+\xi=2). So r0r_0 is the crossover from strong to weak clustering — an amplitude marker, not a physical clustering size.

Step 3. Convert to pair counts. In a shell of radius rr the expected neighbour count is dN=nˉ 4πr2[1+ξ(r)] dr\mathrm{d}N=\bar n\,4\pi r^{2}[1+\xi(r)]\,\mathrm{d}r. At 1 h−1Mpc1\,h^{-1}\mathrm{Mpc} the bracket is 1+18.1=19.11+18.1=\mathbf{19.1} — nineteen times the random count of close neighbours; at 10 h−1Mpc10\,h^{-1}\mathrm{Mpc} it is only 1.2871.287, a 29%29\% excess.

Step 4. Infer the spectrum. The relation γ=n+3\gamma=n+3 inverts to n=1.8−3=−1.2n=1.8-3=\mathbf{-1.2}: locally P(k)∝k−1.2P(k)\propto k^{-1.2}, far steeper than the primordial ns≈0.96n_s\approx0.96 — the fingerprint of nonlinear collapse piling power onto small scales. ξ(1) = 18.1, ξ(10) = 0.287, r₀ marks ξ = 1, effective n = −1.2.

2 · The baryon acoustic oscillation as a standard ruler (worked with numbers)

The comoving sound horizon at the drag epoch is rd=147.1 Mpcr_{\mathrm{d}}=147.1\,\mathrm{Mpc} (Planck). Convert it to h−1Mpch^{-1}\mathrm{Mpc} with h=0.674h=0.674, then find the redshift width Δz\Delta z of the radial BAO feature at z=0.35z=0.35 for a flat ΛCDM cosmology with Ωm=0.315\Omega_m=0.315.

Step 1. Convert units: rd=147.1 Mpc=147.1×0.674 h−1Mpc=99.1 h−1Mpcr_{\mathrm{d}}=147.1\,\mathrm{Mpc}=147.1\times 0.674\,h^{-1}\mathrm{Mpc}=\mathbf{99.1\,h^{-1}\mathrm{Mpc}} — exactly the ∼100 h−1Mpc\sim100\,h^{-1}\mathrm{Mpc} bump in ξ(r)\xi(r). The acoustic feature lands there because it is the sound horizon frozen into the baryons at recombination, carried outward and imprinted on the galaxy distribution.

Step 2. Compute the expansion rate at z=0.35z=0.35: H(z)=H0Ωm(1+z)3+ΩΛH(z)=H_0\sqrt{\Omega_m(1+z)^{3}+\Omega_\Lambda} with ΩΛ=0.685\Omega_\Lambda=0.685. Then (1.35)3=2.460(1.35)^{3}=2.460, so Ωm(1+z)3=0.775\Omega_m(1+z)^3=0.775, and H(z)=67.40.775+0.685=67.4×1.208=81.4 km s−1Mpc−1H(z)=67.4\sqrt{0.775+0.685}=67.4\times1.208=\mathbf{81.4\ km\,s^{-1}Mpc^{-1}}.

Step 3. The radial feature spans Δz=rd H(z)/c=(147.1)(81.4)/(2.998×105)=0.040\Delta z=r_{\mathrm{d}}\,H(z)/c=(147.1)(81.4)/(2.998\times10^{5})=\mathbf{0.040}. Measuring this Δz\Delta z at known rdr_{\mathrm{d}} delivers H(z)H(z); the transverse extent delivers the angular-diameter distance DA(z)D_A(z) — a geometric probe of dark energy independent of supernovae. r_d = 99.1 h⁻¹ Mpc, H(0.35) = 81.4 km s⁻¹ Mpc⁻¹, radial Δz ≈ 0.040.

3 · Galaxy bias from clustering amplitude (worked)

A luminous red galaxy sample has a measured clustering amplitude equivalent to σ8,gal=1.6\sigma_{8,\mathrm{gal}}=1.6, while the ΛCDM matter value is σ8,m=0.811\sigma_{8,\mathrm{m}}=0.811. Estimate the linear bias and compare ξgal\xi_{\mathrm{gal}} with ξm\xi_{\mathrm{m}}.

Step 1. On large scales the amplitudes relate by the bias, σ8,gal=b σ8,m\sigma_{8,\mathrm{gal}}=b\,\sigma_{8,\mathrm{m}}, so b=1.6/0.811=1.97≈2b=1.6/0.811=\mathbf{1.97}\approx2 — just what is expected for luminous red galaxies in massive halos at high density peaks.

Step 2. Since ξ\xi is quadratic in the contrast, ξgal=b2 ξm=(1.97)2 ξm=3.9 ξm\xi_{\mathrm{gal}}=b^{2}\,\xi_{\mathrm{m}}=(1.97)^{2}\,\xi_{\mathrm{m}}=\mathbf{3.9}\,\xi_{\mathrm{m}}: these galaxies cluster nearly four times more strongly than the dark matter, an offset that must be divided out before ξgal\xi_{\mathrm{gal}} probes the mass. b ≈ 1.97, so ξ_gal ≈ 3.9 ξ_m on linear scales.

Physics of the portrait

The correlation function is where a borrowed tool from statistical physics became the central instrument of a science: an excess-pair statistic from turbulence theory that one physicist recognized as the right lens for the clumping of galaxies, and built into a quantitative theory of how structure grows.

Origins · 1 minP. J. E. Peebles and how this was discovered
PJ
P. J. E. Peebles
1935– · physical cosmology and large-scale structure

Phillip James Edwin Peebles, born in Winnipeg in 1935, spent his career at Princeton turning cosmology from a speculative pursuit into a precision science, and the two-point correlation function was one of his sharpest tools. In a landmark series of papers through the late 1960s and 1970s, and in his 1980 monograph The Large-Scale Structure of the Universe, he imported the correlation function from statistical physics and made ξ(r)\xi(r) the standard descriptor of galaxy clustering, establishing the power-law form ξ∝r−1.8\xi\propto r^{-1.8} and connecting it, through the Fourier relation, to the matter power spectrum and the gravitational-instability picture of structure formation. He recognized that a nearly Gaussian primordial field renders the two-point function nearly complete, that galaxies must be biased tracers of an unseen mass, and that clustering statistics would become the arena in which the initial conditions of the early Universe are read off empirically — the same work anticipating the acoustic features later measured as the BAO ruler. For this foundational role Peebles shared the 2019 Nobel Prize in Physics.

History
1934
— Hubble finds that galaxy counts depart from a uniform random distribution
1969
— Totsuji and Kihara fit the power law ξ(r) ∝ r⁻¹·⁸ to galaxy catalogs
1980
— Peebles publishes The Large-Scale Structure of the Universe
1993
— Landy and Szalay introduce the minimum-variance correlation estimator
2005
— Eisenstein and collaborators detect the BAO peak in the SDSS correlation function

Going further

Texts. Peebles, The Large-Scale Structure of the Universe (1980) — the founding treatment of ξ(r)\xi(r) · Dodelson, Modern Cosmology, 2nd ed., Ch. 11 (correlation functions and power spectra) · Mo, van den Bosch & White, Galaxy Formation and Evolution, Ch. 6 (clustering statistics and bias) · Ryden, Introduction to Cosmology, 2nd ed., Ch. 12 · Carroll & Ostlie, An Introduction to Modern Astrophysics, 2nd ed., Ch. 30 · Binney & Tremaine, Galactic Dynamics, 2nd ed. · Landy & Szalay, ApJ 412, 64 (1993) for the estimator; Eisenstein et al., ApJ 633, 560 (2005) for the BAO detection. The key numerical anchors are r0≈5 h−1Mpcr_0\approx5\,h^{-1}\mathrm{Mpc}, γ≈1.8\gamma\approx1.8, and rBAO≈100 h−1Mpcr_{\mathrm{BAO}}\approx100\,h^{-1}\mathrm{Mpc}; the sound-horizon conversion rd=147.1 Mpc⇒99.1 h−1Mpcr_{\mathrm{d}}=147.1\,\mathrm{Mpc}\Rightarrow99.1\,h^{-1}\mathrm{Mpc} and the radial Δz=0.040\Delta z=0.040 at z=0.35z=0.35 use Planck 2018 parameters (h=0.674h=0.674, Ωm=0.315\Omega_m=0.315).

Prerequisites. Density Perturbations and the Density Contrast · Transfer Function and the Matter Power Spectrum · The Primordial Power Spectrum · Linear Growth of Structure. Leads to. Baryon Acoustic Oscillations · The Sound Horizon and the Standard Ruler · The ΛCDM Concordance Model · Dark Energy and the Cosmological Constant.

Ask the tutor: Ask me to walk through the sinc-kernel derivation, work a BAO distance measurement at your favourite redshift, or show exactly how redshift-space multipoles pull the growth rate fσ8f\sigma_8 out of ξ(s,μ)\xi(s,\mu).

Can the correlation function really be negative, and what would that mean

Yes, and it must be somewhere. A negative ξ\xi means anticorrelation — fewer pairs than random — which is what you find on very large scales just beyond the acoustic peak. It is not an error but a requirement: because the total number of galaxies is fixed, the strong positive clustering at small separations must be balanced so that ∫0∞ξ(r) 4πr2 dr=0\int_0^\infty\xi(r)\,4\pi r^{2}\,\mathrm{d}r=0. The matter piled into overdense filaments was drawn out of the surrounding voids, so an excess of close pairs is paid for by a deficit of distant ones.

Why is ξ(r) a power law on small scales when the primordial spectrum is nearly scale-invariant

The primordial spectrum has ns≈0.96n_s\approx0.96, which by γ=n+3\gamma=n+3 would give a very steep γ≈4\gamma\approx4 — nothing like the observed 1.81.8. The resolution is nonlinear gravitational collapse: as overdense regions virialize into halos, they transfer power to small scales, steepening the effective slope of P(k)P(k) there to n≈−1.2n\approx-1.2 and yielding γ≈1.8\gamma\approx1.8. The clean power law is thus not primordial but an attractor of nonlinear clustering — a striking regularity and a stringent test of NN-body simulations.

What is galaxy bias, and why does it enter ξ as b squared

Galaxies do not form everywhere the dark matter is dense — they light up only where the field crests high enough to collapse and cool gas into a halo, so they mark high peaks and clump more strongly than the mass. Linear bias captures this as δgal=b δm\delta_{\mathrm{gal}}=b\,\delta_{\mathrm{m}}. Because the correlation function is quadratic in the contrast, ξ=⟨δδ⟩\xi=\langle\delta\delta\rangle, the bias appears squared: ξgal=b2ξm\xi_{\mathrm{gal}}=b^{2}\xi_{\mathrm{m}}. For bright red galaxies b∼2b\sim2, so ξgal∼4 ξm\xi_{\mathrm{gal}}\sim4\,\xi_{\mathrm{m}}. On small scales bias becomes scale-dependent and stochastic, a central modeling challenge in precision clustering cosmology.

How is ξ(r) actually estimated from a galaxy survey

You compare the survey to synthetic randomness. Build a random catalog with the same sky footprint, radial selection, and completeness as the data but no clustering, then count pairs at each separation: data–data (DDDD), data–random (DRDR), and random–random (RRRR). The Landy–Szalay estimator ξ^=(DD−2DR+RR)/RR\hat\xi=(DD-2DR+RR)/RR combines them into a nearly unbiased, minimum-variance estimate whose −2DR+RR-2DR+RR structure cancels survey edges and holes to first order. The random catalog must be large — tens of times the data — to keep shot noise in RRRR negligible.

Why distinguish redshift-space ξ(s) from real-space ξ(r)

Because we infer distances from redshifts, which mix the smooth Hubble flow with peculiar velocities. On large scales, coherent infall toward overdensities compresses structures along the line of sight and boosts the apparent clustering (the Kaiser effect); on small scales, random virial motions stretch them into radial "fingers of God". Both distort the map from true rr to observed ss, so ξ(s)≠ξ(r)\xi(s)\neq\xi(r) and the redshift-space function is anisotropic. That anisotropy is not a nuisance but signal: its Legendre multipoles ξℓ(s,μ)\xi_\ell(s,\mu) isolate the growth rate fσ8f\sigma_8, a direct test of gravity.

Does the BAO bump give the same information as the acoustic wiggles in P(k)

Yes — they are the same feature in two Fourier-conjugate domains. The baryon acoustic oscillations are a genuine oscillation in P(k)P(k), a series of decaying wiggles spaced by the sound horizon; transforming to real space, that oscillation localizes into a single broad bump in ξ(r)\xi(r) near 100 h−1Mpc100\,h^{-1}\mathrm{Mpc}, because the sinc kernel maps a periodic feature to a localized one. Analysts use both: ξ(r)\xi(r) shows the ruler as one isolated peak easy to model against a smooth background, while P(k)P(k) spreads the signal across many wiggles convenient for the survey window. Fitting either yields the same distance–redshift constraints.

Before this

Where this leads