Literature cutoff: 11 July 2026.

On 28–29 June 2023, four collaboration papers appeared together: NANOGrav (North American Nanohertz Observatory for Gravitational Waves), the EPTA (European Pulsar Timing Array) with InPTA (Indian Pulsar Timing Array) data, the PPTA (Parkes Pulsar Timing Array, Australian), and the CPTA (Chinese Pulsar Timing Array). They reported evidence for a stochastic gravitational-wave background at nanohertz frequencies. The four analyses used different data sets and pipelines; the EPTA and InPTA result was a joint analysis, not a fifth independent measurement.

The detector for all four collaborations was not a machine. It was the Milky Way — specifically, a collection of millisecond pulsars scattered across the galaxy, separated by thousands of light-years. The arms of this instrument make LIGO look microscopic by comparison, and they are pointed in every direction at once.

There is a professional habit, well established in physics, of finding the scale of things either comforting or terrifying. This one is genuinely awe-inspiring.

The scale problem

There is no general rule that a gravitational-wave detector must be as large as the wavelength it observes. LIGO’s 4-km arms are much shorter than the wavelengths it detects; its sensitivity comes from interferometry, long effective optical paths, and extraordinary displacement precision. For a pulsar timing array, the observing span sets the lowest frequencies that can be distinguished, while Earth–pulsar baselines provide the galaxy-scale paths on which a wave imprints a timing shift.

LIGO detects gravitational waves in the frequency range of roughly 10–1000 Hz. A wave at 100 Hz has a wavelength of

$$\lambda = \frac{c}{f} = \frac{3 \times 10^8 \text{ m/s}}{100 \text{ Hz}} \approx 3000 \text{ km},$$

about a quarter of the Earth’s diameter. LIGO’s 4-km Fabry-Pérot arm cavities are short compared to the wavelength, but the cavities fold the light hundreds of times to achieve an effective path length of about 1600 km, and the interferometer is exquisitely sensitive to differential arm-length changes of order $10^{-19}$ m — a fraction of the proton radius. The physics is heroic. But it works at the scale of stellar-mass binary mergers: neutron star–neutron star, neutron star–black hole, stellar-mass black hole–black hole. Those systems radiate at audible frequencies. LIGO famously converts its strain signals to sound, and the chirps do sound like something from a science fiction film.

Now consider pulsar timing arrays. The targets are nanohertz gravitational waves: frequencies of order $f \sim 10^{-9}$ to $10^{-8}$ Hz. At $10^{-9}$ Hz, the wavelength is

$$\lambda = \frac{c}{f} \approx \frac{3 \times 10^8 \text{ m/s}}{10^{-9} \text{ Hz}} = 3 \times 10^{17} \text{ m} \approx 10 \text{ pc} \approx 33 \text{ light-years}.$$

To detect oscillations at these frequencies, you need arms measured in light-years. Tens to thousands of light-years. The nearest millisecond pulsars are a few hundred light-years away; the most distant ones used in PTA arrays are several thousand light-years distant. By accident of the galaxy’s size and the distribution of recycled pulsars within it, we happen to live inside an instrument of approximately the right dimensions to detect nanohertz gravitational waves. The galaxy did not plan this. We got lucky.

Millisecond pulsars as cosmic clocks

A gravitational wave detector is only as good as its clock. LIGO measures differential length changes in its arm cavities using the interference of laser light — effectively using the constancy of the speed of light as a metronome. Pulsar timing arrays use a different clock: the rotation of a neutron star.

Ordinary pulsars are the collapsed remnants of massive stars, born in supernova explosions. A neutron star of roughly 1.4 solar masses is compressed into a sphere about 10 km across, rotating rapidly and radiating beams of radio waves from near its magnetic poles. The rotation periods of young pulsars are typically of order seconds, and they spin down on timescales of millions of years as they lose rotational energy to magnetic dipole radiation. They are not particularly stable clocks — the spin-down is uneven, and they occasionally “glitch,” suddenly spinning up by a small amount.

Millisecond pulsars (MSPs) are a different beast entirely. They have been recycled: spun up to near-millisecond rotation periods by accreting matter from a binary companion star. The accretion process deposits angular momentum onto the neutron star, increasing its spin rate dramatically, while simultaneously burying and weakening its magnetic field. The typical MSP has a surface magnetic field of order $10^4$–$10^5$ G, four to five orders of magnitude weaker than a young pulsar’s $10^9$–$10^{12}$ G. Since the spin-down torque scales as $B^2$, the weaker field means the MSP loses rotational energy far more slowly. Once the accretion stops — when the companion has exhausted its transferable mass — the MSP is left spinning rapidly and stably, with a spin-down rate of order $10^{-20}$ s/s.

