Pair Correlation and Random Matrices

In the spring of 1972, a young number theorist named Hugh Montgomery walked into the common room at the Institute for Advanced Study in Princeton to take tea. He had been studying how the non-trivial zeros of the zeta function are spaced — not just that they sit on the critical line, but how they cluster and repel each other going up it. He had found a formula for the statistical pattern of the gaps. Over tea, the physicist Freeman Dyson asked what he was working on. Montgomery wrote down his function. Dyson looked at it and said: that is the pair correlation of the eigenvalues of a random Hermitian matrix.

Neither man had gone looking for the other's subject. Montgomery was doing pure number theory; Dyson had spent the 1960s modelling the energy levels of heavy atomic nuclei with random matrices. And yet the same formula governed both. This is one of the most astonishing coincidences in all of mathematics — a bridge, still only partly built, between the primes and quantum physics. This page is about that bridge.

First, unfold the zeros

Write the non-trivial zeros as \tfrac12 + i\gamma_n, with heights 0 < \gamma_1 \le \gamma_2 \le \cdots marching up the critical line. To talk about their spacing we hit a snag: the zeros get denser higher up. The number of zeros with height up to T grows like

N(T) \sim \frac{T}{2\pi}\log\frac{T}{2\pi},

so near height \gamma the local density of zeros is about \tfrac{1}{2\pi}\log\tfrac{\gamma}{2\pi} and the average gap shrinks to zero as you climb. Comparing raw gaps at different heights would be comparing apples to oranges. The fix is unfolding: rescale each zero so that the average gap becomes exactly 1 everywhere. Define

\tilde\gamma_n = \frac{\gamma_n}{2\pi}\log\frac{\gamma_n}{2\pi}.

Now the rescaled heights \tilde\gamma_n have unit mean spacing at every height. Only after this normalisation does it make sense to ask what the gaps look like — and only the local statistics of the unfolded zeros are what match random matrices.

Montgomery's pair correlation conjecture

The pair correlation asks: given a zero, how likely is another zero to sit a (normalised) distance u away? Averaged over the whole line, this density is captured by a single function.

The term \big(\tfrac{\sin\pi u}{\pi u}\big)^2 is the square of the sinc function. Everything interesting about the zeros' fine structure is hiding in that one subtracted bump — so let us look at it directly.

Seeing the level repulsion

Below is Montgomery's density R_2(u) = 1 - \big(\tfrac{\sin\pi u}{\pi u}\big)^2. The horizontal axis is the normalised separation u between two zeros; the vertical axis is how densely pairs occur at that separation, relative to a purely random (Poisson) scatter for which the density would be a flat 1.

Two features jump out. Near u = 0 the curve is pushed down to zero: the zeros repel each other, almost never crowding together. Far out, the curve flattens to 1: widely separated zeros behave as if independent. Between them, a gentle overshoot. This shape is not what random dots would give — random dots (a Poisson process) have pair correlation identically 1, with no repulsion at all. The zeros are more ordered than random.

Worked example — the two ends of the curve

Near u = 0 (repulsion). Use the Taylor expansion \sin x = x - \tfrac{x^3}{6} + \cdots with x = \pi u:

\frac{\sin \pi u}{\pi u} = 1 - \frac{(\pi u)^2}{6} + O(u^4), \qquad \left(\frac{\sin \pi u}{\pi u}\right)^2 = 1 - \frac{(\pi u)^2}{3} + O(u^4).

Subtracting from 1,

R_2(u) = 1 - \left(\frac{\sin\pi u}{\pi u}\right)^2 = \frac{\pi^2 u^2}{3} + O(u^4) \;\longrightarrow\; 0 \quad\text{as } u\to 0.

So the density doesn't just dip — it vanishes quadratically. The chance of finding two zeros a tiny distance apart falls off like u^2. That is level repulsion: the zeros actively avoid one another.

For large u (independence). Since |\sin \pi u| \le 1, the sinc term is squeezed:

\left(\frac{\sin\pi u}{\pi u}\right)^2 \le \frac{1}{\pi^2 u^2} \;\longrightarrow\; 0, \qquad\text{so}\qquad R_2(u) \longrightarrow 1.

The correction decays like 1/u^2 (with a tiny ripple from the \sin^2), and the density relaxes to the "uncorrelated" value 1. Distant zeros carry no memory of one another.

Enter the Gaussian Unitary Ensemble

What Dyson recognised is that 1 - \big(\tfrac{\sin\pi u}{\pi u}\big)^2 is exactly the pair-correlation function for the eigenvalues of a large random matrix drawn from the Gaussian Unitary Ensemble (GUE). A GUE matrix is a big N\times N Hermitian matrix whose entries are independent Gaussians (real on the diagonal, complex off it), chosen so the whole distribution is invariant under unitary change of basis. Such a matrix has real eigenvalues, and — after the same kind of unfolding to unit mean spacing — those eigenvalues repel each other with precisely Montgomery's density.

Dyson had built the GUE (and its cousins the GOE and GSE) to model the energy levels of heavy nuclei, where the exact Hamiltonian is hopelessly complicated but its statistics are universal. That the same universal law should describe the zeros of a function built from the prime numbers is the deep mystery — and the strongest hint yet at why the Riemann Hypothesis might be true.

