The Dawn of Gravitational Wave Astronomy at Light-year Wavelengths: Insights from Pulsar Timing Arrays

Abstract

Arrays of precisely-timed millisecond pulsars are used to search for gravitational waves with periods of months to decades. Gravitational waves affect the path of radio pulses propagating from a pulsar to Earth, causing the arrival times of those pulses to deviate from expectations based on the physical characteristics of the pulsar system. By correlating these timing residuals in a pulsar timing array (PTA), one can search for a statistically isotropic background of gravitational waves by revealing evidence for a distinctive pattern predicted by General Relativity, known as the Hellings \& Downs curve. On June 29 2023, five regional PTA collaborations announced the first evidence for GWs at light-year wavelengths, predicated on support for this correlation pattern with statistical significances ranging from $\sim\!2-4σ$. The amplitude and shape of the recovered GW spectrum has also allowed many investigations of the expected source characteristics, ranging from a cosmic population of supermassive binary black holes to numerous processes in the early Universe. In the future, we expect to resolve signals from individual binary systems of supermassive black holes, and probe fundamental assumptions about the background, including its polarization, anisotropy, Gaussianity, and stationarity, all of which will aid efforts to discriminate its origin. In tandem with new facilities like DSA-2000 and the SKA, fueling further observations by regional PTAs and the International Pulsar Timing Array, PTAs have extraordinary potential to be engines of nanohertz GW discovery.

Figures (5)

• Figure 1: A gravitational-wave of extragalactic origin---here shown coming from a black-hole binary---enters the Milky Way and passes through an Earth-pulsar system. The train of radio pulses that we measure from the pulsar, and which propagate through the perturbed spacetime represented by the GW wavefronts in blue, are represented in red. With Earth (or the Solar System barycenter) as the origin, the pulsar is at distance $L$ in direction $\hat{p}$, while the GW with wavelength $\lambda$ originates at distance $d$ from direction $-\hat{\Omega}$. The hierarchy of distances is such that $\lambda\ll L\ll d$, where $\lambda$ is of order lightyears (or a parsec), $L$ of order $10^2-10^3$ parsecs, and $d$ is likely greater than a Megaparsec.
• Figure 2: A visual primer on the Hellings & Downs curve. For a given pair of pulsars (green stars), the Hellings & Downs coefficient is given by the correlated response functions of the pulsars, averaged over the (isotropic) distribution of GWB power on the sky. Each panel showing a map displays the correlated responses of a pair of pulsars, which when integrated over the sky, gives the Hellings & Downs coefficient with the matching color in the bottom right panel.
• Figure 3: A summary of current evidence for a nanohertz-frequency gravitational-wave background, assessed by consistency of pulsar-pulsar correlations with the Hellings & Downs curve. Results from NANOGrav 2023ApJ...951L...8A, the EPTA+InPTA 2023AA...678A..50E, PPTA 2023ApJ...951L...6R, and CPTA 2023RAA....23g5024X were all released in June 2023, while the MPTA results 2025MNRAS.536.1489M were released in December 2024.
• Figure 4: The Bayesian spectral recovery of the GWB signal from the EPTA+InPTA, NANOGrav, and the PPTA. On the left, the power spectral density of induced timing residuals is converted into an RMS timing noise excess, and shown using "violins" that represent the posterior spread at each frequency. On the right, these constraints are cast into the parameter space of a power-law spectrum. The dashed vertical then represents the fiducial power-law exponent for the ensemble-averaged power spectral density of timing residuals due to GWs from a population of circular supermassive black-hole binaries, inspiraling due to GW emission. This is $\gamma=13/3$, which is equivalent $h_c\propto f^{-2/3}$. In both panels, a simple product of all PTA constraints approximates the joint spectrum in black dashed on the left, and the black power-law parameters contours on the right. Figure reproduced from 2024ApJ...966..105A.
• Figure 5: Posterior credible contours from a power-law spectral model of the GWB are contrasted with power-law fits to the strain spectra from many synthesized SMBHB populations, using the Holodeck software package (https://github.com/nanograv/holodeck). The exponent of the ensemble-averaged power-law spectrum is shown as a black dashed line ($\gamma=13/3$). Spectral fluctuations from population discreteness could be a significant mitigating factor when comparing the shallowness of the measured GWB spectrum to the ensemble-averaged spectrum prediction. Figure reproduced from 2023ApJ...951L...8A.