The result is rotational stability competitive with terrestrial atomic clocks. PSR J0437$-$4715, one of the best-timed MSPs, has a rotation period of $P \approx 5.76$ ms and a timing residual — the scatter of individual pulse arrival times around the best-fit timing model — of order 100 ns over decades of observation. For a pulsar completing about 174 rotations per second, a residual of 100 ns over a baseline of years is remarkable. The fractional frequency stability is $\delta P / P \sim 10^{-14}$ or better. These are not merely good clocks; they are among the most stable periodic phenomena known to physics.

The timing model accounts for everything we know about the pulsar’s environment: its spin period and spin-down rate, its proper motion across the sky, parallax (from which we get the distance), Shapiro delay from any companions, dispersion measure variations in the interstellar medium, and more. After subtracting all modelled effects, what remains are timing residuals — small, unexplained deviations in the pulse arrival times. If a gravitational wave passes through, it will appear in those residuals.

How a gravitational wave shifts a pulse arrival time

A gravitational wave is a propagating perturbation of the spacetime metric. In the transverse-traceless (TT) gauge, a wave propagating along the $z$-direction perturbs the metric as

$$ds^2 = -c^2 dt^2 + \bigl(1 + h_+\bigr)dx^2 + \bigl(1 - h_+\bigr)dy^2 + 2h_\times \, dx \, dy + dz^2,$$

where $h_+$ and $h_\times$ are the two polarisation amplitudes, both functions of $t - z/c$. As the wave passes, the proper distances in the $x$- and $y$-directions oscillate: space itself is being stretched and squeezed at the wave frequency.

A radio pulse travelling from a pulsar to Earth traverses this oscillating spacetime. The proper path length changes, and so does the travel time. The timing residual induced in a pulsar in direction $\hat{n}$ by a gravitational wave with wavevector $\hat{k}$ is

$$R(t) \propto \int_0^t dt' \, \frac{\Delta\nu(t')}{\nu} \propto \int_0^t dt' \, h(t', \hat{n}, \hat{k}),$$

where $h$ is an appropriate contraction of the metric perturbation with the geometry of the Earth-pulsar baseline. The key point: the timing residual is the time integral of the fractional frequency shift. For a sinusoidal contribution of strain amplitude $h$ at frequency $f$, its residual scale is of order $h/(2\pi f)$, with geometry and the pulsar term included in the full response. The observing span determines which low frequencies can be separated from fitted timing-model trends; it is not an extra factor in this single-wave amplitude.

For a single pulsar, a timing residual tells you that something disturbed the spacetime between Earth and the pulsar — but it could be a systematic in the timing model, interstellar medium fluctuations, or intrinsic pulsar noise. You cannot claim a gravitational wave detection from one pulsar alone. What you need is the correlation between many pulsars.

The Hellings-Downs curve

Here is the central idea of pulsar timing array science. Consider two pulsars in directions $\hat{n}_a$ and $\hat{n}_b$, separated on the sky by an angle $\theta$ such that $\cos\theta = \hat{n}_a \cdot \hat{n}_b$. Both are embedded in the same stochastic gravitational wave background — a superposition of waves arriving from all directions, at all frequencies in the nanohertz band, with random phases and amplitudes. The timing residuals of both pulsars will be perturbed by this background. The question is: what is the expected cross-correlation between their residuals?

Hellings and Downs (1983) computed this, assuming the background is isotropic (equal power from all directions), unpolarised, and stationary. The answer is now called the Hellings-Downs (HD) curve:

$$\Gamma(\theta) = \frac{3}{2} \left(\frac{1-\cos\theta}{2}\right) \ln\!\left(\frac{1-\cos\theta}{2}\right) - \frac{1}{4}\left(\frac{1-\cos\theta}{2}\right) + \frac{1}{2}.$$

Let me unpack the features of this function:

  • At $\theta = 0$ for two distinct pulsars in the same direction: $\Gamma(0) = 1/2$, the maximum cross-correlation. For a pulsar correlated with itself, the conventional autocorrelation includes an additional $1/2$ term and is normalised to one.
  • At $\theta \approx 50°$: the curve crosses zero and turns negative.
  • At $\theta = \pi/2$ (pulsars at right angles): $\Gamma(\pi/2) \approx -0.15$, near the minimum of the curve.
  • At $\theta = \pi$ (antipodal pulsars, opposite directions on the sky): $\Gamma(\pi) = 1/4$. Positive correlation even for pulsars in opposite directions — counterintuitive but correct.
  • In between, the curve dips negative (anticorrelated) for angles roughly $50°$–$120°$.