Odlyzko's numerical bombshell

A conjecture this strange demands evidence, and in the 1980s Andrew Odlyzko supplied it on a breathtaking scale. Using the Riemann–Siegel formula he computed zeros not near the bottom of the critical line but enormously high up — around the 10^{20}\text{th} zero and beyond — and gathered millions of consecutive ones. Then he formed the histogram of their normalised gaps and overlaid the GUE prediction.

The match was, in his own word, stunning. The empirical nearest-neighbour spacing curve lay on top of the GUE curve so tightly that the two were nearly indistinguishable, right down to the small wiggles. There was the level repulsion at short range, the correct hump, the correct tail — all of it. No theorem forced this agreement; it simply appeared, out to as many decimal places as the computation could resolve. The pictures from the "Odlyzko–Schönhage" computations are now the canonical evidence that the zeta zeros really do live in the random-matrix universality class.

Hilbert–Pólya: could the zeros be eigenvalues?

Why would prime-built zeros know about random matrices at all? The oldest guess is the Hilbert–Pólya conjecture: perhaps there is a self-adjoint (Hermitian) operator H whose eigenvalues are exactly the imaginary parts \gamma_n of the zeros.

This would turn RH from a statement about a function into a statement about spectra: find the right H and the hypothesis falls out for free. Nobody has found it. But the random-matrix statistics are exactly what you would expect if such an H existed, which is why the coincidence feels less like an accident and more like a clue.

Berry–Keating and quantum chaos

Michael Berry and Jonathan Keating sharpened the guess. In quantum physics, when the underlying classical system is chaotic and has no time-reversal symmetry, the energy levels of the quantised system follow GUE statistics — this is the Bohigas–Giannoni–Schmit universality of quantum chaos. So the sought-after H should be the quantisation of a chaotic classical system. Their candidate is startlingly simple: a Hamiltonian like

H = xp \quad(\text{position times momentum}),

whose classical flow is unstable (chaotic), and whose semiclassical level counting reproduces the smooth term \tfrac{T}{2\pi}\log\tfrac{T}{2\pi} - \tfrac{T}{2\pi} in N(T). The missing ingredient — the classical orbits that would pin down the operator — correspond mysteriously to the primes: in these heuristics the primes play the role of periodic orbits, and the zeros the role of energy levels. It is a dictionary between number theory and quantum chaos that no one has yet made rigorous, but that keeps predicting the right answers.

Consequences: moments and the Keating–Snaith conjecture

If zeta really is a random-matrix object, you can turn the analogy into predictions. The hardest open problems about the moments of zeta on the critical line,

M_k(T) = \frac{1}{T}\int_0^T \left|\zeta\!\left(\tfrac12 + it\right)\right|^{2k}\,dt,

had been solved only for k = 1 (Hardy–Littlewood, \sim \log T) and k = 2 (Ingham, \sim \tfrac{1}{2\pi^2}(\log T)^4) for most of a century. Keating and Snaith computed the corresponding moments of the characteristic polynomial of a random unitary matrix (the CUE), matched the size of the matrix to the density of zeros, and read off a prediction for every k:

These moment predictions were completely out of reach of classical analytic number theory, and yet the random-matrix model handed them over almost for free. That the model predicts new arithmetic, rather than merely matching old numerics, is why it is taken so seriously.

Two cautions, both easy to overstate in the excitement.

It is not proven. Montgomery proved his pair-correlation formula only in a restricted form (tested against smooth weights whose "frequencies" stay in a bounded range), and even that assumed the Riemann Hypothesis. The full pair correlation, the higher correlations, the GUE gap distribution, the Keating–Snaith moments — these are supported by beautiful heuristics and overwhelming numerics (Odlyzko), but they remain conjectures. No self-adjoint Hilbert–Pólya operator has ever been exhibited.

It is a statement about LOCAL statistics, after unfolding. Random matrices do not predict where any individual zero is, nor the global count N(T) (that is fixed by the smooth main term). What matches GUE is the fine-grained pattern of spacings once you have rescaled to unit mean gap. Say "the local correlations of the unfolded zeros look like GUE", not "the zeros are eigenvalues of a random matrix" — the second sentence throws away every caveat that makes the first one true.

The joint probability density of the eigenvalues \lambda_1,\dots,\lambda_N of a GUE matrix contains the factor \prod_{i. Look at what happens when two eigenvalues try to coincide: as \lambda_i \to \lambda_j that factor goes to zero, so configurations with near-collisions are strongly suppressed. The eigenvalues behave like charged particles that hate being close. The exponent 2 is special to the unitary ensemble (it is the "\beta = 2" case); the orthogonal and symplectic ensembles use exponents 1 and 4 and repel more weakly or more strongly. The zeta zeros pick out \beta = 2 — the same class as chaotic systems with broken time-reversal symmetry.

The Montgomery–Dyson conversation really did happen over afternoon tea at the Institute for Advanced Study, reportedly arranged by the number theorist Sarvadaman Chowla, who insisted the two be introduced. Montgomery, then a graduate student, had no reason to know that Dyson's random-matrix formula from nuclear physics would be the same one he had just derived for the primes. The story is now a favourite parable about why mathematicians and physicists should keep talking to one another — and about how much can turn on a single cup of tea.