Pulsar timing arrays: the emerging gravitational-wave landscape

Abstract

Pulsar Timing Array (PTA) experiments have entered a new era with evidence for a nanoHertz gravitational wave background (GWB). This review describes the physics of detection, detailing the noise models and cross-correlation techniques required to isolate the Hellings-Downs curve. We discuss astrophysical implications, arguing that the perceived tension between current amplitudes and standard merger models is largely resolved by new insights into supermassive black hole binary populations. Beyond the stochastic background, we review the framework for multi-messenger continuous gravitational-wave searches, highlighting targeted search campaigns and rigorous detection protocols. We also examine the potential to probe New Physics, including cosmic strings and ultralight dark matter. Critical challenges are addressed, including small-scale leakage bias in anisotropy searches and the separation of deterministic signals from the GWB and pulsar noise. Finally, we outline the field's future, from rapid data combination strategies to the sensitivity gains expected from the Square Kilometre Array Observatory (SKAO) and DSA-2000.

Figures (23)

• Figure 1: Timeline of key milestones in PTAs. From the first proposals for GW detection via spacecraft Doppler tracking (1975) and pulsar timing (1978--79), through the near-simultaneous development of the Hellings--Downs correlation framework (1983) and the discovery of the first millisecond pulsar (1982), to the construction of the first PTA by Foster & Backer (1990) FosterBacker1990 , the first stringent GWB upper limit by Kaspi et al. (1994) Kaspi1994, Phinney's GWB strain spectrum (2001) phinney_practical_2001, the formation of PTA collaborations, the independent detections of a common red-noise process (2020--2022) NG12p5_gwbGoncharov2021_CRNChen2021Antoniadis2022, and the landmark 2023 evidence for a GWB.
• Figure 2: The GWB in the PTA frequency band. The characteristic strain spectrum of the GWB from a cosmic population of SMBHBs. At low frequencies the background is well approximated by a power law $h_c \propto f^{-2/3}$, representing the Phinney phinney_practical_2001 ensemble average (i.e. the mean over many realizations). Since we observe a single realization of the Universe, the measured spectrum will be jagged, with excursions above and below this mean; the median spectrum (from discrete realizations of the GWB) may better represent what any one realization looks like, see e.g. NG15_discreteness. The number of binaries per frequency bin, $\Delta N(f)$, drops steeply with frequency (eq. \ref{['eq:delta_N']}): from $\sim 10^6$ sources at 2 nHz to $\sim 10^3$ at 20 nHz, because higher-frequency binaries are more massive and rarer, and they evolve through each bin more quickly ($\dot{f} \propto f^{11/3}$). This steep decline marks the transition from a confusion-dominated stochastic regime to one where individual binaries become resolvable as CW sources (sect. \ref{['sec:CW']}) and the background develops measurable anisotropy (sect. \ref{['sec:anisotropy']}).
• Figure 3: Computational frame used to define the ORF(s). Pulsar $a$ lies at distance $L_a$ on the $+\hat{\mathbf{z}}$ axis and pulsar $b$ lies at distance $L_b$ in the $x$--$z$ plane at an angular separation $\zeta$. The GW propagation direction $\hat{\boldsymbol{\Omega}}$ is specified by polar angle $\theta$ and azimuthal angle $\phi$ with respect to $+\hat{\mathbf{z}}$. The principal polarization axes $\hat{m}$ and $\hat{n}$ are perpendicular to $\hat{\boldsymbol{\Omega}}$, with $\hat{m} \times \hat{n} = \hat{\boldsymbol{\Omega}}$, and define the $+$ and $\times$ polarization tensors that enter the antenna beam pattern $F_a^{A}(\hat{\boldsymbol{\Omega}})$ (eq. \ref{['eq:beam_pattern']}). The polarization angle $\psi$ rotates $(\hat{m}, \hat{n})$ about $\hat{\boldsymbol{\Omega}}$ (sect. \ref{['sec:signal_model']}). This geometry determines the antenna pattern functions $F_a^{A}$ and $F_b^{A}$ and hence any spatial correlation.
• Figure 4: Left: The ORF for an isotropic stochastic GWB, called the Hellings and Downs curve. Right: The optimal statistic anholm2009Chamberlin2015 result for the GWB search from Agazie et al. NG15_gwb. The quantity actually plotted on the $y$-axis is $\hat{A}^2 \Gamma_{ab}$, where $\hat{A}^2 \sim 10^{-30}$ is the maximum-likelihood estimator of the squared GWB amplitude; the axis label shows only $\Gamma(\zeta_{ab})$ for presentation clarity. The binning of pulsar pairs by angular separation is for visual clarity only and has no bearing on the result.
• Figure 5: Predicted GWB characteristic strain amplitudes at $f = 1\,\mathrm{yr}^{-1}$ from 33 astrophysical models spanning three decades of work Rajagopal1995Jaffe2003Wyithe2003Sesana2004Enoki2004Sesana2008PTAoccupancySVV2009KocsisSesana2011Barausse2012Sesana13McWilliams2014Kulier2015Ravi2015RosadoSesana2015Roebber2016sesana2016Kelley2017Dvorkin2017Rasskazov2017Bonetti2018Ryu2018chen_constraining_2019Zhu2019Chen2020Siwek2020Izquierdo2021Simon2023Barausse2023Curylo2024SatoPolito2023LiepoldMa2024Chen2025_ASTRIDTillman2026. *CPTA upper limits span three order of magnitude in amplitude and has been truncated for clarity.