The shape is quadrupolar. Under the usual assumptions of an isotropic, unpolarised tensor background in general relativity, it arises from spin-2 gravitational waves. A clock error common to all pulsars would produce a flat, angle-independent correlation; an error in Earth’s ephemeris would produce an approximately dipolar, $\cos\theta$ pattern. These are useful competing models, which is why PTA analyses fit them rather than merely looking for low-frequency noise.

This is why the HD curve is the central discriminant. Cross-correlations with this specific angular dependence are evidence for the tensor-background model, but their strength and robustness still depend on noise modelling, sampling, and the alternatives tested.

The difficulty is statistical. You need enough pulsars, spread over a wide range of sky angles, each timed with sufficient precision, over a long enough baseline, that you can measure this correlation function with confidence. That has been the programme of the PTA collaborations since the 1990s. In June 2023, they had enough.

The 2023 evidence

The NANOGrav 15-year dataset (Agazie et al., 2023) comprises 68 millisecond pulsars observed for up to 15 years, with an average of roughly 2200 timing observations per pulsar. The dataset represents an enormous investment of telescope time — primarily at the Arecibo Observatory (until its collapse in December 2020) and the Green Bank Telescope.

The analysis found an excess of low-frequency noise common to many pulsars, consistent with a power-law spectral shape. More importantly, when the cross-correlations between all pairs of pulsars were computed and binned by angular separation, the result was consistent with the Hellings-Downs curve. The statistical significance of the HD correlation — that is, the evidence that the spatial correlation pattern matches the quadrupolar prediction rather than some isotropic or zero-correlation model — was 3–4$\sigma$ depending on the analysis method and prior assumptions. The collaboration carefully described this as “evidence for” rather than “detection of” a gravitational wave background, following community conventions (detection would conventionally require $\geq 5\sigma$).

Simultaneously, the EPTA published its second data release (Antoniadis et al., 2023), the PPTA published its results (Reardon et al., 2023), and InPTA contributed its analysis. A fifth collaboration, the Chinese Pulsar Timing Array (CPTA), also published consistent results (Xu et al., 2023). Each group used independent datasets, different telescopes, different software pipelines, different statistical methodologies. All found the same thing.

The agreement among the four contemporaneous analyses is important, but it does not turn their evidence into a single independent five-way measurement: their data volumes, pulsar samples, telescopes, and analysis choices differ, and the EPTA–InPTA analysis is joint. The result is therefore best described, as the collaborations did, as convergent evidence for a background rather than a completed proof of its origin.

The combined interpretation is that the timing data contain a common low-frequency signal with spatial correlations consistent with a tensor gravitational-wave background. The Milky Way is a gravitational-wave detector, and it has measured something.

What is making the noise?

The million-dollar question. The observed signal has a characteristic spectral shape: the power spectral density of the timing residuals scales approximately as $S(f) \propto f^{-13/3}$, or equivalently, the gravitational wave energy density spectrum

$$\Omega_{\rm GW}(f) \propto f^{2/3}.$$

This $f^{2/3}$ scaling is the expected spectrum for an ensemble of circular binary systems driven purely by gravitational wave emission — specifically, a population of supermassive black hole binaries (SMBHBs) distributed across a wide range of masses and redshifts.

Here is the astrophysical picture. When two massive galaxies merge — an event that happens billions of times over cosmic history — their central supermassive black holes (each typically $10^7$–$10^{10}$ solar masses) do not immediately merge. First, they sink toward the centre of the merger remnant by dynamical friction against the stellar background, forming a loosely bound binary on scales of a parsec or so. The binary then hardens: three-body interactions with individual stars passing close to the binary extract orbital energy, driving the pair to smaller separations. Eventually, when the binary has hardened to a separation where gravitational wave emission dominates the energy loss, the pair inspirals and merges on a timescale of millions to billions of years.

The incoherent superposition of gravitational wave emission from the many billions of SMBHB systems across the observable universe — at all masses, at all orbital frequencies, at all redshifts — produces a stochastic background. It is, in a sense, cosmic traffic noise: no individual merger is detectable as a discrete event, but the combined hum from all of them is. The predicted spectral amplitude and shape from this population are broadly consistent with the observed signal, though the uncertainties on the astrophysical model are large enough that the agreement is not a precise test.

Alternative sources are possible. A first-order phase transition in the early universe (such as a QCD or electroweak phase transition) would produce a background of gravitational waves with a different spectral shape — more peaked and potentially harder than the SMBHB prediction. A network of cosmic strings, topological defects from symmetry-breaking phase transitions in the early universe, would produce yet another spectrum, approximately flat in $\Omega_{\rm GW}(f)$. Primordial gravitational waves from inflation are expected at far lower frequencies, near the CMB scale, but with a nearly scale-invariant spectrum that could contribute at nanohertz frequencies in some models. The 2023 data are most consistent with the SMBHB interpretation, but cannot rule out contributions from early-universe sources — or combinations of both. Future data with longer baselines and more pulsars will tighten the spectral measurements and may reveal deviations from the SMBHB prediction.

There is also the tantalising prospect of eventually resolving individual SMBHB systems above the stochastic background — the gravitational wave equivalent of resolving individual radio sources above the extragalactic background. No individual SMBHB has yet been identified in PTA data, but the 15-year NANOGrav dataset already places interesting upper limits on the most massive known candidate systems.

Multi-frequency gravitational wave astronomy

We are now, unambiguously, in the era of multi-frequency gravitational wave astronomy. The universe produces gravitational waves across an enormous span of frequencies, and different frequency windows probe fundamentally different source populations:

  • Primordial / CMB B-modes (~$10^{-17}$ Hz): Wavelengths comparable to the Hubble scale. Primordial tensor perturbations from inflation would imprint a distinctive B-mode polarisation pattern in the cosmic microwave background. No confirmed detection yet; the BICEP/Keck programme and future CMB experiments (LiteBIRD, CMB-S4) are sensitive to this regime.

  • Pulsar timing arrays (~$10^{-9}$–$10^{-8}$ Hz): The nanohertz band, just accessed by NANOGrav and partners. Sources: SMBHB inspirals, possibly early-universe phase transitions and cosmic strings. Arm length: thousands of light-years.

  • LISA (~$10^{-4}$–$10^{-1}$ Hz): The millihertz band. The LISA space interferometer (planned launch in the 2030s) will have arm lengths of 2.5 million km and will detect SMBHB mergers directly as they happen, stellar-mass compact binary inspirals years before their LIGO-band merger, and possibly signals from extreme mass-ratio inspirals. LISA will see the SMBHB sources that PTAs see as a stochastic background, but in their final years of inspiral.

  • LIGO/Virgo/KAGRA (~10–$10^3$ Hz): The audible band. Stellar-mass black hole and neutron star mergers. Over 90 confirmed events as of the end of O3, with many more candidates in O4. Source masses: $\sim$1–100 $M_\odot$.

This is the gravitational wave analogue of multi-wavelength astronomy. Just as the universe looks completely different in radio, infrared, optical, X-ray, and gamma-ray light — each wavelength band revealing different physical processes and source populations — the universe sounds completely different at each gravitational wave frequency. PTAs hear the rumble of cosmic structure formation; LISA will hear the whisper of the final million years before the most massive black hole mergers; LIGO hears the sharp crack of stellar-mass collisions.

The 2023 announcements represent the opening of the nanohertz window. We have gone from one gravitational wave frequency band (LIGO/Virgo) to two. The next decade, with LISA launching and PTAs continuing to accumulate data, will see the opening of a third.

A brief detour: LIGO’s O4 run and the mass-gap object

While the PTA collaborations were making their June 2023 announcement, LIGO’s fourth observing run (O4, 24 May 2023 to 18 November 2025) was proceeding at an extraordinary rate. The network issued roughly 250 public candidate alerts during the run. This is now a production science rather than an exploration.

Among the most scientifically interesting events was GW230529, detected on 29 May 2023 and published by the LIGO Scientific Collaboration (Abbott et al., 2024). The signal is consistent with the merger of a neutron star with a compact object whose mass was measured to be approximately $2.5$–$4.5\ M_\odot$. This mass range sits squarely in what theorists call the “mass gap” — the range between the heaviest neutron stars ($\lesssim 2.3\ M_\odot$, though this upper limit is uncertain) and the lightest stellar-mass black holes inferred from X-ray binaries ($\gtrsim 5\ M_\odot$).

Whether GW230529’s companion was the heaviest neutron star ever observed, or the lightest black hole, is genuinely unknown. The distinction matters enormously for nuclear physics: if it is a neutron star, it constrains the nuclear equation of state at supranuclear densities. If it is a black hole, it means the mass gap is narrower than X-ray observations suggested, and our understanding of compact object formation needs revision. Gravitational wave observations alone cannot distinguish a rapidly spinning heavy neutron star from a slowly spinning light black hole without additional electromagnetic counterpart observations, and no counterpart was found for GW230529. This question will likely not be settled until we have electromagnetic constraints from similar systems, or until we accumulate enough mass-gap events to understand the population statistically.

The instrument we did not build

I want to return to the pulsar timing array concept, because I think it deserves more than a passing technical description. The idea is this: nature has distributed a set of extremely stable clocks across the galaxy. We did not put them there. We did not design them. We simply discovered that neutron stars, after a particular evolutionary pathway involving mass transfer from a binary companion, achieve a rotational stability that happens to be sufficient to detect perturbations of spacetime at cosmological scales.

The instrument is the galaxy itself — or rather, our ability to model it. We build a timing model for each pulsar: a comprehensive description of every known effect that influences pulse arrival times. We subtract the model. What remains, the residuals, contains the signal we cannot yet explain. We cross-correlate the residuals of 68 (or 25, or 57, depending on the collaboration) pulsars in pairs, compute the correlation as a function of angular separation, and compare to the Hellings-Downs prediction.

The engineering challenge is not building the detector. It is characterising it. Understanding the noise. Modelling the interstellar medium, which disperses radio pulses in a frequency-dependent way and varies as the pulsar moves through clouds of ionised gas. Accounting for clock errors in the terrestrial time standards used to timestamp observations. Dealing with instrumental noise in each of the different radio telescopes that contribute data. Building a Bayesian framework that can simultaneously fit the timing model parameters, the pulsar noise properties, and the GWB parameters for dozens of pulsars.

This is painstaking, years-long work. The 15-year NANOGrav dataset reflects something like 600 pulsar-years of observation. The evidence is earned.

The Hellings-Downs correlation — the specific pattern sought in those residuals, consistent with the quadrupolar fingerprint of general relativity’s spin-2 gravitational waves — is one of the more beautiful results I have seen in recent astrophysics. It is evidence for that tensor background at frequencies far below LIGO’s band, using a detector assembled by the Milky Way over the course of 10 billion years of stellar evolution and galaxy mergers.

We are in an age of gravitational wave astronomy. I find that remarkable.

If you are interested in the broader theme of using astronomical observations as physics experiments rather than just cataloguing the sky, the posts on transit photometry and the gift of transits and on smartphone-based exoplanet observations cover similar ground at different scales — stellar radii measured in units of stellar radii, by timing the dimming of a star as a planet crosses its face. The underlying logic is the same: precision timing plus a physical model plus statistics equals a measurement of something you could not directly touch.

References

  • Hellings, R. W., & Downs, G. S. (1983). Upper limits on the isotropic gravitational radiation background from pulsar timing analysis. The Astrophysical Journal Letters, 265, L39–L42. DOI: 10.1086/183954

  • Agazie, G., et al. (NANOGrav Collaboration). (2023). The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background. The Astrophysical Journal Letters, 951, L8. DOI: 10.3847/2041-8213/acdac6

  • Antoniadis, J., et al. (EPTA Collaboration). (2023). The second data release from the European Pulsar Timing Array — III. Search for gravitational wave signals. Astronomy & Astrophysics, 678, A50. DOI: 10.1051/0004-6361/202346844

  • Reardon, D. J., et al. (PPTA Collaboration). (2023). Search for an isotropic gravitational-wave background with the Parkes Pulsar Timing Array. The Astrophysical Journal Letters, 951, L6. DOI: 10.3847/2041-8213/acdd02

  • Xu, H., et al. (CPTA Collaboration). (2023). Searching for the nano-Hertz stochastic gravitational wave background with the Chinese Pulsar Timing Array. Research in Astronomy and Astrophysics, 23(7), 075024. DOI: 10.1088/1674-4527/acdfa5

  • Abbott, R., et al. (LIGO Scientific Collaboration). (2024). Observation of Gravitational Waves from the Coalescence of a 2.5–4.5 $M_\odot$ Compact Object and a Neutron Star. The Astrophysical Journal Letters, 970, L34. DOI: 10.3847/2041-8213/ad5beb

  • Agazie, G., et al. (NANOGrav Collaboration). (2025). The NANOGrav 15 yr data set: Posterior predictive checks for gravitational-wave detection with pulsar timing arrays. Physical Review D, 111, 024039. DOI: 10.1103/PhysRevD.111.024039

Changelog

  • 2025-12-01: Corrected the summary from “four independent pulsar timing arrays” to “five” — the CPTA (Chinese Pulsar Timing Array) also published consistent results in June 2023 and is counted as the fifth group in the body text.
  • 2026-07-11: Corrected the 2023-publication account: four papers reported results from NANOGrav, EPTA/InPTA jointly, PPTA, and CPTA. Replaced several “detection” and “direct measurement” formulations with the collaborations’ evidence-level language; corrected the single-wave timing-residual scaling and the O4 end date.